add more tests (#61)

* add convergence tests

* set development version to v0.1.6-pre

* test convergence of interpolations

* also test in-place versions

* format

Project.toml

@@ -1,7 +1,7 @@
 name = "PositiveIntegrators"
 uuid = "d1b20bf0-b083-4985-a874-dc5121669aa5"
 authors = ["Stefan Kopecz, Hendrik Ranocha, and contributors"]
-version = "0.1.5"
+version = "0.1.6-pre"
 
 [deps]
 FastBroadcast = "7034ab61-46d4-4ed7-9d0f-46aef9175898"

src/PDSProblemLibrary.jl

@@ -33,6 +33,20 @@ There is one independent linear invariant, e.g. ``u_1+u_2 = 1``.
 prob_pds_linmod = ConservativePDSProblem(P_linmod, u0_linmod, (0.0, 2.0),
                                          analytic = f_linmod_analytic)
 
+function P_linmod!(P, u, p, t)
+    P .= P_linmod(u, p, t)
+    return nothing
+end
+
+"""
+    prob_pds_linmod_inplace
+
+Same as [`prob_pds_linmod`](@ref) but with in-place computation.
+"""
+prob_pds_linmod_inplace = ConservativePDSProblem(P_linmod!, Array(u0_linmod),
+                                                 (0.0, 2.0),
+                                                 analytic = f_linmod_analytic)
+
 # nonlinear model problem
 function P_nonlinmod(u, p, t)
     @SMatrix [0.0 0.0 0.0; u[2] * u[1]/(u[1] + 1.0) 0.0 0.0; 0.0 0.3*u[2] 0.0]
@@ -179,7 +193,7 @@ There is one independent linear invariant, e.g. ``u_1+u_2+u_3 = 3.0``.
 
 ## References
 
-- Enrico Bertolazzi. 
+- Enrico Bertolazzi.
   "Positive and conservative schemes for mass action kinetics."
   Computers and Mathematics with Applications 32 (1996): 29-43.
   [DOI: 10.1016/0898-1221(96)00142-3](https://doi.org/10.1016/0898-1221(96)00142-3)

src/PositiveIntegrators.jl

@@ -37,16 +37,17 @@ import OrdinaryDiffEq: @cache,
                        OrdinaryDiffEqMutableCache, OrdinaryDiffEqConstantCache,
                        alg_cache, initialize!, perform_step!,
                        recursivefill!, _vec, DEFAULT_PRECS, wrapprecs,
-                       dolinsolve
+                       dolinsolve, alg_order
 
 # 2. Export functionality defining the public API
 export PDSFunction, PDSProblem
 export ConservativePDSFunction, ConservativePDSProblem
 
 export MPE, MPRK22
 
-export prob_pds_linmod, prob_pds_nonlinmod, prob_pds_robertson, prob_pds_brusselator,
-       prob_pds_sir, prob_pds_bertolazzi, prob_pds_npzd, prob_pds_stratreac
+export prob_pds_linmod, prob_pds_linmod_inplace, prob_pds_nonlinmod,
+       prob_pds_robertson, prob_pds_brusselator, prob_pds_sir,
+       prob_pds_bertolazzi, prob_pds_npzd, prob_pds_stratreac
 
 # 3. Load source code
 

src/mprk.jl

@@ -221,7 +221,7 @@ function Base.getproperty(alg::MPE, f::Symbol)
     end
 end
 
-OrdinaryDiffEq.alg_order(alg::MPE) = 1
+alg_order(alg::MPE) = 1
 
 struct MPECache{uType, rateType, PType, F, uNoUnitsType} <: OrdinaryDiffEqMutableCache
     u::uType
@@ -407,7 +407,7 @@ function Base.getproperty(alg::MPRK22, f::Symbol)
     end
 end
 
-OrdinaryDiffEq.alg_order(alg::MPRK22) = 2
+alg_order(alg::MPRK22) = 2
 
 struct MPRK22Cache{uType, rateType, PType, tabType, Thread, F, uNoUnitsType} <:
        OrdinaryDiffEqMutableCache

test/Project.toml

@@ -4,6 +4,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]

test/runtests.jl

@@ -1,6 +1,7 @@
 using Test
 using LinearAlgebra
 using SparseArrays
+using Statistics: mean
 
 using OrdinaryDiffEq
 using PositiveIntegrators
@@ -9,6 +10,66 @@ using LinearSolve: RFLUFactorization, LUFactorization
 
 using Aqua: Aqua
 
+"""
+    experimental_order_of_convergence(prob, alg, dts, test_times)
+
+Solve `prob` with `alg` and fixed time steps taken from `dts`, and compute
+the mean error at the times `test_times`.
+Return the associated experimental order of convergence.
+"""
+function experimental_order_of_convergence(prob, alg, dts, test_times)
+    @assert length(dts) > 1
+    errors = zeros(eltype(dts), length(dts))
+    analytic = t -> prob.f.analytic(prob.u0, prob.p, t)
+
+    for (i, dt) in enumerate(dts)
+        sol = solve(prob, alg; dt = dt, adaptive = false)
+        errors[i] = mean(test_times) do t
+            norm(sol(t) - analytic(t))
+        end
+    end
+
+    return experimental_order_of_convergence(errors, dts)
+end
+
+"""
+    experimental_order_of_convergence(prob, alg, dts)
+
+Solve `prob` with `alg` and fixed time steps taken from `dts`, and compute
+the mean error at the final time.
+Return the associated experimental order of convergence.
+"""
+function experimental_order_of_convergence(prob, alg, dts)
+    @assert length(dts) > 1
+    errors = zeros(eltype(dts), length(dts))
+    analytic = t -> prob.f.analytic(prob.u0, prob.p, t)
+
+    for (i, dt) in enumerate(dts)
+        sol = solve(prob, alg; dt = dt, adaptive = false, save_everystep = false)
+        errors[i] = norm(sol.u[end] - analytic(sol.t[end]))
+    end
+
+    return experimental_order_of_convergence(errors, dts)
+end
+
+"""
+    experimental_order_of_convergence(errors, dts)
+
+Compute the experimental order of convergence for given `errors` and
+time step sizes `dts`.
+"""
+function experimental_order_of_convergence(errors, dts)
+    Base.require_one_based_indexing(errors, dts)
+    @assert length(errors) == length(dts)
+    orders = zeros(eltype(errors), length(errors) - 1)
+
+    for i in eachindex(orders)
+        orders[i] = log(errors[i] / errors[i + 1]) / log(dts[i] / dts[i + 1])
+    end
+
+    return mean(orders)
+end
+
 @testset "PositiveIntegrators.jl tests" begin
     @testset "Aqua.jl" begin
         # We do not test ambiguities since we get a lot of
@@ -295,6 +356,23 @@ using Aqua: Aqua
                 @test sol_tridiagonal_ip.u ≈ sol_sparse_ip.u
             end
         end
+
+        @testset "Convergence tests" begin
+            alg = MPE()
+            dts = 0.5 .^ (5:10)
+            problems = (prob_pds_linmod, prob_pds_linmod_inplace)
+            for prob in problems
+                eoc = experimental_order_of_convergence(prob, alg, dts)
+                @test isapprox(eoc, PositiveIntegrators.alg_order(alg); atol = 0.2)
+
+                test_times = [
+                    0.123456789, 1 / pi, exp(-1),
+                    1.23456789, 1 + 1 / pi, 1 + exp(-1),
+                ]
+                eoc = experimental_order_of_convergence(prob, alg, dts, test_times)
+                @test isapprox(eoc, PositiveIntegrators.alg_order(alg); atol = 0.2)
+            end
+        end
     end
 
     @testset "MPRK22" begin
@@ -375,6 +453,23 @@ using Aqua: Aqua
                 @test sol_tridiagonal_ip.u ≈ sol_sparse_ip.u
             end
         end
+
+        @testset "Convergence tests" begin
+            dts = 0.5 .^ (5:10)
+            problems = (prob_pds_linmod, prob_pds_linmod_inplace)
+            for alpha in (0.5, 1.0, 2.0), prob in problems
+                alg = MPRK22(alpha)
+                eoc = experimental_order_of_convergence(prob, alg, dts)
+                @test isapprox(eoc, PositiveIntegrators.alg_order(alg); atol = 0.2)
+
+                test_times = [
+                    0.123456789, 1 / pi, exp(-1),
+                    1.23456789, 1 + 1 / pi, 1 + exp(-1),
+                ]
+                eoc = experimental_order_of_convergence(prob, alg, dts, test_times)
+                @test isapprox(eoc, PositiveIntegrators.alg_order(alg); atol = 0.2)
+            end
+        end
     end
 
     # TODO: Do we want to keep the examples and test them or do we want