WIP: Some benchmarks for the documentation

Project.toml

@@ -9,11 +9,13 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
 MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 SimpleUnPack = "ce78b400-467f-4804-87d8-8f486da07d0a"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5"
 
 [compat]
@@ -22,10 +24,12 @@ LinearAlgebra = "1.7"
 LinearSolve = "2.21"
 MuladdMacro = "0.2.1"
 OrdinaryDiffEq = "6.59"
+RecipesBase = "1.3"
 Reexport = "1"
 SciMLBase = "2"
 SimpleUnPack = "1"
 SparseArrays = "1.7"
 StaticArrays = "1.5"
+Statistics = "1"
 SymbolicIndexingInterface = "0.2, 0.3"
 julia = "1.9"

docs/Project.toml

@@ -1,23 +1,27 @@
 [deps]
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
+DiffEqDevTools = "f3b72e0c-5b89-59e1-b016-84e28bfd966d"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
 [compat]
 BenchmarkTools = "1"
+DiffEqDevTools = "2.44"
 Documenter = "1"
 InteractiveUtils = "1"
 LinearAlgebra = "1.7"
 LinearSolve = "2.21"
 OrdinaryDiffEq = "6.59"
 Pkg = "1"
 Plots = "1"
+PrettyTables = "2.3"
 SparseArrays = "1.7"
 StaticArrays = "1.5"

docs/make.jl

@@ -81,6 +81,12 @@ makedocs(modules = [PositiveIntegrators],
                  "Heat Equation, Neumann BCs" => "heat_equation_neumann.md",
                  "Heat Equation, Dirichlet BCs" => "heat_equation_dirichlet.md",
              ],
+             "Benchmarks" => [
+                 "Experimental order of convergence" => "convergence.md",
+                 "NPZD model" => "npzd_model_benchmark.md",
+                 "Robertson problem" => "robertson_benchmark.md",
+                 "Stratospheric reaction problem" => "stratospheric_reaction_benchmark.md",
+             ],
              "Troubleshooting, FAQ" => "faq.md",
              "API reference" => "api_reference.md",
              "Contributing" => "contributing.md",

docs/src/api_reference.md

@@ -42,4 +42,12 @@ SSPMPRK43
 ```@docs
 isnegative
 isnonnegative
+rel_max_error_tend
+rel_max_error_overall
+rel_l1_error_tend
+rel_l2_error_tend
+work_precision_adaptive
+work_precision_adaptive!
+work_precision_fixed
+work_precision_fixed!
 ```

docs/src/convergence.md

@@ -0,0 +1,180 @@
+# [Experimental convergence order of MPRK schemes](@id convergence_mprk)
+
+In this tutorial we check that the implemented MPRK schemes have the expected order of convergence. 
+
+## Conservative production-destruction systems
+
+First, we consider conservative production-destruction systems (PDS). To investigate the convergence order we define the non-autonomous test problem 
+
+```math
+\begin{aligned}
+u_1' &= \cos(\pi t)^2 u_2 - \sin(2\pi t)^2 u_1,\\
+u_2' & = \sin(2\pi t)^2 u_1 - \cos(\pi t)^2 u_2.
+\end{aligned}
+```
+
+The PDS is conservative, since the sum of the right hand side terms equals zero. 
+An implementation of the problem is given next.
+
+
+```@example eoc
+using PositiveIntegrators
+
+# define problem
+P(u, p, t) = [0.0 cos.(π * t) .^ 2 * u[2]; sin.(2 * π * t) .^ 2 * u[1] 0.0]
+prob = ConservativePDSProblem(P, [0.9; 0.1], (0.0, 1.0))
+
+nothing # hide output
+```
+
+To use `analyticless_test_convergence` from [`DiffEqDevTools`](https://github.com/SciML/DiffEqDevTools.jl) we need to pick a solver to compute the reference solution and specify tolerances. Since the problem is not stiff we use the high order explicit solver `Vern9()` from [`OrdinaryDiffEq`](https://docs.sciml.ai/OrdinaryDiffEq/stable/). Moreover, we need to choose the different time step sizes which are used to investigate the convergence behavior. 
+
+```@example eoc
+using OrdinaryDiffEq
+using DiffEqDevTools # load analyticless_test_convergence
+
+# solver and tolerances to compute reference solution
+test_setup = Dict(:alg => Vern9(), :reltol => 1e-14, :abstol => 1e-14)
+
+# choose step sizes
+dts = 0.5 .^ (5:10)
+
+nothing # hide
+```
+
+### Second order MPRK schemes
+
+First, we test several second order MPRK schemes.
+
+```@example eoc
+# select schemes
+algs2 = [MPRK22(0.5); MPRK22(2.0 / 3.0); MPRK22(1.0); SSPMPRK22(0.5, 1.0)]
+labels2 = ["MPRK22(0.5)"; "MPRK22(2.0/3.0)"; "MPRK22(1.0)"; "SSPMPRK22(0.5, 1.0)"]
+
+#compute errors and experimental order of convergence
+err_eoc = []
+for i in eachindex(algs2)
+     sim = analyticless_test_convergence(dts, prob, algs2[i], test_setup)
+
+     err = sim.errors[:l∞]
+     eoc = [NaN; -log2.(err[2:end] ./ err[1:(end - 1)])]
+
+     push!(err_eoc, tuple.(err, eoc))
+end
+```
+
+Next, we print a table with the computed data. The table lists the errors obtained with the respective time step size ``Δ t`` as well as the estimated order of convergence in parenthesis.
+
+```@example eoc
+using Printf # load @sprintf
+using PrettyTables # load pretty_table
+
+# gather data for table
+data = hcat(dts, reduce(hcat,err_eoc))
+
+# print table
+formatter = (v, i, j) ->  (j>1) ? (@sprintf "%5.2e (%4.2f) " v[1] v[2]) : (@sprintf "%5.2e " v)
+pretty_table(data, formatters = formatter, header = ["Δt"; labels2])                  
+```
+
+The table shows that all schemes converge as expected.
+
+### Third order MPRK schemes
+
+In this section, we proceed as above, but consider third order MPRK schemes instead.
+
+```@example eoc
+# select 3rd order schemes
+algs3 = [MPRK43I(1.0, 0.5); MPRK43I(0.5, 0.75); MPRK43II(0.5); MPRK43II(2.0 / 3.0); 
+         SSPMPRK43()]
+labels3 = ["MPRK43I(1.0,0.5)"; "MPRK43I(0.5, 0.75)"; "MPRK43II(0.5)"; "MPRK43II(2.0/3.0)";
+          "SSPMPRK43()"]
+
+#compute errors and experimental order of convergence
+err_eoc = []
+for i in eachindex(algs3)
+     sim = analyticless_test_convergence(dts, prob, algs3[i], test_setup)
+
+     err = sim.errors[:l∞]
+     eoc = [NaN; -log2.(err[2:end] ./ err[1:(end - 1)])]
+
+     push!(err_eoc, tuple.(err, eoc))
+end
+
+# gather data for table
+data = hcat(dts, reduce(hcat,err_eoc))
+
+# print table
+formatter = (v, i, j) ->  (j>1) ? (@sprintf "%5.2e (%4.2f) " v[1] v[2]) : (@sprintf "%5.2e " v)
+pretty_table(data, formatters = formatter, header = ["Δt"; labels3])  
+```
+
+As above, the table shows that all schemes converge as expected.
+
+## Non-conservative PDS
+
+In this section we consider the non-autonomous but non-conservative test problem 
+
+```math
+\begin{aligned}
+u_1' &= \cos(\pi t)^2 u_2 - \sin(2\pi t)^2 u_1 - \cos(2\pi t)^2 u_1,\\
+u_2' & = \sin(2\pi t)^2 u_1 - \cos(\pi t)^2 u_2 - \sin(\pi t)^2 u_2.
+\end{aligned}
+```
+
+Since the sum of the right hand side terms don't cancel, the PDS is indeed non-conservative. Hence, we need to use `PDSProblem` for its implementation.
+
+```@example eoc
+# choose problem
+P(u, p, t) = [0.0 cos.(π * t) .^ 2 * u[2]; sin.(2 * π * t) .^ 2 * u[1] 0.0]
+D(u, p, t) = [cos.(2 * π * t) .^ 2 * u[1]; sin.(π * t) .^ 2 * u[2]]
+prob = PDSProblem(P, D, [0.9; 0.1], (0.0, 1.0))
+
+nothing # hide
+```
+
+The following sections will show that the selected MPRK schemes show the expected convergence order also for this non-conservative PDS.
+
+### Second order MPRK schemes
+
+```@example eoc
+#compute errors and experimental order of convergence
+err_eoc = []
+for i in eachindex(algs2)
+     sim = analyticless_test_convergence(dts, prob, algs2[i], test_setup)
+
+     err = sim.errors[:l∞]
+     eoc = [NaN; -log2.(err[2:end] ./ err[1:(end - 1)])]
+
+     push!(err_eoc, tuple.(err, eoc))
+end
+
+# gather data for table
+data = hcat(dts, reduce(hcat,err_eoc))
+
+# print table
+formatter = (v, i, j) ->  (j>1) ? (@sprintf "%5.2e (%4.2f) " v[1] v[2]) : (@sprintf "%5.2e " v)
+pretty_table(data, formatters = formatter, header = ["Δt"; labels2])                  
+```
+
+### Third order MPRK schemes
+
+```@example eoc
+#compute errors and experimental order of convergence
+err_eoc = []
+for i in eachindex(algs3)
+     sim = analyticless_test_convergence(dts, prob, algs3[i], test_setup)
+
+     err = sim.errors[:l∞]
+     eoc = [NaN; -log2.(err[2:end] ./ err[1:(end - 1)])]
+
+     push!(err_eoc, tuple.(err, eoc))
+end
+
+# gather data for table
+data = hcat(dts, reduce(hcat,err_eoc))
+
+# print table
+formatter = (v, i, j) ->  (j>1) ? (@sprintf "%5.2e (%4.2f) " v[1] v[2]) : (@sprintf "%5.2e " v)
+pretty_table(data, formatters = formatter, header = ["Δt"; labels3])  
+```