CVODE_BDF outperforms Julia solvers for stiff system biology model?

Specific Domains Modelling & Simulations
1

For medium sized stiff system biology ODE models I have recently done two benchmarks that suggest that CVODE_BDF outperform the Julia solvers (in particular the multi-step bdf solvers).

Below is an with 25 ODEs and 39 parameters. The model is imported with SBMLImporter.jl and I used PEtab.jl to set the parameters to the values in the original publication. The model files can be found here. I ran the benchmarks using Julia 1.10.3

using DiffEqBase, OrdinaryDiffEq, Catalyst, SBMLImporter, Sundials, Plots, DiffEqDevTools,
      TimerOutputs, LinearAlgebra, ModelingToolkit, BenchmarkTools, PEtab

# Read as ODESystem
prnbng, cb = load_SBML(joinpath(@__DIR__, "Bachmann_MSB2011", "model_Bachmann_MSB2011.xml"),
                       check_massaction=false)
osys = convert(ODESystem, prnbng.rn)

# Read as PEtabODEProblem and extract parameters for first simulations condition
petab_model = PEtabModel(joinpath(@__DIR__, "Bachmann_MSB2011", "Bachmann_MSB2011.yaml"), verbose=false)
petab_prob = PEtabODEProblem(petab_model, verbose=false)
_p = PEtab.get_ps(petab_prob.θ_nominalT, petab_prob)
_u0 = PEtab.get_u0(petab_prob.θ_nominalT, petab_prob)

tf = 100.0
tspan = (0.0, tf)
oprob = ODEProblem{true, SciMLBase.FullSpecialize}(osys, _p, tspan, _u0)
oprob_jac = ODEProblem{true, SciMLBase.FullSpecialize}(osys, _p, tspan, _u0, jac=true)
oprob_sparse = ODEProblem{true, SciMLBase.FullSpecialize}(osys, _p, tspan, _u0, jac=true, sparse=true)

sol = solve(oprob, CVODE_BDF(), abstol=1/10^14, reltol=1/10^14)
test_sol  = TestSolution(sol)

abstols = 1.0 ./ 10.0 .^ (6:10)
reltols = 1.0 ./ 10.0 .^ (6:10)

setups = [
          Dict(:alg=>CVODE_BDF()),
          Dict(:alg=>TRBDF2()),
          Dict(:alg=>QNDF()),
          Dict(:alg=>FBDF()),
          Dict(:alg=>KenCarp4()),
          Dict(:alg=>Rodas4()),
          Dict(:alg=>Rodas5P())
          ];

wp = WorkPrecisionSet(oprob, abstols, reltols, setups; error_estimate=:l2,
                      saveat=tf/10., appxsol=test_sol, maxiters=Int(1e9), numruns=100)
plot(wp)
title!("No Jacobian")

wp_jac = WorkPrecisionSet(oprob_jac, abstols, reltols, setups; error_estimate=:l2,
                          saveat=tf/10., appxsol=test_sol, maxiters=Int(1e9), numruns=100)
plot(wp_jac)
title!("With Jacobian")

setups[1] = Dict(:alg=>CVODE_BDF(linear_solver=:KLU))
wp_sparse = WorkPrecisionSet(oprob_sparse, abstols, reltols, setups; error_estimate=:l2,
                             saveat=tf/10., appxsol=test_sol, maxiters=Int(1e9), numruns=100)
plot(wp_sparse)
title!("Sparse Jacobian")

With output:
Plot1

Plot1

Plot3

Is there any additional performance option I am missing for the Julia solvers?

2

In this intermediate size the extrapolation methods will dominate. They just need a better interpolation function though. What I mean by this is easy to see with the following:

setups = [
          Dict(:alg=>CVODE_BDF()),
          Dict(:alg=>QNDF()),
          Dict(:alg=>FBDF()),
          Dict(:alg=>KenCarp4()),
          Dict(:alg=>RadauIIA5()),
          Dict(:alg=>ImplicitHairerWannerExtrapolation(min_order = 5, init_order = 3, threading = OrdinaryDiffEq.PolyesterThreads())),
          Dict(:alg=>Rodas5P()),
          ];

wp = WorkPrecisionSet(oprob, abstols, reltols, setups;
          saveat=tf/10., appxsol=test_sol, maxiters=Int(1e9), numruns=100)
p1 = plot(wp, legend=false);

wp = WorkPrecisionSet(oprob, abstols, reltols, setups; error_estimate=:l2,
                      saveat=tf/10., tstops=10:10:100, appxsol=test_sol, maxiters=Int(1e9), numruns=100)
p2 = plot(wp, legend=false)

wp = WorkPrecisionSet(oprob, abstols, reltols, setups;  error_estimate=:l2,
                      saveat=tf/10., appxsol=test_sol, maxiters=Int(1e9), numruns=100)
p3 = plot(wp, legend=false)

wp = WorkPrecisionSet(oprob, abstols, reltols, setups;
                      dt = 0.5, saveat=tf/10., appxsol=test_sol, maxiters=Int(1e9), numruns=100)
p4 = plot(wp)
plot(p1,p2,p3,p4)
savefig("plot1.png")

wp = WorkPrecisionSet(oprob, abstols, reltols, setups;
                      dt = 0.5, saveat=tf/10., appxsol=test_sol, maxiters=Int(1e9), numruns=100)
p5 = plot(wp)
savefig("plot2.png")

plot1
plot2

So using an implicit extrapolation method and pushing it to be aggressive will pretty much always be best in this 20-200 stiff ODE range. I’ll need to derive a better interpolation and then it should shine better.

As for what’s going on with the nonlinear solver though, that’s a bit of a mystery right now.

julia> sol = solve(oprob_jac, FBDF(), abstol=1e-6, reltol=1e-6).stats
SciMLBase.DEStats
Number of function 1 evaluations:                  634
Number of function 2 evaluations:                  0
Number of W matrix evaluations:                    34
Number of linear solves:                           436
Number of Jacobians created:                       2
Number of nonlinear solver iterations:             426
Number of nonlinear solver convergence failures:   0
Number of fixed-point solver iterations:                     0
Number of fixed-point solver convergence failures:           0
Number of rootfind condition calls:                0
Number of accepted steps:                          192
Number of rejected steps:                          3

julia> sol = solve(oprob_jac, CVODE_BDF(), abstol=1e-6, reltol=1e-6).stats
SciMLBase.DEStats
Number of function 1 evaluations:                  227
Number of function 2 evaluations:                  0
Number of W matrix evaluations:                    24
Number of linear solves:                           0
Number of Jacobians created:                       4
Number of nonlinear solver iterations:             224
Number of nonlinear solver convergence failures:   0
Number of fixed-point solver iterations:                     0
Number of fixed-point solver convergence failures:           0
Number of rootfind condition calls:                0
Number of accepted steps:                          183
Number of rejected steps:                          3

The nonlinear solver tolerance seems a bit too strict, not sure why.

But note right now is probably the worst time to benchmark since the NonlinearSolve swap is going and bound to invalidate all benchmarks in the next few weeks.

3

Can you just PR to add the whole set to SciMLBenchmarks? Focusing on just two examples is always a bit misleading (though those 2 show that there is something interesting, possibly a regression), and we might as well let the benchmark machine generate the whole set for a better view of it all.

1 Like
4

I can make a PR with the stiff models in the PEtab benchmark collection, to see if there is any trend.

5

We found what CVODE is doing differently here. Turns out that on the chosen problem the nonlinear solvers seem to always converge in one iteration, while our method always did at least two iterations to check the convergence rate. We just need to add a fast path for that behavior.

That should merge by the end of the week.

2 Likes
6

Sounds great! Looking forward to then test on some other problems when the PR is merged.