diff --git a/NEWS.md b/NEWS.md index f923a648..7d3a7db9 100644 --- a/NEWS.md +++ b/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` + + diff --git a/Project.toml b/Project.toml index fb03dd67..fdc71584 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.1.16" +version = "0.2.0" [deps] FastBroadcast = "7034ab61-46d4-4ed7-9d0f-46aef9175898" diff --git a/src/mprk.jl b/src/mprk.jl index fafda842..6e54e4b2 100644 --- a/src/mprk.jl +++ b/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, diff --git a/src/proddest.jl b/src/proddest.jl index 58e32c3a..0cdf7dd3 100644 --- a/src/proddest.jl +++ b/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,7 +81,6 @@ end # (arbitrary functions) function PDSProblem{iip}(P, D, u0, tspan, p = NullParameters(); p_prototype = nothing, - d_prototype = nothing, analytic = nothing, kwargs...) where {iip} @@ -94,10 +88,9 @@ function PDSProblem{iip}(P, D, u0, tspan, p = NullParameters(); 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)`. diff --git a/src/sspmprk.jl b/src/sspmprk.jl index 89dc362b..c25c79a3 100644 --- a/src/sspmprk.jl +++ b/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, diff --git a/test/runtests.jl b/test/runtests.jl index 1e13ea52..0148c6d6 100644 --- a/test/runtests.jl +++ b/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