Remove d_prototype & create type tests for p_prototype

NEWS.md

@@ -5,6 +5,11 @@ PositiveIntegrators.jl.jl follows the interpretation of
 used in the Julia ecosystem. Notable changes will be documented in this file
 for human readability.
 
+## Changes when updating to v0.2 from v0.1.x
+
+#### Removed
+
+- The optional keyword argument `d_prototype` has been removed from `PDSProblem`
 
 ## Changes in the v0.1 lifecycle
 
@@ -21,3 +26,5 @@ for human readability.
   `prob_pds_linmod`, `prob_pds_nonlinmod`, `prob_pds_npzd`, `prob_pds_robertson`,
   `prob_pds_sir`, `prob_pds_stratreac`
 - Modified Patankar methods `MPE` and `MPRK22`
+
+

Project.toml

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

src/mprk.jl

@@ -363,7 +363,7 @@ function alg_cache(alg::MPE, u, rate_prototype, ::Type{uEltypeNoUnits},
         linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true,
                         assumptions = LinearSolve.OperatorAssumptions(true))
 
-        MPECache(P, zero(u), σ, tab, linsolve_rhs, linsolve)
+        MPECache(P, similar(u), σ, tab, linsolve_rhs, linsolve)
     else
         throw(ArgumentError("MPE can only be applied to production-destruction systems"))
     end
@@ -670,8 +670,8 @@ function alg_cache(alg::MPRK22, u, rate_prototype, ::Type{uEltypeNoUnits},
                         assumptions = LinearSolve.OperatorAssumptions(true))
 
         MPRK22Cache(tmp, P, P2,
-                    zero(u), # D
-                    zero(u), # D2
+                    similar(u), # D
+                    similar(u), # D2
                     σ,
                     tab, #MPRK22ConstantCache
                     linsolve)
@@ -1246,9 +1246,9 @@ function alg_cache(alg::Union{MPRK43I, MPRK43II}, u, rate_prototype, ::Type{uElt
                         assumptions = LinearSolve.OperatorAssumptions(true))
         MPRK43ConservativeCache(tmp, tmp2, P, P2, P3, σ, tab, linsolve)
     elseif f isa PDSFunction
-        D = zero(u)
-        D2 = zero(u)
-        D3 = zero(u)
+        D = similar(u)
+        D2 = similar(u)
+        D3 = similar(u)
 
         linprob = LinearProblem(P3, _vec(tmp))
         linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true,

src/proddest.jl

@@ -3,9 +3,8 @@ abstract type AbstractPDSProblem end
 
 """
     PDSProblem(P, D, u0, tspan, p = NullParameters();
-               p_prototype = nothing,
-               d_prototype = nothing,
-               analytic = nothing)
+                       p_prototype = nothing,
+                       analytic = nothing)
 
 A structure describing a system of ordinary differential equations in form of a production-destruction system (PDS).
 `P` denotes the production matrix.
@@ -20,13 +19,9 @@ The functions `P` and `D` can be used either in the out-of-place form with signa
 
 ### Keyword arguments: ###
 
-- `p_prototype`: If `P` is given in in-place form, `p_prototype` is used to store evaluations of `P`.
+- `p_prototype`: If `P` is given in in-place form, `p_prototype` or copies thereof are used to store evaluations of `P`.
     If `p_prototype` is not specified explicitly and `P` is in-place, then `p_prototype` will be internally
   set to `zeros(eltype(u0), (length(u0), length(u0)))`.
-- `d_prototype`: If `D` is given in in-place form, `d_prototype` is used to store evaluations of `D`.
-  If `d_prototype` is not specified explicitly and `D` is in-place, then `d_prototype` will be internally
-set to `zeros(eltype(u0), (length(u0),))`.
-
 - `analytic`: The analytic solution of a PDS must be given in the form `f(u0,p,t)`.
     Specifying the analytic solution can be useful for plotting and convergence tests.
 
@@ -86,18 +81,16 @@ end
 # (arbitrary functions)
 function PDSProblem{iip}(P, D, u0, tspan, p = NullParameters();
                          p_prototype = nothing,
-                         d_prototype = nothing,
                          analytic = nothing,
                          kwargs...) where {iip}
 
     # 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)))
     end
-    # d_prototype is used to store evaluations of D, if D is in-place.
-    if isnothing(d_prototype) && iip
-        d_prototype = zeros(eltype(u0), (length(u0),))
-    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)
 
     PD = PDSFunction{iip}(P, D; p_prototype = p_prototype, d_prototype = d_prototype,
                           analytic = analytic)
@@ -177,7 +170,7 @@ The function `P` can be given either in the out-of-place form with signature
 
 ### Keyword arguments: ###
 
-- `p_prototype`: If `P` is given in in-place form, `p_prototype` is used to store evaluations of `P`.
+- `p_prototype`: If `P` is given in in-place form, `p_prototype` or copies thereof are used to store evaluations of `P`.
     If `p_prototype` is not specified explicitly and `P` is in-place, then `p_prototype` will be internally
   set to `zeros(eltype(u0), (length(u0), length(u0)))`.
 - `analytic`: The analytic solution of a PDS must be given in the form `f(u0,p,t)`.

src/sspmprk.jl

@@ -252,8 +252,8 @@ function alg_cache(alg::SSPMPRK22, u, rate_prototype, ::Type{uEltypeNoUnits},
                         assumptions = LinearSolve.OperatorAssumptions(true))
 
         SSPMPRK22Cache(tmp, P, P2,
-                       zero(u), # D
-                       zero(u), # D2
+                       similar(u), # D
+                       similar(u), # D2
                        σ,
                        tab, #MPRK22ConstantCache
                        linsolve)
@@ -760,9 +760,9 @@ function alg_cache(alg::SSPMPRK43, u, rate_prototype, ::Type{uEltypeNoUnits},
                         assumptions = LinearSolve.OperatorAssumptions(true))
         SSPMPRK43ConservativeCache(tmp, tmp2, P, P2, P3, σ, ρ, tab, linsolve)
     elseif f isa PDSFunction
-        D = zero(u)
-        D2 = zero(u)
-        D3 = zero(u)
+        D = similar(u)
+        D2 = similar(u)
+        D3 = similar(u)
 
         linprob = LinearProblem(P3, _vec(tmp))
         linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true,

test/runtests.jl

@@ -495,17 +495,14 @@ end
             # problem with sparse matrices
             p_prototype = spdiagm(-1 => ones(eltype(u0), N - 1),
                                   N - 1 => ones(eltype(u0), 1))
-            d_prototype = zero(u0)
             linear_advection_fd_upwind_PDS_sparse = PDSProblem(linear_advection_fd_upwind_P!,
                                                                linear_advection_fd_upwind_D!,
                                                                u0, tspan;
-                                                               p_prototype = p_prototype,
-                                                               d_prototype = d_prototype)
+                                                               p_prototype = p_prototype)
             linear_advection_fd_upwind_PDS_sparse_2 = PDSProblem{true}(linear_advection_fd_upwind_P!,
                                                                        linear_advection_fd_upwind_D!,
                                                                        u0, tspan;
-                                                                       p_prototype = p_prototype,
-                                                                       d_prototype = d_prototype)
+                                                                       p_prototype = p_prototype)
             linear_advection_fd_upwind_ConsPDS_sparse = ConservativePDSProblem(linear_advection_fd_upwind_P!,
                                                                                u0, tspan;
                                                                                p_prototype = p_prototype)
@@ -1201,6 +1198,75 @@ end
             end
         end
 
+        # Here we check that the type of p_prototype actually 
+        # defines the types of the Ps inside the algorithm caches.
+        # We test sparse, tridiagonal, and dense matrices.
+        @testset "Prototype type check" begin
+            #prod and dest functions
+            prod_inner! = (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_sparse! = (P, u, p, t) -> begin
+                @test P isa SparseMatrixCSC
+                prod_inner!(P, u, p, t)
+                return nothing
+            end
+            prod_tridiagonal! = (P, u, p, t) -> begin
+                @test P isa Tridiagonal
+                prod_inner!(P, u, p, t)
+                return nothing
+            end
+            prod_dense! = (P, u, p, t) -> begin
+                @test P isa Matrix
+                prod_inner!(P, u, p, t)
+                return nothing
+            end
+            dest! = (D, u, p, t) -> begin
+                fill!(D, zero(eltype(D)))
+            end
+            #prototypes
+            P_tridiagonal = Tridiagonal([0.1, 0.2, 0.3],
+                                        [0.0, 0.0, 0.0, 0.0],
+                                        [0.4, 0.5, 0.6])
+            P_dense = Matrix(P_tridiagonal)
+            P_sparse = sparse(P_tridiagonal)
+            # problem definition
+            u0 = [1.0, 1.5, 2.0, 2.5]
+            tspan = (0.0, 1.0)
+            dt = 0.5
+            ## conservative PDS
+            prob_default = ConservativePDSProblem(prod_dense!, u0, tspan)
+            prob_tridiagonal = ConservativePDSProblem(prod_tridiagonal!, u0, tspan;
+                                                      p_prototype = P_tridiagonal)
+            prob_dense = ConservativePDSProblem(prod_dense!, u0, tspan;
+                                                p_prototype = P_dense)
+            prob_sparse = ConservativePDSProblem(prod_sparse!, u0, tspan;
+                                                 p_prototype = P_sparse)
+            ## nonconservative PDS                                         
+            prob_default2 = PDSProblem(prod_dense!, dest!, u0, tspan)
+            prob_tridiagonal2 = PDSProblem(prod_tridiagonal!, dest!, u0, tspan;
+                                           p_prototype = P_tridiagonal)
+            prob_dense2 = PDSProblem(prod_dense!, dest!, u0, tspan;
+                                     p_prototype = P_dense)
+            prob_sparse2 = PDSProblem(prod_sparse!, dest!, u0, tspan;
+                                      p_prototype = P_sparse)
+            #solve and test
+            for alg in (MPE(), MPRK22(0.5), MPRK22(1.0), 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 prob in (prob_default, prob_tridiagonal, prob_dense, prob_sparse,
+                             prob_default2,
+                             prob_tridiagonal2, prob_dense2, prob_sparse2)
+                    solve(prob, alg; dt, adaptive = false)
+                end
+            end
+        end
+
         # Here we check the convergence order of pth-order schemes for which
         # also an interpolation of order p is available
         @testset "Convergence tests (conservative)" begin