# CVODE\_BDF outperforms Julia solvers for stiff system biology model?

**URL:** <https://discourse.julialang.org/t/cvode-bdf-outperforms-julia-solvers-for-stiff-system-biology-model/113936>\
**Category:** Modelling & Simulations\
**Tags:** question\
**Created:** [May 7, 2024, 9:19am UTC](https://discourse.julialang.org/t/cvode-bdf-outperforms-julia-solvers-for-stiff-system-biology-model/113936 "2024-05-07T09:19:00Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![sebapersson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sebapersson/32/47456_2.png) [@sebapersson](https://discourse.julialang.org/u/sebapersson)\
**Post date:** [May 7, 2024, 9:19am UTC](https://discourse.julialang.org/t/cvode-bdf-outperforms-julia-solvers-for-stiff-system-biology-model/113936/1 "2024-05-07T09:19:00Z")

</div>

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](https://www.embopress.org/doi/full/10.1038/msb.2011.50). The model files can be found [here](https://github.com/Benchmarking-Initiative/Benchmark-Models-PEtab/tree/master/Benchmark-Models/Bachmann_MSB2011). I ran the benchmarks using Julia 1.10.3

```julia
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](https://global.discourse-cdn.com/julialang/original/3X/2/2/22eca0a484e7ae15fef82ec8893266111326ef96.png)

![Plot1](https://global.discourse-cdn.com/julialang/original/3X/5/1/51099a1ae94e38de3d8c1e21eeb0286e50ae2c45.png)

![Plot3](https://global.discourse-cdn.com/julialang/original/3X/1/b/1b6405d4d81a3430bf2953e805bc2810fc3ab752.png)

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

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [May 7, 2024, 11:01am UTC](https://discourse.julialang.org/t/cvode-bdf-outperforms-julia-solvers-for-stiff-system-biology-model/113936/2 "2024-05-07T11:01:05Z")

</div>

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:

```julia
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](https://global.discourse-cdn.com/julialang/original/3X/f/2/f27db4a408701feef813da04bc5b6aead9a8dcd4.png)  
 ![plot2](https://global.discourse-cdn.com/julialang/original/3X/6/2/626abc90e039246ad3bfdd99ade762d9751bf9dc.png)

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
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.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [May 7, 2024, 11:03am UTC](https://discourse.julialang.org/t/cvode-bdf-outperforms-julia-solvers-for-stiff-system-biology-model/113936/3 "2024-05-07T11:03:52Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![sebapersson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sebapersson/32/47456_2.png) [@sebapersson](https://discourse.julialang.org/u/sebapersson)\
**Post date:** [May 7, 2024, 11:18am UTC](https://discourse.julialang.org/t/cvode-bdf-outperforms-julia-solvers-for-stiff-system-biology-model/113936/4 "2024-05-07T11:18:52Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [May 8, 2024, 3:05pm UTC](https://discourse.julialang.org/t/cvode-bdf-outperforms-julia-solvers-for-stiff-system-biology-model/113936/5 "2024-05-08T15:05:31Z")

</div>

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.

> <https://github.com/SciML/OrdinaryDiffEq.jl/pull/2183>
>
> CVODE estimates the convergence of the first iteration by reusing the
> convergen…ce rate of the previous nonlinear solver iteration.
> 
> Reference: https://github.com/LLNL/sundials/blob/2abd63bd6cbc354fb4861bba8e98d0b95d65e24a/src/cvodes/cvodes\_nls.c#L325-L331
> 
> MWE: https://discourse.julialang.org/t/cvode-bdf-outperforms-julia-solvers-for-stiff-system-biology-model/113936
> 
> Master branch:
> \`\`\`julia
> julia\> sol = solve(oprob\_jac, FBDF()).stats
> SciMLBase.DEStats
> Number of function 1 evaluations: 311
> Number of function 2 evaluations: 0
> Number of W matrix evaluations: 31
> Number of linear solves: 222
> Number of Jacobians created: 2
> Number of nonlinear solver iterations: 212
> 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: 84
> Number of rejected steps: 2
> \`\`\`
> 
> This PR:
> \`\`\`julia
> julia\> sol = solve(oprob\_jac, FBDF()).stats
> SciMLBase.DEStats
> Number of function 1 evaluations: 251
> Number of function 2 evaluations: 0
> Number of W matrix evaluations: 32
> Number of linear solves: 159
> Number of Jacobians created: 2
> Number of nonlinear solver iterations: 149
> 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: 86
> Number of rejected steps: 3
> \`\`\`

That should merge by the end of the week.

---

<div class="post-metadata">

**Author:** ![sebapersson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sebapersson/32/47456_2.png) [@sebapersson](https://discourse.julialang.org/u/sebapersson)\
**Post date:** [May 8, 2024, 3:34pm UTC](https://discourse.julialang.org/t/cvode-bdf-outperforms-julia-solvers-for-stiff-system-biology-model/113936/6 "2024-05-08T15:34:19Z")

</div>

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