Use MPDeC in benchmarks (#160)

* no runtime-dispatch for MPDeC + StaticArrays and first benchmarks npzd model

* order reduction in convergence benchmark

* MPDeC in NPZD benchmark

* MPDeC in NPZD benchmark

* bugfix

* mpded in npzd and robertson benchmark

* increased MPDeC default small_constant to run stratospheric reaction benchmark with constant time steps

* bugfix

* mpdec in stratospheric reaction benchmark

* minor changes

* Apply suggestions from code review

Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com>

* compat entry for DoubleFloats

* removed function sqmtr

* format

---------

Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com>

Author: SKopecz

docs/Project.toml

@@ -4,6 +4,7 @@ DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
 DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def"
 DiffEqDevTools = "f3b72e0c-5b89-59e1-b016-84e28bfd966d"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+DoubleFloats = "497a8b3b-efae-58df-a0af-a86822472b78"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
@@ -25,6 +26,7 @@ DiffEqBase = "6.160"
 DiffEqCallbacks = "4"
 DiffEqDevTools = "2.45.1"
 Documenter = "1"
+DoubleFloats = "1.4"
 InteractiveUtils = "1"
 LinearAlgebra = "1"
 LinearSolve = "3.7.1"

docs/src/convergence.md

@@ -1,6 +1,7 @@
 # [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. 
+In this tutorial, we check that all implemented MPRK schemes can achieve their expected order of convergence.
+We also address the issue that some methods suffer from order reduction when the solution gets too close to zero.
 
 ## Conservative production-destruction systems
 
@@ -14,7 +15,7 @@ u_2' & = \sin(2\pi t)^2 u_1 - \cos(\pi t)^2 u_2, & u_2(0)&=0.1,
 ```
 for ``0≤ t≤ 1``.
 The PDS is conservative since the sum of the right-hand side terms equals zero. 
-An implementation of the problem is given next.
+An implementation of this problem is given next.
 
 
 ```@example eoc
@@ -29,94 +30,102 @@ nothing # hide
 
 To use `analyticless_test_convergence` from [DiffEqDevTools.jl](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 [OrdinaryDiffEqVerner.jl](https://docs.sciml.ai/OrdinaryDiffEq/stable/).
-Moreover, we choose time step sizes to investigate the convergence behavior. 
-
 ```@example eoc
 using OrdinaryDiffEqVerner
 using DiffEqDevTools: 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.
+To keep the code short, we also define an auxiliary function that outputs a convergence table, which lists the errors obtained with the respective time step size ``Δ t`` as well as the estimated order of convergence in parentheses.
 
 ```@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)"]
+using Printf: @sprintf
+using PrettyTables: pretty_table
 
-# compute errors and experimental order of convergence
-err_eoc = []
-for i in eachindex(algs2)
-     sim = analyticless_test_convergence(dts, prob, algs2[i], test_setup)
+# auxiliary function
+function convergence_table(dts, prob, algs, labels, test_setup)
+    # compute errors and estimated convergence orders
+    err_eoc = []
+    for i in eachindex(algs)
+        sim = analyticless_test_convergence(dts, prob, algs[i], test_setup)
+
+        err = sim.errors[:l∞]
+        eoc = [NaN; -log2.(err[2:end] ./ err[1:(end - 1)])]
+
+        push!(err_eoc, tuple.(err, eoc))
+    end
 
-     err = sim.errors[:l∞]
-     eoc = [NaN; -log2.(err[2:end] ./ err[1:(end - 1)])]
+    # gather data for table
+    data = hcat(dts, reduce(hcat,err_eoc))
 
-     push!(err_eoc, tuple.(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"; labels]) 
 end
+
+nothing # hide
 ```
 
-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 parentheses.
+### Second-order and third-order MPRK schemes
+
+First, we test several second-order and third-order MPRK schemes.
 
 ```@example eoc
-using Printf: @sprintf
-using PrettyTables: pretty_table
+# choose step sizes
+dts = 0.5 .^ (5:10)
+
+# select 2nd order schemes
+algs2 = [MPRK22(0.5); MPRK22(2.0 / 3.0); MPRK22(1.0); SSPMPRK22(0.5, 1.0); MPDeC(2)]
+labels2 = ["MPRK22(0.5)"; "MPRK22(2.0/3.0)"; "MPRK22(1.0)"; "SSPMPRK22(0.5, 1.0)"; "MPDeC(2)"]
+
+# select 3rd order schemes
+algs3 = [MPRK43I(1.0, 0.5); MPRK43I(0.5, 0.75); MPRK43II(0.5); MPRK43II(2.0 / 3.0); 
+         SSPMPRK43(); MPDeC(3)]
+labels3 = ["MPRK43I(1.0,0.5)"; "MPRK43I(0.5, 0.75)"; "MPRK43II(0.5)"; "MPRK43II(2.0/3.0)";
+          "SSPMPRK43()"; "MPDeC(3)"]
 
-# gather data for table
-data = hcat(dts, reduce(hcat,err_eoc))
+convergence_table(dts, prob, algs2, labels2, test_setup)
 
-# 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])                  
+convergence_table(dts, prob, algs3, labels3, test_setup)
 ```
 
 The table shows that all schemes converge as expected.
 
-### Third-order MPRK schemes
+### Higher-order MPRK schemes
 
-In this section, we proceed as above, but consider third-order MPRK schemes instead.
+To actually see the order of higher-order methods we need to use more accurate floating-point numbers. Here, we use [`DoubleFloats`](https://github.com/JuliaMath/DoubleFloats.jl).
 
 ```@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()"]
+using DoubleFloats 
 
-# compute errors and experimental order of convergence
-err_eoc = []
-for i in eachindex(algs3)
-     sim = analyticless_test_convergence(dts, prob, algs3[i], test_setup)
+# define problem using Double64
+P(u, p, t) = [0 cospi(t)^2 * u[2]; sinpi(2 * t)^2 * u[1] 0]
+u0 = [Double64(9) / 10; Double64(1) / 10]
+tspan = (Double64(0), Double64(1))
+prob_d64 = ConservativePDSProblem(P, u0, tspan)
 
-     err = sim.errors[:l∞]
-     eoc = [NaN; -log2.(err[2:end] ./ err[1:(end - 1)])]
+# choose step sizes
+dts_d64 = Double64(1/2) .^ (5:9)
 
-     push!(err_eoc, tuple.(err, eoc))
-end
+# select higher-order schemes
+algs4 = [MPDeC(4); MPDeC(5); MPDeC(6); MPDeC(7); MPDeC(8); MPDeC(9); MPDeC(10)]
+labels4 = ["MPDeC(4)"; "MPDeC(5)"; "MPDeC(6)"; "MPDeC(7)"; "MPDeC(8)"; "MPDeC(9)"; "MPDeC(10)"]
 
-# gather data for table
-data = hcat(dts, reduce(hcat,err_eoc))
+# solver and tolerances to compute reference solution
+test_setup_d64 = Dict(:alg => Vern9(), :reltol => 1e-30, :abstol => 1e-30)
 
-# 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])  
+# compute errors and experimental order of convergence
+convergence_table(dts_d64, prob_d64, algs4, labels4, test_setup_d64)
 ```
 
-As above, the table shows that all schemes converge as expected.
+Again, all schemes show the expected converge order.
 
 ## Non-conservative PDS
 
-In this section we consider the non-autonomous but non-conservative test problem 
+Next, we consider the non-autonomous but also non-conservative test problem 
 
 ```math
 \begin{aligned}
@@ -126,60 +135,77 @@ u_2' & = \sin(2\pi t)^2 u_1 - \cos(\pi t)^2 u_2 - \sin(\pi t)^2 u_2, & u_2(0)&=0
 ```
 
 for ``0≤ t≤ 1``.
-Since the sum of the right-hand side terms does not cancel, the PDS is indeed non-conservative.
+Since the sum of the right-hand side terms does not vanish, the PDS is indeed non-conservative.
 Hence, we need to use [`PDSProblem`](@ref) 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]]
+# PDS
+P(u, p, t) = [0.0 cospi(t)^2 * u[2]; sinpi(2 * t)^2 * u[1] 0.0]
+D(u, p, t) = [cospi(2 * t)^2 * u[1]; sinpi(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.
+The following tables demonstrate that the chosen MPRK schemes converge as expected also for this non-conservative PDS.
 
-### Second-order MPRK schemes
+### Second-order and third-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)
+convergence_table(dts, prob, algs2, labels2, test_setup)    
 
-     err = sim.errors[:l∞]
-     eoc = [NaN; -log2.(err[2:end] ./ err[1:(end - 1)])]
+convergence_table(dts, prob, algs3, labels3, test_setup)
+```
 
-     push!(err_eoc, tuple.(err, eoc))
-end
+### Higher-order MPRK schemes
 
-# gather data for table
-data = hcat(dts, reduce(hcat,err_eoc))
+```@example eoc
+# problem implementation using DoubleFloats
+P(u, p, t) = [0 cospi(t)^2 * u[2]; sinpi(2 * t)^2 * u[1] 0]
+D(u, p, t) = [cospi(2 * t)^2 * u[1]; sinpi(t)^2 * u[2]]
+prob_d64 = PDSProblem(P, D, [Double64(9)/10; Double64(1)/10], (Double64(0), Double64(1)))
 
-# 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])                  
+convergence_table(dts_d64, prob_d64, algs4, labels4, test_setup_d64)
 ```
 
-### Third-order MPRK schemes
+## Order reduction
+
+It was shown in [Torlo, Öffner, Ranocha: Issues with positivity-preserving Patankar-type schemes with positivity-preserving Patankar-type schemes](https://doi.org/10.1016/j.apnum.2022.07.014) that some MPRK methods 
+suffer from order reduction if the solution of the PDS is too close to zero.
+We demonstrate this by solving a problem where one component of the initial condition is equal to zero. 
+
+The problem is
+
+```math
+\begin{aligned}
+u_1' &= -u_1, & u_1(0)&=1, \\
+u_2' & = u_1, & u_2(0)&=0,
+\end{aligned}
+```
+
+for ``0≤ t≤ 1`` and can be implemented as follows.
+
 
 ```@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)
+# PDS
+P(u, p, t) = [0 0; u[1] 0]
+prob = ConservativePDSProblem(P, [1.0; 0.0], (0.0, 1.0))
+nothing # hide
+```
 
-     err = sim.errors[:l∞]
-     eoc = [NaN; -log2.(err[2:end] ./ err[1:(end - 1)])]
+Next, we generate the corresponding convergence tables as in the sections above.
 
-     push!(err_eoc, tuple.(err, eoc))
-end
+```@example eoc
+test_setup = Dict(:alg => Vern9(), :reltol => 1e-14, :abstol => 1e-14)
+
+dts = 0.5 .^ (6:12)
 
-# gather data for table
-data = hcat(dts, reduce(hcat,err_eoc))
+convergence_table(dts, prob, algs2, labels2, test_setup) 
+convergence_table(dts, prob, algs3, labels3, test_setup) 
+convergence_table(dts, prob, algs4, labels4, test_setup) 
+
+nothing # hide
+```
 
-# 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])  
-```
+We find that all methods apart from MPDeC(``K``) methods with ``K ≥ 3`` converge as expected.
+The MPDeC(``K``) methods with ``K ≥ 3`` suffer from order reduction and show convergence order 2 instead of ``K``.