Opcje rozwiązywania systemów ODE na GPU?

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?

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

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.