BCR Work-Precision Diagrams

Chris Rackauckas

The following benchmark is of a 1122 ODEs with 24388 terms that describe a stiff chemical reaction network.

using ReactionNetworkImporters, OrdinaryDiffEq, DiffEqBiological,
      Sundials, Plots, DiffEqDevTools, ODEInterface, ODEInterfaceDiffEq,
      LSODA, TimerOutputs

gr()
prnbng = loadrxnetwork(BNGNetwork(), "BNGRepressilator",
                       joinpath(dirname(pathof(ReactionNetworkImporters)),"..","data","bcr","bcr.net"))
Parsing parameters...done
Adding parameters...done
Parsing species...done
Adding species...done
Parsing and adding reactions...done
Parsing groups...done
rn = deepcopy(prnbng.rn)
addodes!(rn; build_jac=false, build_symfuncs=false, build_paramjac=false)
tf = 100000.0
oprob = ODEProblem(rn, prnbng.u₀, (0.,tf), prnbng.p);

densejac_rn = deepcopy(prnbng.rn)
# zeroout_jac=true is needed to keep the generated expressions from being too big for the compiler
addodes!(densejac_rn; build_jac=true, zeroout_jac = true, sparse_jac = false, build_symfuncs=false, build_paramjac=false)
densejacprob = ODEProblem(densejac_rn, prnbng.u₀, (0.,tf), prnbng.p);

sparsejac_rn = deepcopy(prnbng.rn)
addodes!(sparsejac_rn; build_jac=true, sparse_jac = true, build_symfuncs=false, build_paramjac=false)
sparsejacprob = ODEProblem(sparsejac_rn, prnbng.u₀, (0.,tf), prnbng.p);
@show numspecies(rn) # Number of ODEs
numspecies(rn) = 1122
@show numreactions(rn) # Apprx. number of terms in the ODE
numreactions(rn) = 24388
@show numparams(rn) # Number of Parameters
numparams(rn) = 128
128

Time ODE derivative function compilation

As compiling the ODE derivative functions has in the past taken longer than running a simulation, we first force compilation by evaluating these functions one time.

const to = TimerOutput()
u₀ = prnbng.u₀
u = copy(u₀);
du = similar(u);
p = prnbng.p
@timeit to "ODERHS Eval1" rn.f(du,u,p,0.)
@timeit to "ODERHS Eval2" rn.f(du,u,p,0.)
sparsejac_rn.f(du,u,p,0.)

J = zeros(length(u),length(u))
@timeit to "DenseJac Eval1" densejac_rn.jac(J,u,p,0.)
@timeit to "DenseJac Eval2" densejac_rn.jac(J,u,p,0.)

Js = similar(sparsejac_rn.odefun.jac_prototype)
@timeit to "SparseJac Eval1" sparsejac_rn.jac(Js,u,p,0.)
@timeit to "SparseJac Eval2" sparsejac_rn.jac(Js,u,p,0.)
show(to)
──────────────────────────────────────────────────────────────────────────
                                   Time                   Allocations      
                           ──────────────────────   ───────────────────────
     Tot / % measured:           129s / 84.8%           17.2GiB / 80.2%    

 Section           ncalls     time   %tot     avg     alloc   %tot      avg
 ──────────────────────────────────────────────────────────────────────────
 DenseJac Eval1         1    53.8s  49.3%   53.8s   5.29GiB  38.3%  5.29GiB
 SparseJac Eval1        1    36.1s  33.1%   36.1s   5.23GiB  37.8%  5.23GiB
 ODERHS Eval1           1    19.3s  17.7%   19.3s   3.31GiB  23.9%  3.31GiB
 DenseJac Eval2         1    976μs  0.00%   976μs     32.0B  0.00%    32.0B
 ODERHS Eval2           1   95.0μs  0.00%  95.0μs     32.0B  0.00%    32.0B
 SparseJac Eval2        1   53.0μs  0.00%  53.0μs     32.0B  0.00%    32.0B
 ──────────────────────────────────────────────────────────────────────────

Picture of the solution

sol = solve(oprob, CVODE_BDF(), saveat=tf/1000., reltol=1e-5, abstol=1e-5)
plot(sol,legend=false, fmt=:png)