bench_ODEs.jl

using OrdinaryDiffEq, DiffEqDevTools, Plots, ParameterizedFunctions, Sundials, ODEInterfaceDiffEq
using LinearAlgebra
BLAS.set_num_threads(1)
postfix = "PR"
setups = [
          Dict(:alg=>Kvaerno5()),
          Dict(:alg=>KenCarp4()),
          Dict(:alg=>radau5()),
          Dict(:alg=>Rodas4()),
          Dict(:alg=>Rodas5()),
          Dict(:alg=>Rosenbrock23()),
          Dict(:alg=>TRBDF2()),
          Dict(:alg=>CVODE_BDF()),
          Dict(:alg=>RadauIIA5()),
         ]
names = ["Kvaerno5" "KenCarp4" "radau5" "Rodas4" "Rodas5" "Rosenbrock23" "TRBDF2" "CVODE_BDF" "RadauIIA5"]
hires = @ode_def Hires begin
  dy1 = -1.71*y1 + 0.43*y2 + 8.32*y3 + 0.0007
  dy2 = 1.71*y1 - 8.75*y2
  dy3 = -10.03*y3 + 0.43*y4 + 0.035*y5
  dy4 = 8.32*y2 + 1.71*y3 - 1.12*y4
  dy5 = -1.745*y5 + 0.43*y6 + 0.43*y7
  dy6 = -280.0*y6*y8 + 0.69*y4 + 1.71*y5 -
           0.43*y6 + 0.69*y7
  dy7 = 280.0*y6*y8 - 1.81*y7
  dy8 = -280.0*y6*y8 + 1.81*y7
end

u0 = zeros(8)
u0[1] = 1
u0[8] = 0.0057
prob = ODEProblem(hires,u0,(0.0,321.8122))

sol = solve(prob,Rodas4(),abstol=1/10^14,reltol=1/10^14)
test_sol = TestSolution(sol)
abstols = 1 ./10 .^(4:11)
reltols = 1 ./10 .^(1:8);

wp = WorkPrecisionSet(prob,abstols,reltols,setups;names=names,
                  save_everystep=false,appxsol=test_sol,maxiters=Int(1e5), seconds=5)
plot(wp, linewidth=1.5, markersize=2, xlims=(1e-12, 1.), dpi=200, title="Hires")
savefig("$(homedir())/hires_$postfix.png")

rober = @ode_def_all Rober begin
  dy₁ = -k₁*y₁+k₃*y₂*y₃
  dy₂ =  k₁*y₁-k₂*y₂^2-k₃*y₂*y₃
  dy₃ =  k₂*y₂^2
end k₁ k₂ k₃
prob = ODEProblem(rober,[1.0;0.0;0.0],(0.0,1e5),(0.04,3e7,1e4))

# compute test solution
test_sol = TestSolution(solve(prob,CVODE_BDF(),abstol=1/10^14,reltol=1/10^14))

abstols = @. 1/10^(6:12)
reltols = @. 1/10^(2:8)

wp = WorkPrecisionSet(prob,abstols,reltols,setups;names=names,
                  save_everystep=false,appxsol=test_sol,maxiters=Int(1e5), seconds=5)
plot(wp, linewidth=1.5, markersize=2, xlims=(1e-12, 1.), dpi=200, title="Rober")
savefig("$(homedir())/rober_$postfix.png")

pollution = @ode_def POL begin
  dy1  = -k1 *y1-k10*y11*y1-k14*y1*y6-k23*y1*y4-k24*y19*y1+
        k2 *y2*y4+k3 *y5*y2+k9 *y11*y2+k11*y13+k12*y10*y2+k22*y19+k25*y20
  dy2  = -k2 *y2*y4-k3 *y5*y2-k9 *y11*y2-k12*y10*y2+k1 *y1+k21*y19
  dy3  = -k15*y3+k1 *y1+k17*y4+k19*y16+k22*y19
  dy4  = -k2 *y2*y4-k16*y4-k17*y4-k23*y1*y4+k15*y3
  dy5  = -k3 *y5*y2+k4 *y7+k4 *y7+k6 *y7*y6+k7 *y9+k13*y14+k20*y17*y6
  dy6  = -k6 *y7*y6-k8 *y9*y6-k14*y1*y6-k20*y17*y6+k3 *y5*y2+k18*y16+k18*y16
  dy7  = -k4 *y7-k5 *y7-k6 *y7*y6+k13*y14
  dy8  = k4 *y7+k5 *y7+k6 *y7*y6+k7 *y9
  dy9  = -k7 *y9-k8 *y9*y6
  dy10 = -k12*y10*y2+k7 *y9+k9 *y11*y2
  dy11 = -k9 *y11*y2-k10*y11*y1+k8 *y9*y6+k11*y13
  dy12 = k9 *y11*y2
  dy13 = -k11*y13+k10*y11*y1
  dy14 = -k13*y14+k12*y10*y2
  dy15 = k14*y1*y6
  dy16 = -k18*y16-k19*y16+k16*y4
  dy17 = -k20*y17*y6
  dy18 = k20*y17*y6
  dy19 = -k21*y19-k22*y19-k24*y19*y1+k23*y1*y4+k25*y20
  dy20 = -k25*y20+k24*y19*y1
end k1  k2  k3  k4  k5  k6  k7  k8  k9  k10  k11  k12  k13  k14  k15  k16  k17  k18  k19  k20  k21  k22  k23  k24  k25
u0 = zeros(20)
p = [.35e0, .266e2, .123e5, .86e-3, .82e-3, .15e5, .13e-3, .24e5, .165e5, .9e4, .22e-1, .12e5, .188e1, .163e5, .48e7, .35e-3, .175e-1, .1e9, .444e12, .124e4, .21e1, .578e1, .474e-1, .178e4, .312e1]
u0[2]  = 0.2
u0[4]  = 0.04
u0[7]  = 0.1
u0[8]  = 0.3
u0[9]  = 0.01
u0[17] = 0.007
prob = ODEProblem(pollution, u0, (0., 60.), p)

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

abstols = @. 1/10^(6:12)
reltols = @. 1/10^(2:8)

wp = WorkPrecisionSet(prob,abstols,reltols,setups;names=names,
                  save_everystep=false,appxsol=test_sol,maxiters=Int(1e5), seconds=5)
plot(wp, linewidth=1.5, markersize=2, xlims=(1e-12, 1.), dpi=200, title="Pollution")
savefig("$(homedir())/pollution_$postfix.png")