Skip to content

Commit

Permalink
revised npzd benchmark
Browse files Browse the repository at this point in the history
  • Loading branch information
@SKopecz
SKopecz committed Sep 10, 2024
1 parent f6f145b commit cf38f90
Show file tree
Hide file tree
Showing 2 changed files with 119 additions and 142 deletions.
22 changes: 11 additions & 11 deletions docs/make.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -73,19 +73,19 @@ makedocs(modules = [PositiveIntegrators],
# Explicitly specify documentation structure
pages = [
"Home" => "index.md",
#"Tutorials" => [
# "NPZD model" => "npzd_model.md",
# "Robertson problem" => "robertson.md",
# "Stratospheric reaction problem" => "stratospheric_reaction.md",
# "Linear Advection" => "linear_advection.md",
# "Heat Equation, Neumann BCs" => "heat_equation_neumann.md",
# "Heat Equation, Dirichlet BCs" => "heat_equation_dirichlet.md",
#],
"Tutorials" => [
"NPZD model" => "npzd_model.md",
"Robertson problem" => "robertson.md",
"Stratospheric reaction problem" => "stratospheric_reaction.md",
"Linear Advection" => "linear_advection.md",
"Heat Equation, Neumann BCs" => "heat_equation_neumann.md",
"Heat Equation, Dirichlet BCs" => "heat_equation_dirichlet.md",
],
"Benchmarks" => [
#"Experimental order of convergence" => "convergence.md",
"Experimental order of convergence" => "convergence.md",
"NPZD model" => "npzd_model_benchmark.md",
#"Robertson problem" => "robertson_benchmark.md",
#"Stratospheric reaction problem" => "stratospheric_reaction_benchmark.md",
"Robertson problem" => "robertson_benchmark.md",
"Stratospheric reaction problem" => "stratospheric_reaction_benchmark.md",
],
"Troubleshooting, FAQ" => "faq.md",
"API reference" => "api_reference.md",
Expand Down
239 changes: 108 additions & 131 deletions docs/src/npzd_model_benchmark.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,7 +11,7 @@ using Plots
# select problem
prob = prob_pds_npzd
# compute reference solution (standard tolerances are too low)
# compute reference solution
ref_sol = solve(prob, Vern7(); abstol = 1e-14, reltol = 1e-13)
# compute solutions with lose tolerances
Expand All @@ -21,12 +21,12 @@ sol_Ros23 = solve(prob, Rosenbrock23(); abstol, reltol)
sol_MPRK = solve(prob, MPRK22(1.0); abstol, reltol)
# plot solutions
p1 = plot(ref_sol, linestyle = :dash, label = "", legend = :right)
p1 = plot(ref_sol, linestyle = :dash, label = "", legend = :right);
plot!(p1, sol_Ros23; denseplot = false, markers = :circle, ylims = (-1.0, 10.0),
title = "Rosenbrock23", label = ["N" "P" "Z" "D"])
p2 = plot(ref_sol, linestyle = :dash, label = "", legend = :right)
title = "Rosenbrock23", label = ["N" "P" "Z" "D"]);
p2 = plot(ref_sol, linestyle = :dash, label = "", legend = :right);
plot!(p2, sol_MPRK; denseplot = false, markers = true, ylims = (-1.0, 10.0),
title = "MPRK22(1.0)", label = ["N" "P" "Z" "D"])
title = "MPRK22(1.0)", label = ["N" "P" "Z" "D"]);
plot(p1, p2)
```

Expand All @@ -36,7 +36,7 @@ Nevertheless, [OrdinaryDiffEq.jl](https://docs.sciml.ai/OrdinaryDiffEq/stable/)
sol_Ros23 = solve(prob, Rosenbrock23(); abstol, reltol,
isoutofdomain = isnegative) #reject negative solutions
plot(ref_sol, linestyle = :dash, label = "", legend = :right)
plot(ref_sol, linestyle = :dash, label = "", legend = :right);
plot!(sol_Ros23; denseplot = false, markers = :circle, ylims = (-1.0, 10.0),
title = "Rosenbrock23", label = ["N" "P" "Z" "D"])
```
Expand All @@ -52,63 +52,19 @@ sol_ref = sol_ref.u[end]
# define error functions
l2_error(sol, sol_ref) = sqrt(sum(((sol .- sol_ref) ./ sol_ref) .^ 2) / length(sol_ref))
l∞_error(sol, sol_ref) = maximum(abs.((sol .- sol_ref) ./ sol_ref))
nothing #hide output
```

### Fixed time steps sizes
```@example NPZD
# set time step sizes
dts = 1.0 ./ 2.0 .^ (0:1:12)
```

#### L∞
```@example NPZD
# choose methods to compare
algs = [MPRK22(0.5)
MPRK22(2.0 / 3.0)
MPRK22(1.0)
SSPMPRK22(0.5, 1.0)
MPRK43I(1.0, 0.5)
MPRK43I(0.5, 0.75)
MPRK43II(0.5)
MPRK43II(2.0 / 3.0)]

names = ["MPRK22(0.5)"
"MPPRK22(2/3)"
"MPRK22(1.0)"
"SSPMPRK22(0.5,1.0)"
"MPRK43I(1.0,0.5)"
"MPRK43I(0.5,0.75)"
"MPRK43II(0.5)"
"MPRK43II(2.0/3.0)"]
### Adaptive schemes

# compute work-precision
wp_l∞ = workprecision_fixed(prob, algs, names, sol_ref, dts;
compute_error = l∞_error)
First we compare different adaptive MPRK schemes described in the literature.

plot(wp_l∞, names; title = "NPZD benchmark (l∞)", legend = :bottomleft,
color = permutedims([repeat([1], 3)..., 2, repeat([3], 2)..., repeat([4], 2)...]),
xlims = (10^-10, 5*10^-1), xticks = 10.0 .^ (-5:1:0),
ylims = (5*10^-6, 10^-1), yticks = 10.0 .^ (-6:1:0), minorticks = 10
)
```
#### L2
```@example NPZD
# compute work-precision
wp_l2 = workprecision_fixed(prob, algs, names, sol_ref, dts;
compute_error = l2_error)
plot(wp_l2, names; title = "NPZD benchmark (l2)", legend = :topright,
color = permutedims([repeat([1], 3)..., 2, repeat([3], 2)..., repeat([4], 2)...]),
xlims = (10^-10, 5*10^-1), xticks = 10.0 .^ (-10:1:0),
ylims = (5*10^-6, 10^-1), yticks = 10.0 .^ (-5:1:0), minorticks = 10)
```
### Adaptive schemes
First we compare different (adaptive) MPRK schemes described in the literature. The chosen `l∞` error computes the maximum of the absolute values of the difference between the numerical solution and the reference solution over all components and all time steps.
```@example NPZD
# set tolerances
abstols = 1.0 ./ 10.0 .^ (2:1:8)
reltols = abstols ./ 10.0
nothing # hide output
```

#### L∞ error
Expand Down Expand Up @@ -143,30 +99,28 @@ plot(wp_l∞, names; title = "NPZD benchmark (l∞)", legend = :topright,
ylims = (10^-5, 10^0), yticks = 10.0 .^ (-5:1:0), minorticks = 10,)
```

The second- and third-order methods behave very similarly. For comparisons with other schemes from [OrdinaryDiffEq.jl](https://docs.sciml.ai/OrdinaryDiffEq/stable/) we choose the schemes with the smallest error for the initial tolerances, respectively. These are `SSPMPRK22(0.5, 1.0)` and `MPRK43I(1.0, 0.5)`.
The second- and third-order methods behave very similarly. For comparisons with other schemes from [OrdinaryDiffEq.jl](https://docs.sciml.ai/OrdinaryDiffEq/stable/) we choose `MPRK22(1.0)` and `MPRK43I(0.5, 0.75)`.

```@example NPZD
sol_SSPMPRK22 = solve(prob, SSPMPRK22(0.5, 1.0); abstol, reltol)
sol_MPRK43 = solve(prob, MPRK43I(1.0, 0.5); abstol, reltol)
sol_MPRK22 = solve(prob, MPRK22(1.0); abstol, reltol)
sol_MPRK43 = solve(prob, MPRK43I(0.5, 0.75); abstol, reltol)
p1 = plot(ref_sol, linestyle = :dash, label = "", legend = :right)
plot!(p1, sol_SSPMPRK22; denseplot = false, markers = :circle, ylims = (-1.0, 10.0),
title = "SSPMPRK22(0.5, 1.0)", label = ["N" "P" "Z" "D"])
p2 = plot(ref_sol, linestyle = :dash, label = "", legend = :right)
p1 = plot(ref_sol, linestyle = :dash, label = "", legend = :right);
plot!(p1, sol_MPRK22; denseplot = false, markers = :circle, ylims = (-1.0, 10.0),
title = "MPRK22(1.0)", label = ["N" "P" "Z" "D"]);
p2 = plot(ref_sol, linestyle = :dash, label = "", legend = :right);
plot!(p2, sol_MPRK43; denseplot = false, markers = true, ylims = (-1.0, 10.0),
title = "MPRK43I(1.0, 0.5)", label = ["N" "P" "Z" "D"])
title = "MPRK43I(0.5, 0.75)", label = ["N" "P" "Z" "D"]);
plot(p1, p2)
```

Although the SSPMPRK22 solution seems to be more accurate at first glance, the `l∞`-error of the SSPMPRK22 scheme is 0.506359, whereas the `l∞`-error of the MPRK43 scheme is 0.413915. Both errors occurs at approximately $t=2$, where there is a sharp kink in the first component.


Next we compare `SSPMPRK22(0.5, 1.0)` and `MPRK43I(1.0, 0.5)` with some second and third order methods from [OrdinaryDiffEq.jl](https://docs.sciml.ai/OrdinaryDiffEq/stable/). To guarantee positive solutions with these methods, we must select the solver option `isoutofdomain = isnegative`.
Next we compare `MPRK22(1.0)` and `MPRK43I(0.5, 0.75)` with some explicit and implicit methods of second and third order from [OrdinaryDiffEq.jl](https://docs.sciml.ai/OrdinaryDiffEq/stable/). To guarantee positive solutions of these methods, we must select the solver option `isoutofdomain = isnegative`.

```@example NPZD
# select methods
algs1 = [SSPMPRK22(0.5, 1.0)
MPRK43I(1.0, 0.5)]
algs1 = [MPRK22(1.0)
MPRK43I(0.5, 0.75)]
algs2 = [Midpoint()
Heun()
Expand All @@ -180,8 +134,8 @@ algs2 = [Midpoint()
ROS3()
Rosenbrock23()]
names1 = ["SSPMPRK22(0.5,1.0)"
"MPRK43I(1.0,0.5)"]
names1 = ["MPRK22(1.0)"
"MPRK43I(0.5,0.75)"]
names2 = ["Midpoint"
"Heun"
Expand All @@ -202,15 +156,16 @@ workprecision_adaptive!(wp_l∞, prob, algs2, names2, sol_ref, abstols, reltols;
compute_error, isoutofdomain=isnegative)
plot(wp_l∞, [names1; names2]; title = "NPZD benchmark (l∞)", legend = :topright,
color = permutedims([2, 3, repeat([4], 3)..., repeat([5], 4)..., repeat([6], 4)...]),
color = permutedims([1, 3, repeat([4], 3)..., repeat([5], 4)..., repeat([6], 4)...]),
xlims = (10^-8, 10^-1), xticks = 10.0 .^ (-8:1:0),
ylims = (10^-5, 10^0), yticks = 10.0 .^ (-5:1:0), minorticks = 10)
```

Comparison to recommend solvers.

```@example NPZD
algs2 = [Tsit5(),
algs3 = [Tsit5(),
BS3(),
Vern6(),
Vern7(),
Expand All @@ -220,7 +175,7 @@ algs2 = [Tsit5(),
Rodas5P(),
Rodas4P()]
names2 = ["Tsit6"
names3 = ["Tsit6"
"BS3"
"Vern6"
"Vern7"
Expand All @@ -233,17 +188,18 @@ names2 = ["Tsit6"
compute_error = l∞_error
wp_l∞ = workprecision_adaptive(prob, algs1, names1, sol_ref, abstols, reltols;
compute_error)
workprecision_adaptive!(wp_l∞, prob, algs2, names2, sol_ref, abstols, reltols;
workprecision_adaptive!(wp_l∞, prob, algs3, names3, sol_ref, abstols, reltols;
compute_error, isoutofdomain=isnegative)
plot(wp_l∞, [names1; names2]; title = "NPZD benchmark (l∞)", legend = :topright,
color = permutedims([2, 3, repeat([4], 3)..., repeat([5], 4)..., repeat([6], 4)...]),
plot(wp_l∞, [names1; names3]; title = "NPZD benchmark (l∞)", legend = :topright,
color = permutedims([1, 3, repeat([4], 5)...,5, repeat([6], 1)...,repeat([7],2)...]),
xlims = (10^-11, 10^-1), xticks = 10.0 .^ (-11:1:0),
ylims = (10^-5, 10^0), yticks = 10.0 .^ (-5:1:0), minorticks = 10)
```

#### L2 error

```@example NPZD
wp_l2 = workprecision_adaptive(prob, algs, names, sol_ref, abstols, reltols,
compute_error = l2_error)
Expand All @@ -253,84 +209,105 @@ plot(wp_l2, names; title = "NPZD benchmark (l2)", legend = :topright,
xlims = (1 * 10^-8, 10^-1), xticks = 10.0 .^ (-8:1:0),
ylims = (10^-5, 10^0), yticks = 10.0 .^ (-5:1:0), minorticks = 10,)
```
```@example NPZD
# select methods
algs1 = [SSPMPRK22(0.5, 1.0)
MPRK43I(1.0, 0.5)]
algs2 = [Midpoint()
Heun()
Ralston()
TRBDF2()
SDIRK2()
Kvaerno3()
KenCarp3()
Rodas3()
ROS2()
ROS3()
Rosenbrock23()]
names1 = ["SSPMPRK22(0.5,1.0)"
"MPRK43I(1.0,0.5)"]
names2 = ["Midpoint"
"Heun"
"Ralston"
"TRBDF2"
"SDIRK2"
"Kvearno3"
"KenCarp3"
"Rodas3"
"ROS2"
"ROS3"
"Rosenbrock23"]

```@example NPZD
compute_error = l2_error
wp_l2 = workprecision_adaptive(prob, algs1, names1, sol_ref, abstols, reltols;
compute_error)
workprecision_adaptive!(wp_l2, prob, algs2, names2, sol_ref, abstols, reltols;
compute_error, isoutofdomain=isnegative)
plot(wp_l2, [names1; names2]; title = "NPZD benchmark (l2)", legend = :topright,
color = permutedims([2, 3, repeat([4], 3)..., repeat([5], 4)..., repeat([6], 4)...]),
color = permutedims([1, 3, repeat([4], 3)..., repeat([5], 4)..., repeat([6], 4)...]),
xlims = (10^-8, 10^-1), xticks = 10.0 .^ (-8:1:0),
ylims = (10^-5, 10^0), yticks = 10.0 .^ (-5:1:0), minorticks = 10)
```
```@example NPZD
algs2 = [Tsit5(),
BS3(),
Vern6(),
Vern7(),
Vern8(),
TRBDF2(),
Rosenbrock32(),
Rodas5P(),
Rodas4P()]
names2 = ["Tsit6"
"BS3"
"Vern6"
"Vern7"
"Vern8"
"TRBDF2"
"Rosenbrock23"
"Rodas5P"
"Rodas4P"]

```@example NPZD
compute_error = l2_error
wp_l2 = workprecision_adaptive(prob, algs1, names1, sol_ref, abstols, reltols;
compute_error)
workprecision_adaptive!(wp_l2, prob, algs2, names2, sol_ref, abstols, reltols;
workprecision_adaptive!(wp_l2, prob, algs3, names3, sol_ref, abstols, reltols;
compute_error, isoutofdomain=isnegative)
plot(wp_l2, [names1; names2]; title = "NPZD benchmark (l2)", legend = :topright,
color = permutedims([2, 3, repeat([4], 3)..., repeat([5], 4)..., repeat([6], 4)...]),
plot(wp_l2, [names1; names3]; title = "NPZD benchmark (l2)", legend = :topright,
color = permutedims([1, 3, repeat([4], 5)...,5, repeat([6], 1)...,repeat([7],2)...]),
xlims = (10^-11, 10^-1), xticks = 10.0 .^ (-11:1:0),
ylims = (10^-5, 10^0), yticks = 10.0 .^ (-5:1:0), minorticks = 10)
```
### Fixed time steps sizes
```@example NPZD
# set time step sizes
dts = 1.0 ./ 2.0 .^ (0:1:12)
nothing #hide output
```

#### L∞
```@example NPZD
# choose methods to compare
algs = [MPE()
MPRK22(0.5)
MPRK22(2.0 / 3.0)
MPRK22(1.0)
SSPMPRK22(0.5, 1.0)
MPRK43I(1.0, 0.5)
MPRK43I(0.5, 0.75)
MPRK43II(0.5)
MPRK43II(2.0 / 3.0)
SSPMPRK43()]
names = ["MPE()"
"MPRK22(0.5)"
"MPPRK22(2/3)"
"MPRK22(1.0)"
"SSPMPRK22(0.5,1.0)"
"MPRK43I(1.0,0.5)"
"MPRK43I(0.5,0.75)"
"MPRK43II(0.5)"
"MPRK43II(2.0/3.0)"
"SSPMPRK43"]
# compute work-precision
wp_l∞ = workprecision_fixed(prob, algs, names, sol_ref, dts;
compute_error = l∞_error)
plot(wp_l∞, names; title = "NPZD benchmark (l∞)", legend = :bottomleft,
color = permutedims([5,repeat([1], 3)..., 2, repeat([3], 2)..., repeat([4], 2)...,6]),
xlims = (10^-10, 5*10^-1), xticks = 10.0 .^ (-10:1:0),
ylims = (5*10^-6, 10^-1), yticks = 10.0 .^ (-6:1:0), minorticks = 10
)
```

```@example NPZD
algs = [MPRK22(1.0)
MPRK43I(1.0, 0.5)
ROS3()]
names = ["MPRK22(1.0)"
"MPRK43I(1.0,0.5)"]
wp_l∞ = workprecision_fixed(prob, algs, names, sol_ref, dts;
compute_error = l∞_error)
plot(wp_l∞, [names1; names2]; title = "NPZD benchmark (l∞)", legend = :topright,
color = permutedims([1, 3, repeat([4], 3)..., repeat([5], 4)..., repeat([6], 4)...]),
xlims = (10^-8, 10^-1), xticks = 10.0 .^ (-8:1:0),
ylims = (10^-5, 10^0), yticks = 10.0 .^ (-5:1:0), minorticks = 10)
```

#### L2
```@example NPZD
# compute work-precision
wp_l2 = workprecision_fixed(prob, algs, names, sol_ref, dts;
compute_error = l2_error)
plot(wp_l2, names; title = "NPZD benchmark (l2)", legend = :bottomleft,
color = permutedims([5,repeat([1], 3)..., 2, repeat([3], 2)..., repeat([4], 2)...,6]),
xlims = (10^-10, 5*10^-1), xticks = 10.0 .^ (-10:1:0),
ylims = (5*10^-6, 10^-1), yticks = 10.0 .^ (-6:1:0), minorticks = 10)
```

## Literature
- Kopecz, Meister 2nd order
Expand Down

0 comments on commit cf38f90

Please sign in to comment.