Opcje rozwiązywania systemów ODE na GPU?

16

Chciałbym wdrożyć systemy rozwiązywania ODE na GPU w „trywialnie równoległym” ustawieniu. Na przykład, wykonując analizę wrażliwości z 512 różnymi zestawami parametrów.

Idealnie chciałbym wykonać rozwiązywanie ODE za pomocą inteligentnego adaptacyjnego solvera pomiaru czasu, takiego jak CVODE, zamiast stałego pomiaru czasu, takiego jak Forward Euler, ale działającego na GPU NVIDIA zamiast na procesorze.

Czy ktoś to zrobił? Czy są do tego biblioteki?

miramy
źródło
Stąd nawiasy! Rozważam technikę opartą na podziale operatora (symulacje elektrofizjologii serca), w której rozwiązujesz ODE w węzłach, aby uzyskać termin źródłowy dla PDE, a następnie zmieniasz parametr ODE do następnej iteracji.
mirams,
1
Ważne jest, aby sprecyzować, czy chcesz używać tego samego odstępu czasu dla każdej ODE, czy nie.
Christian Clason
Ponadto, jeśli jesteś szczególnie zainteresowany równaniami bidomain (lub monodomen), możesz rzucić okiem na to, jak robi to CARP .
Christian Clason
Różne przebiegi czasowe, jeśli metoda jest adaptacyjna, będzie potrzebować ich dla różnych zestawów parametrów ... dzięki za link do tego, co robi CARP - ustalony pomiar czasu solver Rush Larsena ODE, jeśli poprawnie go odczytam.
mirams

Odpowiedzi:

6

Możesz zajrzeć do biblioteki Odeint Boost i Thrust . Można je łączyć, jak omówiono tutaj .

Juan M. Bello-Rivas
źródło
To wydaje się być nieco inne - równoległe rozwiązywanie ogromnych systemów ODE na GPU (z komunikacją). Ten link mówi: „Przekonaliśmy się, że rozmiar wektora, który jest równoległy, powinien być rzędu 10 ^ 6, aby w pełni wykorzystać procesor graficzny”. Szukam fajnego sposobu
wyhodowania
Czy myślałeś o pisaniu bezpośrednio w cuda lub openCL? Jeśli dobrze zrozumiem, to, co robisz, to iteracja po pewnym równaniu ODE w każdym wątku osobno, nie powinno być trudno napisać go bezpośrednio.
Hydro Guy
Wyobrażam sobie, że można kodować Forward Euler lub inną stałą metodę pomiaru czasu, w którym każdy proces GPU używa tego samego pomiaru czasu, dość łatwo, chciałbym wiedzieć, czy ktoś zdołał uzyskać adaptacyjne odmierzanie czasu, takie jak CVODE, lub czy to nie można uczynić wydajnym GPGPU?
mirams,
Problem z GPU polega na tym, że musisz pisać kod równoległy do ​​danych. Jeśli napiszesz tę samą procedurę adaptacyjną, ale absorbujesz całą elastyczność na wartościach niektórych parametrów, prawdopodobnie możliwe jest wydajne kodowanie jej na gpu. Oznacza to również, że nie możesz używać rozgałęzień na instrukcjach, co prawdopodobnie uważasz, że uniemożliwiłoby to.
Hydro Guy
1
@mirams istnieje przykład odeint, który obejmuje dokładnie to, czego szukasz: boost.org/doc/libs/1_59_0/libs/numeric/odeint/doc/html/… , patrz także github.com/boostorg/odeint/blob/ master / przykłady / thrust /… . Ponadto odeint obsługuje adaptacyjne metody wieloetapowe, jak w CVODE: github.com/boostorg/odeint/blob/master/examples/…
Christian Clason
13

Biblioteka DifferentialEquations.jl jest biblioteką dla języka wysokiego poziomu (Julia), która ma narzędzia do automatycznej transformacji systemu ODE w zoptymalizowaną wersję dla równoległego rozwiązania na GPU. Istnieją dwie formy równoległości, które można zastosować: równoległość tablicowa dla dużych systemów ODE i równoległość parametrów do badań parametrów na stosunkowo małych (<100) systemach ODE. Obsługuje niejawne i jawne metody wyższego rzędu oraz rutynowo osiąga lepsze wyniki lub dopasowuje inne systemy w testach porównawczych (przynajmniej otacza inne systemy , dzięki czemu można je łatwo sprawdzić i użyć!)

Aby zobaczyć tę konkretną funkcjonalność, możesz zajrzeć na DiffEqGPU.jl, który jest modułem automatycznej równoległości parametrów. Biblioteka DifferentialEquations.jl ma funkcjonalność dla równoległych badań parametrów , a ten moduł rozszerza istniejące konfiguracje, aby badanie odbywało się automatycznie równolegle. To, co robi, przekształca ich istniejące ODEProblem(lub inne DEProblempodobne SDEProblem) w EnsembleProblemi określa, w prob_funcjaki sposób inne problemy są generowane z prototypu. Poniższe rozwiązuje 10 000 trajektorii równania Lorenza na GPU za pomocą jawnej metody adaptacyjnej wysokiego rzędu:

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)

Zauważ, że użytkownik nie musi pisać kodu GPU, a przy pojedynczym RTX 2080 ten test porównawczy stanowi 5-krotną poprawę w porównaniu z użyciem 16-rdzeniowej maszyny Xeon z wielowątkową równoległością. Następnie można sprawdzić plik README, aby dowiedzieć się, jak korzystać z wielu procesorów graficznych i robić procesory wieloprocesorowe + procesory graficzne w celu jednoczesnego wykorzystania pełnego klastra procesorów graficznych . Zauważ, że przejście na wielowątkowość zamiast GPU to zmiana jednej linii: EnsembleThreads()zamiast EnsembleGPUArray().

Następnie dla niejawnych solverów obowiązuje ten sam interfejs. Na przykład, w poniższym przykładzie zastosowano metody Rosenbrocka wysokiego poziomu i niejawne metody Runge-Kutty:

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)

Chociaż ten formularz wymaga podania jakobianu w celu użycia go na GPU (obecnie zostanie to naprawione wkrótce), dokumentacja DifferentialEquations.jl pokazuje, jak wykonywać automatyczne symboliczne obliczenia jakobianowe na funkcjach zdefiniowanych liczbowo , więc nadal nie ma instrukcji poród tutaj. Zdecydowanie poleciłbym te algorytmy, ponieważ logika rozgałęziania metody takiej jak CVODE generalnie powoduje desynchronizację wątków i wydaje się nie działać tak dobrze, jak metoda Rosenbrock w tego typu scenariuszach.

Korzystając z DifferentialEquations.jl, masz również dostęp do pełnej biblioteki, która zawiera funkcje takie jak globalna analiza wrażliwości, która może korzystać z tego przyspieszenia GPU. Jest także kompatybilny z podwójnymi liczbami w celu szybkiej lokalnej analizy czułości . Kod oparty na GPU ma wszystkie funkcje DifferentialEquations.jl, takie jak obsługa zdarzeń i duży zestaw solverów ODE, które są zoptymalizowane pod kątem różnych rodzajów problemów , co oznacza, że ​​nie jest to zwykły jednorazowy solver ODE GPU, ale zamiast tego część w pełni funkcjonalnego systemu, który ma również wydajną obsługę GPU.

Chris Rackauckas
źródło