# Options for solving ODE systems on GPUs?

I would like to farm out solving systems of ODEs onto GPUs, in a 'trivially parallelisable' setting. For example, doing a sensitivity analysis with 512 different parameter sets.

Ideally I want to do ODE solving with a smart adaptive timestep solver like CVODE, rather than a fixed timestep like Forward Euler, but running it on an NVIDIA GPU instead of CPU.

Has anyone done this? Are there libraries for it?

• Hence the brackets! I'm considering an operator-splitting based technique (cardiac electrophysiology simulations), where you solve ODEs at nodes to get a source term for PDE, then change an ODE parameter for next iteration. – mirams Jan 12 '14 at 9:34
• Maybe related? What's the state of the art in parallel ODE methods? – Kirill Sep 24 '15 at 18:52
• It's important to specify whether you want to use the same time-stepping for every ODE or not. – Christian Clason Sep 30 '15 at 16:15
• Also, if you're specifically interested in the bidomain (or monodomain) equations, you might want to take a look at how CARP does it. – Christian Clason Sep 30 '15 at 16:17
• Different timesteps, if the method is adaptive then it will need them for different parameter sets... thanks for the link to what CARP is doing - a fixed timestep Rush Larsen ODE solver if I read it correctly. – mirams Sep 30 '15 at 22:27

You may want to look into Boost's odeint library and Thrust. They can be combined as discussed here.

• This seems to be a bit different - solving massive ODE systems on the GPU in parallel (with communication). That link says "We have experienced that the vector size over which is parallelized should be of the order of 10^6 to make full use of the GPU.". I'm looking for a nice way of farming out O(10) or O(100) vector sized trivially parallelisable ODE solves... – mirams Jan 19 '14 at 13:05
• Have you thought on writing directly in cuda or openCL? If I undertood right, what you are doing is iterating over some ODE equation in each thread separately, it shouldn't be difficult to write it directly. – Hydro Guy Sep 24 '15 at 13:22
• I imagine it would be possible to code a Forward Euler or other fixed timestep method, where every GPU process uses the same timestep, fairly easily, I'd like to know whether anyone has managed to get adaptive timestepping like CVODE working, or whether this is impossible to make efficient on a GPGPU? – mirams Sep 25 '15 at 7:24
• the problem with gpu is that you need to write data-parallel code. If you write the same adaptive routine but absorbing all that flexibility on the values of some parameters, probably it's possible to code it efficiently on gpu. That also means that you can't use branching on instructions, which is probably what you think that would make it impossible to do it. – Hydro Guy Sep 29 '15 at 3:30
• @mirams there's an example for odeint that covers exactly what you are looking for: boost.org/doc/libs/1_59_0/libs/numeric/odeint/doc/html/…, see also github.com/boostorg/odeint/blob/master/examples/thrust/…. Also, odeint supports adaptive multistep methods as in CVODE: github.com/boostorg/odeint/blob/master/examples/… – Christian Clason Sep 30 '15 at 16:06

DifferentialEquations.jl library is a library for a high level language (Julia) which has tools for automatically transforming the ODE system to an optimized version for parallel solution on GPUs. There are two forms of parallelism that can be employed: array-based parallelism for large ODE systems and parameter parallelism for parameter studies on relatively small (<100) ODE systems. It supports high order implicit and explicit methods and routinely outperforms or matches other systems in benchmarks (at the very least, it wraps the others so it's easy to check and use them!)

For this specific functionality, you may want to take a look at DiffEqGPU.jl which is the module for automated parameter parallelism. The DifferentialEquations.jl library has functionality for parallel parameter studies, and this module augments the existing configurations to make the study happen automatically in parallel. What one does is transforms their existing ODEProblem (or other DEProblem like SDEProblem) into an EnsembleProblem and specify with a prob_func how the other problems are generated from the prototype. The following solves 10,000 trajectories of the Lorenz equation on the GPU with a high order explicit adaptive method:

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)


Notice the user needs to write no GPU code, and with a single RTX 2080 this benchmarks as a 5x improvement over using a 16 core Xeon machine with multithreaded parallelism. One can then check out the README for how to do things like utilize multiple GPUs and doing multiprocessing+GPUs for utilizing a full cluster of GPUs simultaneously. Note that switching to multithreading instead of GPUs is one line change: EnsembleThreads() instead of EnsembleGPUArray().

Then for implicit solvers, the same interface holds. For example, the following uses high order Rosenbrock and implicit Runge-Kutta methods:

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

nothing
end