From 6b3e3fbfd6346a9491358e1b373aa4539e3be086 Mon Sep 17 00:00:00 2001
From: Stefan Kopecz <kopecz@mathematik.uni-kassel.de>
Date: Thu, 18 Jul 2024 21:38:32 +0200
Subject: [PATCH 01/42] show u and sigma in MPRK43 during perform_step

---
 src/mprk.jl      |  7 +++++
 test/runtests.jl | 70 ++++++++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 77 insertions(+)

diff --git a/src/mprk.jl b/src/mprk.jl
index 6e54e4b2..92836e7d 100644
--- a/src/mprk.jl
+++ b/src/mprk.jl
@@ -1443,6 +1443,7 @@ end
 
     # avoid division by zero due to zero Patankar weights
     @.. broadcast=false σ=uprev + small_constant
+    @show σ
 
     build_mprk_matrix!(P3, P3, σ, dt)
 
@@ -1451,6 +1452,7 @@ end
     linres = solve!(linsolve)
 
     u .= linres
+    @show u
     if !(q1 ≈ q2)
         tmp2 .= u #u2 in out-of-place version
     end
@@ -1458,6 +1460,7 @@ end
 
     @.. broadcast=false σ=σ^(1 - q1) * u^q1
     @.. broadcast=false σ=σ + small_constant
+    @show σ
 
     f.p(P2, u, p, t + c2 * dt) # evaluate production terms
     if issparse(P)
@@ -1480,11 +1483,13 @@ end
     linres = solve!(linsolve)
 
     u .= linres
+    @show u
     integrator.stats.nsolve += 1
 
     if !(q1 ≈ q2)
         @.. broadcast=false σ=(uprev + small_constant)^(1 - q2) * tmp2^q2
         @.. broadcast=false σ=σ + small_constant
+        @show σ
     end
 
     if issparse(P)
@@ -1507,6 +1512,7 @@ end
     integrator.stats.nsolve += 1
 
     σ .= linres
+    @show σ
     # avoid division by zero due to zero Patankar weights
     @.. broadcast=false σ=σ + small_constant
 
@@ -1531,6 +1537,7 @@ end
     linres = solve!(linsolve)
 
     u .= linres
+    @show u
     integrator.stats.nsolve += 1
 
     # Now tmp stores the error estimate
diff --git a/test/runtests.jl b/test/runtests.jl
index 0148c6d6..031213cb 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -202,6 +202,75 @@ function linear_advection_fd_upwind_D!(D, u, p, t)
     return nothing
 end
 
+
+@testset "Different matrix types (conservative, adaptive)" begin
+    prod_1! = (P, u, p, t) -> begin
+        fill!(P, zero(eltype(P)))
+        for i in 1:(length(u) - 1)
+            P[i, i + 1] = i * u[i]
+        end
+        return nothing
+    end
+
+    prod_2! = (P, u, p, t) -> begin
+        fill!(P, zero(eltype(P)))
+        for i in 1:(length(u) - 1)
+            P[i + 1, i] = i * u[i + 1]
+        end
+        return nothing
+    end
+
+    prod_3! = (P, u, p, t) -> begin
+        fill!(P, zero(eltype(P)))
+        for i in 1:(length(u) - 1)
+            P[i, i + 1] = i * u[i]
+            P[i + 1, i] = i * u[i + 1]
+        end
+        return nothing
+    end
+
+    n = 4
+    P_tridiagonal = Tridiagonal([0.1, 0.2, 0.3],
+                                zeros(n),
+                                [0.4, 0.5, 0.6])
+    P_dense = Matrix(P_tridiagonal)
+    P_sparse = sparse(P_tridiagonal)
+    u0 = [1.0, 1.5, 2.0, 2.5]
+    tspan = (0.0, 1.0)
+    dt = 0.25
+
+    rtol = sqrt(eps(Float32))
+
+    @testset "$alg" for alg in (MPRK43II(2.0 / 3.0), MPRK43II(0.5))
+
+        #=
+        # Something is going wrong on macOS with sparse matrices for MPRK43II, see
+        # https://github.com/SKopecz/PositiveIntegrators.jl/pull/101
+        if Sys.isapple() && (alg isa MPRK43II)
+            @test_skip false
+            continue
+        end
+        =#
+
+        for prod! in (prod_1!, prod_2!, prod_3!)
+            prod = (u, p, t) -> begin
+                P = similar(u, (length(u), length(u)))
+                prod!(P, u, p, t)
+                return P
+            end
+            prob_sparse_ip = ConservativePDSProblem(prod!, u0, tspan;
+                                                    p_prototype = P_sparse)
+            prob_sparse_op = ConservativePDSProblem(prod, u0, tspan;
+                                                    p_prototype = P_sparse)
+            sol_sparse_ip = solve(prob_sparse_ip, alg; dt)
+            sol_sparse_op = solve(prob_sparse_op, alg; dt)
+
+            @test isapprox(sol_sparse_ip.t, sol_sparse_op.t; rtol)
+            @test isapprox(sol_sparse_ip.u, sol_sparse_op.u; rtol)
+        end
+    end
+end
+#=
 @testset "PositiveIntegrators.jl tests" begin
     @testset "Aqua.jl" begin
         # We do not test ambiguities since we get a lot of
@@ -1735,3 +1804,4 @@ end
         @test_nowarn plot(sol)
     end
 end;
+=#
\ No newline at end of file

From c26a4471890a42534f5caa781eb650419d7be7ce Mon Sep 17 00:00:00 2001
From: Stefan Kopecz <kopecz@mathematik.uni-kassel.de>
Date: Thu, 18 Jul 2024 21:47:38 +0200
Subject: [PATCH 02/42] more output

---
 src/mprk.jl | 14 +++++++-------
 1 file changed, 7 insertions(+), 7 deletions(-)

diff --git a/src/mprk.jl b/src/mprk.jl
index 92836e7d..754db078 100644
--- a/src/mprk.jl
+++ b/src/mprk.jl
@@ -1443,7 +1443,7 @@ end
 
     # avoid division by zero due to zero Patankar weights
     @.. broadcast=false σ=uprev + small_constant
-    @show σ
+    @show ("1",σ)
 
     build_mprk_matrix!(P3, P3, σ, dt)
 
@@ -1452,7 +1452,7 @@ end
     linres = solve!(linsolve)
 
     u .= linres
-    @show u
+    @show ("2",u)
     if !(q1 ≈ q2)
         tmp2 .= u #u2 in out-of-place version
     end
@@ -1460,7 +1460,7 @@ end
 
     @.. broadcast=false σ=σ^(1 - q1) * u^q1
     @.. broadcast=false σ=σ + small_constant
-    @show σ
+    @show ("3",σ)
 
     f.p(P2, u, p, t + c2 * dt) # evaluate production terms
     if issparse(P)
@@ -1483,13 +1483,13 @@ end
     linres = solve!(linsolve)
 
     u .= linres
-    @show u
+    @show ("4",u)
     integrator.stats.nsolve += 1
 
     if !(q1 ≈ q2)
         @.. broadcast=false σ=(uprev + small_constant)^(1 - q2) * tmp2^q2
         @.. broadcast=false σ=σ + small_constant
-        @show σ
+        @show ("5",σ)
     end
 
     if issparse(P)
@@ -1512,7 +1512,7 @@ end
     integrator.stats.nsolve += 1
 
     σ .= linres
-    @show σ
+    @show ("6",σ)
     # avoid division by zero due to zero Patankar weights
     @.. broadcast=false σ=σ + small_constant
 
@@ -1537,7 +1537,7 @@ end
     linres = solve!(linsolve)
 
     u .= linres
-    @show u
+    @show ("7",u)
     integrator.stats.nsolve += 1
 
     # Now tmp stores the error estimate

From c08dccb66b2b95fd242cc28dc583eb48fbc64d1e Mon Sep 17 00:00:00 2001
From: Stefan Kopecz <kopecz@mathematik.uni-kassel.de>
Date: Thu, 18 Jul 2024 21:53:54 +0200
Subject: [PATCH 03/42] more output

---
 src/mprk.jl | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/src/mprk.jl b/src/mprk.jl
index 754db078..7685acd0 100644
--- a/src/mprk.jl
+++ b/src/mprk.jl
@@ -1508,7 +1508,10 @@ end
 
     # Same as linres = P3 \ tmp
     linsolve.A = P3
+    @show linsolve.A
+    @show linsolve.b
     linres = solve!(linsolve)
+    @show linres
     integrator.stats.nsolve += 1
 
     σ .= linres

From 5ee3c3608334e317ff6b3b6b7974b9f483226e1d Mon Sep 17 00:00:00 2001
From: Stefan Kopecz <kopecz@mathematik.uni-kassel.de>
Date: Fri, 19 Jul 2024 08:35:09 +0200
Subject: [PATCH 04/42] more output

---
 src/mprk.jl | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/src/mprk.jl b/src/mprk.jl
index 7685acd0..e6df3315 100644
--- a/src/mprk.jl
+++ b/src/mprk.jl
@@ -1510,8 +1510,11 @@ end
     linsolve.A = P3
     @show linsolve.A
     @show linsolve.b
+    A_tmp = Matrix(linsolve.A)
+    b_tmp = linsolve.b
     linres = solve!(linsolve)
     @show linres
+    @show A_tmp\b_tmp
     integrator.stats.nsolve += 1
 
     σ .= linres

From eb15c45df420834c72c24c3b335dae1b94964b43 Mon Sep 17 00:00:00 2001
From: Stefan Kopecz <kopecz@mathematik.uni-kassel.de>
Date: Fri, 19 Jul 2024 16:18:54 +0200
Subject: [PATCH 05/42] changed MPRK43II small_constant

---
 src/mprk.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/mprk.jl b/src/mprk.jl
index e6df3315..e70e13aa 100644
--- a/src/mprk.jl
+++ b/src/mprk.jl
@@ -998,7 +998,7 @@ struct MPRK43II{T, F, T2} <: OrdinaryDiffEqAdaptiveAlgorithm
     small_constant_function::T2
 end
 
-function MPRK43II(gamma; linsolve = LUFactorization(), small_constant = nothing)
+function MPRK43II(gamma; linsolve = LUFactorization(), small_constant = 1e-50)
     if isnothing(small_constant)
         small_constant_function = floatmin
     elseif small_constant isa Number

From aa4b1656b69a88097dc9af1a47caa66a94078ea3 Mon Sep 17 00:00:00 2001
From: Stefan Kopecz <kopecz@mathematik.uni-kassel.de>
Date: Fri, 19 Jul 2024 17:18:05 +0200
Subject: [PATCH 06/42] removed addtional output

---
 src/mprk.jl      | 13 --------
 test/runtests.jl | 79 +-----------------------------------------------
 2 files changed, 1 insertion(+), 91 deletions(-)

diff --git a/src/mprk.jl b/src/mprk.jl
index e70e13aa..fbb0cf46 100644
--- a/src/mprk.jl
+++ b/src/mprk.jl
@@ -1443,7 +1443,6 @@ end
 
     # avoid division by zero due to zero Patankar weights
     @.. broadcast=false σ=uprev + small_constant
-    @show ("1",σ)
 
     build_mprk_matrix!(P3, P3, σ, dt)
 
@@ -1452,7 +1451,6 @@ end
     linres = solve!(linsolve)
 
     u .= linres
-    @show ("2",u)
     if !(q1 ≈ q2)
         tmp2 .= u #u2 in out-of-place version
     end
@@ -1460,7 +1458,6 @@ end
 
     @.. broadcast=false σ=σ^(1 - q1) * u^q1
     @.. broadcast=false σ=σ + small_constant
-    @show ("3",σ)
 
     f.p(P2, u, p, t + c2 * dt) # evaluate production terms
     if issparse(P)
@@ -1483,13 +1480,11 @@ end
     linres = solve!(linsolve)
 
     u .= linres
-    @show ("4",u)
     integrator.stats.nsolve += 1
 
     if !(q1 ≈ q2)
         @.. broadcast=false σ=(uprev + small_constant)^(1 - q2) * tmp2^q2
         @.. broadcast=false σ=σ + small_constant
-        @show ("5",σ)
     end
 
     if issparse(P)
@@ -1508,17 +1503,10 @@ end
 
     # Same as linres = P3 \ tmp
     linsolve.A = P3
-    @show linsolve.A
-    @show linsolve.b
-    A_tmp = Matrix(linsolve.A)
-    b_tmp = linsolve.b
     linres = solve!(linsolve)
-    @show linres
-    @show A_tmp\b_tmp
     integrator.stats.nsolve += 1
 
     σ .= linres
-    @show ("6",σ)
     # avoid division by zero due to zero Patankar weights
     @.. broadcast=false σ=σ + small_constant
 
@@ -1543,7 +1531,6 @@ end
     linres = solve!(linsolve)
 
     u .= linres
-    @show ("7",u)
     integrator.stats.nsolve += 1
 
     # Now tmp stores the error estimate
diff --git a/test/runtests.jl b/test/runtests.jl
index 031213cb..8e658a84 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -202,75 +202,6 @@ function linear_advection_fd_upwind_D!(D, u, p, t)
     return nothing
 end
 
-
-@testset "Different matrix types (conservative, adaptive)" begin
-    prod_1! = (P, u, p, t) -> begin
-        fill!(P, zero(eltype(P)))
-        for i in 1:(length(u) - 1)
-            P[i, i + 1] = i * u[i]
-        end
-        return nothing
-    end
-
-    prod_2! = (P, u, p, t) -> begin
-        fill!(P, zero(eltype(P)))
-        for i in 1:(length(u) - 1)
-            P[i + 1, i] = i * u[i + 1]
-        end
-        return nothing
-    end
-
-    prod_3! = (P, u, p, t) -> begin
-        fill!(P, zero(eltype(P)))
-        for i in 1:(length(u) - 1)
-            P[i, i + 1] = i * u[i]
-            P[i + 1, i] = i * u[i + 1]
-        end
-        return nothing
-    end
-
-    n = 4
-    P_tridiagonal = Tridiagonal([0.1, 0.2, 0.3],
-                                zeros(n),
-                                [0.4, 0.5, 0.6])
-    P_dense = Matrix(P_tridiagonal)
-    P_sparse = sparse(P_tridiagonal)
-    u0 = [1.0, 1.5, 2.0, 2.5]
-    tspan = (0.0, 1.0)
-    dt = 0.25
-
-    rtol = sqrt(eps(Float32))
-
-    @testset "$alg" for alg in (MPRK43II(2.0 / 3.0), MPRK43II(0.5))
-
-        #=
-        # Something is going wrong on macOS with sparse matrices for MPRK43II, see
-        # https://github.com/SKopecz/PositiveIntegrators.jl/pull/101
-        if Sys.isapple() && (alg isa MPRK43II)
-            @test_skip false
-            continue
-        end
-        =#
-
-        for prod! in (prod_1!, prod_2!, prod_3!)
-            prod = (u, p, t) -> begin
-                P = similar(u, (length(u), length(u)))
-                prod!(P, u, p, t)
-                return P
-            end
-            prob_sparse_ip = ConservativePDSProblem(prod!, u0, tspan;
-                                                    p_prototype = P_sparse)
-            prob_sparse_op = ConservativePDSProblem(prod, u0, tspan;
-                                                    p_prototype = P_sparse)
-            sol_sparse_ip = solve(prob_sparse_ip, alg; dt)
-            sol_sparse_op = solve(prob_sparse_op, alg; dt)
-
-            @test isapprox(sol_sparse_ip.t, sol_sparse_op.t; rtol)
-            @test isapprox(sol_sparse_ip.u, sol_sparse_op.u; rtol)
-        end
-    end
-end
-#=
 @testset "PositiveIntegrators.jl tests" begin
     @testset "Aqua.jl" begin
         # We do not test ambiguities since we get a lot of
@@ -964,13 +895,6 @@ end
                                         MPRK43II(2.0 / 3.0), MPRK43II(0.5),
                                         SSPMPRK22(0.5, 1.0), SSPMPRK43())
 
-                # Something is going wrong on macOS with sparse matrices for MPRK43II, see
-                # https://github.com/SKopecz/PositiveIntegrators.jl/pull/101
-                if Sys.isapple() && (alg isa MPRK43II)
-                    @test_skip false
-                    continue
-                end
-
                 for prod! in (prod_1!, prod_2!, prod_3!)
                     prod = (u, p, t) -> begin
                         P = similar(u, (length(u), length(u)))
@@ -1803,5 +1727,4 @@ end
         sol = solve(prob_pds_linmod, MPRK22(1.0))
         @test_nowarn plot(sol)
     end
-end;
-=#
\ No newline at end of file
+end;
\ No newline at end of file

From e7d37311ef5a3c5d261a4975af426adab4566463 Mon Sep 17 00:00:00 2001
From: Stefan Kopecz <kopecz@mathematik.uni-kassel.de>
Date: Fri, 19 Jul 2024 17:22:18 +0200
Subject: [PATCH 07/42] format

---
 test/runtests.jl | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/test/runtests.jl b/test/runtests.jl
index 8e658a84..1469db90 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -894,7 +894,6 @@ end
                                         MPRK43I(1.0, 0.5), MPRK43I(0.5, 0.75),
                                         MPRK43II(2.0 / 3.0), MPRK43II(0.5),
                                         SSPMPRK22(0.5, 1.0), SSPMPRK43())
-
                 for prod! in (prod_1!, prod_2!, prod_3!)
                     prod = (u, p, t) -> begin
                         P = similar(u, (length(u), length(u)))
@@ -1727,4 +1726,4 @@ end
         sol = solve(prob_pds_linmod, MPRK22(1.0))
         @test_nowarn plot(sol)
     end
-end;
\ No newline at end of file
+end;

From 5b0cac1934b0323dfa4c58124778294310d8f839 Mon Sep 17 00:00:00 2001
From: Stefan Kopecz <kopecz@mathematik.uni-kassel.de>
Date: Fri, 19 Jul 2024 17:45:55 +0200
Subject: [PATCH 08/42] updated doc strings

---
 src/mprk.jl    | 2 +-
 src/sspmprk.jl | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/mprk.jl b/src/mprk.jl
index fbb0cf46..433b3084 100644
--- a/src/mprk.jl
+++ b/src/mprk.jl
@@ -982,7 +982,7 @@ algorithm from [LinearSolve.jl](https://github.com/SciML/LinearSolve.jl)
 as keyword argument `linsolve`.
 You can also choose the parameter `small_constant` which is added to all Patankar-weight denominators
 to avoid divisions by zero. You can pass a value explicitly, otherwise `small_constant` is set to
-`floatmin` of the floating point type used.
+`1e-50`.
 
 ## References
 
diff --git a/src/sspmprk.jl b/src/sspmprk.jl
index c25c79a3..ef0dcf57 100644
--- a/src/sspmprk.jl
+++ b/src/sspmprk.jl
@@ -467,7 +467,7 @@ algorithm from [LinearSolve.jl](https://github.com/SciML/LinearSolve.jl)
 as keyword argument `linsolve`.
 You can also choose the parameter `small_constant` which is added to all Patankar-weight denominators
 to avoid divisions by zero. You can pass a value explicitly, otherwise `small_constant` is set to
-`floatmin` of the floating point type used.
+`1e-50`.
 
 The current implementation only supports fixed time steps.
 

From a0dbe12e39373638a7beb1e016ed4cc024cb15d6 Mon Sep 17 00:00:00 2001
From: Stefan Kopecz <kopecz@mathematik.uni-kassel.de>
Date: Wed, 21 Aug 2024 14:34:34 +0200
Subject: [PATCH 09/42] added conversion to Float32 in SSPMPRK43

---
 src/sspmprk.jl   | 32 +++++++++++++++++++++++++-------
 test/runtests.jl | 18 ++++++++++++++++++
 2 files changed, 43 insertions(+), 7 deletions(-)

diff --git a/src/sspmprk.jl b/src/sspmprk.jl
index ef0dcf57..dedc39e9 100644
--- a/src/sspmprk.jl
+++ b/src/sspmprk.jl
@@ -488,7 +488,25 @@ struct SSPMPRK43{F, T} <: OrdinaryDiffEqAlgorithm
     small_constant_function::T
 end
 
-function SSPMPRK43(; linsolve = LUFactorization(), small_constant = 1e-50)
+function small_constant_function_SSPMPRK43(type)
+    if type == Float64
+        small_constant = 1e-50
+    elseif type == Float32
+        # small_constant_function is chosen such that the problem below
+        # (zero initial condition) can be solved
+        # P_linmod(u, p, t) = [0 u[2]; 5*u[1] 0]
+        # u0 = [1.0f0, 0.0f0]
+        # prob = ConservativePDSProblem(P_linmod, u0, (0.0f0, 2.0f0))
+        # sol = solve(prob, SSPMPRK43(small_constant=1f-8); dt=0.1f0)
+        small_constant = 1.0f-8
+    else
+        small_constant = floatmin(type)
+    end
+    return small_constant
+end
+
+function SSPMPRK43(; linsolve = LUFactorization(),
+                   small_constant = small_constant_function_SSPMPRK43)
     if isnothing(small_constant)
         small_constant_function = floatmin
     elseif small_constant isa Number
@@ -534,8 +552,7 @@ function get_constant_parameters(alg::SSPMPRK43)
     c3 = β20 + α21 * β10 + β21
 
     return n1, n2, z, η1, η2, η3, η4, η5, η6, s, α10, α20, α21, α30, α31, α32, β10, β20,
-           β21, β30,
-           β31, β32, c3
+           β21, β30, β31, β32, c3
 end
 
 struct SSPMPRK43ConstantCache{T} <: OrdinaryDiffEqConstantCache
@@ -573,11 +590,12 @@ function alg_cache(alg::SSPMPRK43, u, rate_prototype, ::Type{uEltypeNoUnits},
     if !(f isa PDSFunction || f isa ConservativePDSFunction)
         throw(ArgumentError("SSPMPRK43 can only be applied to production-destruction systems"))
     end
-    n1, n2, z, η1, η2, η3, η4, η5, η6, s, α10, α20, α21, α30, α31, α32, β10, β20, β21, β30, β31, β32, c3 = get_constant_parameters(alg)
+    const_param = get_constant_parameters(alg)
+    const_param = convert.(uEltypeNoUnits, const_param)
+    n1, n2, z, η1, η2, η3, η4, η5, η6, s, α10, α20, α21, α30, α31, α32, β10, β20, β21, β30, β31, β32, c3 = const_param
+    small_constant = alg.small_constant_function(uEltypeNoUnits)
     SSPMPRK43ConstantCache(n1, n2, z, η1, η2, η3, η4, η5, η6, s, α10, α20, α21, α30, α31,
-                           α32, β10,
-                           β20, β21, β30,
-                           β31, β32, c3, alg.small_constant_function(uEltypeNoUnits))
+                           α32, β10, β20, β21, β30, β31, β32, c3, small_constant)
 end
 
 function initialize!(integrator, cache::SSPMPRK43ConstantCache)
diff --git a/test/runtests.jl b/test/runtests.jl
index 1469db90..346cbbcf 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -633,6 +633,24 @@ end
                                                                                                  dt = 0.1)
         end
 
+        # Here we check that algorithms with input parameters return constants 
+        # of the same type as the inputs
+        @testset "Constant types" begin
+            algs = (MPRK22(0.5f0), MPRK22(1.0f0), MPRK22(2.0f0), MPRK43I(1.0f0, 0.5f0),
+                    MPRK43I(0.5f0, 0.75f0), MPRK43II(0.5f0), MPRK43II(2.0f0 / 3.0f0),
+                    SSPMPRK22(0.5f0, 1.0f0))
+            @testset "$i" for (i, alg) in enumerate(algs)
+                @test eltype(PositiveIntegrators.get_constant_parameters(alg)) == Float32
+            end
+
+            algs = (MPRK22(0.5), MPRK22(1.0), MPRK22(2.0), MPRK43I(1.0, 0.5),
+                    MPRK43I(0.5, 0.75), MPRK43II(0.5), MPRK43II(2.0 / 3.0),
+                    SSPMPRK22(0.5, 1.0))
+            @testset "$i" for (i, alg) in enumerate(algs)
+                @test eltype(PositiveIntegrators.get_constant_parameters(alg)) == Float64
+            end
+        end
+
         # Here we check that MPE equals implicit Euler (IE) for a linear PDS
         @testset "Linear model problem: MPE = IE?" begin
             # problem data

From b042cbf45a1c553c0d0fd9fb4adfd8b978fcb8dc Mon Sep 17 00:00:00 2001
From: Stefan Kopecz <kopecz@mathematik.uni-kassel.de>
Date: Wed, 21 Aug 2024 16:28:10 +0200
Subject: [PATCH 10/42] Added small_constant_function_MPRK43II

---
 src/mprk.jl    | 14 +++++++++++---
 src/sspmprk.jl |  7 +++----
 2 files changed, 14 insertions(+), 7 deletions(-)

diff --git a/src/mprk.jl b/src/mprk.jl
index 433b3084..5ca1248c 100644
--- a/src/mprk.jl
+++ b/src/mprk.jl
@@ -981,8 +981,7 @@ You can optionally choose the linear solver to be used by passing an
 algorithm from [LinearSolve.jl](https://github.com/SciML/LinearSolve.jl)
 as keyword argument `linsolve`.
 You can also choose the parameter `small_constant` which is added to all Patankar-weight denominators
-to avoid divisions by zero. You can pass a value explicitly, otherwise `small_constant` is set to
-`1e-50`.
+to avoid divisions by zero. The default for `Float64` computations is `1e-50`.
 
 ## References
 
@@ -998,7 +997,16 @@ struct MPRK43II{T, F, T2} <: OrdinaryDiffEqAdaptiveAlgorithm
     small_constant_function::T2
 end
 
-function MPRK43II(gamma; linsolve = LUFactorization(), small_constant = 1e-50)
+function small_constant_function_MPRK43II(type)
+    if type == Float64
+        small_constant = 1e-50
+    else
+        small_constant = floatmin(type)
+    end
+    return small_constant
+end
+
+function MPRK43II(gamma; linsolve = LUFactorization(), small_constant = small_constant_function_MPRK43II)
     if isnothing(small_constant)
         small_constant_function = floatmin
     elseif small_constant isa Number
diff --git a/src/sspmprk.jl b/src/sspmprk.jl
index dedc39e9..5562b13c 100644
--- a/src/sspmprk.jl
+++ b/src/sspmprk.jl
@@ -466,8 +466,7 @@ You can optionally choose the linear solver to be used by passing an
 algorithm from [LinearSolve.jl](https://github.com/SciML/LinearSolve.jl)
 as keyword argument `linsolve`.
 You can also choose the parameter `small_constant` which is added to all Patankar-weight denominators
-to avoid divisions by zero. You can pass a value explicitly, otherwise `small_constant` is set to
-`1e-50`.
+to avoid divisions by zero. The default for `Float64` computations is `1e-50`.
 
 The current implementation only supports fixed time steps.
 
@@ -492,12 +491,12 @@ function small_constant_function_SSPMPRK43(type)
     if type == Float64
         small_constant = 1e-50
     elseif type == Float32
-        # small_constant_function is chosen such that the problem below
+        # small_constant is chosen such that the problem below
         # (zero initial condition) can be solved
         # P_linmod(u, p, t) = [0 u[2]; 5*u[1] 0]
         # u0 = [1.0f0, 0.0f0]
         # prob = ConservativePDSProblem(P_linmod, u0, (0.0f0, 2.0f0))
-        # sol = solve(prob, SSPMPRK43(small_constant=1f-8); dt=0.1f0)
+        # sol = solve(prob, SSPMPRK43(); dt=0.1f0)
         small_constant = 1.0f-8
     else
         small_constant = floatmin(type)

From 9bd98430489407be1156a4c8d2997b1041c17594 Mon Sep 17 00:00:00 2001
From: Stefan Kopecz <kopecz@mathematik.uni-kassel.de>
Date: Wed, 21 Aug 2024 16:32:55 +0200
Subject: [PATCH 11/42] revised docstrings

---
 src/mprk.jl    | 3 ++-
 src/sspmprk.jl | 3 ++-
 2 files changed, 4 insertions(+), 2 deletions(-)

diff --git a/src/mprk.jl b/src/mprk.jl
index 5ca1248c..b3087712 100644
--- a/src/mprk.jl
+++ b/src/mprk.jl
@@ -981,7 +981,8 @@ You can optionally choose the linear solver to be used by passing an
 algorithm from [LinearSolve.jl](https://github.com/SciML/LinearSolve.jl)
 as keyword argument `linsolve`.
 You can also choose the parameter `small_constant` which is added to all Patankar-weight denominators
-to avoid divisions by zero. The default for `Float64` computations is `1e-50`.
+to avoid divisions by zero. To display the default value for data type `type` evaluate 
+`PositiveIntegrators.small_constant_function_MPRK43II(type)`.
 
 ## References
 
diff --git a/src/sspmprk.jl b/src/sspmprk.jl
index 5562b13c..23df69fd 100644
--- a/src/sspmprk.jl
+++ b/src/sspmprk.jl
@@ -466,7 +466,8 @@ You can optionally choose the linear solver to be used by passing an
 algorithm from [LinearSolve.jl](https://github.com/SciML/LinearSolve.jl)
 as keyword argument `linsolve`.
 You can also choose the parameter `small_constant` which is added to all Patankar-weight denominators
-to avoid divisions by zero. The default for `Float64` computations is `1e-50`.
+to avoid divisions by zero. To display the default value for data type `type` evaluate 
+`PositiveIntegrators.small_constant_function_SSPMPRK43(type)`.
 
 The current implementation only supports fixed time steps.
 

From e4d18c9690a6e567ffaba0dff0d6b9e7a813f1be Mon Sep 17 00:00:00 2001
From: Stefan Kopecz <kopecz@mathematik.uni-kassel.de>
Date: Wed, 21 Aug 2024 16:34:28 +0200
Subject: [PATCH 12/42] format

---
 src/mprk.jl | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/mprk.jl b/src/mprk.jl
index b3087712..cc6b80a8 100644
--- a/src/mprk.jl
+++ b/src/mprk.jl
@@ -1007,7 +1007,8 @@ function small_constant_function_MPRK43II(type)
     return small_constant
 end
 
-function MPRK43II(gamma; linsolve = LUFactorization(), small_constant = small_constant_function_MPRK43II)
+function MPRK43II(gamma; linsolve = LUFactorization(),
+                  small_constant = small_constant_function_MPRK43II)
     if isnothing(small_constant)
         small_constant_function = floatmin
     elseif small_constant isa Number

From 9a9a9f733d81e6513bd470ede571a465b0a4acaa Mon Sep 17 00:00:00 2001
From: Stefan Kopecz <kopecz@mathematik.uni-kassel.de>
Date: Wed, 21 Aug 2024 16:43:41 +0200
Subject: [PATCH 13/42] typo

---
 test/runtests.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/runtests.jl b/test/runtests.jl
index 346cbbcf..fab4f007 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -633,7 +633,7 @@ end
                                                                                                  dt = 0.1)
         end
 
-        # Here we check that algorithms with input parameters return constants 
+        # Here we check that algorithms which accept input parameters return constants 
         # of the same type as the inputs
         @testset "Constant types" begin
             algs = (MPRK22(0.5f0), MPRK22(1.0f0), MPRK22(2.0f0), MPRK43I(1.0f0, 0.5f0),

From 2b62b7a088d279d889dcddc9c3224387faa9aecf Mon Sep 17 00:00:00 2001
From: Stefan Kopecz <kopecz@mathematik.uni-kassel.de>
Date: Fri, 23 Aug 2024 13:46:24 +0200
Subject: [PATCH 14/42] Evaluation of PDSFunction and ConservativePDSFunction
 can handle units

---
 src/proddest.jl  |  35 ++++++---
 test/runtests.jl | 180 ++++++++++++++++++++++++++++++++++++-----------
 2 files changed, 163 insertions(+), 52 deletions(-)

diff --git a/src/proddest.jl b/src/proddest.jl
index 0cdf7dd3..9295a7bd 100644
--- a/src/proddest.jl
+++ b/src/proddest.jl
@@ -130,9 +130,15 @@ function (PD::PDSFunction)(du, u, p, t)
     PD.p(PD.p_prototype, u, p, t)
 
     if PD.p_prototype isa AbstractSparseMatrix
-        # Same result but more efficient - at least currently for SparseMatrixCSC
-        fill!(PD.d_prototype, one(eltype(PD.d_prototype)))
-        mul!(vec(du), PD.p_prototype, PD.d_prototype)
+        # Within mul! PD.d_prototype must not have units, in sum! PD.d_prototype needs units.   
+        # To avoid an additional auxiliary vector we use sum! also to compute 
+        # the row sum.
+
+        # row sum coded as matrix-vector product 
+        # fill!(PD.d_prototype, one(eltype(PD.d_prototype)))
+        # mul!(vec(du), PD.p_prototype, PD.d_prototype)
+        sum!(vec(du), PD.p_prototype)
+
         for i in 1:length(u)  #vec(du) .+= diag(PD.p_prototype)
             du[i] += PD.p_prototype[i, i]
         end
@@ -267,7 +273,7 @@ function (PD::ConservativePDSFunction)(u, p, t)
     #vec(sum(PD.p(u, p, t), dims = 2)) - vec(sum(PD.p(u, p, t), dims = 1))
     P = PD.p(u, p, t)
 
-    f = zero(u)
+    f = zeros(eltype(P), size(u))
     @fastmath @inbounds @simd for I in CartesianIndices(P)
         if !iszero(P[I])
             f[I[1]] += P[I]
@@ -280,8 +286,8 @@ end
 function (PD::ConservativePDSFunction)(u::SVector, p, t)
     P = PD.p(u, p, t)
 
-    f = similar(u) #constructs MVector
-    zeroT = zero(eltype(u))
+    f = similar(P[:, 1]) #constructs MVector
+    zeroT = zero(eltype(P))
     for i in eachindex(f)
         f[i] = zeroT
     end
@@ -306,9 +312,9 @@ end
 # Generic fallback (for dense arrays)
 # This implementation does not need any auxiliary vectors
 @inline function sum_terms!(du, tmp, P)
-    for i in 1:length(du)
+    for i in eachindex(du)
         du[i] = zero(eltype(du))
-        for j in 1:length(du)
+        for j in eachindex(du)
             du[i] += P[i, j] - P[j, i]
         end
     end
@@ -317,9 +323,18 @@ end
 
 # Same result but more efficient - at least currently for SparseMatrixCSC
 @inline function sum_terms!(du, tmp, P::AbstractSparseMatrix)
-    fill!(tmp, one(eltype(tmp)))
-    mul!(vec(du), P, tmp)
+    # Within mul! tmp must not have units, in sum! tmp needs units.   
+    # To avoid an additional auxiliary vector we use sum! also to compute 
+    # the row sum.
+
+    # # row sum coded as matrix vector product
+    # fill!(tmp, one(eltype(tmp)))
+    # mul!(vec(du), P, tmp)
+    sum!(vec(du), P)
+
+    #column sum
     sum!(tmp', P)
+
     vec(du) .-= tmp
     return nothing
 end
diff --git a/test/runtests.jl b/test/runtests.jl
index fab4f007..5468b64f 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -210,56 +210,152 @@ end
                       ambiguities = false,)
     end
 
-    @testset "ConservativePDSFunction" begin
-        prod_1! = (P, u, p, t) -> begin
-            fill!(P, zero(eltype(P)))
-            for i in 1:(length(u) - 1)
-                P[i, i + 1] = i * u[i]
+    @testset "Test for (Conservative)PDSFunctions" begin
+        @testset "ConservativePDSFunction" begin
+            prod_1! = (P, u, p, t) -> begin
+                fill!(P, zero(eltype(P)))
+                for i in 1:(length(u) - 1)
+                    P[i, i + 1] = i * u[i]
+                end
+                return nothing
             end
-            return nothing
-        end
-        prod_2! = (P, u, p, t) -> begin
-            fill!(P, zero(eltype(P)))
-            for i in 1:(length(u) - 1)
-                P[i + 1, i] = i * u[i + 1]
+            prod_2! = (P, u, p, t) -> begin
+                fill!(P, zero(eltype(P)))
+                for i in 1:(length(u) - 1)
+                    P[i + 1, i] = i * u[i + 1]
+                end
+                return nothing
             end
-            return nothing
-        end
-        prod_3! = (P, u, p, t) -> begin
-            fill!(P, zero(eltype(P)))
-            for i in 1:(length(u) - 1)
-                P[i, i + 1] = i * u[i]
-                P[i + 1, i] = i * u[i + 1]
+            prod_3! = (P, u, p, t) -> begin
+                fill!(P, zero(eltype(P)))
+                for i in 1:(length(u) - 1)
+                    P[i, i + 1] = i * u[i]
+                    P[i + 1, i] = i * u[i + 1]
+                end
+                return nothing
             end
-            return nothing
-        end
 
-        n = 10
-        P_tridiagonal = Tridiagonal(rand(n - 1), zeros(n), rand(n - 1))
-        P_dense = Matrix(P_tridiagonal)
-        P_sparse = sparse(P_tridiagonal)
-        u0 = rand(n)
-        tspan = (0.0, 1.0)
+            n = 10
+            P_tridiagonal = Tridiagonal(rand(n - 1), zeros(n), rand(n - 1))
+            P_dense = Matrix(P_tridiagonal)
+            P_sparse = sparse(P_tridiagonal)
+            u0 = rand(n)
+            tspan = (0.0, 1.0)
 
-        du_tridiagonal = similar(u0)
-        du_dense = similar(u0)
-        du_sparse = similar(u0)
+            du_tridiagonal = similar(u0)
+            du_dense = similar(u0)
+            du_sparse = similar(u0)
 
-        for prod! in (prod_1!, prod_2!, prod_3!)
-            prob_tridiagonal = ConservativePDSProblem(prod!, u0, tspan;
-                                                      p_prototype = P_tridiagonal)
-            prob_dense = ConservativePDSProblem(prod!, u0, tspan;
-                                                p_prototype = P_dense)
-            prob_sparse = ConservativePDSProblem(prod!, u0, tspan;
-                                                 p_prototype = P_sparse)
+            for prod! in (prod_1!, prod_2!, prod_3!)
+                prob_tridiagonal = ConservativePDSProblem(prod!, u0, tspan;
+                                                          p_prototype = P_tridiagonal)
+                prob_dense = ConservativePDSProblem(prod!, u0, tspan;
+                                                    p_prototype = P_dense)
+                prob_sparse = ConservativePDSProblem(prod!, u0, tspan;
+                                                     p_prototype = P_sparse)
 
-            prob_tridiagonal.f(du_tridiagonal, u0, nothing, 0.0)
-            prob_dense.f(du_dense, u0, nothing, 0.0)
-            prob_sparse.f(du_sparse, u0, nothing, 0.0)
+                prob_tridiagonal.f(du_tridiagonal, u0, nothing, 0.0)
+                prob_dense.f(du_dense, u0, nothing, 0.0)
+                prob_sparse.f(du_sparse, u0, nothing, 0.0)
 
-            @test du_tridiagonal ≈ du_dense
-            @test du_tridiagonal ≈ du_sparse
-            @test du_dense ≈ du_sparse
+                @test du_tridiagonal ≈ du_dense
+                @test du_tridiagonal ≈ du_sparse
+                @test du_dense ≈ du_sparse
+            end
+        end
+
+        # Here we check that PDSFunctions and ConservativePDSFunctions can be evaluated
+        # at (u, p, t) also in cases where u and t carry units.
+        # This is required when solvers form OrdinarDiffEq are used to solve a PDS.
+        @testset "(Conservative)PDSFunction and units" begin
+            @testset "(Conservative)PDSFunction and units (out-of-place)" begin
+                u0 = [0.9u"N", 0.1u"N"]
+                u0_static = @SVector [0.9u"N", 0.1u"N"]
+                t0 = 0.0u"s"
+                tspan = (t0, 2.0u"s")
+
+                f(u, p, t) = [u[2] - 5 * u[1]; -u[2] + 5 * u[1]] / u"s"
+                P(u, p, t) = [0u"N/s" u[2]/u"s"; 5 * u[1]/u"s" 0u"N/s"]
+                D(u, p, t) = [0.0; 0.0]u"N/s"
+                g = PositiveIntegrators.ConservativePDSFunction{false}(P)
+                h = PositiveIntegrators.PDSFunction{false}(P, D)
+
+                f_static(u, p, t) = SA[(u[2] - 5 * u[1]) / u"s"; (-u[2] + 5 * u[1]) / u"s"]
+                P_static(u, p, t) = SA[0u"N/s" u[2]/u"s"; 5 * u[1]/u"s" 0u"N/s"]
+                D_static(u, p, t) = SA[0.0u"N/s"; 0.0u"N/s"]
+                g_static = PositiveIntegrators.ConservativePDSFunction{false}(P_static)
+                h_static = PositiveIntegrators.PDSFunction{false}(P_static, D_static)
+
+                f1 = f(u0, nothing, t0)
+                g1 = g(u0, nothing, t0)
+                h1 = h(u0, nothing, t0)
+                f2 = f_static(u0_static, nothing, t0)
+                g2 = g_static(u0_static, nothing, t0)
+                h2 = h_static(u0_static, nothing, t0)
+
+                @test f1 == g1 == h1 == f2 == g2 == h2
+            end
+            @testset "(Conservative)PDSFunction and units (out-of-place)" begin
+                u0 = [1.0, 1.5, 2.0, 2.5]u"N"
+                t0 = 0.0u"s"
+
+                function f!(du, u, p, t)
+                    fill!(du, zero(eltype(du)))
+
+                    du[1] = (u[1] - u[2]) / u"s"
+                    du[2] = (3 * u[2] - 2 * u[3] - u[1]) / u"s"
+                    du[3] = (-2 * u[2] + 5 * u[3] - 3 * u[4]) / u"s"
+                    du[4] = (3 * u[4] - 3 * u[3]) / u"s"
+
+                    return nothing
+                end
+                function P!(P, u, p, t)
+                    fill!(P, zero(eltype(P)))
+                    for i in 1:(length(u) - 1)
+                        P[i, i + 1] = i * u[i] / u"s"
+                        P[i + 1, i] = i * u[i + 1] / u"s"
+                    end
+                    return nothing
+                end
+                function D!(D, u, p, t)
+                    fill!(D, zero(eltype(D)))
+                    return nothing
+                end
+
+                P_tridiagonal = Tridiagonal([0.1, 0.2, 0.3], zeros(4),
+                                            [0.4, 0.5, 0.6])u"N/s"
+                P_dense = Matrix(P_tridiagonal)
+                P_sparse = sparse(P_tridiagonal)
+                D_dense = zeros(4)u"N/s"
+
+                f_dense! = PositiveIntegrators.ConservativePDSFunction{true}(P!,
+                                                                             p_prototype = P_dense)
+                f_tridiagonal! = PositiveIntegrators.ConservativePDSFunction{true}(P!,
+                                                                                   p_prototype = P_tridiagonal)
+                f_sparse! = PositiveIntegrators.ConservativePDSFunction{true}(P!,
+                                                                              p_prototype = P_sparse)
+                g_dense! = PositiveIntegrators.PDSFunction{true}(P!, D!,
+                                                                 p_prototype = P_dense,
+                                                                 d_prototype = D_dense)
+                g_tridiagonal! = PositiveIntegrators.PDSFunction{true}(P!, D!,
+                                                                       p_prototype = P_tridiagonal,
+                                                                       d_prototype = D_dense)
+                g_sparse! = PositiveIntegrators.PDSFunction{true}(P!, D!,
+                                                                  p_prototype = P_sparse,
+                                                                  d_prototype = D_dense)
+
+                du1 = du2 = du3 = du4 = du5 = du6 = du7 = zeros(4)u"N/s"
+
+                f!(du1, u0, nothing, t0)
+                f_dense!(du2, u0, nothing, t0)
+                f_tridiagonal!(du3, u0, nothing, t0)
+                f_sparse!(du4, u0, nothing, t0)
+                g_dense!(du5, u0, nothing, t0)
+                g_tridiagonal!(du6, u0, nothing, t0)
+                g_sparse!(du7, u0, nothing, t0)
+
+                @test du1 == du2 == du3 == du4 == du5 == du6 == du7
+            end
         end
     end
 

From 3c98f617ccf9dda68e848a9325f0203fe780f4c1 Mon Sep 17 00:00:00 2001
From: Stefan Kopecz <kopecz@mathematik.uni-kassel.de>
Date: Fri, 23 Aug 2024 13:55:40 +0200
Subject: [PATCH 15/42] structure runtest

---
 test/runtests.jl | 318 +++++++++++++++++++++++++++--------------------
 1 file changed, 181 insertions(+), 137 deletions(-)

diff --git a/test/runtests.jl b/test/runtests.jl
index 5468b64f..55d2541d 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -210,152 +210,56 @@ end
                       ambiguities = false,)
     end
 
-    @testset "Test for (Conservative)PDSFunctions" begin
-        @testset "ConservativePDSFunction" begin
-            prod_1! = (P, u, p, t) -> begin
-                fill!(P, zero(eltype(P)))
-                for i in 1:(length(u) - 1)
-                    P[i, i + 1] = i * u[i]
-                end
-                return nothing
+    @testset "ConservativePDSFunction" begin
+        prod_1! = (P, u, p, t) -> begin
+            fill!(P, zero(eltype(P)))
+            for i in 1:(length(u) - 1)
+                P[i, i + 1] = i * u[i]
             end
-            prod_2! = (P, u, p, t) -> begin
-                fill!(P, zero(eltype(P)))
-                for i in 1:(length(u) - 1)
-                    P[i + 1, i] = i * u[i + 1]
-                end
-                return nothing
-            end
-            prod_3! = (P, u, p, t) -> begin
-                fill!(P, zero(eltype(P)))
-                for i in 1:(length(u) - 1)
-                    P[i, i + 1] = i * u[i]
-                    P[i + 1, i] = i * u[i + 1]
-                end
-                return nothing
-            end
-
-            n = 10
-            P_tridiagonal = Tridiagonal(rand(n - 1), zeros(n), rand(n - 1))
-            P_dense = Matrix(P_tridiagonal)
-            P_sparse = sparse(P_tridiagonal)
-            u0 = rand(n)
-            tspan = (0.0, 1.0)
-
-            du_tridiagonal = similar(u0)
-            du_dense = similar(u0)
-            du_sparse = similar(u0)
-
-            for prod! in (prod_1!, prod_2!, prod_3!)
-                prob_tridiagonal = ConservativePDSProblem(prod!, u0, tspan;
-                                                          p_prototype = P_tridiagonal)
-                prob_dense = ConservativePDSProblem(prod!, u0, tspan;
-                                                    p_prototype = P_dense)
-                prob_sparse = ConservativePDSProblem(prod!, u0, tspan;
-                                                     p_prototype = P_sparse)
-
-                prob_tridiagonal.f(du_tridiagonal, u0, nothing, 0.0)
-                prob_dense.f(du_dense, u0, nothing, 0.0)
-                prob_sparse.f(du_sparse, u0, nothing, 0.0)
-
-                @test du_tridiagonal ≈ du_dense
-                @test du_tridiagonal ≈ du_sparse
-                @test du_dense ≈ du_sparse
+            return nothing
+        end
+        prod_2! = (P, u, p, t) -> begin
+            fill!(P, zero(eltype(P)))
+            for i in 1:(length(u) - 1)
+                P[i + 1, i] = i * u[i + 1]
             end
+            return nothing
         end
-
-        # Here we check that PDSFunctions and ConservativePDSFunctions can be evaluated
-        # at (u, p, t) also in cases where u and t carry units.
-        # This is required when solvers form OrdinarDiffEq are used to solve a PDS.
-        @testset "(Conservative)PDSFunction and units" begin
-            @testset "(Conservative)PDSFunction and units (out-of-place)" begin
-                u0 = [0.9u"N", 0.1u"N"]
-                u0_static = @SVector [0.9u"N", 0.1u"N"]
-                t0 = 0.0u"s"
-                tspan = (t0, 2.0u"s")
-
-                f(u, p, t) = [u[2] - 5 * u[1]; -u[2] + 5 * u[1]] / u"s"
-                P(u, p, t) = [0u"N/s" u[2]/u"s"; 5 * u[1]/u"s" 0u"N/s"]
-                D(u, p, t) = [0.0; 0.0]u"N/s"
-                g = PositiveIntegrators.ConservativePDSFunction{false}(P)
-                h = PositiveIntegrators.PDSFunction{false}(P, D)
-
-                f_static(u, p, t) = SA[(u[2] - 5 * u[1]) / u"s"; (-u[2] + 5 * u[1]) / u"s"]
-                P_static(u, p, t) = SA[0u"N/s" u[2]/u"s"; 5 * u[1]/u"s" 0u"N/s"]
-                D_static(u, p, t) = SA[0.0u"N/s"; 0.0u"N/s"]
-                g_static = PositiveIntegrators.ConservativePDSFunction{false}(P_static)
-                h_static = PositiveIntegrators.PDSFunction{false}(P_static, D_static)
-
-                f1 = f(u0, nothing, t0)
-                g1 = g(u0, nothing, t0)
-                h1 = h(u0, nothing, t0)
-                f2 = f_static(u0_static, nothing, t0)
-                g2 = g_static(u0_static, nothing, t0)
-                h2 = h_static(u0_static, nothing, t0)
-
-                @test f1 == g1 == h1 == f2 == g2 == h2
+        prod_3! = (P, u, p, t) -> begin
+            fill!(P, zero(eltype(P)))
+            for i in 1:(length(u) - 1)
+                P[i, i + 1] = i * u[i]
+                P[i + 1, i] = i * u[i + 1]
             end
-            @testset "(Conservative)PDSFunction and units (out-of-place)" begin
-                u0 = [1.0, 1.5, 2.0, 2.5]u"N"
-                t0 = 0.0u"s"
+            return nothing
+        end
 
-                function f!(du, u, p, t)
-                    fill!(du, zero(eltype(du)))
+        n = 10
+        P_tridiagonal = Tridiagonal(rand(n - 1), zeros(n), rand(n - 1))
+        P_dense = Matrix(P_tridiagonal)
+        P_sparse = sparse(P_tridiagonal)
+        u0 = rand(n)
+        tspan = (0.0, 1.0)
 
-                    du[1] = (u[1] - u[2]) / u"s"
-                    du[2] = (3 * u[2] - 2 * u[3] - u[1]) / u"s"
-                    du[3] = (-2 * u[2] + 5 * u[3] - 3 * u[4]) / u"s"
-                    du[4] = (3 * u[4] - 3 * u[3]) / u"s"
+        du_tridiagonal = similar(u0)
+        du_dense = similar(u0)
+        du_sparse = similar(u0)
 
-                    return nothing
-                end
-                function P!(P, u, p, t)
-                    fill!(P, zero(eltype(P)))
-                    for i in 1:(length(u) - 1)
-                        P[i, i + 1] = i * u[i] / u"s"
-                        P[i + 1, i] = i * u[i + 1] / u"s"
-                    end
-                    return nothing
-                end
-                function D!(D, u, p, t)
-                    fill!(D, zero(eltype(D)))
-                    return nothing
-                end
-
-                P_tridiagonal = Tridiagonal([0.1, 0.2, 0.3], zeros(4),
-                                            [0.4, 0.5, 0.6])u"N/s"
-                P_dense = Matrix(P_tridiagonal)
-                P_sparse = sparse(P_tridiagonal)
-                D_dense = zeros(4)u"N/s"
-
-                f_dense! = PositiveIntegrators.ConservativePDSFunction{true}(P!,
-                                                                             p_prototype = P_dense)
-                f_tridiagonal! = PositiveIntegrators.ConservativePDSFunction{true}(P!,
-                                                                                   p_prototype = P_tridiagonal)
-                f_sparse! = PositiveIntegrators.ConservativePDSFunction{true}(P!,
-                                                                              p_prototype = P_sparse)
-                g_dense! = PositiveIntegrators.PDSFunction{true}(P!, D!,
-                                                                 p_prototype = P_dense,
-                                                                 d_prototype = D_dense)
-                g_tridiagonal! = PositiveIntegrators.PDSFunction{true}(P!, D!,
-                                                                       p_prototype = P_tridiagonal,
-                                                                       d_prototype = D_dense)
-                g_sparse! = PositiveIntegrators.PDSFunction{true}(P!, D!,
-                                                                  p_prototype = P_sparse,
-                                                                  d_prototype = D_dense)
+        for prod! in (prod_1!, prod_2!, prod_3!)
+            prob_tridiagonal = ConservativePDSProblem(prod!, u0, tspan;
+                                                      p_prototype = P_tridiagonal)
+            prob_dense = ConservativePDSProblem(prod!, u0, tspan;
+                                                p_prototype = P_dense)
+            prob_sparse = ConservativePDSProblem(prod!, u0, tspan;
+                                                 p_prototype = P_sparse)
 
-                du1 = du2 = du3 = du4 = du5 = du6 = du7 = zeros(4)u"N/s"
+            prob_tridiagonal.f(du_tridiagonal, u0, nothing, 0.0)
+            prob_dense.f(du_dense, u0, nothing, 0.0)
+            prob_sparse.f(du_sparse, u0, nothing, 0.0)
 
-                f!(du1, u0, nothing, t0)
-                f_dense!(du2, u0, nothing, t0)
-                f_tridiagonal!(du3, u0, nothing, t0)
-                f_sparse!(du4, u0, nothing, t0)
-                g_dense!(du5, u0, nothing, t0)
-                g_tridiagonal!(du6, u0, nothing, t0)
-                g_sparse!(du7, u0, nothing, t0)
-
-                @test du1 == du2 == du3 == du4 == du5 == du6 == du7
-            end
+            @test du_tridiagonal ≈ du_dense
+            @test du_tridiagonal ≈ du_sparse
+            @test du_dense ≈ du_sparse
         end
     end
 
@@ -648,6 +552,146 @@ end
             @test 0.95 < alloc1 / alloc5 < 1.05
             @test 0.95 < alloc1 / alloc6 < 1.05
         end
+
+        # Here we check that PDSFunctions and ConservativePDSFunctions can be evaluated
+        # at (u, p, t) also in cases where u and t carry units.
+        # This is required when solvers form OrdinarDiffEq are used to solve a PDS
+        # and std_rhs is not specified.
+        @testset "(Conservative)PDSFunction and units" begin
+            @testset "(Conservative)PDSFunction and units (out-of-place)" begin
+                u0 = [0.9u"N", 0.1u"N"]
+                u0_static = @SVector [0.9u"N", 0.1u"N"]
+                t0 = 0.0u"s"
+                tspan = (t0, 2.0u"s")
+
+                f(u, p, t) = [u[2] - 5 * u[1]; -u[2] + 5 * u[1]] / u"s"
+                P(u, p, t) = [0u"N/s" u[2]/u"s"; 5 * u[1]/u"s" 0u"N/s"]
+                D(u, p, t) = [0.0; 0.0]u"N/s"
+                g = PositiveIntegrators.ConservativePDSFunction{false}(P)
+                h = PositiveIntegrators.PDSFunction{false}(P, D)
+
+                f_static(u, p, t) = SA[(u[2] - 5 * u[1]) / u"s"; (-u[2] + 5 * u[1]) / u"s"]
+                P_static(u, p, t) = SA[0u"N/s" u[2]/u"s"; 5 * u[1]/u"s" 0u"N/s"]
+                D_static(u, p, t) = SA[0.0u"N/s"; 0.0u"N/s"]
+                g_static = PositiveIntegrators.ConservativePDSFunction{false}(P_static)
+                h_static = PositiveIntegrators.PDSFunction{false}(P_static, D_static)
+
+                f1 = f(u0, nothing, t0)
+                g1 = g(u0, nothing, t0)
+                h1 = h(u0, nothing, t0)
+                f2 = f_static(u0_static, nothing, t0)
+                g2 = g_static(u0_static, nothing, t0)
+                h2 = h_static(u0_static, nothing, t0)
+
+                @test f1 == g1 == h1 == f2 == g2 == h2
+            end
+            @testset "(Conservative)PDSFunction and units (out-of-place)" begin
+                u0 = [1.0, 1.5, 2.0, 2.5]u"N"
+                t0 = 0.0u"s"
+
+                function f!(du, u, p, t)
+                    fill!(du, zero(eltype(du)))
+
+                    du[1] = (u[1] - u[2]) / u"s"
+                    du[2] = (3 * u[2] - 2 * u[3] - u[1]) / u"s"
+                    du[3] = (-2 * u[2] + 5 * u[3] - 3 * u[4]) / u"s"
+                    du[4] = (3 * u[4] - 3 * u[3]) / u"s"
+
+                    return nothing
+                end
+                function P!(P, u, p, t)
+                    fill!(P, zero(eltype(P)))
+                    for i in 1:(length(u) - 1)
+                        P[i, i + 1] = i * u[i] / u"s"
+                        P[i + 1, i] = i * u[i + 1] / u"s"
+                    end
+                    return nothing
+                end
+                function D!(D, u, p, t)
+                    fill!(D, zero(eltype(D)))
+                    return nothing
+                end
+
+                P_tridiagonal = Tridiagonal([0.1, 0.2, 0.3], zeros(4),
+                                            [0.4, 0.5, 0.6])u"N/s"
+                P_dense = Matrix(P_tridiagonal)
+                P_sparse = sparse(P_tridiagonal)
+                D_dense = zeros(4)u"N/s"
+
+                f_dense! = PositiveIntegrators.ConservativePDSFunction{true}(P!,
+                                                                             p_prototype = P_dense)
+                f_tridiagonal! = PositiveIntegrators.ConservativePDSFunction{true}(P!,
+                                                                                   p_prototype = P_tridiagonal)
+                f_sparse! = PositiveIntegrators.ConservativePDSFunction{true}(P!,
+                                                                              p_prototype = P_sparse)
+                g_dense! = PositiveIntegrators.PDSFunction{true}(P!, D!,
+                                                                 p_prototype = P_dense,
+                                                                 d_prototype = D_dense)
+                g_tridiagonal! = PositiveIntegrators.PDSFunction{true}(P!, D!,
+                                                                       p_prototype = P_tridiagonal,
+                                                                       d_prototype = D_dense)
+                g_sparse! = PositiveIntegrators.PDSFunction{true}(P!, D!,
+                                                                  p_prototype = P_sparse,
+                                                                  d_prototype = D_dense)
+
+                du1 = du2 = du3 = du4 = du5 = du6 = du7 = zeros(4)u"N/s"
+
+                f!(du1, u0, nothing, t0)
+                f_dense!(du2, u0, nothing, t0)
+                f_tridiagonal!(du3, u0, nothing, t0)
+                f_sparse!(du4, u0, nothing, t0)
+                g_dense!(du5, u0, nothing, t0)
+                g_tridiagonal!(du6, u0, nothing, t0)
+                g_sparse!(du7, u0, nothing, t0)
+
+                @test du1 == du2 == du3 == du4 == du5 == du6 == du7
+            end
+        end
+
+        @testset "Compatibility of OrdinarDiffEq solvers with (Conservative)PDSProblems and units" begin
+            @testset "out-of-place" begin
+                u0 = [0.9u"N", 0.1u"N"]
+                u0_static = @SVector [0.9u"N", 0.1u"N"]
+                tspan = (0.0u"s", 2.0u"s")
+
+                f(u, p, t) = [u[2] - 5 * u[1]; -u[2] + 5 * u[1]] / u"s"
+                P(u, p, t) = [0u"N/s" u[2]/u"s"; 5 * u[1]/u"s" 0u"N/s"]
+                D(u, p, t) = [0.0; 0.0]u"N/s"
+                f_static(u, p, t) = SA[(u[2] - 5 * u[1]) / u"s"; (-u[2] + 5 * u[1]) / u"s"]
+                P_static(u, p, t) = SA[0u"N/s" u[2]/u"s"; 5 * u[1]/u"s" 0u"N/s"]
+                D_static(u, p, t) = SA[0.0u"N/s"; 0.0u"N/s"]
+
+                prob_ode = ODEProblem(f, u0, tspan)
+                prob_cpds = ConservativePDSProblem(P, u0, tspan)
+                prob_pds = PDSProblem(P, D, u0, tspan)
+                #prob_cpds_rhs = ConservativePDSProblem(P, u0, tspan; std_rhs = f)
+                #prob_pds_rhs = PDSProblem(P, D, u0, tspan; std_rhs = f)
+                prob_ode_static = ODEProblem(f, u0_static, tspan)
+                prob_cpds_static = ConservativePDSProblem(P, u0_static, tspan)
+                prob_pds_static = PDSProblem(P, D, u0_static, tspan)
+                #prob_cpds_static_rhs = ConservativePDSProblem(P, u0_static, tspan; std_rhs = f)
+                #prob_pds_static_rhs = PDSProblem(P, D, u0_static, tspan; std_rhs = f)
+
+                algs = (Euler(), Tsit5())
+
+                for alg in algs
+                    alg = Euler()
+                    sol_ode = solve(prob_ode, alg; dt = 0.1u"s")
+                    sol_cpds = solve(prob_cpds, alg; dt = 0.1u"s")
+                    sol_pds = solve(prob_pds, alg; dt = 0.1u"s")
+                    #sol_cpds_rhs = solve(prob_cpds_rhs, alg; dt = 0.1u"s")
+                    #sol_pds_rhs = solve(prob_pds_rhs, alg; dt = 0.1u"s")
+                    sol_ode_static = solve(prob_ode, alg; dt = 0.1u"s")
+                    sol_cpds_static = solve(prob_cpds, alg; dt = 0.1u"s")
+                    sol_pds_static = solve(prob_pds, alg; dt = 0.1u"s")
+                    #sol_cpds_static_rhs = solve(prob_cpds_rhs, alg; dt = 0.1u"s")
+                    #sol_pds_static_rhs = solve(prob_pds_rhs, alg; dt = 0.1u"s")        
+
+                    @test sol_ode ≈ sol_cpds ≈ sol_pds ≈ sol_ode_static ≈ sol_cpds_static ≈
+                          sol_pds_static
+                end
+            end
+        end
     end
 
     @testset "PDS Solvers" begin

From 27adf7a238e084ac655c1dae8e33a74b75ec6c32 Mon Sep 17 00:00:00 2001
From: Stefan Kopecz <kopecz@mathematik.uni-kassel.de>
Date: Fri, 23 Aug 2024 14:01:23 +0200
Subject: [PATCH 16/42] bugfix runtest

---
 test/runtests.jl | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/test/runtests.jl b/test/runtests.jl
index d19ebc29..432e097f 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -802,7 +802,9 @@ end
                     @test sol_ode ≈ sol_cpds ≈ sol_pds ≈ sol_ode_static ≈ sol_cpds_static ≈
                           sol_pds_static
                 end
-            
+            end
+        end
+
         # Here we check that production-destruction form and standard ODE form fit together in predefinded problems,
         # i.e. standard solvers using std_rhs should generate results that are equal to those without specifying std_rhs
         @testset "Check that production-destruction form and standard ODE fit together in predefinded problems" begin
@@ -888,7 +890,6 @@ end
                 sol3 = solve(ODEProblem(prob.f.std_rhs, prob.u0, prob.tspan), alg; dt) # use f to create ODEProblem
                 @test sol.t ≈ sol2.t ≈ sol3.t
                 @test sol.u ≈ sol2.u ≈ sol3.u
-
             end
         end
     end

From 9df751d5658fa01b84ccb3101b9ed0aa64a0fb91 Mon Sep 17 00:00:00 2001
From: Stefan Kopecz <kopecz@mathematik.uni-kassel.de>
Date: Fri, 23 Aug 2024 14:06:37 +0200
Subject: [PATCH 17/42] additoinal tests

---
 test/runtests.jl | 22 ++++++++++++----------
 1 file changed, 12 insertions(+), 10 deletions(-)

diff --git a/test/runtests.jl b/test/runtests.jl
index 432e097f..dbc16e90 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -776,13 +776,14 @@ end
                 prob_ode = ODEProblem(f, u0, tspan)
                 prob_cpds = ConservativePDSProblem(P, u0, tspan)
                 prob_pds = PDSProblem(P, D, u0, tspan)
-                #prob_cpds_rhs = ConservativePDSProblem(P, u0, tspan; std_rhs = f)
-                #prob_pds_rhs = PDSProblem(P, D, u0, tspan; std_rhs = f)
+                prob_cpds_rhs = ConservativePDSProblem(P, u0, tspan; std_rhs = f)
+                prob_pds_rhs = PDSProblem(P, D, u0, tspan; std_rhs = f)
                 prob_ode_static = ODEProblem(f, u0_static, tspan)
                 prob_cpds_static = ConservativePDSProblem(P, u0_static, tspan)
                 prob_pds_static = PDSProblem(P, D, u0_static, tspan)
-                #prob_cpds_static_rhs = ConservativePDSProblem(P, u0_static, tspan; std_rhs = f)
-                #prob_pds_static_rhs = PDSProblem(P, D, u0_static, tspan; std_rhs = f)
+                prob_cpds_static_rhs = ConservativePDSProblem(P, u0_static, tspan;
+                                                              std_rhs = f_static)
+                prob_pds_static_rhs = PDSProblem(P, D, u0_static, tspan; std_rhs = f_static)
 
                 algs = (Euler(), Tsit5())
 
@@ -791,16 +792,17 @@ end
                     sol_ode = solve(prob_ode, alg; dt = 0.1u"s")
                     sol_cpds = solve(prob_cpds, alg; dt = 0.1u"s")
                     sol_pds = solve(prob_pds, alg; dt = 0.1u"s")
-                    #sol_cpds_rhs = solve(prob_cpds_rhs, alg; dt = 0.1u"s")
-                    #sol_pds_rhs = solve(prob_pds_rhs, alg; dt = 0.1u"s")
+                    sol_cpds_rhs = solve(prob_cpds_rhs, alg; dt = 0.1u"s")
+                    sol_pds_rhs = solve(prob_pds_rhs, alg; dt = 0.1u"s")
                     sol_ode_static = solve(prob_ode, alg; dt = 0.1u"s")
                     sol_cpds_static = solve(prob_cpds, alg; dt = 0.1u"s")
                     sol_pds_static = solve(prob_pds, alg; dt = 0.1u"s")
-                    #sol_cpds_static_rhs = solve(prob_cpds_rhs, alg; dt = 0.1u"s")
-                    #sol_pds_static_rhs = solve(prob_pds_rhs, alg; dt = 0.1u"s")        
+                    sol_cpds_static_rhs = solve(prob_cpds_rhs, alg; dt = 0.1u"s")
+                    sol_pds_static_rhs = solve(prob_pds_rhs, alg; dt = 0.1u"s")
 
-                    @test sol_ode ≈ sol_cpds ≈ sol_pds ≈ sol_ode_static ≈ sol_cpds_static ≈
-                          sol_pds_static
+                    @test sol_ode ≈ sol_cpds ≈ sol_pds ≈ sol_cpds_rhs ≈ sol_pds_rhs ≈
+                          sol_ode_static ≈ sol_cpds_static ≈
+                          sol_pds_static ≈ sol_cpds_static_rhs ≈ sol_pds_static_rhs
                 end
             end
         end

From 28f762835e46bcaca1770090dacde73ebbe8968f Mon Sep 17 00:00:00 2001
From: Stefan Kopecz <kopecz@mathematik.uni-kassel.de>
Date: Fri, 23 Aug 2024 14:15:30 +0200
Subject: [PATCH 18/42] Use Unitful in tests

---
 test/Project.toml | 2 ++
 test/runtests.jl  | 2 ++
 2 files changed, 4 insertions(+)

diff --git a/test/Project.toml b/test/Project.toml
index 8855f505..5182e22f 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -9,6 +9,7 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [compat]
 Aqua = "0.7, 0.8"
@@ -18,3 +19,4 @@ OrdinaryDiffEq = "6.59"
 Plots = "1.11.1"
 StaticArrays = "1.5"
 Statistics = "1"
+Unitful = "1.21"
diff --git a/test/runtests.jl b/test/runtests.jl
index dbc16e90..7625e9ba 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -5,6 +5,8 @@ using Statistics: mean
 
 using StaticArrays: MVector
 
+using Unitful
+
 using OrdinaryDiffEq
 using PositiveIntegrators
 

From be16982d1b5cd89da46dcf011c64b6f3c40b6ed1 Mon Sep 17 00:00:00 2001
From: Stefan Kopecz <kopecz@mathematik.uni-kassel.de>
Date: Fri, 23 Aug 2024 14:39:30 +0200
Subject: [PATCH 19/42] More functionality from StaticArrays

---
 test/runtests.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/runtests.jl b/test/runtests.jl
index 7625e9ba..04d986e7 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -3,7 +3,7 @@ using LinearAlgebra
 using SparseArrays
 using Statistics: mean
 
-using StaticArrays: MVector
+using StaticArrays: MVector, @SVector, SA
 
 using Unitful
 

From d41c143850f226f6ea676e578bef3dae7a1cfb4b Mon Sep 17 00:00:00 2001
From: Stefan Kopecz <kopecz@mathematik.uni-kassel.de>
Date: Fri, 23 Aug 2024 15:29:07 +0200
Subject: [PATCH 20/42] bugfixes

---
 src/PositiveIntegrators.jl | 2 +-
 src/proddest.jl            | 2 +-
 test/runtests.jl           | 2 +-
 3 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/PositiveIntegrators.jl b/src/PositiveIntegrators.jl
index 82d54784..979e8927 100644
--- a/src/PositiveIntegrators.jl
+++ b/src/PositiveIntegrators.jl
@@ -1,7 +1,7 @@
 module PositiveIntegrators
 
 # 1. Load dependencies
-using LinearAlgebra: LinearAlgebra, Tridiagonal, I, diag, mul!
+using LinearAlgebra: LinearAlgebra, Tridiagonal, I, diag#, mul!
 using SparseArrays: SparseArrays, AbstractSparseMatrix,
                     issparse, nonzeros, nzrange, rowvals, spdiagm
 using StaticArrays: SVector, SMatrix, StaticArray, @SVector, @SMatrix
diff --git a/src/proddest.jl b/src/proddest.jl
index a90caf15..ba2e0df4 100644
--- a/src/proddest.jl
+++ b/src/proddest.jl
@@ -121,7 +121,7 @@ end
 # Most specific constructor for PDSFunction
 function PDSFunction{iip, FullSpecialize}(P, D;
                                           p_prototype = nothing,
-                                          d_prototype,
+                                          d_prototype = nothing,
                                           analytic = nothing,
                                           std_rhs = nothing) where {iip}
     if std_rhs === nothing
diff --git a/test/runtests.jl b/test/runtests.jl
index 04d986e7..559180bb 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -699,7 +699,7 @@ end
 
                 @test f1 == g1 == h1 == f2 == g2 == h2
             end
-            @testset "(Conservative)PDSFunction and units (out-of-place)" begin
+            @testset "(Conservative)PDSFunction and units (in-place)" begin
                 u0 = [1.0, 1.5, 2.0, 2.5]u"N"
                 t0 = 0.0u"s"
 

From faa1c340811c5c4a92377f9cd4ff1c841473fb49 Mon Sep 17 00:00:00 2001
From: Stefan Kopecz <kopecz@mathematik.uni-kassel.de>
Date: Fri, 23 Aug 2024 16:18:44 +0200
Subject: [PATCH 21/42] Reset evaluation of PDSFunctions with sparse prototyp

---
 src/PositiveIntegrators.jl |  2 +-
 src/proddest.jl            | 21 ++++++---------------
 2 files changed, 7 insertions(+), 16 deletions(-)

diff --git a/src/PositiveIntegrators.jl b/src/PositiveIntegrators.jl
index 979e8927..82d54784 100644
--- a/src/PositiveIntegrators.jl
+++ b/src/PositiveIntegrators.jl
@@ -1,7 +1,7 @@
 module PositiveIntegrators
 
 # 1. Load dependencies
-using LinearAlgebra: LinearAlgebra, Tridiagonal, I, diag#, mul!
+using LinearAlgebra: LinearAlgebra, Tridiagonal, I, diag, mul!
 using SparseArrays: SparseArrays, AbstractSparseMatrix,
                     issparse, nonzeros, nzrange, rowvals, spdiagm
 using StaticArrays: SVector, SMatrix, StaticArray, @SVector, @SMatrix
diff --git a/src/proddest.jl b/src/proddest.jl
index ba2e0df4..3a30c4d9 100644
--- a/src/proddest.jl
+++ b/src/proddest.jl
@@ -163,14 +163,9 @@ function (PD::PDSStdRHS)(du, u, p, t)
     PD.p(PD.p_prototype, u, p, t)
 
     if PD.p_prototype isa AbstractSparseMatrix
-        # Within mul! PD.d_prototype must not have units, in sum! PD.d_prototype needs units.   
-        # To avoid an additional auxiliary vector we use sum! also to compute 
-        # the row sum.
-
         # row sum coded as matrix-vector product 
-        # fill!(PD.d_prototype, one(eltype(PD.d_prototype)))
-        # mul!(vec(du), PD.p_prototype, PD.d_prototype)
-        sum!(vec(du), PD.p_prototype)
+        fill!(PD.d_prototype, one(eltype(PD.tmp)))
+        mul!(vec(du), PD.p_prototype, PD.tmp)
 
         for i in 1:length(u)  #vec(du) .+= diag(PD.p_prototype)
             du[i] += PD.p_prototype[i, i]
@@ -390,14 +385,10 @@ end
 
 # Same result but more efficient - at least currently for SparseMatrixCSC
 @inline function sum_terms!(du, tmp, P::AbstractSparseMatrix)
-    # Within mul! tmp must not have units, in sum! tmp needs units.   
-    # To avoid an additional auxiliary vector we use sum! also to compute 
-    # the row sum.
-
-    # # row sum coded as matrix vector product
-    # fill!(tmp, one(eltype(tmp)))
-    # mul!(vec(du), P, tmp)
-    sum!(vec(du), P)
+    # row sum coded as matrix vector product
+    fill!(tmp, one(eltype(tmp)))
+    mul!(vec(du), P, tmp)
+    #sum!(vec(du), P)
 
     #column sum
     sum!(tmp', P)

From 58f234fe5bd94b470632e54bea62f3ed61297422 Mon Sep 17 00:00:00 2001
From: Stefan Kopecz <kopecz@mathematik.uni-kassel.de>
Date: Fri, 23 Aug 2024 16:52:10 +0200
Subject: [PATCH 22/42] bugfix

---
 src/proddest.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/proddest.jl b/src/proddest.jl
index 3a30c4d9..e626d88e 100644
--- a/src/proddest.jl
+++ b/src/proddest.jl
@@ -164,8 +164,8 @@ function (PD::PDSStdRHS)(du, u, p, t)
 
     if PD.p_prototype isa AbstractSparseMatrix
         # row sum coded as matrix-vector product 
-        fill!(PD.d_prototype, one(eltype(PD.tmp)))
-        mul!(vec(du), PD.p_prototype, PD.tmp)
+        fill!(PD.d_prototype, one(eltype(PD.d_prototype)))
+        mul!(vec(du), PD.p_prototype, PD.d_prototype)
 
         for i in 1:length(u)  #vec(du) .+= diag(PD.p_prototype)
             du[i] += PD.p_prototype[i, i]

From fd6806590133cd0728edb12f5762fd9b273b8bb3 Mon Sep 17 00:00:00 2001
From: Stefan Kopecz <kopecz@mathematik.uni-kassel.de>
Date: Fri, 23 Aug 2024 17:47:31 +0200
Subject: [PATCH 23/42] Evaluation of PDSFunctions with sparse Prototype and
 units is working

---
 src/proddest.jl | 35 ++++++++++++++++++++++++-----------
 1 file changed, 24 insertions(+), 11 deletions(-)

diff --git a/src/proddest.jl b/src/proddest.jl
index e626d88e..8fd5df13 100644
--- a/src/proddest.jl
+++ b/src/proddest.jl
@@ -143,11 +143,21 @@ function (PD::PDSFunction)(du, u, p, t)
 end
 
 # Default implementation of the standard right-hand side evaluation function
-struct PDSStdRHS{P, D, PrototypeP, PrototypeD} <: Function
+struct PDSStdRHS{P, D, PrototypeP, PrototypeD, TMP} <: Function
     p::P
     d::D
     p_prototype::PrototypeP
     d_prototype::PrototypeD
+    tmp::TMP
+end
+
+function PDSStdRHS(P, D, p_prototype, d_prototype)
+    if p_prototype isa AbstractSparseMatrix
+        tmp = zeros(eltype(p_prototype), (size(p_prototype, 1),))/oneunit(first(p_prototype)) # drop units
+    else
+        tmp = nothing
+    end
+    PDSStdRHS(P, D, p_prototype, d_prototype, tmp)
 end
 
 # Evaluation of a PDSStdRHS (out-of-place)
@@ -164,8 +174,8 @@ function (PD::PDSStdRHS)(du, u, p, t)
 
     if PD.p_prototype isa AbstractSparseMatrix
         # row sum coded as matrix-vector product 
-        fill!(PD.d_prototype, one(eltype(PD.d_prototype)))
-        mul!(vec(du), PD.p_prototype, PD.d_prototype)
+        fill!(PD.tmp, one(eltype(PD.tmp)))
+        mul!(vec(du), PD.p_prototype, PD.tmp)
 
         for i in 1:length(u)  #vec(du) .+= diag(PD.p_prototype)
             du[i] += PD.p_prototype[i, i]
@@ -315,19 +325,22 @@ function (PD::ConservativePDSFunction)(du, u, p, t)
 end
 
 # Default implementation of the standard right-hand side evaluation function
-struct ConservativePDSStdRHS{P, PrototypeP, TMP} <: Function
+struct ConservativePDSStdRHS{P, PrototypeP, TMP, TMP2} <: Function
     p::P
     p_prototype::PrototypeP
     tmp::TMP
+    tmp2::TMP2
 end
 
 function ConservativePDSStdRHS(P, p_prototype)
     if p_prototype isa AbstractSparseMatrix
         tmp = zeros(eltype(p_prototype), (size(p_prototype, 1),))
+        tmp2 = tmp/oneunit(first(tmp)) # drop units
     else
         tmp = nothing
+        tmp2 = nothing
     end
-    ConservativePDSStdRHS(P, p_prototype, tmp)
+    ConservativePDSStdRHS(P, p_prototype, tmp, tmp2)
 end
 
 # Evaluation of a ConservativePDSStdRHS (out-of-place)
@@ -367,13 +380,13 @@ end
 # Evaluation of a ConservativePDSStdRHS (in-place)
 function (PD::ConservativePDSStdRHS)(du, u, p, t)
     PD.p(PD.p_prototype, u, p, t)
-    sum_terms!(du, PD.tmp, PD.p_prototype)
+    sum_terms!(du, PD.tmp, PD.tmp2, PD.p_prototype)
     return nothing
 end
 
 # Generic fallback (for dense arrays)
 # This implementation does not need any auxiliary vectors
-@inline function sum_terms!(du, tmp, P)
+@inline function sum_terms!(du, tmp, tmp2, P)
     for i in eachindex(du)
         du[i] = zero(eltype(du))
         for j in eachindex(du)
@@ -384,10 +397,10 @@ end
 end
 
 # Same result but more efficient - at least currently for SparseMatrixCSC
-@inline function sum_terms!(du, tmp, P::AbstractSparseMatrix)
+@inline function sum_terms!(du, tmp, tmp2, P::AbstractSparseMatrix)
     # row sum coded as matrix vector product
-    fill!(tmp, one(eltype(tmp)))
-    mul!(vec(du), P, tmp)
+    fill!(tmp2, one(eltype(tmp2)))
+    mul!(vec(du), P, tmp2)
     #sum!(vec(du), P)
 
     #column sum
@@ -397,7 +410,7 @@ end
     return nothing
 end
 
-@inline function sum_terms!(du, tmp, P::Tridiagonal)
+@inline function sum_terms!(du, tmp, tmp2, P::Tridiagonal)
     Base.require_one_based_indexing(du, P.dl, P.du)
     @assert length(du) == length(P.dl) + 1 == length(P.du) + 1
 

From 4af97d88290e7c0aec1786edd12906742e5513bf Mon Sep 17 00:00:00 2001
From: Stefan Kopecz <kopecz@mathematik.uni-kassel.de>
Date: Fri, 23 Aug 2024 17:49:20 +0200
Subject: [PATCH 24/42] format

---
 src/proddest.jl | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/src/proddest.jl b/src/proddest.jl
index 8fd5df13..7305cd86 100644
--- a/src/proddest.jl
+++ b/src/proddest.jl
@@ -153,7 +153,8 @@ end
 
 function PDSStdRHS(P, D, p_prototype, d_prototype)
     if p_prototype isa AbstractSparseMatrix
-        tmp = zeros(eltype(p_prototype), (size(p_prototype, 1),))/oneunit(first(p_prototype)) # drop units
+        tmp = zeros(eltype(p_prototype), (size(p_prototype, 1),)) /
+              oneunit(first(p_prototype)) # drop units
     else
         tmp = nothing
     end
@@ -335,7 +336,7 @@ end
 function ConservativePDSStdRHS(P, p_prototype)
     if p_prototype isa AbstractSparseMatrix
         tmp = zeros(eltype(p_prototype), (size(p_prototype, 1),))
-        tmp2 = tmp/oneunit(first(tmp)) # drop units
+        tmp2 = tmp / oneunit(first(tmp)) # drop units
     else
         tmp = nothing
         tmp2 = nothing

From bc6bfdbb405821cda2e8de496f782031d9699677 Mon Sep 17 00:00:00 2001
From: Stefan Kopecz <kopecz@mathematik.uni-kassel.de>
Date: Sun, 25 Aug 2024 21:21:28 +0200
Subject: [PATCH 25/42] Enable solution of PDS with unity by OrdinaryDiffEq
 solvers

---
 src/proddest.jl  |  8 +++---
 test/runtests.jl | 68 ++++++++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 72 insertions(+), 4 deletions(-)

diff --git a/src/proddest.jl b/src/proddest.jl
index 7305cd86..0cb88e7d 100644
--- a/src/proddest.jl
+++ b/src/proddest.jl
@@ -95,11 +95,11 @@ function PDSProblem{iip}(P, D, u0, tspan, p = NullParameters();
 
     # p_prototype is used to store evaluations of P, if P is in-place.
     if isnothing(p_prototype) && iip
-        p_prototype = zeros(eltype(u0), (length(u0), length(u0)))
+        p_prototype = zeros(eltype(u0), (length(u0), length(u0))) / oneunit(first(tspan))
     end
     # If a PDSFunction is to be evaluated and D is in-place, then d_prototype is used to store
     # evaluations of D.
-    d_prototype = similar(u0)
+    d_prototype = similar(u0) / oneunit(first(tspan))
 
     PD = PDSFunction{iip}(P, D; p_prototype, d_prototype,
                           analytic, std_rhs)
@@ -199,7 +199,7 @@ function (PD::PDSStdRHS)(du, u, p, t)
 end
 
 # New problem type ConservativePDSProblem
-"""
+"""ConservativePDSProblem
     ConservativePDSProblem(P, u0, tspan, p = NullParameters();
                            p_prototype = nothing,
                            analytic = nothing,
@@ -284,7 +284,7 @@ function ConservativePDSProblem{iip}(P, u0, tspan, p = NullParameters();
 
     # p_prototype is used to store evaluations of P, if P is in-place.
     if isnothing(p_prototype) && iip
-        p_prototype = zeros(eltype(u0), (length(u0), length(u0)))
+        p_prototype = zeros(eltype(u0), (length(u0), length(u0))) / oneunit(first(tspan))
     end
 
     PD = ConservativePDSFunction{iip}(P; p_prototype, analytic, std_rhs)
diff --git a/test/runtests.jl b/test/runtests.jl
index 559180bb..3beb6167 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -807,6 +807,74 @@ end
                           sol_pds_static ≈ sol_cpds_static_rhs ≈ sol_pds_static_rhs
                 end
             end
+            @testset "in-place" begin
+                u0 = [1.0, 1.5, 2.0, 2.5]u"N"
+                tspan = (0.0u"s", 1.0u"s")
+
+                function f!(du, u, p, t)
+                    fill!(du, zero(eltype(du)))
+
+                    du[1] = (u[1] - u[2]) / u"s"
+                    du[2] = (3 * u[2] - 2 * u[3] - u[1]) / u"s"
+                    du[3] = (-2 * u[2] + 5 * u[3] - 3 * u[4]) / u"s"
+                    du[4] = (3 * u[4] - 3 * u[3]) / u"s"
+
+                    return nothing
+                end
+                function P!(P, u, p, t)
+                    fill!(P, zero(eltype(P)))
+                    for i in 1:(length(u) - 1)
+                        P[i, i + 1] = i * u[i] / u"s"
+                        P[i + 1, i] = i * u[i + 1] / u"s"
+                    end
+                    return nothing
+                end
+                function D!(D, u, p, t)
+                    fill!(D, zero(eltype(D)))
+                    return nothing
+                end
+
+                P_tridiagonal = Tridiagonal([0.1, 0.2, 0.3], zeros(4),
+                                            [0.4, 0.5, 0.6])u"N/s"
+                P_dense = Matrix(P_tridiagonal)
+                P_sparse = sparse(P_tridiagonal)
+                D_dense = zeros(4)u"N/s"
+
+                alg = Euler()
+
+                prob0 = ODEProblem(f!, u0, tspan)
+                sol0 = solve(prob0, alg; dt = 0.2u"s")
+
+                probs = Array{Any}(undef, 14)
+                probs[1] = ConservativePDSProblem(P!, u0, tspan)
+                probs[2] = ConservativePDSProblem(P!, u0, tspan; p_prototype = P_dense)
+                probs[3] = ConservativePDSProblem(P!, u0, tspan; p_prototype = P_dense,
+                                                  std_rhs = f!)
+                probs[4] = ConservativePDSProblem(P!, u0, tspan;
+                                                  p_prototype = P_tridiagonal)
+                probs[5] = ConservativePDSProblem(P!, u0, tspan;
+                                                  p_prototype = P_tridiagonal, std_rhs = f!)
+                probs[6] = ConservativePDSProblem(P!, u0, tspan; p_prototype = P_sparse)
+                probs[7] = ConservativePDSProblem(P!, u0, tspan; p_prototype = P_sparse,
+                                                  std_rhs = f!)
+                probs[8] = PDSProblem(P!, D!, u0, tspan)
+                probs[9] = PDSProblem(P!, D!, u0, tspan; p_prototype = P_dense)
+                probs[10] = PDSProblem(P!, D!, u0, tspan; p_prototype = P_dense,
+                                       std_rhs = f!)
+                probs[11] = PDSProblem(P!, D!, u0, tspan; p_prototype = P_tridiagonal)
+                probs[12] = PDSProblem(P!, D!, u0, tspan; p_prototype = P_tridiagonal,
+                                       std_rhs = f!)
+                probs[13] = PDSProblem(P!, D!, u0, tspan; p_prototype = P_sparse)
+                probs[14] = PDSProblem(P!, D!, u0, tspan; p_prototype = P_sparse,
+                                       std_rhs = f!)
+
+                for prob in probs
+                    sol = solve(prob, alg; dt = 0.2u"s")
+
+                    @test sol0.t ≈ sol.t
+                    @test ustrip.(sol0.u) ≈ ustrip.(sol.u)
+                end
+            end
         end
 
         # Here we check that production-destruction form and standard ODE form fit together in predefinded problems,

From b6c35576c850c74883eb87d70839693e1a8d8019 Mon Sep 17 00:00:00 2001
From: Stefan Kopecz <kopecz@mathematik.uni-kassel.de>
Date: Tue, 27 Aug 2024 15:20:54 +0200
Subject: [PATCH 26/42] test more OrdinaryDiffEq algorithms solving problems
 with units

---
 test/runtests.jl | 18 +++++++++++-------
 1 file changed, 11 insertions(+), 7 deletions(-)

diff --git a/test/runtests.jl b/test/runtests.jl
index 3beb6167..87e6f5f0 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -787,7 +787,8 @@ end
                                                               std_rhs = f_static)
                 prob_pds_static_rhs = PDSProblem(P, D, u0_static, tspan; std_rhs = f_static)
 
-                algs = (Euler(), Tsit5())
+                algs = (Euler(), ImplicitEuler(), Tsit5(), Rosenbrock23(), SDIRK2(),
+                        TRBDF2())
 
                 for alg in algs
                     alg = Euler()
@@ -840,10 +841,10 @@ end
                 P_sparse = sparse(P_tridiagonal)
                 D_dense = zeros(4)u"N/s"
 
-                alg = Euler()
+                algs = (Euler(), ImplicitEuler(), Tsit5(), Rosenbrock23(), SDIRK2(),
+                        TRBDF2())
 
                 prob0 = ODEProblem(f!, u0, tspan)
-                sol0 = solve(prob0, alg; dt = 0.2u"s")
 
                 probs = Array{Any}(undef, 14)
                 probs[1] = ConservativePDSProblem(P!, u0, tspan)
@@ -868,11 +869,14 @@ end
                 probs[14] = PDSProblem(P!, D!, u0, tspan; p_prototype = P_sparse,
                                        std_rhs = f!)
 
-                for prob in probs
-                    sol = solve(prob, alg; dt = 0.2u"s")
+                for alg in algs
+                    sol0 = solve(prob0, alg; dt = 0.2u"s")
+                    for prob in probs
+                        sol = solve(prob, alg; dt = 0.2u"s")
 
-                    @test sol0.t ≈ sol.t
-                    @test ustrip.(sol0.u) ≈ ustrip.(sol.u)
+                        @test sol0.t ≈ sol.t
+                        @test ustrip.(sol0.u) ≈ ustrip.(sol.u)
+                    end
                 end
             end
         end

From 80ad9b9cdd33ce4158b773727480a8fa09a10718 Mon Sep 17 00:00:00 2001
From: Stefan Kopecz <kopecz@mathematik.uni-kassel.de>
Date: Tue, 27 Aug 2024 15:36:40 +0200
Subject: [PATCH 27/42] Float32 test

---
 test/runtests.jl | 17 +++++++++++++++++
 1 file changed, 17 insertions(+)

diff --git a/test/runtests.jl b/test/runtests.jl
index 87e6f5f0..9a07af89 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1067,6 +1067,23 @@ end
             end
         end
 
+        # We check that Float32 inputs result in Float32 outputs
+        # TODO: Use a tool such as LIKWID to check that no double-precision calculations are performed internally
+        @testset "Float32 input data" begin
+            P_linmod(u, p, t) = [0 u[2]; 5*u[1] 0]
+            u0 = [0.9f0, 0.1f0]
+            prob = ConservativePDSProblem(P_linmod, u0, (0.0f0, 2.0f0))
+
+            algs = (MPRK22(0.5f0), MPRK22(1.0f0), MPRK22(2.0f0), MPRK43I(1.0f0, 0.5f0),
+                    MPRK43I(0.5f0, 0.75f0), MPRK43II(0.5f0), MPRK43II(2.0f0 / 3.0f0),
+                    SSPMPRK22(0.5f0, 1.0f0))
+            for alg in algs
+                sol = solve(prob, alg; dt = 0.1f0, save_everystep = false)
+                @test sol.t isa Vector{Float32}
+                @test sol.u isa Vector{Vector{Float32}}
+            end
+        end
+
         # Here we check that MPE equals implicit Euler (IE) for a linear PDS
         @testset "Linear model problem: MPE = IE?" begin
             # problem data

From 98cc818acd5b58a0ac6edafcedad116a88db628b Mon Sep 17 00:00:00 2001
From: Stefan Kopecz <kopecz@mathematik.uni-kassel.de>
Date: Tue, 27 Aug 2024 16:05:40 +0200
Subject: [PATCH 28/42] removed tests for solving PDSProblems with units by
 implicit solvers

---
 test/runtests.jl | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/test/runtests.jl b/test/runtests.jl
index 9a07af89..230ea0d9 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -841,8 +841,9 @@ end
                 P_sparse = sparse(P_tridiagonal)
                 D_dense = zeros(4)u"N/s"
 
-                algs = (Euler(), ImplicitEuler(), Tsit5(), Rosenbrock23(), SDIRK2(),
-                        TRBDF2())
+                # implicit solvers and units don't work
+                #algs = (Euler(), ImplicitEuler(), Tsit5(), Rosenbrock23(), SDIRK2(), TRBDF2())
+                algs = (Euler(), Tsit5(), Vern9())
 
                 prob0 = ODEProblem(f!, u0, tspan)
 

From a7756c9fcfea380fe26ea537d4f006a63e39b5e9 Mon Sep 17 00:00:00 2001
From: Stefan Kopecz <kopecz@mathematik.uni-kassel.de>
Date: Tue, 27 Aug 2024 16:31:26 +0200
Subject: [PATCH 29/42] bugfix test Compatibility of OrdinarDiffEq solvers with
 (Conservative)PDSProblems and units

---
 test/runtests.jl | 1 -
 1 file changed, 1 deletion(-)

diff --git a/test/runtests.jl b/test/runtests.jl
index 230ea0d9..ce4b9632 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -791,7 +791,6 @@ end
                         TRBDF2())
 
                 for alg in algs
-                    alg = Euler()
                     sol_ode = solve(prob_ode, alg; dt = 0.1u"s")
                     sol_cpds = solve(prob_cpds, alg; dt = 0.1u"s")
                     sol_pds = solve(prob_pds, alg; dt = 0.1u"s")

From 16a513cd99798df88e0a94eaa234f53ee4a5da06 Mon Sep 17 00:00:00 2001
From: Stefan Kopecz <kopecz@mathematik.uni-kassel.de>
Date: Tue, 27 Aug 2024 16:52:48 +0200
Subject: [PATCH 30/42] Don't use implicit OrdinaryDiffEq solvers with units in
 tests

---
 test/runtests.jl | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/test/runtests.jl b/test/runtests.jl
index ce4b9632..983ab74a 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -787,8 +787,9 @@ end
                                                               std_rhs = f_static)
                 prob_pds_static_rhs = PDSProblem(P, D, u0_static, tspan; std_rhs = f_static)
 
-                algs = (Euler(), ImplicitEuler(), Tsit5(), Rosenbrock23(), SDIRK2(),
-                        TRBDF2())
+                # implicit solvers and units don't work
+                # algs = (Euler(), ImplicitEuler(), Tsit5(), Rosenbrock23(), SDIRK2(),TRBDF2())
+                algs = (Euler(), Tsit5(), Vern9())
 
                 for alg in algs
                     sol_ode = solve(prob_ode, alg; dt = 0.1u"s")

From c035e486567b0ce64eba3271714f2af7b2f732c5 Mon Sep 17 00:00:00 2001
From: Stefan Kopecz <SKopecz@users.noreply.github.com>
Date: Wed, 28 Aug 2024 11:58:52 +0200
Subject: [PATCH 31/42] Update src/mprk.jl

Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com>
---
 src/mprk.jl | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/mprk.jl b/src/mprk.jl
index 705aa277..131cfaaa 100644
--- a/src/mprk.jl
+++ b/src/mprk.jl
@@ -988,7 +988,8 @@ algorithm from [LinearSolve.jl](https://github.com/SciML/LinearSolve.jl)
 as keyword argument `linsolve`.
 You can also choose the parameter `small_constant` which is added to all Patankar-weight denominators
 to avoid divisions by zero. To display the default value for data type `type` evaluate 
-`PositiveIntegrators.small_constant_function_MPRK43II(type)`.
+`MPRK43II(gamma).small_constant_function(type)`, where `type` can be, e.g.,
+`Float64`.
 
 ## References
 

From 93b0d8a52eaadf379f30bbf079db5dfc5b806339 Mon Sep 17 00:00:00 2001
From: Stefan Kopecz <SKopecz@users.noreply.github.com>
Date: Wed, 28 Aug 2024 12:00:09 +0200
Subject: [PATCH 32/42] Update src/mprk.jl

Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com>
---
 src/mprk.jl | 4 +---
 1 file changed, 1 insertion(+), 3 deletions(-)

diff --git a/src/mprk.jl b/src/mprk.jl
index 131cfaaa..cd369eea 100644
--- a/src/mprk.jl
+++ b/src/mprk.jl
@@ -1016,9 +1016,7 @@ end
 
 function MPRK43II(gamma; linsolve = LUFactorization(),
                   small_constant = small_constant_function_MPRK43II)
-    if isnothing(small_constant)
-        small_constant_function = floatmin
-    elseif small_constant isa Number
+    if small_constant isa Number
         small_constant_function = Returns(small_constant)
     else # assume small_constant isa Function
         small_constant_function = small_constant

From e35ba2eeab2dc6d08b37bf5c7daf53b861e394f8 Mon Sep 17 00:00:00 2001
From: Stefan Kopecz <SKopecz@users.noreply.github.com>
Date: Wed, 28 Aug 2024 12:00:39 +0200
Subject: [PATCH 33/42] Update src/proddest.jl

Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com>
---
 src/proddest.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/proddest.jl b/src/proddest.jl
index 0cb88e7d..fb1bdfe0 100644
--- a/src/proddest.jl
+++ b/src/proddest.jl
@@ -199,7 +199,7 @@ function (PD::PDSStdRHS)(du, u, p, t)
 end
 
 # New problem type ConservativePDSProblem
-"""ConservativePDSProblem
+"""
     ConservativePDSProblem(P, u0, tspan, p = NullParameters();
                            p_prototype = nothing,
                            analytic = nothing,

From 713dc751ecdecb464436590e7dad9c50262f7e98 Mon Sep 17 00:00:00 2001
From: Stefan Kopecz <SKopecz@users.noreply.github.com>
Date: Wed, 28 Aug 2024 12:01:24 +0200
Subject: [PATCH 34/42] Update src/sspmprk.jl

Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com>
---
 src/sspmprk.jl | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/sspmprk.jl b/src/sspmprk.jl
index faa337a7..bd2dfdf5 100644
--- a/src/sspmprk.jl
+++ b/src/sspmprk.jl
@@ -467,7 +467,8 @@ algorithm from [LinearSolve.jl](https://github.com/SciML/LinearSolve.jl)
 as keyword argument `linsolve`.
 You can also choose the parameter `small_constant` which is added to all Patankar-weight denominators
 to avoid divisions by zero. To display the default value for data type `type` evaluate 
-`PositiveIntegrators.small_constant_function_SSPMPRK43(type)`.
+`SSPMPRK43. small_constant_function(type)`, where `type` can be, e.g.,
+`Float64`.
 
 The current implementation only supports fixed time steps.
 

From a5f123da2ec6098660f220429d155c672aa368a1 Mon Sep 17 00:00:00 2001
From: Stefan Kopecz <SKopecz@users.noreply.github.com>
Date: Wed, 28 Aug 2024 12:01:45 +0200
Subject: [PATCH 35/42] Update test/runtests.jl

Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com>
---
 test/runtests.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/runtests.jl b/test/runtests.jl
index 983ab74a..1a9328c6 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -5,7 +5,7 @@ using Statistics: mean
 
 using StaticArrays: MVector, @SVector, SA
 
-using Unitful
+using Unitful: @u_str
 
 using OrdinaryDiffEq
 using PositiveIntegrators

From a91cf205d13708a2a9208b6e9631a30889f41d11 Mon Sep 17 00:00:00 2001
From: Stefan Kopecz <kopecz@mathematik.uni-kassel.de>
Date: Wed, 28 Aug 2024 12:58:42 +0200
Subject: [PATCH 36/42] use unitful.jl's ustirp in tests

---
 test/runtests.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/runtests.jl b/test/runtests.jl
index 1a9328c6..b15a5a6e 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -5,7 +5,7 @@ using Statistics: mean
 
 using StaticArrays: MVector, @SVector, SA
 
-using Unitful: @u_str
+using Unitful: @u_str, ustrip
 
 using OrdinaryDiffEq
 using PositiveIntegrators

From 6fb3b4cc45ba05d8c910716468c30064b8eb2051 Mon Sep 17 00:00:00 2001
From: Stefan Kopecz <kopecz@mathematik.uni-kassel.de>
Date: Thu, 29 Aug 2024 10:38:11 +0200
Subject: [PATCH 37/42] Compare unitful data without ustrip

---
 test/runtests.jl | 11 +++++++++--
 1 file changed, 9 insertions(+), 2 deletions(-)

diff --git a/test/runtests.jl b/test/runtests.jl
index b15a5a6e..bb2d5234 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -875,8 +875,15 @@ end
                     for prob in probs
                         sol = solve(prob, alg; dt = 0.2u"s")
 
-                        @test sol0.t ≈ sol.t
-                        @test ustrip.(sol0.u) ≈ ustrip.(sol.u)
+                        #@test sol0.t ≈ sol.t
+                        # isapprox(sol0.u,sol.u) tries to check norm(sol0.u-sol.u) < 0.0,
+                        # which fails, because we have units on the left, but not on the right
+                        # In addition, abstol with units is not allowed.
+                        # Hence, we strip the units.
+                        # For Vector{Quantity{}} like sol.t there is a isapprox method in quantities.jl, which
+                        # takes care of the stripping.
+                        #@test ustrip.(sol0.u) ≈ ustrip.(sol.u)
+                        @test sol0 ≈ sol
                     end
                 end
             end

From eff9d3f44d8acb0ec1489a2cf2be6cdef8850d99 Mon Sep 17 00:00:00 2001
From: Stefan Kopecz <kopecz@mathematik.uni-kassel.de>
Date: Thu, 29 Aug 2024 11:19:57 +0200
Subject: [PATCH 38/42] Use ustrip before comparing solutions with units in
 tests

---
 test/runtests.jl | 35 ++++++++++++++++++++++++++++-------
 1 file changed, 28 insertions(+), 7 deletions(-)

diff --git a/test/runtests.jl b/test/runtests.jl
index bb2d5234..76a60978 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -803,9 +803,31 @@ end
                     sol_cpds_static_rhs = solve(prob_cpds_rhs, alg; dt = 0.1u"s")
                     sol_pds_static_rhs = solve(prob_pds_rhs, alg; dt = 0.1u"s")
 
-                    @test sol_ode ≈ sol_cpds ≈ sol_pds ≈ sol_cpds_rhs ≈ sol_pds_rhs ≈
-                          sol_ode_static ≈ sol_cpds_static ≈
-                          sol_pds_static ≈ sol_cpds_static_rhs ≈ sol_pds_static_rhs
+                    @test sol_ode.t ≈ sol_cpds.t
+                    @test sol_ode.t ≈ sol_pds.t
+                    @test sol_ode.t ≈ sol_cpds_rhs.t
+                    @test sol_ode.t ≈ sol_pds_rhs.t
+                    @test sol_ode.t ≈ sol_ode_static.t
+                    @test sol_ode.t ≈ sol_cpds_static.t
+                    @test sol_ode.t ≈ sol_pds_static.t
+                    @test sol_ode.t ≈ sol_cpds_static_rhs.t
+                    @test sol_ode.t ≈ sol_pds_static_rhs.t
+
+                    # isapprox(sol0.u,sol.u) tries to check norm(sol0.u-sol.u) < 0.0,
+                    # which fails, because we have units on the left, but not on the right
+                    # In addition, abstol with units is not allowed.
+                    # Hence, we strip the units.
+                    # For Vector{Quantity{}} like sol.t there is an isapprox method in quantities.jl, which
+                    # takes care of the stripping.
+                    @test ustrip.(sol_ode.u) ≈ ustrip.(sol_cpds.u)
+                    @test ustrip.(sol_ode.u) ≈ ustrip.(sol_pds.u)
+                    @test ustrip.(sol_ode.u) ≈ ustrip.(sol_cpds_rhs.u)
+                    @test ustrip.(sol_ode.u) ≈ ustrip.(sol_pds_rhs.u)
+                    @test ustrip.(sol_ode.u) ≈ ustrip.(sol_ode_static.u)
+                    @test ustrip.(sol_ode.u) ≈ ustrip.(sol_cpds_static.u)
+                    @test ustrip.(sol_ode.u) ≈ ustrip.(sol_pds_static.u)
+                    @test ustrip.(sol_ode.u) ≈ ustrip.(sol_cpds_static_rhs.u)
+                    @test ustrip.(sol_ode.u) ≈ ustrip.(sol_pds_static_rhs.u)
                 end
             end
             @testset "in-place" begin
@@ -875,15 +897,14 @@ end
                     for prob in probs
                         sol = solve(prob, alg; dt = 0.2u"s")
 
-                        #@test sol0.t ≈ sol.t
+                        @test sol0.t ≈ sol.t
                         # isapprox(sol0.u,sol.u) tries to check norm(sol0.u-sol.u) < 0.0,
                         # which fails, because we have units on the left, but not on the right
                         # In addition, abstol with units is not allowed.
                         # Hence, we strip the units.
-                        # For Vector{Quantity{}} like sol.t there is a isapprox method in quantities.jl, which
+                        # For Vector{Quantity{}} like sol.t there is an isapprox method in quantities.jl, which
                         # takes care of the stripping.
-                        #@test ustrip.(sol0.u) ≈ ustrip.(sol.u)
-                        @test sol0 ≈ sol
+                        @test ustrip.(sol0.u) ≈ ustrip.(sol.u)
                     end
                 end
             end

From 4842b5079b0601a1b8eb2d3f65cccb0d0a3c8581 Mon Sep 17 00:00:00 2001
From: Stefan Kopecz <kopecz@mathematik.uni-kassel.de>
Date: Thu, 29 Aug 2024 14:13:57 +0200
Subject: [PATCH 39/42] Added missing test (SSPMPRK43 and Float32)

---
 src/sspmprk.jl   | 4 +---
 test/runtests.jl | 2 +-
 2 files changed, 2 insertions(+), 4 deletions(-)

diff --git a/src/sspmprk.jl b/src/sspmprk.jl
index bd2dfdf5..447eb0d4 100644
--- a/src/sspmprk.jl
+++ b/src/sspmprk.jl
@@ -508,9 +508,7 @@ end
 
 function SSPMPRK43(; linsolve = LUFactorization(),
                    small_constant = small_constant_function_SSPMPRK43)
-    if isnothing(small_constant)
-        small_constant_function = floatmin
-    elseif small_constant isa Number
+    if small_constant isa Number
         small_constant_function = Returns(small_constant)
     else # assume small_constant isa Function
         small_constant_function = small_constant
diff --git a/test/runtests.jl b/test/runtests.jl
index 76a60978..3c0b1829 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1105,7 +1105,7 @@ end
 
             algs = (MPRK22(0.5f0), MPRK22(1.0f0), MPRK22(2.0f0), MPRK43I(1.0f0, 0.5f0),
                     MPRK43I(0.5f0, 0.75f0), MPRK43II(0.5f0), MPRK43II(2.0f0 / 3.0f0),
-                    SSPMPRK22(0.5f0, 1.0f0))
+                    SSPMPRK22(0.5f0, 1.0f0), SSPMPRK43())
             for alg in algs
                 sol = solve(prob, alg; dt = 0.1f0, save_everystep = false)
                 @test sol.t isa Vector{Float32}

From f4bdcb41701a97ebaefb008c98358e30e792935e Mon Sep 17 00:00:00 2001
From: Stefan Kopecz <kopecz@mathematik.uni-kassel.de>
Date: Thu, 29 Aug 2024 14:50:24 +0200
Subject: [PATCH 40/42] similar and oneunit

---
 src/proddest.jl | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/proddest.jl b/src/proddest.jl
index fb1bdfe0..8fefd61e 100644
--- a/src/proddest.jl
+++ b/src/proddest.jl
@@ -99,7 +99,7 @@ function PDSProblem{iip}(P, D, u0, tspan, p = NullParameters();
     end
     # If a PDSFunction is to be evaluated and D is in-place, then d_prototype is used to store
     # evaluations of D.
-    d_prototype = similar(u0) / oneunit(first(tspan))
+    d_prototype = similar(u0 ./ oneunit(first(tspan)))
 
     PD = PDSFunction{iip}(P, D; p_prototype, d_prototype,
                           analytic, std_rhs)
@@ -362,8 +362,8 @@ end
 function (PD::ConservativePDSStdRHS)(u::SVector, p, t)
     P = PD.p(u, p, t)
 
-    f = similar(P[:, 1]) #constructs MVector
-    zeroT = zero(eltype(P))
+    f = similar(u ./ oneunit(t)) #constructs MVector
+    zeroT = zero(eltype(f))
     for i in eachindex(f)
         f[i] = zeroT
     end

From cea2491b952d44977b74e381d45bb2e3fa37e9ed Mon Sep 17 00:00:00 2001
From: Stefan Kopecz <kopecz@mathematik.uni-kassel.de>
Date: Thu, 29 Aug 2024 18:31:04 +0200
Subject: [PATCH 41/42] similar(P[:,1])

---
 src/proddest.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/proddest.jl b/src/proddest.jl
index 8fefd61e..646f466f 100644
--- a/src/proddest.jl
+++ b/src/proddest.jl
@@ -362,8 +362,8 @@ end
 function (PD::ConservativePDSStdRHS)(u::SVector, p, t)
     P = PD.p(u, p, t)
 
-    f = similar(u ./ oneunit(t)) #constructs MVector
-    zeroT = zero(eltype(f))
+    f = similar(P[:, 1]) #constructs MVector
+    zeroT = zero(eltype(P))
     for i in eachindex(f)
         f[i] = zeroT
     end

From cba3acbde5fa8f9d63e82b229cf5c069747b9262 Mon Sep 17 00:00:00 2001
From: Hendrik Ranocha <ranocha@users.noreply.github.com>
Date: Fri, 30 Aug 2024 08:37:50 +0200
Subject: [PATCH 42/42] set version to v0.2.6

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index 510c40f5..48f7eec0 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "PositiveIntegrators"
 uuid = "d1b20bf0-b083-4985-a874-dc5121669aa5"
 authors = ["Stefan Kopecz, Hendrik Ranocha, and contributors"]
-version = "0.2.5"
+version = "0.2.6"
 
 [deps]
 FastBroadcast = "7034ab61-46d4-4ed7-9d0f-46aef9175898"