Spatial ODE simulations

Our introduction to spatial lattice simulations has already provided an extensive description of how to simulate LatticeReactionSystems using ODEs. Further tutorials have also shown how to retrieve values from simulations and or how to plot them. Here we will build on this, primarily discussing strategies for increasing ODE simulation performance. This is especially important for spatial simulations, as these typically are more computationally demanding as compared to non-spatial ones. While focusing on non-spatial simulations, this ODE performance tutorial is also be useful to read.

Solver selection for spatial ODE simulations

Previously we have described how to select ODE solvers, and how this can impact simulation performance. This is especially relevant for spatial simulations. For stiff problems, FBDF is a good first solver to try. For non-stiff problems, ROCK2 is instead a good first alternative. However, it is still worthwhile to explore a range of alternative solvers.

Jacobian options for spatial ODE simulations

We have previously described how, when implicit solvers are used to solve stiff ODEs, the strategy for computing the system Jacobian is important. This is especially the case for spatial simulations, where the Jacobian often is large and highly sparse. Catalyst implements special methods for spatial Jacobians. To utilise these, provide the jac = true argument to your ODEProblem when it is created (if jac = false, which is the default, automatic differentiation will be used for Jacobian computation). Here we simulate a Brusselator while designating to use Catalyst's computed Jacobian:

using Catalyst, OrdinaryDiffEqBDF
brusselator = @reaction_network begin
    A, ∅ --> X
    1, 2X + Y --> 3X
    B, X --> Y
    1, X --> ∅
end
diffusion_rx = @transport_reaction D X
lattice = CartesianGrid((20,20))
lrs = LatticeReactionSystem(brusselator, [diffusion_rx], lattice)

u0 = [:X => rand(20, 20), :Y => 10.0]
tspan = (0.0, 1.0)
ps = [:A => 1.0, :B => 4.0, :D => 0.2]
oprob = ODEProblem(lrs, u0, tspan, ps; jac = true)
sol = solve(oprob, FBDF())

For large systems, building a dense Jacobian can be problematic, in which case a sparse Jacobian can be designated using sparse = true:

oprob = ODEProblem(lrs, u0, tspan, ps; jac = true, sparse = true)
sol = solve(oprob, FBDF())

It is possible to use sparse = true while jac = false, in which case a sparse Jacobian is computed using automatic differentiation.