Use sparse Jacobians for large ODE sizes #9
Description
When the size of the (perturbation) ODE system increases, solving the linear system at every time step becomes a bottleneck of the implicit ODE solver. The system size increases quickly particularly when
- the multipole cutoff
$l_\text{max}$ increases (for e.g. photons, massless neutrinos and massive neutrinos), - the number of momentum bins increases (e.g. for massive neutrinos).
It should be possible to speed up the linear solves by exploiting the sparsity pattern of the Jacobian (which contains mostly zeros) when the system size exceeds some threshold under which dense methods are faster:
using SymBoltz, ModelingToolkit, Plots
M = ΛCDM(lmax = 50)
Ms = mtkcompile(M)
sparsity = ModelingToolkit.jacobian_sparsity(Ms)
nonzerofrac = length(findall(!iszero, sparsity)) / length(sparsity)
spy(sparsity; markersize = 5, title = "$(nonzerofrac*100) % nonzeros")
Some resources: