Opsi untuk menyelesaikan sistem ODE pada GPU?

15

Saya ingin menerapkan sistem penyelesaian ODE ke GPU, dalam pengaturan yang 'sepele parallelisable'. Misalnya, melakukan analisis sensitivitas dengan 512 set parameter yang berbeda.

Idealnya saya ingin melakukan penyelesaian ODE dengan solver pencatat waktu adaptif yang cerdas seperti CVODE, daripada pencatat waktu yang tetap seperti Forward Euler, tetapi menjalankannya pada GPU NVIDIA alih-alih CPU.

Adakah yang melakukan ini? Apakah ada perpustakaan untuk itu?

keajaiban
sumber
Karena itu kurung! Saya sedang mempertimbangkan teknik berbasis operator-splitting (simulasi electrophysiology jantung), di mana Anda memecahkan ODE pada node untuk mendapatkan istilah sumber untuk PDE, kemudian mengubah parameter ODE untuk iterasi berikutnya.
mirams
1
Sangat penting untuk menentukan apakah Anda ingin menggunakan loncatan waktu yang sama untuk setiap ODE atau tidak.
Christian Clason
Juga, jika Anda secara khusus tertarik pada persamaan bidomain (atau monodomain), Anda mungkin ingin melihat bagaimana CARP melakukannya .
Christian Clason
Stempel waktu yang berbeda, jika metode ini adaptif maka akan membutuhkannya untuk set parameter yang berbeda ... terima kasih atas tautannya ke apa yang dilakukan CARP - pemecah stempel waktu Rush Larsen ODE jika saya membacanya dengan benar.
mirams

Jawaban:

6

Anda mungkin ingin melihat ke perpustakaan odeint Boost dan Thrust . Mereka dapat digabungkan sebagaimana dibahas di sini .

Juan M. Bello-Rivas
sumber
Ini tampaknya sedikit berbeda - memecahkan sistem ODE besar pada GPU secara paralel (dengan komunikasi). Tautan itu mengatakan, "Kami telah mengalami bahwa ukuran vektor yang diparalelkan harus dalam urutan 10 ^ 6 untuk memanfaatkan GPU secara penuh." Saya sedang mencari cara yang bagus untuk keluar dari O (10) atau O (100) berukuran vektor yang dapat dipecahkan secara parallelisable ODE ...
mirams
Pernahkah Anda berpikir untuk menulis langsung di cuda atau openCL? Jika saya benar, apa yang Anda lakukan adalah mengulangi beberapa persamaan ODE di setiap utas secara terpisah, seharusnya tidak sulit untuk menulisnya secara langsung.
Hydro Guy
Saya membayangkan akan mungkin untuk membuat kode Forward Euler atau metode cap waktu tetap lainnya, di mana setiap proses GPU menggunakan catatan waktu yang sama, cukup mudah, saya ingin tahu apakah ada yang berhasil mendapatkan catatan waktu adaptif seperti CVODE berfungsi, atau apakah ini tidak mungkin membuat efisien pada GPGPU?
mirams
masalah dengan GPU adalah bahwa Anda perlu menulis kode data-paralel. Jika Anda menulis rutin adaptif yang sama tetapi menyerap semua fleksibilitas pada nilai-nilai beberapa parameter, mungkin dimungkinkan untuk mengkodekannya secara efisien pada GPU. Itu juga berarti bahwa Anda tidak dapat menggunakan percabangan pada instruksi, yang mungkin menurut Anda tidak memungkinkan untuk melakukannya.
Hydro Guy
1
@mirams ada contoh untuk odeint yang mencakup persis apa yang Anda cari: boost.org/doc/libs/1_59_0/libs/numeric/odeint/doc/html/… , lihat juga github.com/boostorg/odeint/blob/ master / contoh / dorong / ... . Juga, odeint mendukung metode multistep adaptif seperti pada CVODE: github.com/boostorg/odeint/blob/master/examples/…
Christian Clason
12

Perpustakaan DifferentialEquations.jl adalah perpustakaan untuk bahasa tingkat tinggi (Julia) yang memiliki alat untuk secara otomatis mengubah sistem ODE ke versi yang dioptimalkan untuk solusi paralel pada GPU. Ada dua bentuk paralelisme yang dapat digunakan: paralelisme berbasis array untuk sistem ODE besar dan paralelisme parameter untuk studi parameter pada sistem ODE yang relatif kecil (<100). Ini mendukung metode tersirat dan eksplisit urutan tinggi dan secara rutin mengungguli atau mencocokkan sistem lain dalam tolok ukur (setidaknya, itu membungkus yang lain sehingga mudah untuk memeriksa dan menggunakannya!)

Untuk fungsi khusus ini, Anda mungkin ingin melihat DiffEqGPU.jl yang merupakan modul untuk paralelisme parameter otomatis. Pustaka DifferentialEquations.jl memiliki fungsi untuk studi parameter paralel , dan modul ini menambah konfigurasi yang ada untuk membuat studi terjadi secara otomatis secara paralel. Apa yang dilakukan adalah mengubah yang sudah ada ODEProblem(atau DEProblemsejenisnya SDEProblem) menjadi EnsembleProblemdan tentukan dengan prob_funcbagaimana masalah lain dihasilkan dari prototipe. Berikut ini memecahkan 10.000 lintasan persamaan Lorenz pada GPU dengan metode adaptif eksplisit orde tinggi:

using OrdinaryDiffEq, DiffEqGPU
function lorenz(du,u,p,t)
 @inbounds begin
     du[1] = p[1]*(u[2]-u[1])
     du[2] = u[1]*(p[2]-u[3]) - u[2]
     du[3] = u[1]*u[2] - p[3]*u[3]
 end
 nothing
end

u0 = Float32[1.0;0.0;0.0]
tspan = (0.0f0,100.0f0)
p = (10.0f0,28.0f0,8/3f0)
prob = ODEProblem(lorenz,u0,tspan,p)
prob_func = (prob,i,repeat) -> remake(prob,p=rand(Float32,3).*p)
monteprob = EnsembleProblem(prob, prob_func = prob_func)
@time sol = solve(monteprob,Tsit5(),EnsembleGPUArray(),trajectories=10_000,saveat=1.0f0)

Perhatikan bahwa pengguna tidak perlu menulis kode GPU, dan dengan satu RTX 2080 tolok ukur ini sebagai peningkatan 5x dibandingkan menggunakan mesin 16 inti Xeon dengan paralelisme multithreaded. Satu kemudian dapat memeriksa README untuk bagaimana melakukan hal-hal seperti memanfaatkan beberapa GPU dan melakukan multiprocessing + GPU untuk memanfaatkan sekelompok penuh GPU secara bersamaan . Perhatikan bahwa beralih ke multithreading alih-alih GPU adalah perubahan satu baris: EnsembleThreads()alih-alih EnsembleGPUArray().

Kemudian untuk pemecah implisit, antarmuka yang sama berlaku. Misalnya, berikut ini menggunakan Rosenbrock orde tinggi dan metode Runge-Kutta implisit:

function lorenz_jac(J,u,p,t)
 @inbounds begin
     σ = p[1]
     ρ = p[2]
     β = p[3]
     x = u[1]
     y = u[2]
     z = u[3]
     J[1,1] = -σ
     J[2,1] = ρ - z
     J[3,1] = y
     J[1,2] = σ
     J[2,2] = -1
     J[3,2] = x
     J[1,3] = 0
     J[2,3] = -x
     J[3,3] = -β
 end
 nothing
end

function lorenz_tgrad(J,u,p,t)
 nothing
end

func = ODEFunction(lorenz,jac=lorenz_jac,tgrad=lorenz_tgrad)
prob_jac = ODEProblem(func,u0,tspan,p)
monteprob_jac = EnsembleProblem(prob_jac, prob_func = prob_func)

@time solve(monteprob_jac,Rodas5(linsolve=LinSolveGPUSplitFactorize()),EnsembleGPUArray(),dt=0.1,trajectories=10_000,saveat=1.0f0)
@time solve(monteprob_jac,TRBDF2(linsolve=LinSolveGPUSplitFactorize()),EnsembleGPUArray(),dt=0.1,trajectories=10_000,saveat=1.0f0)

Meskipun formulir ini mengharuskan Anda memberikan Jacobian agar dapat digunakan pada GPU (saat ini, akan segera diperbaiki), dokumentasi DifferentialEquations.jl menunjukkan cara melakukan perhitungan Jacobian simbolis simbolis pada fungsi yang ditentukan secara numerik , sehingga masih belum ada manual kerja di sini. Saya akan sangat merekomendasikan algoritma ini karena logika percabangan dari metode seperti CVODE umumnya menyebabkan desync thread dan tampaknya tidak berkinerja sebaik metode Rosenbrock dalam jenis skenario ini.

Dengan menggunakan DifferentialEquations.jl, Anda juga mendapatkan akses ke perpustakaan lengkap, yang mencakup fungsionalitas seperti analisis sensitivitas global yang dapat memanfaatkan akselerasi GPU ini. Ini juga kompatibel dengan angka ganda untuk analisis sensitivitas lokal cepat . Kode berbasis GPU mendapatkan semua fitur DifferentialEquations.jl, seperti penanganan acara dan serangkaian besar pemecah ODE yang dioptimalkan untuk berbagai jenis masalah , yang berarti bukan hanya pemecah ODE GPU satu-kali yang sederhana tetapi sebagai bagian dari sistem berfitur lengkap yang kebetulan juga memiliki dukungan GPU yang efisien.

Chris Rackauckas
sumber