Skip to content

Commit 3d42af4

Browse files
ChrisRackauckasChrisRackauckas
committed
fix up e5
1 parent 840c155 commit 3d42af4

File tree

1 file changed

+29
-99
lines changed

1 file changed

+29
-99
lines changed

‎benchmarks/StiffODE/E5.jmd‎

Lines changed: 29 additions & 99 deletions
Original file line numberDiff line numberDiff line change
@@ -6,6 +6,7 @@ author: Chris Rackauckas
66
```julia
77
using OrdinaryDiffEq, DiffEqDevTools, Sundials, Plots, ODEInterfaceDiffEq, LSODA
88
using LinearAlgebra, StaticArrays, RecursiveFactorization
9+
using OrdinaryDiffEqFIRK
910
gr()
1011

1112
# E5 Problem from Stiff Test Set
@@ -61,127 +62,58 @@ end
6162
u0 = [1.76e-3, 0.0, 0.0, 0.0]
6263

6364
# Time span - long enough to see both fast and slow time scales
64-
tspan = (0.0, 1e6)
65+
tspan = (0.0, 1e10)
6566

6667
prob = ODEProblem{true, SciMLBase.FullSpecialize}(e5!, u0, tspan)
6768
probstatic = ODEProblem{false}(e5, SVector{4}(u0), tspan)
69+
probbig = ODEProblem{true, SciMLBase.FullSpecialize}(e5!, big.(u0), big.(tspan))
6870

6971
# Generate reference solution
70-
sol = solve(prob, CVODE_BDF(), abstol=1/10^14, reltol=1/10^14)
71-
sol2 = solve(probstatic, Rodas5P(), abstol=1/10^14, reltol=1/10^14)
72+
sol = solve(probbig, Rodas5P(), abstol=1/10^60, reltol=1/10^30)
7273
probs = [prob, probstatic]
73-
test_sol = [sol, sol2]
74+
test_sol = [sol, sol]
7475

75-
abstols = 1.0 ./ 10.0 .^ (4:11)
76-
reltols = 1.0 ./ 10.0 .^ (1:8)
76+
abstols = [1.7e-28 for i in 1:2]
77+
reltols = 10.0 .^ -(6 .+ (1:2)./4)
7778
```
7879

7980
```julia
80-
plot(sol, labels=["A" "B" "C" "D"])
81+
ylabels = 10.0 .^ (-30:2:0)
82+
plot(sol, xscale=:log10, yscale=:log10, tspan = (1e-5, 1e10), labels=["A" "B" "C" "D"], yticks = ylabels)
8183
```
8284

83-
```julia
84-
plot(sol, xscale=:log10, labels=["A" "B" "C" "D"])
85-
```
86-
87-
## High Tolerances
88-
89-
This is the speed when you just want the answer.
85+
## Benchmarks
9086

9187
```julia
92-
abstols = 1.0 ./ 10.0 .^ (5:8)
93-
reltols = 1.0 ./ 10.0 .^ (1:4)
9488
setups = [Dict(:alg=>Rosenbrock23()),
9589
Dict(:alg=>Rosenbrock23(), :prob_choice => 2),
96-
Dict(:alg=>FBDF()),
97-
Dict(:alg=>QNDF()),
98-
Dict(:alg=>TRBDF2()),
99-
Dict(:alg=>CVODE_BDF()),
100-
Dict(:alg=>rodas()),
101-
Dict(:alg=>radau()),
102-
Dict(:alg=>RadauIIA5()),
103-
Dict(:alg=>ROS34PW1a()),
10490
Dict(:alg=>lsoda()),
105-
]
106-
wp = WorkPrecisionSet(probs,abstols,reltols,setups; verbose=false, dense=false,
107-
save_everystep=false,appxsol=test_sol,maxiters=Int(1e5),numruns=10)
108-
plot(wp)
109-
```
110-
111-
```julia
112-
setups = [Dict(:alg=>Rosenbrock23()),
113-
Dict(:alg=>Rosenbrock23(), :prob_choice => 2),
114-
Dict(:alg=>Kvaerno3()),
115-
Dict(:alg=>KenCarp4()),
116-
Dict(:alg=>TRBDF2()),
117-
Dict(:alg=>KenCarp3()),
118-
Dict(:alg=>lsoda()),
119-
Dict(:alg=>radau())]
120-
names = ["Rosenbrock23" "Rosenbrock23 Static" "Kvaerno3" "KenCarp4" "TRBDF2" "KenCarp3" "lsoda" "radau"]
121-
wp = WorkPrecisionSet(probs,abstols,reltols,setups;names=names, verbose=false, dense=false,
122-
save_everystep=false,appxsol=test_sol,maxiters=Int(1e5),numruns=10)
123-
plot(wp)
124-
```
125-
126-
## Medium Tolerances
127-
128-
```julia
129-
abstols = 1.0 ./ 10.0 .^ (6:9)
130-
reltols = 1.0 ./ 10.0 .^ (3:6)
131-
132-
setups = [Dict(:alg=>Rosenbrock23()),
133-
Dict(:alg=>Rosenbrock23(), :prob_choice => 2),
134-
Dict(:alg=>KenCarp5()),
135-
Dict(:alg=>KenCarp4()),
136-
Dict(:alg=>KenCarp4(), :prob_choice => 2),
137-
Dict(:alg=>KenCarp3()),
138-
Dict(:alg=>TRBDF2()),
139-
Dict(:alg=>CVODE_BDF()),
140-
Dict(:alg=>rodas()),
141-
Dict(:alg=>radau()),
142-
]
143-
wp = WorkPrecisionSet(probs,abstols,reltols,setups; verbose=false, dense=false,
144-
save_everystep=false,appxsol=test_sol,maxiters=Int(1e5),numruns=10)
145-
plot(wp)
146-
```
147-
148-
### Timeseries Errors
149-
150-
```julia
151-
abstols = 1.0 ./ 10.0 .^ (5:8)
152-
reltols = 1.0 ./ 10.0 .^ (1:4)
153-
setups = [Dict(:alg=>Rosenbrock23()),
154-
Dict(:alg=>Rosenbrock23(), :prob_choice => 2),
91+
Dict(:alg=>Rodas5P()),
92+
Dict(:alg=>Rodas5P(), :prob_choice => 2),
15593
Dict(:alg=>FBDF()),
15694
Dict(:alg=>QNDF()),
15795
Dict(:alg=>TRBDF2()),
158-
Dict(:alg=>rodas()),
159-
Dict(:alg=>lsoda()),
160-
Dict(:alg=>radau()),
96+
Dict(:alg=>CVODE_BDF()),
97+
#Dict(:alg=>rodas()),
98+
#Dict(:alg=>radau()),
16199
Dict(:alg=>RadauIIA5()),
162100
Dict(:alg=>ROS34PW1a()),
101+
163102
]
164-
wp = WorkPrecisionSet(probs,abstols,reltols,setups;error_estimate=:l2, verbose=false, dense=false,
103+
wp = WorkPrecisionSet(probs,abstols,reltols,setups; error_estimate = :l2, verbose=false, dense=false,
165104
save_everystep=false,appxsol=test_sol,maxiters=Int(1e5),numruns=10)
166105
plot(wp)
167106
```
168107

169-
## Low Tolerances
170-
171-
This is the speed at lower tolerances, measuring what's good when accuracy is needed.
172-
173108
```julia
174-
abstols = 1.0 ./ 10.0 .^ (7:12)
175-
reltols = 1.0 ./ 10.0 .^ (4:9)
176-
177109
setups = [Dict(:alg=>FBDF()),
178110
Dict(:alg=>QNDF()),
179111
Dict(:alg=>CVODE_BDF()),
180-
Dict(:alg=>ddebdf()),
181-
Dict(:alg=>Rodas4()),
182-
Dict(:alg=>Rodas4(), :prob_choice => 2),
183-
Dict(:alg=>rodas()),
184112
Dict(:alg=>lsoda()),
113+
Dict(:alg=>ddebdf()),
114+
Dict(:alg=>Rodas5P()),
115+
Dict(:alg=>Rodas5P(), :prob_choice => 2),
116+
Dict(:alg=>rodas()),
185117
Dict(:alg=>radau()),
186118
Dict(:alg=>RadauIIA5()),
187119
]
@@ -198,23 +130,21 @@ setups = [Dict(:alg=>Kvaerno4()),
198130
Dict(:alg=>KenCarp47()),
199131
Dict(:alg=>KenCarp47(), :prob_choice => 2),
200132
Dict(:alg=>KenCarp5()),
201-
Dict(:alg=>Rodas4()),
202-
Dict(:alg=>Rodas4(), :prob_choice => 2),
203133
Dict(:alg=>lsoda()),
204-
Dict(:alg=>radau()),
205-
Dict(:alg=>ImplicitEulerExtrapolation()),
206-
Dict(:alg=>ImplicitEulerBarycentricExtrapolation()),
207-
Dict(:alg=>ImplicitHairerWannerExtrapolation()),
134+
Dict(:alg=>Rodas5P()),
135+
Dict(:alg=>Rodas5P(), :prob_choice => 2),
136+
Dict(:alg=>ImplicitEulerExtrapolation(min_order = 4, init_order = 11,threading = OrdinaryDiffEq.PolyesterThreads())),
137+
Dict(:alg=>ImplicitEulerExtrapolation(min_order = 4, init_order = 11,threading = false)),
138+
Dict(:alg=>ImplicitEulerBarycentricExtrapolation(min_order = 4, init_order = 11, threading = OrdinaryDiffEq.PolyesterThreads())),
139+
Dict(:alg=>ImplicitEulerBarycentricExtrapolation(min_order = 4, init_order = 11, threading = false)),
140+
Dict(:alg=>ImplicitHairerWannerExtrapolation(min_order = 3, init_order = 111,threading = OrdinaryDiffEq.PolyesterThreads())),
141+
Dict(:alg=>ImplicitHairerWannerExtrapolation(min_order = 3, init_order = 11,threading = false)),
208142
]
209143
wp = WorkPrecisionSet(probs,abstols,reltols,setups; verbose=false, dense=false,
210144
save_everystep=false,appxsol=test_sol,maxiters=Int(1e5),numruns=10)
211145
plot(wp)
212146
```
213147

214-
### Conclusion
215-
216-
The E5 problem demonstrates a classic stiff chemical pyrolysis system with widely separated time scales. The slow formation of intermediate B is followed by rapid conversion to C and D. At high tolerances, `Rosenbrock23` performs well, while `Rodas4` and `KenCarp4` excel at medium tolerances. For high precision requirements, `radau` and the Rodas family show excellent efficiency.
217-
218148
```julia, echo = false
219149
using SciMLBenchmarks
220150
SciMLBenchmarks.bench_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])

0 commit comments

Comments
 (0)