From 8b1f5d1502256a9ce5fcdc82724257b033489960 Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Wed, 17 Jul 2024 21:54:56 +0200 Subject: [PATCH 1/5] Removed d_prototype --- src/mprk.jl | 12 ++++----- src/proddest.jl | 17 ++++-------- src/sspmprk.jl | 10 +++---- test/runtests.jl | 70 ++++++++++++++++++++++++++++++++++++++++++++++++ 4 files changed, 86 insertions(+), 23 deletions(-) diff --git a/src/mprk.jl b/src/mprk.jl index 9434cec6..dbc5b7f1 100644 --- a/src/mprk.jl +++ b/src/mprk.jl @@ -260,7 +260,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 @@ -567,8 +567,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) @@ -1107,9 +1107,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 01bb76c3..0cdf7dd3 100644 --- a/src/proddest.jl +++ b/src/proddest.jl @@ -4,7 +4,6 @@ abstract type AbstractPDSProblem end """ PDSProblem(P, D, u0, tspan, p = NullParameters(); p_prototype = nothing, - d_prototype = nothing, analytic = nothing) A structure describing a system of ordinary differential equations in form of a production-destruction system (PDS). @@ -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 bd871d33..10afd92e 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) @@ -724,9 +724,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 3b14f06a..48a6ad0c 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1193,6 +1193,76 @@ end end end + # Here we check that the types of p_prototype and d_prototype actually + # define the types of the Ps and Ds inside the algorithm caches. + # We test sparse, tridiagonal and dense matrices as well as sparse and + # dense vectors + @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 From fee912bc241904bf74f55010fbe3886be37f985a Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Wed, 17 Jul 2024 22:27:34 +0200 Subject: [PATCH 2/5] bug fix: removed d_prototype keywords in tests --- test/runtests.jl | 11 ++++------- 1 file changed, 4 insertions(+), 7 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index 48a6ad0c..fd2af5a3 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) @@ -1193,8 +1190,8 @@ end end end - # Here we check that the types of p_prototype and d_prototype actually - # define the types of the Ps and Ds inside the algorithm caches. + # 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 as well as sparse and # dense vectors @testset "Prototype type check" begin From d3da03d98b551f9d2c20c368880b3c00eb9e2294 Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Thu, 18 Jul 2024 15:22:32 +0200 Subject: [PATCH 3/5] Update test/runtests.jl Co-authored-by: Hendrik Ranocha --- test/runtests.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index fd2af5a3..4842bd65 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1192,8 +1192,7 @@ 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 as well as sparse and - # dense vectors + # We test sparse, tridiagonal, and dense matrices. @testset "Prototype type check" begin #prod and dest functions prod_inner! = (P, u, p, t) -> begin From 9a9b5637772f4bb4d3a68b76817a50a8549fda62 Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Thu, 18 Jul 2024 15:23:10 +0200 Subject: [PATCH 4/5] Set version to v0.2.0 --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 5bb53446..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.15" +version = "0.2.0" [deps] FastBroadcast = "7034ab61-46d4-4ed7-9d0f-46aef9175898" From 3452cc21a83e85c446a9fde9168a82d1929c1981 Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Thu, 18 Jul 2024 15:25:21 +0200 Subject: [PATCH 5/5] Updated News --- NEWS.md | 7 +++++++ 1 file changed, 7 insertions(+) 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` + +