From fff0edbb6c413819f6d792644de3d79a0a4e9244 Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Wed, 26 Jun 2024 11:47:29 +0200 Subject: [PATCH 01/48] Added SSPMPRK22. --- src/PositiveIntegrators.jl | 4 ++++ src/mprk.jl | 22 ++++++++++++++-------- 2 files changed, 18 insertions(+), 8 deletions(-) diff --git a/src/PositiveIntegrators.jl b/src/PositiveIntegrators.jl index 8929959c..64d953ea 100644 --- a/src/PositiveIntegrators.jl +++ b/src/PositiveIntegrators.jl @@ -48,6 +48,7 @@ export PDSFunction, PDSProblem export ConservativePDSFunction, ConservativePDSProblem export MPE, MPRK22, MPRK43I, MPRK43II +export SSPMPRK22 export prob_pds_linmod, prob_pds_linmod_inplace, prob_pds_nonlinmod, prob_pds_robertson, prob_pds_brusselator, prob_pds_sir, @@ -61,6 +62,9 @@ include("proddest.jl") # modified Patankar-Runge-Kutta (MPRK) methods include("mprk.jl") +# modified Patankar-Runge-Kutta based on the SSP formulation of RK methods (SSPMPRK) +include("sspmprk.jl") + # predefined PDS problems include("PDSProblemLibrary.jl") diff --git a/src/mprk.jl b/src/mprk.jl index 8a9e7ca3..5662a9f9 100644 --- a/src/mprk.jl +++ b/src/mprk.jl @@ -553,9 +553,9 @@ function perform_step!(integrator, cache::MPRK22ConstantCache, repeat_step = fal u = sol.u integrator.stats.nsolve += 1 - # copied from perform_step for HeunConstantCache - # If a21 = 1.0, then σ is the MPE approximation and thus suited for stiff problems. - # If a21 ≠ 1.0, σ might be a bad choice to estimate errors. + # If a21 = 1.0, then σ is the MPE approximation, i.e. a first order approximation + # of the solution, and can be used for error estimation. Moreover, MPE is suited for stiff problems. + # TODO: Find first order approximation if a21≠ 1.0. tmp = u - σ atmp = calculate_residuals(tmp, uprev, u, integrator.opts.abstol, integrator.opts.reltol, integrator.opts.internalnorm, t) @@ -693,6 +693,9 @@ function perform_step!(integrator, cache::MPRK22Cache, repeat_step = false) integrator.stats.nsolve += 1 # Now σ stores the error estimate + # If a21 = 1.0, then σ is the MPE approximation, i.e. a first order approximation + # of the solution, and can be used for error estimation. Moreover, MPE is suited for stiff problems. + # TODO: Find first order approximation if a21≠ 1.0. @.. broadcast=false σ=u - σ # Now tmp stores error residuals @@ -750,6 +753,9 @@ function perform_step!(integrator, cache::MPRK22ConservativeCache, repeat_step = integrator.stats.nsolve += 1 # Now σ stores the error estimate + # If a21 = 1.0, then σ is the MPE approximation, i.e. a first order approximation + # of the solution, and can be used for error estimation. Moreover, MPE is suited for stiff problems. + # TODO: Find first order approximation if a21≠ 1.0. @.. broadcast=false σ=u - σ # Now tmp stores error residuals @@ -1343,13 +1349,13 @@ end function interp_summary(::Type{cacheType}, dense::Bool) where { cacheType <: - Union{MPRK22ConstantCache, MPRK22Cache, + Union{MPRK22ConstantCache, MPRK22Cache,SSPMPRK22ConstantCache, SSPMPRK22Cache, MPRK43ConstantCache, MPRK43Cache}} "1st order linear" end function _ode_interpolant(Θ, dt, u0, u1, k, - cache::Union{MPRK22ConstantCache, MPRK22Cache, + cache::Union{MPRK22ConstantCache, MPRK22Cache,SSPMPRK22ConstantCache, SSPMPRK22Cache, MPRK43ConstantCache, MPRK43Cache}, idxs, # Optionally specialize for ::Nothing and others T::Type{Val{0}}, @@ -1358,7 +1364,7 @@ function _ode_interpolant(Θ, dt, u0, u1, k, end function _ode_interpolant!(out, Θ, dt, u0, u1, k, - cache::Union{MPRK22ConstantCache, MPRK22Cache, + cache::Union{MPRK22ConstantCache, MPRK22Cache,SSPMPRK22ConstantCache, SSPMPRK22Cache, MPRK43ConstantCache, MPRK43Cache}, idxs, # Optionally specialize for ::Nothing and others T::Type{Val{0}}, @@ -1367,7 +1373,7 @@ function _ode_interpolant!(out, Θ, dt, u0, u1, k, end function _ode_interpolant(Θ, dt, u0, u1, k, - cache::Union{MPRK22ConstantCache, MPRK22Cache, + cache::Union{MPRK22ConstantCache, MPRK22Cache,SSPMPRK22ConstantCache, SSPMPRK22Cache, MPRK43ConstantCache, MPRK43Cache}, idxs, # Optionally specialize for ::Nothing and others T::Type{Val{1}}, @@ -1376,7 +1382,7 @@ function _ode_interpolant(Θ, dt, u0, u1, k, end function _ode_interpolant!(out, Θ, dt, u0, u1, k, - cache::Union{MPRK22ConstantCache, MPRK22Cache, + cache::Union{MPRK22ConstantCache, MPRK22Cache,SSPMPRK22ConstantCache, SSPMPRK22Cache, MPRK43ConstantCache, MPRK43Cache}, idxs, # Optionally specialize for ::Nothing and others T::Type{Val{1}}, From 91bacdae97af05a470fd5b7aa08d1cbb5e7de01a Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Wed, 26 Jun 2024 11:53:13 +0200 Subject: [PATCH 02/48] format --- src/mprk.jl | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/src/mprk.jl b/src/mprk.jl index 5662a9f9..784d7e4d 100644 --- a/src/mprk.jl +++ b/src/mprk.jl @@ -1349,13 +1349,15 @@ end function interp_summary(::Type{cacheType}, dense::Bool) where { cacheType <: - Union{MPRK22ConstantCache, MPRK22Cache,SSPMPRK22ConstantCache, SSPMPRK22Cache, + Union{MPRK22ConstantCache, MPRK22Cache, + SSPMPRK22ConstantCache, SSPMPRK22Cache, MPRK43ConstantCache, MPRK43Cache}} "1st order linear" end function _ode_interpolant(Θ, dt, u0, u1, k, - cache::Union{MPRK22ConstantCache, MPRK22Cache,SSPMPRK22ConstantCache, SSPMPRK22Cache, + cache::Union{MPRK22ConstantCache, MPRK22Cache, + SSPMPRK22ConstantCache, SSPMPRK22Cache, MPRK43ConstantCache, MPRK43Cache}, idxs, # Optionally specialize for ::Nothing and others T::Type{Val{0}}, @@ -1364,7 +1366,8 @@ function _ode_interpolant(Θ, dt, u0, u1, k, end function _ode_interpolant!(out, Θ, dt, u0, u1, k, - cache::Union{MPRK22ConstantCache, MPRK22Cache,SSPMPRK22ConstantCache, SSPMPRK22Cache, + cache::Union{MPRK22ConstantCache, MPRK22Cache, + SSPMPRK22ConstantCache, SSPMPRK22Cache, MPRK43ConstantCache, MPRK43Cache}, idxs, # Optionally specialize for ::Nothing and others T::Type{Val{0}}, @@ -1373,7 +1376,8 @@ function _ode_interpolant!(out, Θ, dt, u0, u1, k, end function _ode_interpolant(Θ, dt, u0, u1, k, - cache::Union{MPRK22ConstantCache, MPRK22Cache,SSPMPRK22ConstantCache, SSPMPRK22Cache, + cache::Union{MPRK22ConstantCache, MPRK22Cache, + SSPMPRK22ConstantCache, SSPMPRK22Cache, MPRK43ConstantCache, MPRK43Cache}, idxs, # Optionally specialize for ::Nothing and others T::Type{Val{1}}, @@ -1382,7 +1386,8 @@ function _ode_interpolant(Θ, dt, u0, u1, k, end function _ode_interpolant!(out, Θ, dt, u0, u1, k, - cache::Union{MPRK22ConstantCache, MPRK22Cache,SSPMPRK22ConstantCache, SSPMPRK22Cache, + cache::Union{MPRK22ConstantCache, MPRK22Cache, + SSPMPRK22ConstantCache, SSPMPRK22Cache, MPRK43ConstantCache, MPRK43Cache}, idxs, # Optionally specialize for ::Nothing and others T::Type{Val{1}}, From cd6c5b726cf7abb9c10b9d5ed32cb72f1924d9b3 Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Wed, 26 Jun 2024 13:07:59 +0200 Subject: [PATCH 03/48] Added sspmprk.jl file --- src/sspmprk.jl | 372 +++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 372 insertions(+) create mode 100644 src/sspmprk.jl diff --git a/src/sspmprk.jl b/src/sspmprk.jl new file mode 100644 index 00000000..a030a68a --- /dev/null +++ b/src/sspmprk.jl @@ -0,0 +1,372 @@ +### SSPMPRK ##################################################################################### +""" + SSPMPRK22(α, β; [linsolve = ...]) + +A family of second-order modified Patankar-Runge-Kutta algorithms for +production-destruction systems. Each member of this family is an one-step, two-stage method which is +second-order accurate, unconditionally positivity-preserving, and linearly +implicit. The parameters `α` and `β` are described by Huang and Shu (2019) and +studied by Huang, Izgin, Kopecz, Meister and Shu (2023). +The difference to [`MPRK22`](@ref) is that this method is based on the SSP formulation of +an explicit second-order Runge-Kutta method. This family of schemes contains the [`MPRK22`](@ref) family, +where `MPRK22(α) = SSMPRK22(0, α)` applies. + +The scheme was introduced by Huang and Shu for conservative production-destruction systems. +For nonconservative production–destruction systems we use the straight forward extension +analogous to [`MPE`](@ref). + +This modified Patankar-Runge-Kutta method requires the special structure of a +[`PDSProblem`](@ref) or a [`ConservativePDSProblem`](@ref). + +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`. + +## References + +- Juntao Huang and Chi-Wang Shu. + "Positivity-Preserving Time Discretizations for Production–Destruction Equations + with Applications to Non-equilibrium Flows." + Journal of Scientific Computing 78 (2019): 1811–1839 + [DOI: 10.1007/s10915-018-0852-1](https://doi.org/10.1007/s10915-018-0852-1) +- Juntao Huang, Thomas Izgin, Stefan Kopecz, Andreas Meister and Chi-Wang Shu. + "On the stability of strong-stability-preserving modified Patankar-Runge-Kutta schemes." + ESAIM: Mathematical Modelling and Numerical Analysis 57 (2023):1063–1086 + [DOI: 10.1051/m2an/2023005](https://doi.org/10.1051/m2an/2023005) +""" +struct SSPMPRK22{T, F} <: OrdinaryDiffEqAdaptiveAlgorithm + alpha::T + beta::T + linsolve::F +end + +function SSPMPRK22(alpha, beta; linsolve = LUFactorization()) + SSPMPRK22{typeof(alpha), typeof(linsolve)}(alpha, beta, linsolve) +end + +alg_order(::SSPMPRK22) = 2 +isfsal(::SSPMPRK22) = false + +function get_constant_parameters(alg::SSPMPRK22) + if !((0 ≤ alg.alpha ≤ 1) && (alg.beta ≥ 0) && + (alg.alpha * alg.beta + 1 / (2 * alg.beta) ≤ 1)) + throw(ArgumentError("SSPMPRK22 requires 0 ≤ α ≤ 1, β ≥ 0 and αβ + 1/(2β) ≤ 1.")) + end + + a21 = alg.alpha + a10 = one(alg.alpha) + a20 = 1 - a21 + + b10 = alg.beta + b20 = 1 - 1 / (2 * b10) - a21 * b10 + b21 = 1 / (2 * b10) + + s = (b20 + b21 + a21 * b10^2) / (b10 * (b20 + b21)) + + # This should never happen + if !all((a21, a10, a20, b10, b20, b21) .≥ 0) + throw(ArgumentError("SSPMPRK22 requires nonnegative SSP coefficients.")) + end + return a21, a10, a20, b10, b20, b21, s +end + +struct SSPMPRK22ConstantCache{T} <: OrdinaryDiffEqConstantCache + a21::T + a10::T + a20::T + b10::T + b20::T + b21::T + s::T + small_constant::T +end + +# Out-of-place +function alg_cache(alg::SSPMPRK22, u, rate_prototype, ::Type{uEltypeNoUnits}, + ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, + uprev, uprev2, f, t, dt, reltol, p, calck, + ::Val{false}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} + if !(f isa PDSFunction || f isa ConservativePDSFunction) + throw(ArgumentError("SSPMPRK22 can only be applied to production-destruction systems")) + end + + a21, a10, a20, b10, b20, b21, s = get_constant_parameters(alg) + SSPMPRK22ConstantCache(a21, a10, a20, b10, b20, b21, s, floatmin(uEltypeNoUnits)) +end + +function initialize!(integrator, cache::SSPMPRK22ConstantCache) +end + +function perform_step!(integrator, cache::SSPMPRK22ConstantCache, repeat_step = false) + @unpack alg, t, dt, uprev, f, p = integrator + @unpack a21, a10, a20, b10, b20, b21, s, small_constant = cache + + f = integrator.f + + # evaluate production matrix + P = f.p(uprev, p, t) + Ptmp = b10 * P + integrator.stats.nf += 1 + + # avoid division by zero due to zero Patankar weights + σ = add_small_constant(uprev, small_constant) + + # build linear system matrix and rhs + if f isa PDSFunction + d = f.d(uprev, p, t) # evaluate nonconservative destruction terms + dtmp = b10 * d + rhs = a10 * uprev + dt * diag(Ptmp) + M = build_mprk_matrix(Ptmp, σ, dt, dtmp) + else + # f isa ConservativePDSFunction + M = build_mprk_matrix(Ptmp, σ, dt) + rhs = a10 * uprev + end + + # solve linear system + linprob = LinearProblem(M, rhs) + sol = solve(linprob, alg.linsolve) + u = sol.u + integrator.stats.nsolve += 1 + + # compute Patankar weight denominator + if isone(s) + σ = u + else + σ = σ .^ (1 - s) .* u .^ s + end + # avoid division by zero due to zero Patankar weights + σ = add_small_constant(σ, small_constant) + + P2 = f.p(u, p, t + a21 * dt) + Ptmp = b20 * P + b21 * P2 + integrator.stats.nf += 1 + + # build linear system matrix and rhs + if f isa PDSFunction + d2 = f.d(u, p, t + a21 * dt) # evaluate nonconservative destruction terms + dtmp = b20 * d + b21 * d2 + rhs = a20 * uprev + a21 * u + dt * diag(Ptmp) + M = build_mprk_matrix(Ptmp, σ, dt, dtmp) + else + # f isa ConservativePDSFunction + M = build_mprk_matrix(Ptmp, σ, dt) + rhs = a20 * uprev + a21 * u + end + + # solve linear system + linprob = LinearProblem(M, rhs) + sol = solve(linprob, alg.linsolve) + u = sol.u + integrator.stats.nsolve += 1 + + # If a21 = 0 or b10 = 0, then σ is a first order approximation of the solution and + # can be used for error estimation. + # TODO: Find first order approximation, if a21*b10 ≠ 0. + tmp = u - σ + atmp = calculate_residuals(tmp, uprev, u, integrator.opts.abstol, + integrator.opts.reltol, integrator.opts.internalnorm, t) + integrator.EEst = integrator.opts.internalnorm(atmp, t) + + integrator.u = u +end + +struct SSPMPRK22Cache{uType, PType, tabType, F} <: + OrdinaryDiffEqMutableCache + tmp::uType + P::PType + P2::PType + D::uType + D2::uType + σ::uType + tab::tabType + linsolve::F +end + +struct SSPMPRK22ConservativeCache{uType, PType, tabType, F} <: + OrdinaryDiffEqMutableCache + tmp::uType + P::PType + P2::PType + σ::uType + tab::tabType + linsolve::F +end + +# In-place +function alg_cache(alg::SSPMPRK22, u, rate_prototype, ::Type{uEltypeNoUnits}, + ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, + uprev, uprev2, f, t, dt, reltol, p, calck, + ::Val{true}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} + a21, a10, a20, b10, b20, b21, s = get_constant_parameters(alg) + tab = SSPMPRK22ConstantCache(a21, a10, a20, b10, b20, b21, s, floatmin(uEltypeNoUnits)) + tmp = zero(u) + P = p_prototype(u, f) + # We use P2 to store the last evaluation of the PDS + # as well as to store the system matrix of the linear system + P2 = p_prototype(u, f) + σ = zero(u) + + if f isa ConservativePDSFunction + linprob = LinearProblem(P2, _vec(tmp)) + linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true, + assumptions = LinearSolve.OperatorAssumptions(true)) + + SSPMPRK22ConservativeCache(tmp, P, P2, σ, + tab, #MPRK22ConstantCache + linsolve) + elseif f isa PDSFunction + linprob = LinearProblem(P2, _vec(tmp)) + linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true, + assumptions = LinearSolve.OperatorAssumptions(true)) + + SSPMPRK22Cache(tmp, P, P2, + zero(u), # D + zero(u), # D2 + σ, + tab, #MPRK22ConstantCache + linsolve) + else + throw(ArgumentError("MPRK22 can only be applied to production-destruction systems")) + end +end + +function initialize!(integrator, cache::Union{SSPMPRK22Cache, SSPMPRK22ConservativeCache}) +end + +function perform_step!(integrator, cache::SSPMPRK22Cache, repeat_step = false) + @unpack t, dt, uprev, u, f, p = integrator + @unpack tmp, P, P2, D, D2, σ, linsolve = cache + @unpack a21, a10, a20, b10, b20, b21, s, small_constant = cache.tab + + # We use P2 to store the last evaluation of the PDS + # as well as to store the system matrix of the linear system + + f.p(P, uprev, p, t) # evaluate production terms + f.d(D, uprev, p, t) # evaluate nonconservative destruction terms + integrator.stats.nf += 1 + @.. broadcast=false P2=b10 * P + @.. broadcast=false D2=b10 * D + + # avoid division by zero due to zero Patankar weights + @.. broadcast=false σ=uprev + small_constant + + # tmp holds the right hand side of the linear system + @.. broadcast=false tmp=a10 * uprev + @inbounds for i in eachindex(tmp) + tmp[i] += dt * P2[i, i] + end + + build_mprk_matrix!(P2, P2, σ, dt, D2) + + # Same as linres = P2 \ tmp + linsolve.A = P2 + linres = solve!(linsolve) + + u .= linres + integrator.stats.nsolve += 1 + + if isone(s) + σ .= u + else + @.. broadcast=false σ=σ^(1 - s) * u^s + end + @.. broadcast=false σ=σ + small_constant + + f.p(P2, u, p, t + a21 * dt) # evaluate production terms + f.d(D2, u, p, t + a21 * dt) # evaluate nonconservative destruction terms + integrator.stats.nf += 1 + + @.. broadcast=false P2=b20 * P + b21 * P2 + @.. broadcast=false D2=b20 * D + b21 * D2 + + # tmp holds the right hand side of the linear system + @.. broadcast=false tmp=a20 * uprev + a21 * u + @inbounds for i in eachindex(tmp) + tmp[i] += dt * P2[i, i] + end + + build_mprk_matrix!(P2, P2, σ, dt, D2) + + # Same as linres = P2 \ tmp + linsolve.A = P2 + linres = solve!(linsolve) + + u .= linres + integrator.stats.nsolve += 1 + + # If a21 = 0 or b10 = 0, then σ is a first order approximation of the solution and + # can be used for error estimation. + # TODO: Find first order approximation, if a21*b10 ≠ 0. + + # Now σ stores the error estimate + @.. broadcast=false σ=u - σ + + # Now tmp stores error residuals + calculate_residuals!(tmp, σ, uprev, u, integrator.opts.abstol, + integrator.opts.reltol, integrator.opts.internalnorm, t, + False()) + integrator.EEst = integrator.opts.internalnorm(tmp, t) +end + +function perform_step!(integrator, cache::SSPMPRK22ConservativeCache, repeat_step = false) + @unpack t, dt, uprev, u, f, p = integrator + @unpack tmp, P, P2, σ, linsolve = cache + @unpack a21, a10, a20, b10, b20, b21, s, small_constant = cache.tab + + # We use P2 to store the last evaluation of the PDS + # as well as to store the system matrix of the linear system + f.p(P, uprev, p, t) # evaluate production terms + integrator.stats.nf += 1 + @.. broadcast=false P2=b10 * P + + # Avoid division by zero due to zero Patankar weights + @.. broadcast=false σ=uprev + small_constant + + # tmp holds the right hand side of the linear system + @.. broadcast=false tmp=a10 * uprev + + build_mprk_matrix!(P2, P2, σ, dt) + + # Same as linres = P2 \ tmp + linsolve.A = P2 + linres = solve!(linsolve) + + u .= linres + integrator.stats.nsolve += 1 + + if isone(s) + σ .= u + else + @.. broadcast=false σ=σ^(1 - s) * u^s + end + @.. broadcast=false σ=σ + small_constant + + f.p(P2, u, p, t + a21 * dt) # evaluate production terms + integrator.stats.nf += 1 + + @.. broadcast=false P2=b20 * P + b21 * P2 + + @.. broadcast=false tmp=a20 * uprev + a21 * u + + build_mprk_matrix!(P2, P2, σ, dt) + + # Same as linres = P2 \ tmp + linsolve.A = P2 + linres = solve!(linsolve) + + u .= linres + integrator.stats.nsolve += 1 + + # If a21 = 0 or b10 = 0, then σ is a first order approximation of the solution and + # can be used for error estimation. + # TODO: Find first order approximation, if a21*b10 ≠ 0. + # Now σ stores the error estimate + @.. broadcast=false σ=u - σ + + # Now tmp stores error residuals + calculate_residuals!(tmp, σ, uprev, u, integrator.opts.abstol, + integrator.opts.reltol, integrator.opts.internalnorm, t, + False()) + integrator.EEst = integrator.opts.internalnorm(tmp, t) +end From b1c91ea705d4c1f29d59837b9f298fbfb5c4e45b Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Wed, 26 Jun 2024 13:22:37 +0200 Subject: [PATCH 04/48] Include sspmprk.jl before mprk.jl --- src/PositiveIntegrators.jl | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/src/PositiveIntegrators.jl b/src/PositiveIntegrators.jl index 64d953ea..0ed50f58 100644 --- a/src/PositiveIntegrators.jl +++ b/src/PositiveIntegrators.jl @@ -59,12 +59,16 @@ export prob_pds_linmod, prob_pds_linmod_inplace, prob_pds_nonlinmod, # production-destruction systems include("proddest.jl") -# modified Patankar-Runge-Kutta (MPRK) methods -include("mprk.jl") - +#TODO: sspmprk.jl must be included before mprk.jl since the interpolation specializations +# at the end of mprk.jl need to know SSPMPRK22ConstantCache. +# It would probably be a good idea to separate schemes and interpolations. +# # modified Patankar-Runge-Kutta based on the SSP formulation of RK methods (SSPMPRK) include("sspmprk.jl") +# modified Patankar-Runge-Kutta (MPRK) methods +include("mprk.jl") + # predefined PDS problems include("PDSProblemLibrary.jl") From 15ab1fccb72cb06bcd110aee649303a7ff890c16 Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Wed, 26 Jun 2024 14:33:21 +0200 Subject: [PATCH 05/48] Added tests for SSPMPRK22 --- src/mprk.jl | 2 ++ test/runtests.jl | 33 +++++++++++++++++++++++++-------- 2 files changed, 27 insertions(+), 8 deletions(-) diff --git a/src/mprk.jl b/src/mprk.jl index 784d7e4d..51110c1a 100644 --- a/src/mprk.jl +++ b/src/mprk.jl @@ -1345,6 +1345,8 @@ end ######################################################################################## + +####################################################################################### # interpolation specializations function interp_summary(::Type{cacheType}, dense::Bool) where { diff --git a/test/runtests.jl b/test/runtests.jl index bd9bd703..a7ae8f56 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -324,6 +324,9 @@ const prob_pds_linmod_nonconservative_inplace = PDSProblem(linmodP!, linmodD!, [ MPRK43II(1.0)) @test_throws "MPRK43II requires 3/8 ≤ γ ≤ 3/4." solve(prob_pds_linmod, MPRK43II(0.0)) + @test_throws "SSPMPRK22 requires 0 ≤ α ≤ 1, β ≥ 0 and αβ + 1/(2β) ≤ 1." solve(prob_pds_linmod, SSPMPRK22(-1.0, 1.0)) + @test_throws "SSPMPRK22 requires 0 ≤ α ≤ 1, β ≥ 0 and αβ + 1/(2β) ≤ 1." solve(prob_pds_linmod, SSPMPRK22(0.0, -1.0)) + @test_throws "SSPMPRK22 requires 0 ≤ α ≤ 1, β ≥ 0 and αβ + 1/(2β) ≤ 1." solve(prob_pds_linmod, SSPMPRK22(1.0, 10.0)) end # Here we check that MPE equals implicit Euler (IE) for a linear PDS @@ -380,6 +383,17 @@ const prob_pds_linmod_nonconservative_inplace = PDSProblem(linmodP!, linmodD!, [ @test sol_MPE_ip.u ≈ sol_MPE_ip_2.u ≈ sol_IE_ip.u end + # Here we check that MPRK22(α) = SSPMPRK22(0,α) + @testset "MPRK22(α) = SSPMPRK22(0, α)" begin + for α in (0.5, 2.0/3.0, 1.0, 2.0) + sol1 = solve(prob_pds_linmod, MPRK22(α)) + sol2 = solve(prob_pds_linmod, SSPMPRK22(0.0,α)) + sol3 = solve(prob_pds_linmod_inplace, MPRK22(α)) + sol4 = solve(prob_pds_linmod_inplace, SSPMPRK22(0.0,α)) + @test sol1.u ≈ sol2.u ≈ sol3.u ≈ sol4.u + end + end + # Here we check that different linear solvers can be used @testset "Different linear solvers" begin # problem data @@ -417,7 +431,8 @@ const prob_pds_linmod_nonconservative_inplace = PDSProblem(linmodP!, linmodD!, [ (; kwargs...) -> MPRK43I(1.0, 0.5; kwargs...), (; kwargs...) -> MPRK43I(0.5, 0.75; kwargs...), (; kwargs...) -> MPRK43II(0.5; kwargs...), - (; kwargs...) -> MPRK43II(2.0 / 3.0; kwargs...)) + (; kwargs...) -> MPRK43II(2.0 / 3.0; kwargs...), + (; kwargs...) -> SSPMPRK22(0.5, 1.0; kwargs...)) for alg in algs # Check different linear solvers @@ -478,7 +493,7 @@ const prob_pds_linmod_nonconservative_inplace = PDSProblem(linmodP!, linmodD!, [ @testset "$alg" 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)) + MPRK43II(2.0 / 3.0), MPRK43II(0.5), SSPMPRK22(0.5, 1.0)) for prod! in (prod_1!, prod_2!, prod_3!) prod = (u, p, t) -> begin P = similar(u, (length(u), length(u))) @@ -594,7 +609,7 @@ const prob_pds_linmod_nonconservative_inplace = PDSProblem(linmodP!, linmodD!, [ @testset "$alg" 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)) + MPRK43II(2.0 / 3.0), MPRK43II(0.5),SSPMPRK22(0.5, 1.0)) for (prod!, dest!) in zip((prod_1!, prod_2!, prod_3!), (dest_1!, dest_2!, dest_3!)) prod = (u, p, t) -> begin @@ -651,7 +666,7 @@ const prob_pds_linmod_nonconservative_inplace = PDSProblem(linmodP!, linmodD!, [ # 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 - algs = (MPE(), MPRK22(0.5), MPRK22(1.0), MPRK22(2.0)) + algs = (MPE(), MPRK22(0.5), MPRK22(1.0), MPRK22(2.0),SSPMPRK22(0.5, 1.0)) dts = 0.5 .^ (6:11) problems = (prob_pds_linmod, prob_pds_linmod_array, prob_pds_linmod_mvector, prob_pds_linmod_inplace) @@ -677,7 +692,7 @@ const prob_pds_linmod_nonconservative_inplace = PDSProblem(linmodP!, linmodD!, [ # Here we check the convergence order of pth-order schemes for which # also an interpolation of order p is available @testset "Convergence tests (nonconservative)" begin - algs = (MPE(), MPRK22(0.5), MPRK22(1.0), MPRK22(2.0)) + algs = (MPE(), MPRK22(0.5), MPRK22(1.0), MPRK22(2.0),SSPMPRK22(0.5, 1.0)) dts = 0.5 .^ (6:11) problems = (prob_pds_linmod_nonconservative, prob_pds_linmod_nonconservative_inplace) @@ -748,7 +763,7 @@ const prob_pds_linmod_nonconservative_inplace = PDSProblem(linmodP!, linmodD!, [ @testset "Interpolation tests (nonconservative)" begin algs = (MPE(), 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)) + MPRK43I(0.5, 0.75), MPRK43II(0.5), MPRK43II(2.0 / 3.0),SSPMPRK22(0.5, 1.0)) dt = 0.5^6 problems = (prob_pds_linmod_nonconservative, prob_pds_linmod_nonconservative_inplace) @@ -795,7 +810,7 @@ const prob_pds_linmod_nonconservative_inplace = PDSProblem(linmodP!, linmodD!, [ algs = (MPE(), 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)) + MPRK43II(2.0 / 3.0),SSPMPRK22(0.5, 1.0)) for alg in algs sol = solve(prob_ip, alg; dt = dt, adaptive = false) @@ -843,7 +858,7 @@ const prob_pds_linmod_nonconservative_inplace = PDSProblem(linmodP!, linmodD!, [ algs = (MPE(), 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)) + MPRK43II(2.0 / 3.0),SSPMPRK22(0.5, 1.0)) dt = 1e-3 for alg in algs @@ -860,6 +875,7 @@ const prob_pds_linmod_nonconservative_inplace = PDSProblem(linmodP!, linmodD!, [ end end + #= # TODO: Do we want to keep the examples and test them or do we want # to switch to real docs/tutorials instead? @testset "Examples" begin @@ -883,4 +899,5 @@ const prob_pds_linmod_nonconservative_inplace = PDSProblem(linmodP!, linmodD!, [ end end end + =# end; From a5e9ffc1bebee0947180d88e857e5fc5a4ef94eb Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Wed, 26 Jun 2024 14:35:41 +0200 Subject: [PATCH 06/48] format --- src/mprk.jl | 1 - test/runtests.jl | 35 ++++++++++++++++++++++------------- 2 files changed, 22 insertions(+), 14 deletions(-) diff --git a/src/mprk.jl b/src/mprk.jl index 51110c1a..a09b46b0 100644 --- a/src/mprk.jl +++ b/src/mprk.jl @@ -1345,7 +1345,6 @@ end ######################################################################################## - ####################################################################################### # interpolation specializations function interp_summary(::Type{cacheType}, diff --git a/test/runtests.jl b/test/runtests.jl index a7ae8f56..e4f8dbb7 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -324,9 +324,15 @@ const prob_pds_linmod_nonconservative_inplace = PDSProblem(linmodP!, linmodD!, [ MPRK43II(1.0)) @test_throws "MPRK43II requires 3/8 ≤ γ ≤ 3/4." solve(prob_pds_linmod, MPRK43II(0.0)) - @test_throws "SSPMPRK22 requires 0 ≤ α ≤ 1, β ≥ 0 and αβ + 1/(2β) ≤ 1." solve(prob_pds_linmod, SSPMPRK22(-1.0, 1.0)) - @test_throws "SSPMPRK22 requires 0 ≤ α ≤ 1, β ≥ 0 and αβ + 1/(2β) ≤ 1." solve(prob_pds_linmod, SSPMPRK22(0.0, -1.0)) - @test_throws "SSPMPRK22 requires 0 ≤ α ≤ 1, β ≥ 0 and αβ + 1/(2β) ≤ 1." solve(prob_pds_linmod, SSPMPRK22(1.0, 10.0)) + @test_throws "SSPMPRK22 requires 0 ≤ α ≤ 1, β ≥ 0 and αβ + 1/(2β) ≤ 1." solve(prob_pds_linmod, + SSPMPRK22(-1.0, + 1.0)) + @test_throws "SSPMPRK22 requires 0 ≤ α ≤ 1, β ≥ 0 and αβ + 1/(2β) ≤ 1." solve(prob_pds_linmod, + SSPMPRK22(0.0, + -1.0)) + @test_throws "SSPMPRK22 requires 0 ≤ α ≤ 1, β ≥ 0 and αβ + 1/(2β) ≤ 1." solve(prob_pds_linmod, + SSPMPRK22(1.0, + 10.0)) end # Here we check that MPE equals implicit Euler (IE) for a linear PDS @@ -385,11 +391,11 @@ const prob_pds_linmod_nonconservative_inplace = PDSProblem(linmodP!, linmodD!, [ # Here we check that MPRK22(α) = SSPMPRK22(0,α) @testset "MPRK22(α) = SSPMPRK22(0, α)" begin - for α in (0.5, 2.0/3.0, 1.0, 2.0) + for α in (0.5, 2.0 / 3.0, 1.0, 2.0) sol1 = solve(prob_pds_linmod, MPRK22(α)) - sol2 = solve(prob_pds_linmod, SSPMPRK22(0.0,α)) + sol2 = solve(prob_pds_linmod, SSPMPRK22(0.0, α)) sol3 = solve(prob_pds_linmod_inplace, MPRK22(α)) - sol4 = solve(prob_pds_linmod_inplace, SSPMPRK22(0.0,α)) + sol4 = solve(prob_pds_linmod_inplace, SSPMPRK22(0.0, α)) @test sol1.u ≈ sol2.u ≈ sol3.u ≈ sol4.u end end @@ -493,7 +499,8 @@ const prob_pds_linmod_nonconservative_inplace = PDSProblem(linmodP!, linmodD!, [ @testset "$alg" 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)) + MPRK43II(2.0 / 3.0), MPRK43II(0.5), + SSPMPRK22(0.5, 1.0)) for prod! in (prod_1!, prod_2!, prod_3!) prod = (u, p, t) -> begin P = similar(u, (length(u), length(u))) @@ -609,7 +616,8 @@ const prob_pds_linmod_nonconservative_inplace = PDSProblem(linmodP!, linmodD!, [ @testset "$alg" 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)) + MPRK43II(2.0 / 3.0), MPRK43II(0.5), + SSPMPRK22(0.5, 1.0)) for (prod!, dest!) in zip((prod_1!, prod_2!, prod_3!), (dest_1!, dest_2!, dest_3!)) prod = (u, p, t) -> begin @@ -666,7 +674,7 @@ const prob_pds_linmod_nonconservative_inplace = PDSProblem(linmodP!, linmodD!, [ # 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 - algs = (MPE(), MPRK22(0.5), MPRK22(1.0), MPRK22(2.0),SSPMPRK22(0.5, 1.0)) + algs = (MPE(), MPRK22(0.5), MPRK22(1.0), MPRK22(2.0), SSPMPRK22(0.5, 1.0)) dts = 0.5 .^ (6:11) problems = (prob_pds_linmod, prob_pds_linmod_array, prob_pds_linmod_mvector, prob_pds_linmod_inplace) @@ -692,7 +700,7 @@ const prob_pds_linmod_nonconservative_inplace = PDSProblem(linmodP!, linmodD!, [ # Here we check the convergence order of pth-order schemes for which # also an interpolation of order p is available @testset "Convergence tests (nonconservative)" begin - algs = (MPE(), MPRK22(0.5), MPRK22(1.0), MPRK22(2.0),SSPMPRK22(0.5, 1.0)) + algs = (MPE(), MPRK22(0.5), MPRK22(1.0), MPRK22(2.0), SSPMPRK22(0.5, 1.0)) dts = 0.5 .^ (6:11) problems = (prob_pds_linmod_nonconservative, prob_pds_linmod_nonconservative_inplace) @@ -763,7 +771,8 @@ const prob_pds_linmod_nonconservative_inplace = PDSProblem(linmodP!, linmodD!, [ @testset "Interpolation tests (nonconservative)" begin algs = (MPE(), 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)) + MPRK43I(0.5, 0.75), MPRK43II(0.5), MPRK43II(2.0 / 3.0), + SSPMPRK22(0.5, 1.0)) dt = 0.5^6 problems = (prob_pds_linmod_nonconservative, prob_pds_linmod_nonconservative_inplace) @@ -810,7 +819,7 @@ const prob_pds_linmod_nonconservative_inplace = PDSProblem(linmodP!, linmodD!, [ algs = (MPE(), 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)) + MPRK43II(2.0 / 3.0), SSPMPRK22(0.5, 1.0)) for alg in algs sol = solve(prob_ip, alg; dt = dt, adaptive = false) @@ -858,7 +867,7 @@ const prob_pds_linmod_nonconservative_inplace = PDSProblem(linmodP!, linmodD!, [ algs = (MPE(), 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)) + MPRK43II(2.0 / 3.0), SSPMPRK22(0.5, 1.0)) dt = 1e-3 for alg in algs From d179e0180007e0dcbac2ebd035d66026cb8601dc Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Wed, 26 Jun 2024 15:23:26 +0200 Subject: [PATCH 07/48] Additional tests --- test/runtests.jl | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index e4f8dbb7..7f98bfff 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -324,6 +324,10 @@ const prob_pds_linmod_nonconservative_inplace = PDSProblem(linmodP!, linmodD!, [ MPRK43II(1.0)) @test_throws "MPRK43II requires 3/8 ≤ γ ≤ 3/4." solve(prob_pds_linmod, MPRK43II(0.0)) + @test_throws "SSPMPRK22 can only be applied to production-destruction systems" solve(prob_oop, + SSPMPRK22(0.5,1.0)) + @test_throws "SSPMPRK22 can only be applied to production-destruction systems" solve(prob_ip, + SSPMPRK22(0.5,1.0)) @test_throws "SSPMPRK22 requires 0 ≤ α ≤ 1, β ≥ 0 and αβ + 1/(2β) ≤ 1." solve(prob_pds_linmod, SSPMPRK22(-1.0, 1.0)) @@ -884,7 +888,7 @@ const prob_pds_linmod_nonconservative_inplace = PDSProblem(linmodP!, linmodD!, [ end end - #= + # TODO: Do we want to keep the examples and test them or do we want # to switch to real docs/tutorials instead? @testset "Examples" begin @@ -908,5 +912,5 @@ const prob_pds_linmod_nonconservative_inplace = PDSProblem(linmodP!, linmodD!, [ end end end - =# + end; From 5a2de58d20f13d1b3911b1b5fea451ad2f97f77a Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Wed, 26 Jun 2024 15:25:42 +0200 Subject: [PATCH 08/48] format --- test/runtests.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index 7f98bfff..8e575459 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -325,9 +325,11 @@ const prob_pds_linmod_nonconservative_inplace = PDSProblem(linmodP!, linmodD!, [ @test_throws "MPRK43II requires 3/8 ≤ γ ≤ 3/4." solve(prob_pds_linmod, MPRK43II(0.0)) @test_throws "SSPMPRK22 can only be applied to production-destruction systems" solve(prob_oop, - SSPMPRK22(0.5,1.0)) + SSPMPRK22(0.5, + 1.0)) @test_throws "SSPMPRK22 can only be applied to production-destruction systems" solve(prob_ip, - SSPMPRK22(0.5,1.0)) + SSPMPRK22(0.5, + 1.0)) @test_throws "SSPMPRK22 requires 0 ≤ α ≤ 1, β ≥ 0 and αβ + 1/(2β) ≤ 1." solve(prob_pds_linmod, SSPMPRK22(-1.0, 1.0)) @@ -888,7 +890,6 @@ const prob_pds_linmod_nonconservative_inplace = PDSProblem(linmodP!, linmodD!, [ end end - # TODO: Do we want to keep the examples and test them or do we want # to switch to real docs/tutorials instead? @testset "Examples" begin @@ -912,5 +913,4 @@ const prob_pds_linmod_nonconservative_inplace = PDSProblem(linmodP!, linmodD!, [ end end end - end; From ae7dc07954f0f0c205bffca8502621a44209fd90 Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Wed, 26 Jun 2024 16:16:07 +0200 Subject: [PATCH 09/48] bug fix: wrong error message --- src/sspmprk.jl | 512 ++++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 510 insertions(+), 2 deletions(-) diff --git a/src/sspmprk.jl b/src/sspmprk.jl index a030a68a..7b626d75 100644 --- a/src/sspmprk.jl +++ b/src/sspmprk.jl @@ -3,7 +3,7 @@ SSPMPRK22(α, β; [linsolve = ...]) A family of second-order modified Patankar-Runge-Kutta algorithms for -production-destruction systems. Each member of this family is an one-step, two-stage method which is +production-destruction systems. Each member of this family is a one-step, two-stage method which is second-order accurate, unconditionally positivity-preserving, and linearly implicit. The parameters `α` and `β` are described by Huang and Shu (2019) and studied by Huang, Izgin, Kopecz, Meister and Shu (2023). @@ -227,7 +227,7 @@ function alg_cache(alg::SSPMPRK22, u, rate_prototype, ::Type{uEltypeNoUnits}, tab, #MPRK22ConstantCache linsolve) else - throw(ArgumentError("MPRK22 can only be applied to production-destruction systems")) + throw(ArgumentError("SSPMPRK22 can only be applied to production-destruction systems")) end end @@ -370,3 +370,511 @@ function perform_step!(integrator, cache::SSPMPRK22ConservativeCache, repeat_ste False()) integrator.EEst = integrator.opts.internalnorm(tmp, t) end + +""" + SSPMPRK43([linsolve = ...]) + +A third-order modified Patankar-Runge-Kutta algorithm for +production-destruction systems. This scheme is a one-step, two-stage method which is +third-order accurate, unconditionally positivity-preserving, and linearly +implicit. The scheme is described by Huang, Zhao and Shu (2019) and +studied by Huang, Izgin, Kopecz, Meister and Shu (2023). +The difference to [`MPRK43I`](@ref) or [`MPRK43II`](@ref) is that this method is based on the SSP formulation of +an explicit third-order Runge-Kutta method. + +The scheme was introduced by Huang, Zhao and Shu for conservative production-destruction systems. +For nonconservative production–destruction systems we use the straight forward extension +analogous to [`MPE`](@ref). + +This modified Patankar-Runge-Kutta method requires the special structure of a +[`PDSProblem`](@ref) or a [`ConservativePDSProblem`](@ref). + +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`. + +## References + +- Juntao Huang, Weifeng Zhao and Chi-Wang Shu. + "A Third-Order Unconditionally Positivity-Preserving Scheme for + Production–Destruction Equations with Applications to Non-equilibrium Flows." + Journal of Scientific Computing 79 (2019): 1015–1056 + [DOI: 10.1007/s10915-018-0881-9](https://doi.org/10.1007/s10915-018-0881-9) +- Juntao Huang, Thomas Izgin, Stefan Kopecz, Andreas Meister and Chi-Wang Shu. + "On the stability of strong-stability-preserving modified Patankar-Runge-Kutta schemes." + ESAIM: Mathematical Modelling and Numerical Analysis 57 (2023):1063–1086 + [DOI: 10.1051/m2an/2023005](https://doi.org/10.1051/m2an/2023005) +""" +struct SSPMPRK43{F} <: OrdinaryDiffEqAdaptiveAlgorithm + linsolve::F +end + +function SSPMPRK43(; linsolve = LUFactorization()) + SSPMPRK43{typeof(linsolve)}(linsolve) +end + +alg_order(::SSPMPRK43) = 3 +isfsal(::SSPMPRK43) = false + +function get_constant_parameters(alg::SSPMPRK43) + n1 = 2.569046025732011E-01 + n2 = 7.430953974267989E-01 + z = 6.288938077828750E-01 + η1 = 3.777285888379173E-02 + η2 = 1.0 / 3.0 + η3 = 1.868649805549811E-01 + η4 = 2.224876040351123 + s = 5.721964308755304 + + α10 = 1.0 + α20 = 9.2600312554031827E-01 + α21 = 7.3996874459681783E-02 + α30 = 7.0439040373427619E-01 + α31 = 2.0662904223744017E-10 + α32 = 2.9560959605909481E-01 + β10 = 4.7620819268131703E-01 + β20 = 7.7545442722396801E-02 + β21 = 5.9197500149679749E-01 + β30 = 2.0044747790361456E-01 + β31 = 6.8214380786704851E-10 + β32 = 5.9121918658514827E-01 + + c3 = β20 + α21 * β10 + β21 + + return n1, n2, z, η1, η2, η3, η4, s, α10, α20, α21, α30, α31, α32, β10, β20, β21, β30, + β31, β32, c3 +end + +struct SSPMPRK43ConstantCache{T} <: OrdinaryDiffEqConstantCache + n1::T + n2::T + z::T + η1::T + η2::T + η3::T + η4::T + s::T + α10::T + α20::T + α21::T + α30::T + α31::T + α32::T + β10::T + β20::T + β21::T + β30::T + β31::T + β32::T + c3::T + small_constant::T +end + +# Out-of-place +function alg_cache(alg::SSPMPRK43, u, rate_prototype, ::Type{uEltypeNoUnits}, + ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, + uprev, uprev2, f, t, dt, reltol, p, calck, + ::Val{false}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} + 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, s, α10, α20, α21, α30, α31, α32, β10, β20, β21, β30, β31, β32, c3 = get_constant_parameters(alg) + SSPMPRK43ConstantCache(n1, n2, z, η1, η2, η3, η4, s, α10, α20, α21, α30, α31, α32, β10, + β20, β21, β30, + β31, β32, c3, floatmin(uEltypeNoUnits)) +end + +function initialize!(integrator, cache::SSPMPRK43ConstantCache) +end + +function perform_step!(integrator, cache::SSPMPRK43ConstantCache, repeat_step = false) + @unpack alg, t, dt, uprev, f, p = integrator + @unpack n1, n2, z, η1, η2, η3, η4, s, α10, α20, α21, α30, α31, α32, β10, β20, β21, β30, β31, β32, c3, small_constant = cache + + f = integrator.f + + # evaluate production matrix + P = f.p(uprev, p, t) + Ptmp = β10 * P + integrator.stats.nf += 1 + + # avoid division by zero due to zero Patankar weights + σ = add_small_constant(uprev, small_constant) + + # build linear system matrix and rhs + if f isa PDSFunction + d = f.d(uprev, p, t) + dtmp = β10 * d + rhs = α10 * uprev + dt * diag(Ptmp) + M = build_mprk_matrix(Ptmp, σ, dt, dtmp) + else + rhs = α10 * uprev + M = build_mprk_matrix(Ptmp, σ, dt) + end + + # solve linear system + linprob = LinearProblem(M, rhs) + sol = solve(linprob, alg.linsolve) + u2 = sol.u + u = u2 + integrator.stats.nsolve += 1 + + # compute Patankar weight denominator + ρ = n1 * u + n2 * u .^ 2 ./ σ + # avoid division by zero due to zero Patankar weights + ρ = add_small_constant(ρ, small_constant) + + P2 = f.p(u, p, t + β10 * dt) + Ptmp = β20 * P + β21 * P2 + integrator.stats.nf += 1 + + # build linear system matrix and rhs + if f isa PDSFunction + d2 = f.d(u, p, t + β10 * dt) # evaluate nonconservative destruction terms + dtmp = β20 * d + β21 * d2 + rhs = α20 * uprev + α21 * u2 + dt * diag(Ptmp) + M = build_mprk_matrix(Ptmp, ρ, dt, dtmp) + + else + rhs = α20 * uprev + α21 * u2 + M = build_mprk_matrix(Ptmp, ρ, dt) + end + + # solve linear system + linprob = LinearProblem(M, rhs) + sol = solve(linprob, alg.linsolve) + u = sol.u + integrator.stats.nsolve += 1 + + # compute Patankar weight denominator + σ = σ .^ (1 - s) .* u2 .^ s + # avoid division by zero due to zero Patankar weights + σ = add_small_constant(σ, small_constant) + + Ptmp = η3 * P + η4 * P2 + + # build linear system matrix and rhs + if f isa PDSFunction + dtmp = η3 * d + η4 * d2 + rhs = η1 * uprev + η2 * u2 + dt * diag(Ptmp) + M = build_mprk_matrix(Ptmp, σ, dt, dtmp) + else + rhs = η1 * uprev + η2 * u2 + M = build_mprk_matrix(Ptmp, σ, dt) + end + + # solve linear system + linprob = LinearProblem(M, rhs) + sol = solve(linprob, alg.linsolve) + σ = sol.u + integrator.stats.nsolve += 1 + + # compute Patankar weight denominator + σ = σ + z .* uprev .* u ./ rho + # avoid division by zero due to zero Patankar weights + σ = add_small_constant(σ, small_constant) + + P3 = f.p(u, p, t + c3 * dt) + Ptmp = β30 * P + β31 * P2 + β32 * P3 + integrator.stats.nf += 1 + + # build linear system matrix + if f isa PDSFunction + d3 = f.d(u, p, t + c3 * dt) # evaluate nonconservative destruction terms + dtmp = β30 * d + β31 * d2 + β30 * d3 + rhs = α30 * uprev + α31 * u2 + α32 * u + dt * diag(Ptmp) + M = build_mprk_matrix(Ptmp, σ, dt, dtmp) + else + rhs = α30 * uprev + α31 * u2 + α32 * u + M = build_mprk_matrix(Ptmp, σ, dt) + end + + # solve linear system + linprob = LinearProblem(M, rhs) + sol = solve(linprob, alg.linsolve) + u = sol.u + integrator.stats.nsolve += 1 + + #TODO: Figure out if a second order approximation of the solution + # is hidden somewhere. + #= + tmp = u - σ + atmp = calculate_residuals(tmp, uprev, u, integrator.opts.abstol, + integrator.opts.reltol, integrator.opts.internalnorm, t) + integrator.EEst = integrator.opts.internalnorm(atmp, t) + =# + + integrator.u = u +end + +#= +struct MPRK43Cache{uType, PType, tabType, F} <: OrdinaryDiffEqMutableCache + tmp::uType + tmp2::uType + P::PType + P2::PType + P3::PType + D::uType + D2::uType + D3::uType + σ::uType + tab::tabType + linsolve::F +end + +struct MPRK43ConservativeCache{uType, PType, tabType, F} <: OrdinaryDiffEqMutableCache + tmp::uType + tmp2::uType + P::PType + P2::PType + P3::PType + σ::uType + tab::tabType + linsolve::F +end + +# In-place +function alg_cache(alg::Union{MPRK43I, MPRK43II}, u, rate_prototype, ::Type{uEltypeNoUnits}, + ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, + uprev, uprev2, f, t, dt, reltol, p, calck, + ::Val{true}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} + a21, a31, a32, b1, b2, b3, c2, c3, beta1, beta2, q1, q2 = get_constant_parameters(alg) + tab = MPRK43ConstantCache(a21, a31, a32, b1, b2, b3, c2, c3, + beta1, beta2, q1, q2, floatmin(uEltypeNoUnits)) + tmp = zero(u) + tmp2 = zero(u) + P = p_prototype(u, f) + P2 = p_prototype(u, f) + # We use P3 to store the last evaluation of the PDS + # as well as to store the system matrix of the linear system + P3 = p_prototype(u, f) + σ = zero(u) + + if f isa ConservativePDSFunction + # The right hand side of the linear system is always uprev. But using + # tmp instead of uprev for the rhs we allow `alias_b=true`. uprev must + # not be altered, since it is needed to compute the adaptive time step + # size. + linprob = LinearProblem(P3, _vec(tmp)) + linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true, + 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) + + linprob = LinearProblem(P3, _vec(tmp)) + linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true, + assumptions = LinearSolve.OperatorAssumptions(true)) + + MPRK43Cache(tmp, tmp2, P, P2, P3, D, D2, D3, σ, tab, linsolve) + else + throw(ArgumentError("MPRK43 can only be applied to production-destruction systems")) + end +end + +function initialize!(integrator, cache::Union{MPRK43Cache, MPRK43ConservativeCache}) +end + +function perform_step!(integrator, cache::MPRK43Cache, repeat_step = false) + @unpack t, dt, uprev, u, f, p = integrator + @unpack tmp, tmp2, P, P2, P3, D, D2, D3, σ, linsolve = cache + @unpack a21, a31, a32, b1, b2, b3, c2, c3, beta1, beta2, q1, q2, small_constant = cache.tab + + # We use P3 to store the last evaluation of the PDS + # as well as to store the system matrix of the linear system + + f.p(P, uprev, p, t) # evaluate production terms + f.d(D, uprev, p, t) # evaluate nonconservative destruction terms + @.. broadcast=false P3=a21 * P + @.. broadcast=false D3=a21 * D + integrator.stats.nf += 1 + + # avoid division by zero due to zero Patankar weights + @.. broadcast=false σ=uprev + small_constant + + # tmp holds the right hand side of the linear system + tmp .= uprev + @inbounds for i in eachindex(tmp) + tmp[i] += dt * P3[i, i] + end + + build_mprk_matrix!(P3, P3, σ, dt, D3) + + # Same as linres = P3 \ tmp + linsolve.A = P3 + linres = solve!(linsolve) + + u .= linres + if !(q1 ≈ q2) + tmp2 .= u #u2 in out-of-place version + end + integrator.stats.nsolve += 1 + + @.. broadcast=false σ=σ^(1 - q1) * u^q1 + @.. broadcast=false σ=σ + small_constant + + f.p(P2, u, p, t + c2 * dt) # evaluate production terms + f.d(D2, u, p, t + c2 * dt) # evaluate nonconservative destruction terms + @.. broadcast=false P3=a31 * P + a32 * P2 + @.. broadcast=false D3=a31 * D + a32 * D2 + integrator.stats.nf += 1 + + # tmp holds the right hand side of the linear system + tmp .= uprev + @inbounds for i in eachindex(tmp) + tmp[i] += dt * P3[i, i] + end + + build_mprk_matrix!(P3, P3, σ, dt, D3) + + # Same as linres = P3 \ tmp + linsolve.A = P3 + linres = solve!(linsolve) + + u .= linres + integrator.stats.nsolve += 1 + + if !(q1 ≈ q2) + @.. broadcast=false σ=(uprev + small_constant)^(1 - q2) * tmp2^q2 + @.. broadcast=false σ=σ + small_constant + end + + @.. broadcast=false P3=beta1 * P + beta2 * P2 + @.. broadcast=false D3=beta1 * D + beta2 * D2 + + # tmp holds the right hand side of the linear system + tmp .= uprev + @inbounds for i in eachindex(tmp) + tmp[i] += dt * P3[i, i] + end + + build_mprk_matrix!(P3, P3, σ, dt, D3) + + # Same as linres = P3 \ tmp + linsolve.A = P3 + linres = solve!(linsolve) + integrator.stats.nsolve += 1 + + σ .= linres + # avoid division by zero due to zero Patankar weights + @.. broadcast=false σ=σ + small_constant + + f.p(P3, u, p, t + c3 * dt) # evaluate production terms + f.d(D3, u, p, t + c3 * dt) # evaluate nonconservative destruction terms + @.. broadcast=false P3=b1 * P + b2 * P2 + b3 * P3 + @.. broadcast=false D3=b1 * D + b2 * D2 + b3 * D3 + integrator.stats.nf += 1 + + # tmp holds the right hand side of the linear system + tmp .= uprev + @inbounds for i in eachindex(tmp) + tmp[i] += dt * P3[i, i] + end + + build_mprk_matrix!(P3, P3, σ, dt, D3) + + # Same as linres = P3 \ tmp + linsolve.A = P3 + linres = solve!(linsolve) + + u .= linres + integrator.stats.nsolve += 1 + + # Now tmp stores the error estimate + @.. broadcast=false tmp=u - σ + + # Now tmp2 stores error residuals + calculate_residuals!(tmp2, tmp, uprev, u, integrator.opts.abstol, + integrator.opts.reltol, integrator.opts.internalnorm, t, + False()) + integrator.EEst = integrator.opts.internalnorm(tmp2, t) +end + +function perform_step!(integrator, cache::MPRK43ConservativeCache, repeat_step = false) + @unpack t, dt, uprev, u, f, p = integrator + @unpack tmp, tmp2, P, P2, P3, σ, linsolve = cache + @unpack a21, a31, a32, b1, b2, b3, c2, c3, beta1, beta2, q1, q2, small_constant = cache.tab + + # Set right hand side of linear system + tmp .= uprev + + # We use P3 to store the last evaluation of the PDS + # as well as to store the system matrix of the linear system + f.p(P, uprev, p, t) # evaluate production terms + @.. broadcast=false P3=a21 * P + integrator.stats.nf += 1 + + # avoid division by zero due to zero Patankar weights + @.. broadcast=false σ=uprev + small_constant + + build_mprk_matrix!(P3, P3, σ, dt) + + # Same as linres = P3 \ tmp + linsolve.A = P3 + linres = solve!(linsolve) + + u .= linres + if !(q1 ≈ q2) + tmp2 .= u #u2 in out-of-place version + end + integrator.stats.nsolve += 1 + + @.. broadcast=false σ=σ^(1 - q1) * u^q1 + @.. broadcast=false σ=σ + small_constant + + f.p(P2, u, p, t + c2 * dt) # evaluate production terms + @.. broadcast=false P3=a31 * P + a32 * P2 + integrator.stats.nf += 1 + + build_mprk_matrix!(P3, P3, σ, dt) + + # Same as linres = P3 \ tmp + linsolve.A = P3 + linres = solve!(linsolve) + + u .= linres + integrator.stats.nsolve += 1 + + if !(q1 ≈ q2) + @.. broadcast=false σ=(uprev + small_constant)^(1 - q2) * tmp2^q2 + @.. broadcast=false σ=σ + small_constant + end + + @.. broadcast=false P3=beta1 * P + beta2 * P2 + + build_mprk_matrix!(P3, P3, σ, dt) + + # Same as linres = P3 \ tmp + linsolve.A = P3 + linres = solve!(linsolve) + integrator.stats.nsolve += 1 + + σ .= linres + # avoid division by zero due to zero Patankar weights + @.. broadcast=false σ=σ + small_constant + + f.p(P3, u, p, t + c3 * dt) # evaluate production terms + @.. broadcast=false P3=b1 * P + b2 * P2 + b3 * P3 + integrator.stats.nf += 1 + + build_mprk_matrix!(P3, P3, σ, dt) + + # Same as linres = P3 \ tmp + linsolve.A = P3 + linres = solve!(linsolve) + + u .= linres + integrator.stats.nsolve += 1 + + # Now tmp stores the error estimate + @.. broadcast=false tmp=u - σ + + # Now tmp2 stores error residuals + calculate_residuals!(tmp2, tmp, uprev, u, integrator.opts.abstol, + integrator.opts.reltol, integrator.opts.internalnorm, t, + False()) + integrator.EEst = integrator.opts.internalnorm(tmp2, t) +end +=# From f1d1203db86369c42bd25f8b730b1816c5baaa5a Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Thu, 27 Jun 2024 14:01:12 +0200 Subject: [PATCH 10/48] Added SSPMPRK43 --- src/PositiveIntegrators.jl | 2 +- src/sspmprk.jl | 134 +++++++++++++++++++------------------ 2 files changed, 69 insertions(+), 67 deletions(-) diff --git a/src/PositiveIntegrators.jl b/src/PositiveIntegrators.jl index 0ed50f58..5eb8feef 100644 --- a/src/PositiveIntegrators.jl +++ b/src/PositiveIntegrators.jl @@ -48,7 +48,7 @@ export PDSFunction, PDSProblem export ConservativePDSFunction, ConservativePDSProblem export MPE, MPRK22, MPRK43I, MPRK43II -export SSPMPRK22 +export SSPMPRK22, SSPMPRK43 export prob_pds_linmod, prob_pds_linmod_inplace, prob_pds_nonlinmod, prob_pds_robertson, prob_pds_brusselator, prob_pds_sir, diff --git a/src/sspmprk.jl b/src/sspmprk.jl index 7b626d75..8aa4a9ce 100644 --- a/src/sspmprk.jl +++ b/src/sspmprk.jl @@ -405,7 +405,7 @@ as keyword argument `linsolve`. ESAIM: Mathematical Modelling and Numerical Analysis 57 (2023):1063–1086 [DOI: 10.1051/m2an/2023005](https://doi.org/10.1051/m2an/2023005) """ -struct SSPMPRK43{F} <: OrdinaryDiffEqAdaptiveAlgorithm +struct SSPMPRK43{F} <: OrdinaryDiffEqAlgorithm linsolve::F end @@ -570,7 +570,7 @@ function perform_step!(integrator, cache::SSPMPRK43ConstantCache, repeat_step = integrator.stats.nsolve += 1 # compute Patankar weight denominator - σ = σ + z .* uprev .* u ./ rho + σ = σ + z .* uprev .* u ./ ρ # avoid division by zero due to zero Patankar weights σ = add_small_constant(σ, small_constant) @@ -607,8 +607,7 @@ function perform_step!(integrator, cache::SSPMPRK43ConstantCache, repeat_step = integrator.u = u end -#= -struct MPRK43Cache{uType, PType, tabType, F} <: OrdinaryDiffEqMutableCache +struct SSPMPRK43Cache{uType, PType, tabType, F} <: OrdinaryDiffEqMutableCache tmp::uType tmp2::uType P::PType @@ -618,29 +617,31 @@ struct MPRK43Cache{uType, PType, tabType, F} <: OrdinaryDiffEqMutableCache D2::uType D3::uType σ::uType + ρ::uType tab::tabType linsolve::F end -struct MPRK43ConservativeCache{uType, PType, tabType, F} <: OrdinaryDiffEqMutableCache +struct SSPMPRK43ConservativeCache{uType, PType, tabType, F} <: OrdinaryDiffEqMutableCache tmp::uType tmp2::uType P::PType P2::PType P3::PType σ::uType + ρ::uType tab::tabType linsolve::F end # In-place -function alg_cache(alg::Union{MPRK43I, MPRK43II}, u, rate_prototype, ::Type{uEltypeNoUnits}, +function alg_cache(alg::SSPMPRK43, u, rate_prototype, ::Type{uEltypeNoUnits}, ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, dt, reltol, p, calck, ::Val{true}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} - a21, a31, a32, b1, b2, b3, c2, c3, beta1, beta2, q1, q2 = get_constant_parameters(alg) - tab = MPRK43ConstantCache(a21, a31, a32, b1, b2, b3, c2, c3, - beta1, beta2, q1, q2, floatmin(uEltypeNoUnits)) + n1, n2, z, η1, η2, η3, η4, s, α10, α20, α21, α30, α31, α32, β10, β20, β21, β30, β31, β32, c3 = get_constant_parameters(alg) + tab = SSPMPRK43ConstantCache(n1, n2, z, η1, η2, η3, η4, s, α10, α20, α21, α30, α31, α32, + β10, β20, β21, β30, β31, β32, c3, floatmin(uEltypeNoUnits)) tmp = zero(u) tmp2 = zero(u) P = p_prototype(u, f) @@ -649,16 +650,13 @@ function alg_cache(alg::Union{MPRK43I, MPRK43II}, u, rate_prototype, ::Type{uElt # as well as to store the system matrix of the linear system P3 = p_prototype(u, f) σ = zero(u) + ρ = zero(u) if f isa ConservativePDSFunction - # The right hand side of the linear system is always uprev. But using - # tmp instead of uprev for the rhs we allow `alias_b=true`. uprev must - # not be altered, since it is needed to compute the adaptive time step - # size. linprob = LinearProblem(P3, _vec(tmp)) linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true, assumptions = LinearSolve.OperatorAssumptions(true)) - MPRK43ConservativeCache(tmp, tmp2, P, P2, P3, σ, tab, linsolve) + SSPMPRK43ConservativeCache(tmp, tmp2, P, P2, P3, σ, ρ, tab, linsolve) elseif f isa PDSFunction D = zero(u) D2 = zero(u) @@ -668,34 +666,34 @@ function alg_cache(alg::Union{MPRK43I, MPRK43II}, u, rate_prototype, ::Type{uElt linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true, assumptions = LinearSolve.OperatorAssumptions(true)) - MPRK43Cache(tmp, tmp2, P, P2, P3, D, D2, D3, σ, tab, linsolve) + SSPMPRK43Cache(tmp, tmp2, P, P2, P3, D, D2, D3, σ, ρ, tab, linsolve) else - throw(ArgumentError("MPRK43 can only be applied to production-destruction systems")) + throw(ArgumentError("SSPMPRK43 can only be applied to production-destruction systems")) end end -function initialize!(integrator, cache::Union{MPRK43Cache, MPRK43ConservativeCache}) +function initialize!(integrator, cache::Union{SSPMPRK43ConservativeCache, SSPMPRK43Cache}) end -function perform_step!(integrator, cache::MPRK43Cache, repeat_step = false) +function perform_step!(integrator, cache::SSPMPRK43Cache, repeat_step = false) @unpack t, dt, uprev, u, f, p = integrator - @unpack tmp, tmp2, P, P2, P3, D, D2, D3, σ, linsolve = cache - @unpack a21, a31, a32, b1, b2, b3, c2, c3, beta1, beta2, q1, q2, small_constant = cache.tab + @unpack tmp, tmp2, P, P2, P3, D, D2, D3, σ, ρ, linsolve = cache + @unpack n1, n2, z, η1, η2, η3, η4, s, α10, α20, α21, α30, α31, α32, β10, β20, β21, β30, β31, β32, c3, small_constant = cache.tab # We use P3 to store the last evaluation of the PDS # as well as to store the system matrix of the linear system f.p(P, uprev, p, t) # evaluate production terms f.d(D, uprev, p, t) # evaluate nonconservative destruction terms - @.. broadcast=false P3=a21 * P - @.. broadcast=false D3=a21 * D + @.. broadcast=false P3=β10 * P + @.. broadcast=false D3=β10 * D integrator.stats.nf += 1 # avoid division by zero due to zero Patankar weights @.. broadcast=false σ=uprev + small_constant # tmp holds the right hand side of the linear system - tmp .= uprev + @.. broadcast=false tmp=α10 * uprev @inbounds for i in eachindex(tmp) tmp[i] += dt * P3[i, i] end @@ -707,27 +705,25 @@ function perform_step!(integrator, cache::MPRK43Cache, repeat_step = false) linres = solve!(linsolve) u .= linres - if !(q1 ≈ q2) - tmp2 .= u #u2 in out-of-place version - end + tmp2 .= u integrator.stats.nsolve += 1 - @.. broadcast=false σ=σ^(1 - q1) * u^q1 - @.. broadcast=false σ=σ + small_constant + @.. broadcast=false ρ=n1 * u + n2 * u^2 / σ + @.. broadcast=false ρ=ρ + small_constant - f.p(P2, u, p, t + c2 * dt) # evaluate production terms - f.d(D2, u, p, t + c2 * dt) # evaluate nonconservative destruction terms - @.. broadcast=false P3=a31 * P + a32 * P2 - @.. broadcast=false D3=a31 * D + a32 * D2 + f.p(P2, u, p, t + β10 * dt) # evaluate production terms + f.d(D2, u, p, t + β10 * dt) # evaluate nonconservative destruction terms + @.. broadcast=false P3=β20 * P + β21 * P2 + @.. broadcast=false D3=β20 * D + β21 * D2 integrator.stats.nf += 1 # tmp holds the right hand side of the linear system - tmp .= uprev + @.. broadcast=false tmp=α20 * uprev + α21 * tmp2 @inbounds for i in eachindex(tmp) tmp[i] += dt * P3[i, i] end - build_mprk_matrix!(P3, P3, σ, dt, D3) + build_mprk_matrix!(P3, P3, ρ, dt, D3) # Same as linres = P3 \ tmp linsolve.A = P3 @@ -736,16 +732,14 @@ function perform_step!(integrator, cache::MPRK43Cache, repeat_step = false) u .= linres integrator.stats.nsolve += 1 - if !(q1 ≈ q2) - @.. broadcast=false σ=(uprev + small_constant)^(1 - q2) * tmp2^q2 - @.. broadcast=false σ=σ + small_constant - end + @.. broadcast=false σ=σ^(1 - s) * tmp2^s + @.. broadcast=false σ=σ + small_constant - @.. broadcast=false P3=beta1 * P + beta2 * P2 - @.. broadcast=false D3=beta1 * D + beta2 * D2 + @.. broadcast=false P3=η3 * P + η4 * P2 + @.. broadcast=false D3=η3 * D + η4 * D2 # tmp holds the right hand side of the linear system - tmp .= uprev + @.. broadcast=false tmp=η1 * uprev + η2 * tmp2 @inbounds for i in eachindex(tmp) tmp[i] += dt * P3[i, i] end @@ -758,17 +752,18 @@ function perform_step!(integrator, cache::MPRK43Cache, repeat_step = false) integrator.stats.nsolve += 1 σ .= linres + @.. broadcast=false σ=σ + z * uprev * u / ρ # avoid division by zero due to zero Patankar weights @.. broadcast=false σ=σ + small_constant f.p(P3, u, p, t + c3 * dt) # evaluate production terms f.d(D3, u, p, t + c3 * dt) # evaluate nonconservative destruction terms - @.. broadcast=false P3=b1 * P + b2 * P2 + b3 * P3 - @.. broadcast=false D3=b1 * D + b2 * D2 + b3 * D3 + @.. broadcast=false P3=β30 * P + β31 * P2 + β32 * P3 + @.. broadcast=false D3=β30 * D + β31 * D2 + β32 * D3 integrator.stats.nf += 1 # tmp holds the right hand side of the linear system - tmp .= uprev + @.. broadcast=false tmp=α30 * uprev + α31 * tmp2 + α32 * u @inbounds for i in eachindex(tmp) tmp[i] += dt * P3[i, i] end @@ -782,6 +777,9 @@ function perform_step!(integrator, cache::MPRK43Cache, repeat_step = false) u .= linres integrator.stats.nsolve += 1 + #TODO: Figure out if a second order approximation of the solution + # is hidden somewhere. + #= # Now tmp stores the error estimate @.. broadcast=false tmp=u - σ @@ -790,25 +788,26 @@ function perform_step!(integrator, cache::MPRK43Cache, repeat_step = false) integrator.opts.reltol, integrator.opts.internalnorm, t, False()) integrator.EEst = integrator.opts.internalnorm(tmp2, t) + =# end -function perform_step!(integrator, cache::MPRK43ConservativeCache, repeat_step = false) +function perform_step!(integrator, cache::SSPMPRK43ConservativeCache, repeat_step = false) @unpack t, dt, uprev, u, f, p = integrator - @unpack tmp, tmp2, P, P2, P3, σ, linsolve = cache - @unpack a21, a31, a32, b1, b2, b3, c2, c3, beta1, beta2, q1, q2, small_constant = cache.tab - - # Set right hand side of linear system - tmp .= uprev + @unpack tmp, tmp2, P, P2, P3, σ, ρ, linsolve = cache + @unpack n1, n2, z, η1, η2, η3, η4, s, α10, α20, α21, α30, α31, α32, β10, β20, β21, β30, β31, β32, c3, small_constant = cache.tab # We use P3 to store the last evaluation of the PDS # as well as to store the system matrix of the linear system f.p(P, uprev, p, t) # evaluate production terms - @.. broadcast=false P3=a21 * P + @.. broadcast=false P3=β10 * P integrator.stats.nf += 1 # avoid division by zero due to zero Patankar weights @.. broadcast=false σ=uprev + small_constant + # tmp holds the right hand side of the linear system + @.. broadcast=false tmp=α10 * uprev + build_mprk_matrix!(P3, P3, σ, dt) # Same as linres = P3 \ tmp @@ -816,19 +815,18 @@ function perform_step!(integrator, cache::MPRK43ConservativeCache, repeat_step = linres = solve!(linsolve) u .= linres - if !(q1 ≈ q2) - tmp2 .= u #u2 in out-of-place version - end + tmp2 .= u integrator.stats.nsolve += 1 - @.. broadcast=false σ=σ^(1 - q1) * u^q1 - @.. broadcast=false σ=σ + small_constant + @.. broadcast=false ρ=n1 * u + n2 * u^2 / σ + @.. broadcast=false ρ=ρ + small_constant - f.p(P2, u, p, t + c2 * dt) # evaluate production terms - @.. broadcast=false P3=a31 * P + a32 * P2 + f.p(P2, u, p, t + β10 * dt) # evaluate production terms + @.. broadcast=false P3=β20 * P + β21 * P2 integrator.stats.nf += 1 - build_mprk_matrix!(P3, P3, σ, dt) + @.. broadcast=false tmp=α20 * uprev + α21 * tmp2 + build_mprk_matrix!(P3, P3, ρ, dt) # Same as linres = P3 \ tmp linsolve.A = P3 @@ -837,12 +835,11 @@ function perform_step!(integrator, cache::MPRK43ConservativeCache, repeat_step = u .= linres integrator.stats.nsolve += 1 - if !(q1 ≈ q2) - @.. broadcast=false σ=(uprev + small_constant)^(1 - q2) * tmp2^q2 - @.. broadcast=false σ=σ + small_constant - end + @.. broadcast=false σ=σ^(1 - s) * tmp2^s + @.. broadcast=false σ=σ + small_constant - @.. broadcast=false P3=beta1 * P + beta2 * P2 + @.. broadcast=false P3=η3 * P + η4 * P2 + @.. broadcast=false tmp=η1 * uprev + η2 * tmp2 build_mprk_matrix!(P3, P3, σ, dt) @@ -852,13 +849,15 @@ function perform_step!(integrator, cache::MPRK43ConservativeCache, repeat_step = integrator.stats.nsolve += 1 σ .= linres + @.. broadcast=false σ=σ + z * uprev * u / ρ # avoid division by zero due to zero Patankar weights @.. broadcast=false σ=σ + small_constant f.p(P3, u, p, t + c3 * dt) # evaluate production terms - @.. broadcast=false P3=b1 * P + b2 * P2 + b3 * P3 + @.. broadcast=false P3=β30 * P + β31 * P2 + β32 * P3 integrator.stats.nf += 1 + @.. broadcast=false tmp=α30 * uprev + α31 * tmp2 + α32 * u build_mprk_matrix!(P3, P3, σ, dt) # Same as linres = P3 \ tmp @@ -868,6 +867,9 @@ function perform_step!(integrator, cache::MPRK43ConservativeCache, repeat_step = u .= linres integrator.stats.nsolve += 1 + #TODO: Figure out if a second order approximation of the solution + # is hidden somewhere. + #= # Now tmp stores the error estimate @.. broadcast=false tmp=u - σ @@ -876,5 +878,5 @@ function perform_step!(integrator, cache::MPRK43ConservativeCache, repeat_step = integrator.opts.reltol, integrator.opts.internalnorm, t, False()) integrator.EEst = integrator.opts.internalnorm(tmp2, t) + =# end -=# From a5b08c4a20b02f70cff6576a70797e905b3273e6 Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Thu, 27 Jun 2024 14:11:47 +0200 Subject: [PATCH 11/48] Added tests for SSPMPRK43 --- test/runtests.jl | 31 +++++++++++++++++-------------- 1 file changed, 17 insertions(+), 14 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index 8e575459..5b4125e2 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -339,6 +339,10 @@ const prob_pds_linmod_nonconservative_inplace = PDSProblem(linmodP!, linmodD!, [ @test_throws "SSPMPRK22 requires 0 ≤ α ≤ 1, β ≥ 0 and αβ + 1/(2β) ≤ 1." solve(prob_pds_linmod, SSPMPRK22(1.0, 10.0)) + @test_throws "SSPMPRK43 can only be applied to production-destruction systems" solve(prob_oop, + SSPMPRK43()) + @test_throws "SSPMPRK43 can only be applied to production-destruction systems" solve(prob_ip, + SSPMPRK43()) end # Here we check that MPE equals implicit Euler (IE) for a linear PDS @@ -444,7 +448,8 @@ const prob_pds_linmod_nonconservative_inplace = PDSProblem(linmodP!, linmodD!, [ (; kwargs...) -> MPRK43I(0.5, 0.75; kwargs...), (; kwargs...) -> MPRK43II(0.5; kwargs...), (; kwargs...) -> MPRK43II(2.0 / 3.0; kwargs...), - (; kwargs...) -> SSPMPRK22(0.5, 1.0; kwargs...)) + (; kwargs...) -> SSPMPRK22(0.5, 1.0; kwargs...), + (; kwargs...) -> SSPMPRK43(; kwargs...)) for alg in algs # Check different linear solvers @@ -506,7 +511,7 @@ const prob_pds_linmod_nonconservative_inplace = PDSProblem(linmodP!, linmodD!, [ 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)) + 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))) @@ -623,7 +628,7 @@ const prob_pds_linmod_nonconservative_inplace = PDSProblem(linmodP!, linmodD!, [ 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)) + SSPMPRK22(0.5, 1.0), SSPMPRK43()) for (prod!, dest!) in zip((prod_1!, prod_2!, prod_3!), (dest_1!, dest_2!, dest_3!)) prod = (u, p, t) -> begin @@ -734,12 +739,9 @@ const prob_pds_linmod_nonconservative_inplace = PDSProblem(linmodP!, linmodD!, [ dts = 0.5 .^ (6:11) problems = (prob_pds_linmod, prob_pds_linmod_array, prob_pds_linmod_mvector, prob_pds_linmod_inplace) - for alg in [ - MPRK43I(1.0, 0.5), - MPRK43I(0.5, 0.75), - MPRK43II(0.5), - MPRK43II(2.0 / 3.0), - ], prob in problems + algs = (MPRK43I(1.0, 0.5), MPRK43I(0.5, 0.75), + MPRK43II(0.5), MPRK43II(2.0 / 3.0), SSPMPRK43()) + for alg in algs, prob in problems eoc = experimental_order_of_convergence(prob, alg, dts) @test isapprox(eoc, PositiveIntegrators.alg_order(alg); atol = 0.2) end @@ -752,7 +754,7 @@ const prob_pds_linmod_nonconservative_inplace = PDSProblem(linmodP!, linmodD!, [ problems = (prob_pds_linmod_nonconservative, prob_pds_linmod_nonconservative_inplace) algs = (MPRK43I(1.0, 0.5), MPRK43I(0.5, 0.75), - MPRK43II(0.5), MPRK43II(2.0 / 3.0)) + MPRK43II(0.5), MPRK43II(2.0 / 3.0), SSPMPRK43()) for alg in algs, prob in problems eoc = experimental_order_of_convergence(prob, alg, dts) @test isapprox(eoc, PositiveIntegrators.alg_order(alg); atol = 0.2) @@ -761,7 +763,8 @@ const prob_pds_linmod_nonconservative_inplace = PDSProblem(linmodP!, linmodD!, [ @testset "Interpolation tests (conservative)" begin algs = (MPE(), 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)) + MPRK43I(0.5, 0.75), MPRK43II(0.5), MPRK43II(2.0 / 3.0), + SSPMPRK22(0.5, 1.0), SSPMPRK43()) dt = 0.5^6 problems = (prob_pds_linmod, prob_pds_linmod_array, prob_pds_linmod_mvector, prob_pds_linmod_inplace) @@ -778,7 +781,7 @@ const prob_pds_linmod_nonconservative_inplace = PDSProblem(linmodP!, linmodD!, [ @testset "Interpolation tests (nonconservative)" begin algs = (MPE(), 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)) + SSPMPRK22(0.5, 1.0), SSPMPRK43()) dt = 0.5^6 problems = (prob_pds_linmod_nonconservative, prob_pds_linmod_nonconservative_inplace) @@ -825,7 +828,7 @@ const prob_pds_linmod_nonconservative_inplace = PDSProblem(linmodP!, linmodD!, [ algs = (MPE(), 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)) + MPRK43II(2.0 / 3.0), SSPMPRK22(0.5, 1.0), SSPMPRK43()) for alg in algs sol = solve(prob_ip, alg; dt = dt, adaptive = false) @@ -873,7 +876,7 @@ const prob_pds_linmod_nonconservative_inplace = PDSProblem(linmodP!, linmodD!, [ algs = (MPE(), 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)) + MPRK43II(2.0 / 3.0), SSPMPRK22(0.5, 1.0), SSPMPRK43()) dt = 1e-3 for alg in algs From 9fa8ce72c45e82ecf02b62f823806790b1cbd10e Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Thu, 27 Jun 2024 14:17:22 +0200 Subject: [PATCH 12/48] fixed bug in tests --- test/runtests.jl | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index 5b4125e2..c735ff89 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -340,9 +340,11 @@ const prob_pds_linmod_nonconservative_inplace = PDSProblem(linmodP!, linmodD!, [ SSPMPRK22(1.0, 10.0)) @test_throws "SSPMPRK43 can only be applied to production-destruction systems" solve(prob_oop, - SSPMPRK43()) + SSPMPRK43(), + dt = 0.1) @test_throws "SSPMPRK43 can only be applied to production-destruction systems" solve(prob_ip, - SSPMPRK43()) + SSPMPRK43(), + dt = 0.1) end # Here we check that MPE equals implicit Euler (IE) for a linear PDS From b95d44ff67f79596b682af3dc09370b4d743f2bd Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Thu, 27 Jun 2024 15:03:45 +0200 Subject: [PATCH 13/48] Fixed bug in SSPMPRK43: Typo --- src/sspmprk.jl | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/sspmprk.jl b/src/sspmprk.jl index 8aa4a9ce..278f7aec 100644 --- a/src/sspmprk.jl +++ b/src/sspmprk.jl @@ -556,8 +556,11 @@ function perform_step!(integrator, cache::SSPMPRK43ConstantCache, repeat_step = # build linear system matrix and rhs if f isa PDSFunction dtmp = η3 * d + η4 * d2 + rhs = η1 * uprev + η2 * u2 + dt * diag(Ptmp) + M = build_mprk_matrix(Ptmp, σ, dt, dtmp) + else rhs = η1 * uprev + η2 * u2 M = build_mprk_matrix(Ptmp, σ, dt) @@ -581,7 +584,7 @@ function perform_step!(integrator, cache::SSPMPRK43ConstantCache, repeat_step = # build linear system matrix if f isa PDSFunction d3 = f.d(u, p, t + c3 * dt) # evaluate nonconservative destruction terms - dtmp = β30 * d + β31 * d2 + β30 * d3 + dtmp = β30 * d + β31 * d2 + β32 * d3 rhs = α30 * uprev + α31 * u2 + α32 * u + dt * diag(Ptmp) M = build_mprk_matrix(Ptmp, σ, dt, dtmp) else @@ -752,6 +755,7 @@ function perform_step!(integrator, cache::SSPMPRK43Cache, repeat_step = false) integrator.stats.nsolve += 1 σ .= linres + @.. broadcast=false σ=σ + z * uprev * u / ρ # avoid division by zero due to zero Patankar weights @.. broadcast=false σ=σ + small_constant From cf8a475093fdd6ef554c6a354f37d27ec959ef4e Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Mon, 1 Jul 2024 11:07:18 +0200 Subject: [PATCH 14/48] Fixed bug to obtain 3rd order for nonconservative PDS --- src/sspmprk.jl | 31 +++++++++++++++++++++---------- 1 file changed, 21 insertions(+), 10 deletions(-) diff --git a/src/sspmprk.jl b/src/sspmprk.jl index 278f7aec..061fe754 100644 --- a/src/sspmprk.jl +++ b/src/sspmprk.jl @@ -417,6 +417,8 @@ alg_order(::SSPMPRK43) = 3 isfsal(::SSPMPRK43) = false function get_constant_parameters(alg::SSPMPRK43) + # parameters from original paper + n1 = 2.569046025732011E-01 n2 = 7.430953974267989E-01 z = 6.288938077828750E-01 @@ -426,6 +428,9 @@ function get_constant_parameters(alg::SSPMPRK43) η4 = 2.224876040351123 s = 5.721964308755304 + η5 = η3 * (η1 + η2) + η6 = η4 * (η1 + η2) + α10 = 1.0 α20 = 9.2600312554031827E-01 α21 = 7.3996874459681783E-02 @@ -441,7 +446,8 @@ function get_constant_parameters(alg::SSPMPRK43) c3 = β20 + α21 * β10 + β21 - return n1, n2, z, η1, η2, η3, η4, s, α10, α20, α21, α30, α31, α32, β10, β20, β21, β30, + return n1, n2, z, η1, η2, η3, η4, η5, η6, s, α10, α20, α21, α30, α31, α32, β10, β20, + β21, β30, β31, β32, c3 end @@ -453,6 +459,8 @@ struct SSPMPRK43ConstantCache{T} <: OrdinaryDiffEqConstantCache η2::T η3::T η4::T + η5::T + η6::T s::T α10::T α20::T @@ -478,8 +486,9 @@ 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, s, α10, α20, α21, α30, α31, α32, β10, β20, β21, β30, β31, β32, c3 = get_constant_parameters(alg) - SSPMPRK43ConstantCache(n1, n2, z, η1, η2, η3, η4, s, α10, α20, α21, α30, α31, α32, β10, + 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) + SSPMPRK43ConstantCache(n1, n2, z, η1, η2, η3, η4, η5, η6, s, α10, α20, α21, α30, α31, + α32, β10, β20, β21, β30, β31, β32, c3, floatmin(uEltypeNoUnits)) end @@ -489,7 +498,7 @@ end function perform_step!(integrator, cache::SSPMPRK43ConstantCache, repeat_step = false) @unpack alg, t, dt, uprev, f, p = integrator - @unpack n1, n2, z, η1, η2, η3, η4, s, α10, α20, α21, α30, α31, α32, β10, β20, β21, β30, β31, β32, c3, small_constant = cache + @unpack n1, n2, z, η1, η2, η3, η4, η5, η6, s, α10, α20, α21, α30, α31, α32, β10, β20, β21, β30, β31, β32, c3, small_constant = cache f = integrator.f @@ -557,10 +566,10 @@ function perform_step!(integrator, cache::SSPMPRK43ConstantCache, repeat_step = if f isa PDSFunction dtmp = η3 * d + η4 * d2 - rhs = η1 * uprev + η2 * u2 + dt * diag(Ptmp) + # see (3.25 f) in original paper + rhs = η1 * uprev + η2 * u2 + dt * (η5 * diag(P) + η6 * diag(P2)) M = build_mprk_matrix(Ptmp, σ, dt, dtmp) - else rhs = η1 * uprev + η2 * u2 M = build_mprk_matrix(Ptmp, σ, dt) @@ -642,8 +651,9 @@ function alg_cache(alg::SSPMPRK43, u, rate_prototype, ::Type{uEltypeNoUnits}, ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, dt, reltol, p, calck, ::Val{true}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} - n1, n2, z, η1, η2, η3, η4, s, α10, α20, α21, α30, α31, α32, β10, β20, β21, β30, β31, β32, c3 = get_constant_parameters(alg) - tab = SSPMPRK43ConstantCache(n1, n2, z, η1, η2, η3, η4, s, α10, α20, α21, α30, α31, α32, + 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) + tab = SSPMPRK43ConstantCache(n1, n2, z, η1, η2, η3, η4, η5, η6, s, α10, α20, α21, α30, + α31, α32, β10, β20, β21, β30, β31, β32, c3, floatmin(uEltypeNoUnits)) tmp = zero(u) tmp2 = zero(u) @@ -681,7 +691,7 @@ end function perform_step!(integrator, cache::SSPMPRK43Cache, repeat_step = false) @unpack t, dt, uprev, u, f, p = integrator @unpack tmp, tmp2, P, P2, P3, D, D2, D3, σ, ρ, linsolve = cache - @unpack n1, n2, z, η1, η2, η3, η4, s, α10, α20, α21, α30, α31, α32, β10, β20, β21, β30, β31, β32, c3, small_constant = cache.tab + @unpack n1, n2, z, η1, η2, η3, η4, η5, η6, s, α10, α20, α21, α30, α31, α32, β10, β20, β21, β30, β31, β32, c3, small_constant = cache.tab # We use P3 to store the last evaluation of the PDS # as well as to store the system matrix of the linear system @@ -744,7 +754,8 @@ function perform_step!(integrator, cache::SSPMPRK43Cache, repeat_step = false) # tmp holds the right hand side of the linear system @.. broadcast=false tmp=η1 * uprev + η2 * tmp2 @inbounds for i in eachindex(tmp) - tmp[i] += dt * P3[i, i] + # see (3.25 f) in original paper + tmp[i] += dt * (η5 * P[i, i] + η6 * P2[i, i]) end build_mprk_matrix!(P3, P3, σ, dt, D3) From b59d1d41b0d392289714e88c00b867c77cf41035 Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Mon, 1 Jul 2024 11:59:51 +0200 Subject: [PATCH 15/48] Changed computation of experimental order. --- test/runtests.jl | 146 +++++++++++++++++++++++++++++------------------ 1 file changed, 89 insertions(+), 57 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index c735ff89..b4779694 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -13,53 +13,46 @@ using LinearSolve: RFLUFactorization, LUFactorization, KrylovJL_GMRES using Aqua: Aqua """ - experimental_order_of_convergence(prob, alg, dts, test_times; - only_first_index = false) + experimental_order_of_convergence(prob, alg, dts; test_time = nothing) Solve `prob` with `alg` and fixed time steps taken from `dts`, and compute -the mean error at the times `test_times`. -Return the associated experimental order of convergence. +the errors at `test_time`. If`test_time` is not specified the error is computed +at the final time. +Return the associated experimental orders of convergence. If `only_first_index == true`, only the first solution component is used to compute the error. """ -function experimental_order_of_convergence(prob, alg, dts, test_times; +function experimental_order_of_convergence(prob, alg, dts; test_time = nothing, only_first_index = false) @assert length(dts) > 1 errors = zeros(eltype(dts), length(dts)) analytic = t -> prob.f.analytic(prob.u0, prob.p, t) - for (i, dt) in enumerate(dts) - sol = solve(prob, alg; dt = dt, adaptive = false) + if isnothing(test_time) if only_first_index - errors[i] = mean(test_times) do t - norm(sol(t; idxs = 1) - first(analytic(t))) + for (i, dt) in enumerate(dts) + sol = solve(prob, alg; dt = dt, adaptive = false, save_everystep = false) + errors[i] = norm(sol.u[end][1] - analytic(sol.t[end])[1]) end else - errors[i] = mean(test_times) do t - norm(sol(t) - analytic(t)) + for (i, dt) in enumerate(dts) + sol = solve(prob, alg; dt = dt, adaptive = false, save_everystep = false) + errors[i] = norm(sol.u[end] - analytic(sol.t[end])) + end + end + else + if only_first_index + for (i, dt) in enumerate(dts) + sol = solve(prob, alg; dt = dt, adaptive = false) + errors[i] = norm(sol(test_time; idxs = 1) - first(analytic(test_time))) + end + else + for (i, dt) in enumerate(dts) + sol = solve(prob, alg; dt = dt, adaptive = false) + errors[i] = norm(sol(test_time) - analytic(test_time)) end end - end - - return experimental_order_of_convergence(errors, dts) -end - -""" - experimental_order_of_convergence(prob, alg, dts) - -Solve `prob` with `alg` and fixed time steps taken from `dts`, and compute -the mean error at the final time. -Return the associated experimental order of convergence. -""" -function experimental_order_of_convergence(prob, alg, dts) - @assert length(dts) > 1 - errors = zeros(eltype(dts), length(dts)) - analytic = t -> prob.f.analytic(prob.u0, prob.p, t) - - for (i, dt) in enumerate(dts) - sol = solve(prob, alg; dt = dt, adaptive = false, save_everystep = false) - errors[i] = norm(sol.u[end] - analytic(sol.t[end])) end return experimental_order_of_convergence(errors, dts) @@ -68,7 +61,7 @@ end """ experimental_order_of_convergence(errors, dts) -Compute the experimental order of convergence for given `errors` and +Compute the experimental orders of convergence for given `errors` and time step sizes `dts`. """ function experimental_order_of_convergence(errors, dts) @@ -80,7 +73,31 @@ function experimental_order_of_convergence(errors, dts) orders[i] = log(errors[i] / errors[i + 1]) / log(dts[i] / dts[i + 1]) end - return mean(orders) + return orders +end + +""" + check_order(orders, alg_order; N = 3, atol = 0.1) + +Check if `alg_order` can be found approximately in `N` consecutive elements of `orders`. +""" +function check_order(orders, alg_order; N = 3, atol = 0.1) + check = false + # The calculation of the indices of the permissible orders by + # + # indices = findall(x -> isapprox(x, alg_order; atol = atol), orders) + # + # sometimes failed because the experimental order was too good. To accept + # such cases, but avoid to decrease atol, we now use the following asymmetric criterion: + indices = findall(x -> -atol ≤ x - alg_order ≤ 2 * atol, orders) + + for (i, idx) in enumerate(indices) + if i + N - 1 ≤ length(indices) && indices[i:(i + N - 1)] == idx:(idx + N - 1) + check = true + break + end + end + return check end const prob_pds_linmod_array = ConservativePDSProblem(prob_pds_linmod.f, @@ -688,24 +705,32 @@ const prob_pds_linmod_nonconservative_inplace = PDSProblem(linmodP!, linmodD!, [ # also an interpolation of order p is available @testset "Convergence tests (conservative)" begin algs = (MPE(), MPRK22(0.5), MPRK22(1.0), MPRK22(2.0), SSPMPRK22(0.5, 1.0)) - dts = 0.5 .^ (6:11) + dts = 0.5 .^ (4:15) problems = (prob_pds_linmod, prob_pds_linmod_array, prob_pds_linmod_mvector, prob_pds_linmod_inplace) - for alg in algs + @testset "$alg" for alg in algs + alg = MPRK22(1.0) for prob in problems - eoc = experimental_order_of_convergence(prob, alg, dts) - @test isapprox(eoc, PositiveIntegrators.alg_order(alg); atol = 0.2) + prob = problems[1] + orders = experimental_order_of_convergence(prob, alg, dts) + @test check_order(orders, PositiveIntegrators.alg_order(alg)) test_times = [ 0.123456789, 1 / pi, exp(-1), 1.23456789, 1 + 1 / pi, 1 + exp(-1), ] - eoc = experimental_order_of_convergence(prob, alg, dts, test_times) - @test isapprox(eoc, PositiveIntegrators.alg_order(alg); atol = 0.2) - eoc = experimental_order_of_convergence(prob, alg, dts, test_times; - only_first_index = true) - @test isapprox(eoc, PositiveIntegrators.alg_order(alg); atol = 0.2) + for test_time in test_times + orders = experimental_order_of_convergence(prob, alg, dts; + test_time) + @test check_order(orders, PositiveIntegrators.alg_order(alg), + atol = 0.2) + orders = experimental_order_of_convergence(prob, alg, dts; + test_time, + only_first_index = true) + @test check_order(orders, PositiveIntegrators.alg_order(alg), + atol = 0.2) + end end end end @@ -714,23 +739,30 @@ const prob_pds_linmod_nonconservative_inplace = PDSProblem(linmodP!, linmodD!, [ # also an interpolation of order p is available @testset "Convergence tests (nonconservative)" begin algs = (MPE(), MPRK22(0.5), MPRK22(1.0), MPRK22(2.0), SSPMPRK22(0.5, 1.0)) - dts = 0.5 .^ (6:11) + dts = 0.5 .^ (4:15) problems = (prob_pds_linmod_nonconservative, prob_pds_linmod_nonconservative_inplace) - for alg in algs + @testset "$alg" for alg in algs + alg = MPRK22(1.0) for prob in problems - eoc = experimental_order_of_convergence(prob, alg, dts) - @test isapprox(eoc, PositiveIntegrators.alg_order(alg); atol = 0.2) + orders = experimental_order_of_convergence(prob, alg, dts) + @test check_order(orders, PositiveIntegrators.alg_order(alg)) test_times = [ 0.123456789, 1 / pi, exp(-1), 1.23456789, 1 + 1 / pi, 1 + exp(-1), ] - eoc = experimental_order_of_convergence(prob, alg, dts, test_times) - @test isapprox(eoc, PositiveIntegrators.alg_order(alg); atol = 0.2) - eoc = experimental_order_of_convergence(prob, alg, dts, test_times; - only_first_index = true) - @test isapprox(eoc, PositiveIntegrators.alg_order(alg); atol = 0.2) + for test_time in test_times + orders = experimental_order_of_convergence(prob, alg, dts; + test_time) + @test check_order(orders, PositiveIntegrators.alg_order(alg), + atol = 0.2) + orders = experimental_order_of_convergence(prob, alg, dts; + test_time, + only_first_index = true) + @test check_order(orders, PositiveIntegrators.alg_order(alg), + atol = 0.2) + end end end end @@ -738,28 +770,28 @@ const prob_pds_linmod_nonconservative_inplace = PDSProblem(linmodP!, linmodD!, [ # Here we check the convergence order of pth-order schemes for which # no interpolation of order p is available @testset "Convergence tests (conservative)" begin - dts = 0.5 .^ (6:11) + dts = 0.5 .^ (4:12) problems = (prob_pds_linmod, prob_pds_linmod_array, prob_pds_linmod_mvector, prob_pds_linmod_inplace) algs = (MPRK43I(1.0, 0.5), MPRK43I(0.5, 0.75), MPRK43II(0.5), MPRK43II(2.0 / 3.0), SSPMPRK43()) for alg in algs, prob in problems - eoc = experimental_order_of_convergence(prob, alg, dts) - @test isapprox(eoc, PositiveIntegrators.alg_order(alg); atol = 0.2) + orders = experimental_order_of_convergence(prob, alg, dts) + @test check_order(orders, PositiveIntegrators.alg_order(alg), atol = 0.2) end end # Here we check the convergence order of pth-order schemes for which # no interpolation of order p is available @testset "Convergence tests (nonconservative)" begin - dts = 0.5 .^ (6:11) + dts = 0.5 .^ (4:12) problems = (prob_pds_linmod_nonconservative, prob_pds_linmod_nonconservative_inplace) algs = (MPRK43I(1.0, 0.5), MPRK43I(0.5, 0.75), MPRK43II(0.5), MPRK43II(2.0 / 3.0), SSPMPRK43()) for alg in algs, prob in problems - eoc = experimental_order_of_convergence(prob, alg, dts) - @test isapprox(eoc, PositiveIntegrators.alg_order(alg); atol = 0.2) + orders = experimental_order_of_convergence(prob, alg, dts) + @test check_order(orders, PositiveIntegrators.alg_order(alg), atol = 0.2) end end From da78e9c87ce28c011d5bf3bccf0682a42b400464 Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Mon, 1 Jul 2024 13:55:20 +0200 Subject: [PATCH 16/48] Increased small_constant for SSPMPRK43 to 1e-50 and introduced abstract type MPRKCache as super type for all MPRK caches to simplify interpolation decisions --- src/PositiveIntegrators.jl | 2 ++ src/mprk.jl | 34 +++++++++++----------------------- src/sspmprk.jl | 16 ++++++++++------ 3 files changed, 23 insertions(+), 29 deletions(-) diff --git a/src/PositiveIntegrators.jl b/src/PositiveIntegrators.jl index 5eb8feef..4bc31c53 100644 --- a/src/PositiveIntegrators.jl +++ b/src/PositiveIntegrators.jl @@ -59,6 +59,8 @@ export prob_pds_linmod, prob_pds_linmod_inplace, prob_pds_nonlinmod, # production-destruction systems include("proddest.jl") +abstract type MPRKCache <: OrdinaryDiffEqMutableCache end + #TODO: sspmprk.jl must be included before mprk.jl since the interpolation specializations # at the end of mprk.jl need to know SSPMPRK22ConstantCache. # It would probably be a good idea to separate schemes and interpolations. diff --git a/src/mprk.jl b/src/mprk.jl index a09b46b0..348c627f 100644 --- a/src/mprk.jl +++ b/src/mprk.jl @@ -247,7 +247,7 @@ function perform_step!(integrator, cache::MPEConstantCache, repeat_step = false) integrator.u = u end -struct MPECache{PType, uType, tabType, F} <: OrdinaryDiffEqMutableCache +struct MPECache{PType, uType, tabType, F} <: MPRKCache P::PType D::uType σ::uType @@ -256,7 +256,7 @@ struct MPECache{PType, uType, tabType, F} <: OrdinaryDiffEqMutableCache linsolve::F end -struct MPEConservativeCache{PType, uType, tabType, F} <: OrdinaryDiffEqMutableCache +struct MPEConservativeCache{PType, uType, tabType, F} <: MPRKCache P::PType σ::uType tab::tabType @@ -565,7 +565,7 @@ function perform_step!(integrator, cache::MPRK22ConstantCache, repeat_step = fal end struct MPRK22Cache{uType, PType, tabType, F} <: - OrdinaryDiffEqMutableCache + MPRKCache tmp::uType P::PType P2::PType @@ -577,7 +577,7 @@ struct MPRK22Cache{uType, PType, tabType, F} <: end struct MPRK22ConservativeCache{uType, PType, tabType, F} <: - OrdinaryDiffEqMutableCache + MPRKCache tmp::uType P::PType P2::PType @@ -1073,7 +1073,7 @@ function perform_step!(integrator, cache::MPRK43ConstantCache, repeat_step = fal integrator.u = u end -struct MPRK43Cache{uType, PType, tabType, F} <: OrdinaryDiffEqMutableCache +struct MPRK43Cache{uType, PType, tabType, F} <: MPRKCache tmp::uType tmp2::uType P::PType @@ -1087,7 +1087,7 @@ struct MPRK43Cache{uType, PType, tabType, F} <: OrdinaryDiffEqMutableCache linsolve::F end -struct MPRK43ConservativeCache{uType, PType, tabType, F} <: OrdinaryDiffEqMutableCache +struct MPRK43ConservativeCache{uType, PType, tabType, F} <: MPRKCache tmp::uType tmp2::uType P::PType @@ -1348,18 +1348,12 @@ end ####################################################################################### # interpolation specializations function interp_summary(::Type{cacheType}, - dense::Bool) where { - cacheType <: - Union{MPRK22ConstantCache, MPRK22Cache, - SSPMPRK22ConstantCache, SSPMPRK22Cache, - MPRK43ConstantCache, MPRK43Cache}} + dense::Bool) where {cacheType <: MPRKCache} "1st order linear" end function _ode_interpolant(Θ, dt, u0, u1, k, - cache::Union{MPRK22ConstantCache, MPRK22Cache, - SSPMPRK22ConstantCache, SSPMPRK22Cache, - MPRK43ConstantCache, MPRK43Cache}, + cache::MPRKCache, idxs, # Optionally specialize for ::Nothing and others T::Type{Val{0}}, differential_vars::Nothing) @@ -1367,9 +1361,7 @@ function _ode_interpolant(Θ, dt, u0, u1, k, end function _ode_interpolant!(out, Θ, dt, u0, u1, k, - cache::Union{MPRK22ConstantCache, MPRK22Cache, - SSPMPRK22ConstantCache, SSPMPRK22Cache, - MPRK43ConstantCache, MPRK43Cache}, + cache::MPRKCache, idxs, # Optionally specialize for ::Nothing and others T::Type{Val{0}}, differential_vars::Nothing) @@ -1377,9 +1369,7 @@ function _ode_interpolant!(out, Θ, dt, u0, u1, k, end function _ode_interpolant(Θ, dt, u0, u1, k, - cache::Union{MPRK22ConstantCache, MPRK22Cache, - SSPMPRK22ConstantCache, SSPMPRK22Cache, - MPRK43ConstantCache, MPRK43Cache}, + cache::MPRKCache, idxs, # Optionally specialize for ::Nothing and others T::Type{Val{1}}, differential_vars::Nothing) @@ -1387,9 +1377,7 @@ function _ode_interpolant(Θ, dt, u0, u1, k, end function _ode_interpolant!(out, Θ, dt, u0, u1, k, - cache::Union{MPRK22ConstantCache, MPRK22Cache, - SSPMPRK22ConstantCache, SSPMPRK22Cache, - MPRK43ConstantCache, MPRK43Cache}, + cache::MPRKCache, idxs, # Optionally specialize for ::Nothing and others T::Type{Val{1}}, differential_vars::Nothing) diff --git a/src/sspmprk.jl b/src/sspmprk.jl index 061fe754..875ba3f4 100644 --- a/src/sspmprk.jl +++ b/src/sspmprk.jl @@ -172,7 +172,7 @@ function perform_step!(integrator, cache::SSPMPRK22ConstantCache, repeat_step = end struct SSPMPRK22Cache{uType, PType, tabType, F} <: - OrdinaryDiffEqMutableCache + MPRKCache tmp::uType P::PType P2::PType @@ -184,7 +184,7 @@ struct SSPMPRK22Cache{uType, PType, tabType, F} <: end struct SSPMPRK22ConservativeCache{uType, PType, tabType, F} <: - OrdinaryDiffEqMutableCache + MPRKCache tmp::uType P::PType P2::PType @@ -487,10 +487,12 @@ function alg_cache(alg::SSPMPRK43, u, rate_prototype, ::Type{uEltypeNoUnits}, 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) + # small_constant = floatmin(uEltypeNoUnits) + small_constant = 1e-50 SSPMPRK43ConstantCache(n1, n2, z, η1, η2, η3, η4, η5, η6, s, α10, α20, α21, α30, α31, α32, β10, β20, β21, β30, - β31, β32, c3, floatmin(uEltypeNoUnits)) + β31, β32, c3, small_constant) end function initialize!(integrator, cache::SSPMPRK43ConstantCache) @@ -619,7 +621,7 @@ function perform_step!(integrator, cache::SSPMPRK43ConstantCache, repeat_step = integrator.u = u end -struct SSPMPRK43Cache{uType, PType, tabType, F} <: OrdinaryDiffEqMutableCache +struct SSPMPRK43Cache{uType, PType, tabType, F} <: MPRKCache tmp::uType tmp2::uType P::PType @@ -634,7 +636,7 @@ struct SSPMPRK43Cache{uType, PType, tabType, F} <: OrdinaryDiffEqMutableCache linsolve::F end -struct SSPMPRK43ConservativeCache{uType, PType, tabType, F} <: OrdinaryDiffEqMutableCache +struct SSPMPRK43ConservativeCache{uType, PType, tabType, F} <: MPRKCache tmp::uType tmp2::uType P::PType @@ -652,9 +654,11 @@ function alg_cache(alg::SSPMPRK43, u, rate_prototype, ::Type{uEltypeNoUnits}, uprev, uprev2, f, t, dt, reltol, p, calck, ::Val{true}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} 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) + # small_constant = floatmin(uEltypeNoUnits) + small_constant = 1e-50 tab = SSPMPRK43ConstantCache(n1, n2, z, η1, η2, η3, η4, η5, η6, s, α10, α20, α21, α30, α31, α32, - β10, β20, β21, β30, β31, β32, c3, floatmin(uEltypeNoUnits)) + β10, β20, β21, β30, β31, β32, c3, small_constant) tmp = zero(u) tmp2 = zero(u) P = p_prototype(u, f) From 85cee209e268f22fe2a0a8bd26deef1ed6ab07f7 Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Mon, 1 Jul 2024 14:28:02 +0200 Subject: [PATCH 17/48] Added interpolation.jl --- src/PositiveIntegrators.jl | 10 ++- src/interpolation.jl | 80 ++++++++++++++++++++++++ src/mprk.jl | 121 ------------------------------------- 3 files changed, 84 insertions(+), 127 deletions(-) create mode 100644 src/interpolation.jl diff --git a/src/PositiveIntegrators.jl b/src/PositiveIntegrators.jl index 4bc31c53..d710bc37 100644 --- a/src/PositiveIntegrators.jl +++ b/src/PositiveIntegrators.jl @@ -61,15 +61,13 @@ include("proddest.jl") abstract type MPRKCache <: OrdinaryDiffEqMutableCache end -#TODO: sspmprk.jl must be included before mprk.jl since the interpolation specializations -# at the end of mprk.jl need to know SSPMPRK22ConstantCache. -# It would probably be a good idea to separate schemes and interpolations. -# +# modified Patankar-Runge-Kutta (MPRK) methods +include("mprk.jl") + # modified Patankar-Runge-Kutta based on the SSP formulation of RK methods (SSPMPRK) include("sspmprk.jl") -# modified Patankar-Runge-Kutta (MPRK) methods -include("mprk.jl") +include("interpolation.jl") # predefined PDS problems include("PDSProblemLibrary.jl") diff --git a/src/interpolation.jl b/src/interpolation.jl new file mode 100644 index 00000000..7a551804 --- /dev/null +++ b/src/interpolation.jl @@ -0,0 +1,80 @@ +# Linear interpolations +@muladd @inline function linear_interpolant(Θ, dt, u0, u1, idxs::Nothing, T::Type{Val{0}}) + Θm1 = (1 - Θ) + @.. broadcast=false Θm1 * u0+Θ * u1 +end + +@muladd @inline function linear_interpolant(Θ, dt, u0, u1, idxs, T::Type{Val{0}}) + Θm1 = (1 - Θ) + @.. broadcast=false Θm1 * u0[idxs]+Θ * u1[idxs] +end + +@muladd @inline function linear_interpolant!(out, Θ, dt, u0, u1, idxs::Nothing, + T::Type{Val{0}}) + Θm1 = (1 - Θ) + @.. broadcast=false out=Θm1 * u0 + Θ * u1 + out +end + +@muladd @inline function linear_interpolant!(out, Θ, dt, u0, u1, idxs, T::Type{Val{0}}) + Θm1 = (1 - Θ) + @views @.. broadcast=false out=Θm1 * u0[idxs] + Θ * u1[idxs] + out +end + +@inline function linear_interpolant(Θ, dt, u0, u1, idxs::Nothing, T::Type{Val{1}}) + @.. broadcast=false (u1 - u0)/dt +end + +@inline function linear_interpolant(Θ, dt, u0, u1, idxs, T::Type{Val{1}}) + @.. broadcast=false (u1[idxs] - u0[idxs])/dt +end + +@inline function linear_interpolant!(out, Θ, dt, u0, u1, idxs::Nothing, T::Type{Val{1}}) + @.. broadcast=false out=(u1 - u0) / dt + out +end + +@inline function linear_interpolant!(out, Θ, dt, u0, u1, idxs, T::Type{Val{1}}) + @views @.. broadcast=false out=(u1[idxs] - u0[idxs]) / dt + out +end + +####################################################################################### +# interpolation specializations +function interp_summary(::Type{cacheType}, + dense::Bool) where {cacheType <: MPRKCache} + "1st order linear" +end + +function _ode_interpolant(Θ, dt, u0, u1, k, + cache::MPRKCache, + idxs, # Optionally specialize for ::Nothing and others + T::Type{Val{0}}, + differential_vars::Nothing) + linear_interpolant(Θ, dt, u0, u1, idxs, T) +end + +function _ode_interpolant!(out, Θ, dt, u0, u1, k, + cache::MPRKCache, + idxs, # Optionally specialize for ::Nothing and others + T::Type{Val{0}}, + differential_vars::Nothing) + linear_interpolant!(out, Θ, dt, u0, u1, idxs, T) +end + +function _ode_interpolant(Θ, dt, u0, u1, k, + cache::MPRKCache, + idxs, # Optionally specialize for ::Nothing and others + T::Type{Val{1}}, + differential_vars::Nothing) + linear_interpolant(Θ, dt, u0, u1, idxs, T) +end + +function _ode_interpolant!(out, Θ, dt, u0, u1, k, + cache::MPRKCache, + idxs, # Optionally specialize for ::Nothing and others + T::Type{Val{1}}, + differential_vars::Nothing) + linear_interpolant!(out, Θ, dt, u0, u1, idxs, T) +end diff --git a/src/mprk.jl b/src/mprk.jl index 348c627f..a78a8caf 100644 --- a/src/mprk.jl +++ b/src/mprk.jl @@ -114,49 +114,6 @@ function build_mprk_matrix!(M::Tridiagonal, P::Tridiagonal, σ, dt, d = nothing) return nothing end -##################################################################### -# Linear interpolations -@muladd @inline function linear_interpolant(Θ, dt, u0, u1, idxs::Nothing, T::Type{Val{0}}) - Θm1 = (1 - Θ) - @.. broadcast=false Θm1 * u0+Θ * u1 -end - -@muladd @inline function linear_interpolant(Θ, dt, u0, u1, idxs, T::Type{Val{0}}) - Θm1 = (1 - Θ) - @.. broadcast=false Θm1 * u0[idxs]+Θ * u1[idxs] -end - -@muladd @inline function linear_interpolant!(out, Θ, dt, u0, u1, idxs::Nothing, - T::Type{Val{0}}) - Θm1 = (1 - Θ) - @.. broadcast=false out=Θm1 * u0 + Θ * u1 - out -end - -@muladd @inline function linear_interpolant!(out, Θ, dt, u0, u1, idxs, T::Type{Val{0}}) - Θm1 = (1 - Θ) - @views @.. broadcast=false out=Θm1 * u0[idxs] + Θ * u1[idxs] - out -end - -@inline function linear_interpolant(Θ, dt, u0, u1, idxs::Nothing, T::Type{Val{1}}) - @.. broadcast=false (u1 - u0)/dt -end - -@inline function linear_interpolant(Θ, dt, u0, u1, idxs, T::Type{Val{1}}) - @.. broadcast=false (u1[idxs] - u0[idxs])/dt -end - -@inline function linear_interpolant!(out, Θ, dt, u0, u1, idxs::Nothing, T::Type{Val{1}}) - @.. broadcast=false out=(u1 - u0) / dt - out -end - -@inline function linear_interpolant!(out, Θ, dt, u0, u1, idxs, T::Type{Val{1}}) - @views @.. broadcast=false out=(u1[idxs] - u0[idxs]) / dt - out -end - ### MPE ##################################################################################### """ MPE([linsolve = ...]) @@ -356,45 +313,6 @@ function perform_step!(integrator, cache::MPEConservativeCache, repeat_step = fa integrator.stats.nsolve += 1 end -# interpolation specializations -function interp_summary(::Type{cacheType}, - dense::Bool) where { - cacheType <: Union{MPEConstantCache, MPECache}} - "1st order linear" -end - -function _ode_interpolant(Θ, dt, u0, u1, k, - cache::Union{MPEConstantCache, MPECache}, - idxs, # Optionally specialize for ::Nothing and others - T::Type{Val{0}}, - differential_vars::Nothing) - linear_interpolant(Θ, dt, u0, u1, idxs, T) -end - -function _ode_interpolant!(out, Θ, dt, u0, u1, k, - cache::Union{MPEConstantCache, MPECache}, - idxs, # Optionally specialize for ::Nothing and others - T::Type{Val{0}}, - differential_vars::Nothing) - linear_interpolant!(out, Θ, dt, u0, u1, idxs, T) -end - -function _ode_interpolant(Θ, dt, u0, u1, k, - cache::Union{MPEConstantCache, MPECache}, - idxs, # Optionally specialize for ::Nothing and others - T::Type{Val{1}}, - differential_vars::Nothing) - linear_interpolant(Θ, dt, u0, u1, idxs, T) -end - -function _ode_interpolant!(out, Θ, dt, u0, u1, k, - cache::Union{MPEConstantCache, MPECache}, - idxs, # Optionally specialize for ::Nothing and others - T::Type{Val{1}}, - differential_vars::Nothing) - linear_interpolant!(out, Θ, dt, u0, u1, idxs, T) -end - ### MPRK22 ##################################################################################### """ MPRK22(α; [linsolve = ...]) @@ -1344,42 +1262,3 @@ function perform_step!(integrator, cache::MPRK43ConservativeCache, repeat_step = end ######################################################################################## - -####################################################################################### -# interpolation specializations -function interp_summary(::Type{cacheType}, - dense::Bool) where {cacheType <: MPRKCache} - "1st order linear" -end - -function _ode_interpolant(Θ, dt, u0, u1, k, - cache::MPRKCache, - idxs, # Optionally specialize for ::Nothing and others - T::Type{Val{0}}, - differential_vars::Nothing) - linear_interpolant(Θ, dt, u0, u1, idxs, T) -end - -function _ode_interpolant!(out, Θ, dt, u0, u1, k, - cache::MPRKCache, - idxs, # Optionally specialize for ::Nothing and others - T::Type{Val{0}}, - differential_vars::Nothing) - linear_interpolant!(out, Θ, dt, u0, u1, idxs, T) -end - -function _ode_interpolant(Θ, dt, u0, u1, k, - cache::MPRKCache, - idxs, # Optionally specialize for ::Nothing and others - T::Type{Val{1}}, - differential_vars::Nothing) - linear_interpolant(Θ, dt, u0, u1, idxs, T) -end - -function _ode_interpolant!(out, Θ, dt, u0, u1, k, - cache::MPRKCache, - idxs, # Optionally specialize for ::Nothing and others - T::Type{Val{1}}, - differential_vars::Nothing) - linear_interpolant!(out, Θ, dt, u0, u1, idxs, T) -end From fffbbac092ee8d6eec73a9732ec2da1a9ce065c4 Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Mon, 1 Jul 2024 15:41:18 +0200 Subject: [PATCH 18/48] Removed abstract super type. Didn't work as expected. --- src/PositiveIntegrators.jl | 1 + src/interpolation.jl | 32 ++++++++++++++++++++++++++----- src/mprk.jl | 39 +++++++++++++++++++++++--------------- src/sspmprk.jl | 8 ++++---- 4 files changed, 56 insertions(+), 24 deletions(-) diff --git a/src/PositiveIntegrators.jl b/src/PositiveIntegrators.jl index d710bc37..6cdf7ad7 100644 --- a/src/PositiveIntegrators.jl +++ b/src/PositiveIntegrators.jl @@ -60,6 +60,7 @@ export prob_pds_linmod, prob_pds_linmod_inplace, prob_pds_nonlinmod, include("proddest.jl") abstract type MPRKCache <: OrdinaryDiffEqMutableCache end +abstract type MPRKConstantCache <: OrdinaryDiffEqConstantCache end # modified Patankar-Runge-Kutta (MPRK) methods include("mprk.jl") diff --git a/src/interpolation.jl b/src/interpolation.jl index 7a551804..d64656d3 100644 --- a/src/interpolation.jl +++ b/src/interpolation.jl @@ -43,12 +43,22 @@ end ####################################################################################### # interpolation specializations function interp_summary(::Type{cacheType}, - dense::Bool) where {cacheType <: MPRKCache} + dense::Bool) where { + cacheType <: + Union{MPEConstantCache, MPECache, + MPRK22ConstantCache, MPRK22Cache, + MPRK43ConstantCache, MPRK43Cache, + SSPMPRK22ConstantCache, SSPMPRK22Cache, + SSPMPRK43ConstantCache, SSPMPRK43Cache}} "1st order linear" end function _ode_interpolant(Θ, dt, u0, u1, k, - cache::MPRKCache, + cache::Union{MPEConstantCache, MPECache, MPRK22ConstantCache, + MPRK22Cache, + MPRK43ConstantCache, MPRK43Cache, + SSPMPRK22ConstantCache, SSPMPRK22Cache, + SSPMPRK43ConstantCache, SSPMPRK43Cache}, idxs, # Optionally specialize for ::Nothing and others T::Type{Val{0}}, differential_vars::Nothing) @@ -56,7 +66,11 @@ function _ode_interpolant(Θ, dt, u0, u1, k, end function _ode_interpolant!(out, Θ, dt, u0, u1, k, - cache::MPRKCache, + cache::Union{MPEConstantCache, MPECache, MPRK22ConstantCache, + MPRK22Cache, + MPRK43ConstantCache, MPRK43Cache, + SSPMPRK22ConstantCache, SSPMPRK22Cache, + SSPMPRK43ConstantCache, SSPMPRK43Cache}, idxs, # Optionally specialize for ::Nothing and others T::Type{Val{0}}, differential_vars::Nothing) @@ -64,7 +78,11 @@ function _ode_interpolant!(out, Θ, dt, u0, u1, k, end function _ode_interpolant(Θ, dt, u0, u1, k, - cache::MPRKCache, + cache::Union{MPEConstantCache, MPECache, MPRK22ConstantCache, + MPRK22Cache, + MPRK43ConstantCache, MPRK43Cache, + SSPMPRK22ConstantCache, SSPMPRK22Cache, + SSPMPRK43ConstantCache, SSPMPRK43Cache}, idxs, # Optionally specialize for ::Nothing and others T::Type{Val{1}}, differential_vars::Nothing) @@ -72,7 +90,11 @@ function _ode_interpolant(Θ, dt, u0, u1, k, end function _ode_interpolant!(out, Θ, dt, u0, u1, k, - cache::MPRKCache, + cache::Union{MPEConstantCache, MPECache, MPRK22ConstantCache, + MPRK22Cache, + MPRK43ConstantCache, MPRK43Cache, + SSPMPRK22ConstantCache, SSPMPRK22Cache, + SSPMPRK43ConstantCache, SSPMPRK43Cache}, idxs, # Optionally specialize for ::Nothing and others T::Type{Val{1}}, differential_vars::Nothing) diff --git a/src/mprk.jl b/src/mprk.jl index a78a8caf..a2f01186 100644 --- a/src/mprk.jl +++ b/src/mprk.jl @@ -116,7 +116,7 @@ end ### MPE ##################################################################################### """ - MPE([linsolve = ...]) + MPE([linsolve = ..., small_constant = ...]) The first-order modified Patankar-Euler algorithm for production-destruction systems. This one-step, one-stage method is first-order accurate, unconditionally positivity-preserving, and @@ -133,6 +133,8 @@ The modified Patankar-Euler method requires the special structure of a 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. ## References @@ -142,12 +144,13 @@ as keyword argument `linsolve`. Applied Numerical Mathematics 47.1 (2003): 1-30. [DOI: 10.1016/S0168-9274(03)00101-6](https://doi.org/10.1016/S0168-9274(03)00101-6) """ -struct MPE{F} <: OrdinaryDiffEqAlgorithm +struct MPE{F, T} <: OrdinaryDiffEqAlgorithm linsolve::F + small_constant::T end -function MPE(; linsolve = LUFactorization()) - MPE(linsolve) +function MPE(; linsolve = LUFactorization(), small_constant = floatmin()) + MPE(linsolve, small_constant) end alg_order(::MPE) = 1 @@ -165,7 +168,7 @@ function alg_cache(alg::MPE, u, rate_prototype, ::Type{uEltypeNoUnits}, if !(f isa PDSFunction || f isa ConservativePDSFunction) throw(ArgumentError("MPE can only be applied to production-destruction systems")) end - MPEConstantCache(floatmin(uEltypeNoUnits)) + MPEConstantCache(convert(uEltypeNoUnits, alg.small_constant)) end function initialize!(integrator, cache::MPEConstantCache) @@ -204,7 +207,7 @@ function perform_step!(integrator, cache::MPEConstantCache, repeat_step = false) integrator.u = u end -struct MPECache{PType, uType, tabType, F} <: MPRKCache +struct MPECache{PType, uType, tabType, F} <: OrdinaryDiffEqMutableCache P::PType D::uType σ::uType @@ -213,7 +216,7 @@ struct MPECache{PType, uType, tabType, F} <: MPRKCache linsolve::F end -struct MPEConservativeCache{PType, uType, tabType, F} <: MPRKCache +struct MPEConservativeCache{PType, uType, tabType, F} <: OrdinaryDiffEqMutableCache P::PType σ::uType tab::tabType @@ -227,7 +230,7 @@ function alg_cache(alg::MPE, u, rate_prototype, ::Type{uEltypeNoUnits}, ::Val{true}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} P = p_prototype(u, f) σ = zero(u) - tab = MPEConstantCache(floatmin(uEltypeNoUnits)) + tab = MPEConstantCache(alg.small_constant) if f isa ConservativePDSFunction # We use P to store the evaluation of the PDS @@ -315,7 +318,7 @@ end ### MPRK22 ##################################################################################### """ - MPRK22(α; [linsolve = ...]) + MPRK22(α; [linsolve = ..., small_constant = ...]) The second-order modified Patankar-Runge-Kutta algorithm for production-destruction systems. This one-step, two-stage method is @@ -334,6 +337,8 @@ This modified Patankar-Runge-Kutta method requires the special structure of a 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. ## References @@ -483,7 +488,7 @@ function perform_step!(integrator, cache::MPRK22ConstantCache, repeat_step = fal end struct MPRK22Cache{uType, PType, tabType, F} <: - MPRKCache + OrdinaryDiffEqMutableCache tmp::uType P::PType P2::PType @@ -495,7 +500,7 @@ struct MPRK22Cache{uType, PType, tabType, F} <: end struct MPRK22ConservativeCache{uType, PType, tabType, F} <: - MPRKCache + OrdinaryDiffEqMutableCache tmp::uType P::PType P2::PType @@ -685,7 +690,7 @@ end ### MPRK43 ##################################################################################### """ - MPRK43I(α, β; [linsolve = ...]) + MPRK43I(α, β; [linsolve = ..., small_constant = ...]) A family of third-order modified Patankar-Runge-Kutta schemes for (conservative) production-destruction systems, which is based on the two-parameter family of third order explicit Runge--Kutta schemes. @@ -705,6 +710,8 @@ These modified Patankar-Runge-Kutta methods require the special structure of a 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. ## References @@ -770,7 +777,7 @@ function get_constant_parameters(alg::MPRK43I) end """ - MPRK43II(γ; [linsolve = ...]) + MPRK43II(γ; [linsolve = ..., small_constant = ...]) A family of third-order modified Patankar-Runge-Kutta schemes for (conservative) production-destruction systems, which is based on the one-parameter family of third order explicit Runge--Kutta schemes with @@ -791,6 +798,8 @@ These modified Patankar-Runge-Kutta methods require the special structure of a 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. ## References @@ -991,7 +1000,7 @@ function perform_step!(integrator, cache::MPRK43ConstantCache, repeat_step = fal integrator.u = u end -struct MPRK43Cache{uType, PType, tabType, F} <: MPRKCache +struct MPRK43Cache{uType, PType, tabType, F} <: OrdinaryDiffEqMutableCache tmp::uType tmp2::uType P::PType @@ -1005,7 +1014,7 @@ struct MPRK43Cache{uType, PType, tabType, F} <: MPRKCache linsolve::F end -struct MPRK43ConservativeCache{uType, PType, tabType, F} <: MPRKCache +struct MPRK43ConservativeCache{uType, PType, tabType, F} <: OrdinaryDiffEqMutableCache tmp::uType tmp2::uType P::PType diff --git a/src/sspmprk.jl b/src/sspmprk.jl index 875ba3f4..9e7a4124 100644 --- a/src/sspmprk.jl +++ b/src/sspmprk.jl @@ -172,7 +172,7 @@ function perform_step!(integrator, cache::SSPMPRK22ConstantCache, repeat_step = end struct SSPMPRK22Cache{uType, PType, tabType, F} <: - MPRKCache + OrdinaryDiffEqMutableCache tmp::uType P::PType P2::PType @@ -184,7 +184,7 @@ struct SSPMPRK22Cache{uType, PType, tabType, F} <: end struct SSPMPRK22ConservativeCache{uType, PType, tabType, F} <: - MPRKCache + OrdinaryDiffEqMutableCache tmp::uType P::PType P2::PType @@ -621,7 +621,7 @@ function perform_step!(integrator, cache::SSPMPRK43ConstantCache, repeat_step = integrator.u = u end -struct SSPMPRK43Cache{uType, PType, tabType, F} <: MPRKCache +struct SSPMPRK43Cache{uType, PType, tabType, F} <: OrdinaryDiffEqMutableCache tmp::uType tmp2::uType P::PType @@ -636,7 +636,7 @@ struct SSPMPRK43Cache{uType, PType, tabType, F} <: MPRKCache linsolve::F end -struct SSPMPRK43ConservativeCache{uType, PType, tabType, F} <: MPRKCache +struct SSPMPRK43ConservativeCache{uType, PType, tabType, F} <: OrdinaryDiffEqMutableCache tmp::uType tmp2::uType P::PType From 2de46131b97c899e8b6ea42032c62c7c59b34f14 Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Mon, 1 Jul 2024 15:54:06 +0200 Subject: [PATCH 19/48] small_constant is now an optional parameter of the algorithms --- src/mprk.jl | 24 ++++++++++++++++-------- src/sspmprk.jl | 29 ++++++++++++++++++----------- 2 files changed, 34 insertions(+), 19 deletions(-) diff --git a/src/mprk.jl b/src/mprk.jl index a2f01186..0571b838 100644 --- a/src/mprk.jl +++ b/src/mprk.jl @@ -364,10 +364,12 @@ to avoid divisions by zero. struct MPRK22{T, F} <: OrdinaryDiffEqAdaptiveAlgorithm alpha::T linsolve::F + small_constant::T end -function MPRK22(alpha; linsolve = LUFactorization()) - MPRK22{typeof(alpha), typeof(linsolve)}(alpha, linsolve) +function MPRK22(alpha; linsolve = LUFactorization(), + small_constant = floatmin(typeof(alpha))) + MPRK22{typeof(alpha), typeof(linsolve)}(alpha, linsolve, small_constant) end alg_order(::MPRK22) = 2 @@ -406,7 +408,8 @@ function alg_cache(alg::MPRK22, u, rate_prototype, ::Type{uEltypeNoUnits}, end a21, b1, b2 = get_constant_parameters(alg) - MPRK22ConstantCache(a21, b1, b2, floatmin(uEltypeNoUnits)) + small_constant = alg.small_constant + MPRK22ConstantCache(a21, b1, b2, small_constant) end function initialize!(integrator, cache::MPRK22ConstantCache) @@ -515,7 +518,8 @@ function alg_cache(alg::MPRK22, u, rate_prototype, ::Type{uEltypeNoUnits}, uprev, uprev2, f, t, dt, reltol, p, calck, ::Val{true}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} a21, b1, b2 = get_constant_parameters(alg) - tab = MPRK22ConstantCache(a21, b1, b2, floatmin(uEltypeNoUnits)) + small_constant = alg.small_constant + tab = MPRK22ConstantCache(a21, b1, b2, small_constant) tmp = zero(u) P = p_prototype(u, f) # We use P2 to store the last evaluation of the PDS @@ -725,10 +729,12 @@ struct MPRK43I{T, F} <: OrdinaryDiffEqAdaptiveAlgorithm alpha::T beta::T linsolve::F + small_constant::T end -function MPRK43I(alpha, beta; linsolve = LUFactorization()) - MPRK43I{typeof(alpha), typeof(linsolve)}(alpha, beta, linsolve) +function MPRK43I(alpha, beta; linsolve = LUFactorization(), + small_constant = floatmin(typeof(alpha))) + MPRK43I{typeof(alpha), typeof(linsolve)}(alpha, beta, linsolve, small_constant) end alg_order(::MPRK43I) = 3 @@ -873,8 +879,9 @@ function alg_cache(alg::Union{MPRK43I, MPRK43II}, u, rate_prototype, ::Type{uElt throw(ArgumentError("MPRK43 can only be applied to production-destruction systems")) end a21, a31, a32, b1, b2, b3, c2, c3, beta1, beta2, q1, q2 = get_constant_parameters(alg) + small_constant = alg.small_constant MPRK43ConstantCache(a21, a31, a32, b1, b2, b3, c2, c3, - beta1, beta2, q1, q2, floatmin(uEltypeNoUnits)) + beta1, beta2, q1, q2, small_constant) end function initialize!(integrator, cache::MPRK43ConstantCache) @@ -1031,8 +1038,9 @@ function alg_cache(alg::Union{MPRK43I, MPRK43II}, u, rate_prototype, ::Type{uElt uprev, uprev2, f, t, dt, reltol, p, calck, ::Val{true}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} a21, a31, a32, b1, b2, b3, c2, c3, beta1, beta2, q1, q2 = get_constant_parameters(alg) + small_constant = alg.small_constant tab = MPRK43ConstantCache(a21, a31, a32, b1, b2, b3, c2, c3, - beta1, beta2, q1, q2, floatmin(uEltypeNoUnits)) + beta1, beta2, q1, q2, small_constant) tmp = zero(u) tmp2 = zero(u) P = p_prototype(u, f) diff --git a/src/sspmprk.jl b/src/sspmprk.jl index 9e7a4124..a340bc8d 100644 --- a/src/sspmprk.jl +++ b/src/sspmprk.jl @@ -21,6 +21,8 @@ This modified Patankar-Runge-Kutta method requires the special structure of a 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. ## References @@ -38,10 +40,12 @@ struct SSPMPRK22{T, F} <: OrdinaryDiffEqAdaptiveAlgorithm alpha::T beta::T linsolve::F + small_constant::T end -function SSPMPRK22(alpha, beta; linsolve = LUFactorization()) - SSPMPRK22{typeof(alpha), typeof(linsolve)}(alpha, beta, linsolve) +function SSPMPRK22(alpha, beta; linsolve = LUFactorization(), + small_constant = floatmin(alpha)) + SSPMPRK22{typeof(alpha), typeof(linsolve)}(alpha, beta, linsolve, small_constant) end alg_order(::SSPMPRK22) = 2 @@ -91,7 +95,8 @@ function alg_cache(alg::SSPMPRK22, u, rate_prototype, ::Type{uEltypeNoUnits}, end a21, a10, a20, b10, b20, b21, s = get_constant_parameters(alg) - SSPMPRK22ConstantCache(a21, a10, a20, b10, b20, b21, s, floatmin(uEltypeNoUnits)) + small_constant = alg.small_constant + SSPMPRK22ConstantCache(a21, a10, a20, b10, b20, b21, s, small_constant) end function initialize!(integrator, cache::SSPMPRK22ConstantCache) @@ -199,7 +204,8 @@ function alg_cache(alg::SSPMPRK22, u, rate_prototype, ::Type{uEltypeNoUnits}, uprev, uprev2, f, t, dt, reltol, p, calck, ::Val{true}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} a21, a10, a20, b10, b20, b21, s = get_constant_parameters(alg) - tab = SSPMPRK22ConstantCache(a21, a10, a20, b10, b20, b21, s, floatmin(uEltypeNoUnits)) + small_constant = alg.small_constant + tab = SSPMPRK22ConstantCache(a21, a10, a20, b10, b20, b21, s, small_constant) tmp = zero(u) P = p_prototype(u, f) # We use P2 to store the last evaluation of the PDS @@ -392,6 +398,8 @@ This modified Patankar-Runge-Kutta method requires the special structure of a 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. ## References @@ -405,12 +413,13 @@ as keyword argument `linsolve`. ESAIM: Mathematical Modelling and Numerical Analysis 57 (2023):1063–1086 [DOI: 10.1051/m2an/2023005](https://doi.org/10.1051/m2an/2023005) """ -struct SSPMPRK43{F} <: OrdinaryDiffEqAlgorithm +struct SSPMPRK43{F, T} <: OrdinaryDiffEqAlgorithm linsolve::F + small_constant::T end -function SSPMPRK43(; linsolve = LUFactorization()) - SSPMPRK43{typeof(linsolve)}(linsolve) +function SSPMPRK43(; linsolve = LUFactorization(), small_constant = 1e-50) + SSPMPRK43{typeof(linsolve), typeof(small_constant)}(linsolve, small_constant) end alg_order(::SSPMPRK43) = 3 @@ -487,8 +496,7 @@ function alg_cache(alg::SSPMPRK43, u, rate_prototype, ::Type{uEltypeNoUnits}, 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) - # small_constant = floatmin(uEltypeNoUnits) - small_constant = 1e-50 + small_constant = alg.small_constant SSPMPRK43ConstantCache(n1, n2, z, η1, η2, η3, η4, η5, η6, s, α10, α20, α21, α30, α31, α32, β10, β20, β21, β30, @@ -654,8 +662,7 @@ function alg_cache(alg::SSPMPRK43, u, rate_prototype, ::Type{uEltypeNoUnits}, uprev, uprev2, f, t, dt, reltol, p, calck, ::Val{true}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} 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) - # small_constant = floatmin(uEltypeNoUnits) - small_constant = 1e-50 + small_constant = alg.small_constant tab = SSPMPRK43ConstantCache(n1, n2, z, η1, η2, η3, η4, η5, η6, s, α10, α20, α21, α30, α31, α32, β10, β20, β21, β30, β31, β32, c3, small_constant) From ef68a75592e01afeb8e3f3182be28035177edb45 Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Mon, 1 Jul 2024 16:00:33 +0200 Subject: [PATCH 20/48] Removed additional abstract types --- src/PositiveIntegrators.jl | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/PositiveIntegrators.jl b/src/PositiveIntegrators.jl index 6cdf7ad7..f5fd7c26 100644 --- a/src/PositiveIntegrators.jl +++ b/src/PositiveIntegrators.jl @@ -59,9 +59,6 @@ export prob_pds_linmod, prob_pds_linmod_inplace, prob_pds_nonlinmod, # production-destruction systems include("proddest.jl") -abstract type MPRKCache <: OrdinaryDiffEqMutableCache end -abstract type MPRKConstantCache <: OrdinaryDiffEqConstantCache end - # modified Patankar-Runge-Kutta (MPRK) methods include("mprk.jl") From 2ee9a8b9df5af8711ec2c9cde92bea8419b9b855 Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Mon, 1 Jul 2024 16:08:56 +0200 Subject: [PATCH 21/48] Added small_constant to MPRK43II --- src/mprk.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/mprk.jl b/src/mprk.jl index 0571b838..d7e7f99d 100644 --- a/src/mprk.jl +++ b/src/mprk.jl @@ -818,10 +818,11 @@ to avoid divisions by zero. struct MPRK43II{T, F} <: OrdinaryDiffEqAdaptiveAlgorithm gamma::T linsolve::F + small_constant::T end -function MPRK43II(gamma; linsolve = LUFactorization()) - MPRK43II{typeof(gamma), typeof(linsolve)}(gamma, linsolve) +function MPRK43II(gamma; linsolve = LUFactorization(), small_constant = floatmin()) + MPRK43II{typeof(gamma), typeof(linsolve)}(gamma, linsolve, small_constant) end alg_order(::MPRK43II) = 3 From 22943e998639d64036e239bc9a57a95d2d2f32c9 Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Mon, 1 Jul 2024 16:13:59 +0200 Subject: [PATCH 22/48] Increased version number --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 07e341e4..d7d18a7a 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.14" +version = "0.1.15" [deps] FastBroadcast = "7034ab61-46d4-4ed7-9d0f-46aef9175898" From b5bce813d5fb6c3908fc9aa125a8e7dd51258ab3 Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Tue, 2 Jul 2024 09:43:02 +0200 Subject: [PATCH 23/48] Update src/interpolation.jl Co-authored-by: Hendrik Ranocha --- src/interpolation.jl | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/interpolation.jl b/src/interpolation.jl index d64656d3..7e4709ff 100644 --- a/src/interpolation.jl +++ b/src/interpolation.jl @@ -42,6 +42,12 @@ end ####################################################################################### # interpolation specializations +const MPRKCaches = Union{MPEConstantCache, MPECache, MPEConservativeCache, + MPRK22ConstantCache, MPRK22Cache, MPRK22ConservativeCache, + MPRK43ConstantCache, MPRK43Cache, MPRK43ConservativeCache, + SSPMPRK22ConstantCache, SSPMPRK22Cache, SSPMPRK22ConservativeCache, + SSPMPRK43ConstantCache, SSPMPRK43Cache, SSPMPRK43ConservativeCache} + function interp_summary(::Type{cacheType}, dense::Bool) where { cacheType <: From a02dc8b7a55bf04010e088304ee3ac954480481e Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Tue, 2 Jul 2024 09:43:22 +0200 Subject: [PATCH 24/48] Update src/interpolation.jl Co-authored-by: Hendrik Ranocha --- src/interpolation.jl | 8 +------- 1 file changed, 1 insertion(+), 7 deletions(-) diff --git a/src/interpolation.jl b/src/interpolation.jl index 7e4709ff..c3059323 100644 --- a/src/interpolation.jl +++ b/src/interpolation.jl @@ -49,13 +49,7 @@ const MPRKCaches = Union{MPEConstantCache, MPECache, MPEConservativeCache, SSPMPRK43ConstantCache, SSPMPRK43Cache, SSPMPRK43ConservativeCache} function interp_summary(::Type{cacheType}, - dense::Bool) where { - cacheType <: - Union{MPEConstantCache, MPECache, - MPRK22ConstantCache, MPRK22Cache, - MPRK43ConstantCache, MPRK43Cache, - SSPMPRK22ConstantCache, SSPMPRK22Cache, - SSPMPRK43ConstantCache, SSPMPRK43Cache}} + dense::Bool) where {cacheType <: MPRKCaches} "1st order linear" end From 7722762a1c881b0f95af66de2a375edb18724f65 Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Tue, 2 Jul 2024 09:43:41 +0200 Subject: [PATCH 25/48] Update src/interpolation.jl Co-authored-by: Hendrik Ranocha --- src/interpolation.jl | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/src/interpolation.jl b/src/interpolation.jl index c3059323..ca54beff 100644 --- a/src/interpolation.jl +++ b/src/interpolation.jl @@ -54,11 +54,7 @@ function interp_summary(::Type{cacheType}, end function _ode_interpolant(Θ, dt, u0, u1, k, - cache::Union{MPEConstantCache, MPECache, MPRK22ConstantCache, - MPRK22Cache, - MPRK43ConstantCache, MPRK43Cache, - SSPMPRK22ConstantCache, SSPMPRK22Cache, - SSPMPRK43ConstantCache, SSPMPRK43Cache}, + cache:: MPRKCaches, idxs, # Optionally specialize for ::Nothing and others T::Type{Val{0}}, differential_vars::Nothing) From 39d6ab21b0bdc5d47817fe63f4565e5bd5b6bb12 Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Tue, 2 Jul 2024 09:43:53 +0200 Subject: [PATCH 26/48] Update src/interpolation.jl Co-authored-by: Hendrik Ranocha --- src/interpolation.jl | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/src/interpolation.jl b/src/interpolation.jl index ca54beff..170d2267 100644 --- a/src/interpolation.jl +++ b/src/interpolation.jl @@ -62,11 +62,7 @@ function _ode_interpolant(Θ, dt, u0, u1, k, end function _ode_interpolant!(out, Θ, dt, u0, u1, k, - cache::Union{MPEConstantCache, MPECache, MPRK22ConstantCache, - MPRK22Cache, - MPRK43ConstantCache, MPRK43Cache, - SSPMPRK22ConstantCache, SSPMPRK22Cache, - SSPMPRK43ConstantCache, SSPMPRK43Cache}, + cache:: MPRKCaches, idxs, # Optionally specialize for ::Nothing and others T::Type{Val{0}}, differential_vars::Nothing) From f07919e0aedc5e33cb67cff1a057149768640ff2 Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Tue, 2 Jul 2024 09:44:02 +0200 Subject: [PATCH 27/48] Update src/interpolation.jl Co-authored-by: Hendrik Ranocha --- src/interpolation.jl | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/src/interpolation.jl b/src/interpolation.jl index 170d2267..45a4fa9b 100644 --- a/src/interpolation.jl +++ b/src/interpolation.jl @@ -70,11 +70,7 @@ function _ode_interpolant!(out, Θ, dt, u0, u1, k, end function _ode_interpolant(Θ, dt, u0, u1, k, - cache::Union{MPEConstantCache, MPECache, MPRK22ConstantCache, - MPRK22Cache, - MPRK43ConstantCache, MPRK43Cache, - SSPMPRK22ConstantCache, SSPMPRK22Cache, - SSPMPRK43ConstantCache, SSPMPRK43Cache}, + cache:: MPRKCaches, idxs, # Optionally specialize for ::Nothing and others T::Type{Val{1}}, differential_vars::Nothing) From c409ddbe86ea6bcc98122ed189a6e911a0eaa0ec Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Tue, 2 Jul 2024 09:44:37 +0200 Subject: [PATCH 28/48] Update src/interpolation.jl Co-authored-by: Hendrik Ranocha --- src/interpolation.jl | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/src/interpolation.jl b/src/interpolation.jl index 45a4fa9b..1c35e763 100644 --- a/src/interpolation.jl +++ b/src/interpolation.jl @@ -78,11 +78,7 @@ function _ode_interpolant(Θ, dt, u0, u1, k, end function _ode_interpolant!(out, Θ, dt, u0, u1, k, - cache::Union{MPEConstantCache, MPECache, MPRK22ConstantCache, - MPRK22Cache, - MPRK43ConstantCache, MPRK43Cache, - SSPMPRK22ConstantCache, SSPMPRK22Cache, - SSPMPRK43ConstantCache, SSPMPRK43Cache}, + cache:: MPRKCaches, idxs, # Optionally specialize for ::Nothing and others T::Type{Val{1}}, differential_vars::Nothing) From 86a551ffda191009f0957e0c20b4c765489eb2d6 Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Tue, 2 Jul 2024 09:46:37 +0200 Subject: [PATCH 29/48] Update src/sspmprk.jl Co-authored-by: Hendrik Ranocha --- src/sspmprk.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/sspmprk.jl b/src/sspmprk.jl index a340bc8d..219bb7fc 100644 --- a/src/sspmprk.jl +++ b/src/sspmprk.jl @@ -68,7 +68,7 @@ function get_constant_parameters(alg::SSPMPRK22) s = (b20 + b21 + a21 * b10^2) / (b10 * (b20 + b21)) # This should never happen - if !all((a21, a10, a20, b10, b20, b21) .≥ 0) + if any(<(0), (a21, a10, a20, b10, b20, b21)) throw(ArgumentError("SSPMPRK22 requires nonnegative SSP coefficients.")) end return a21, a10, a20, b10, b20, b21, s From 1965ff4aedf16ca5583a1cfcec136cc5fef9ed57 Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Tue, 2 Jul 2024 09:48:29 +0200 Subject: [PATCH 30/48] Update test/runtests.jl Co-authored-by: Hendrik Ranocha --- test/runtests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/runtests.jl b/test/runtests.jl index b4779694..c9299b23 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -64,7 +64,7 @@ end Compute the experimental orders of convergence for given `errors` and time step sizes `dts`. """ -function experimental_order_of_convergence(errors, dts) +function experimental_orders_of_convergence(errors, dts) Base.require_one_based_indexing(errors, dts) @assert length(errors) == length(dts) orders = zeros(eltype(errors), length(errors) - 1) From 25b885d0829fb208109fb7f7c732f82cac6fdb78 Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Tue, 2 Jul 2024 09:49:21 +0200 Subject: [PATCH 31/48] Update test/runtests.jl Co-authored-by: Hendrik Ranocha --- test/runtests.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/test/runtests.jl b/test/runtests.jl index c9299b23..5fc632b2 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -13,7 +13,9 @@ using LinearSolve: RFLUFactorization, LUFactorization, KrylovJL_GMRES using Aqua: Aqua """ - experimental_order_of_convergence(prob, alg, dts; test_time = nothing) + experimental_order_of_convergence(prob, alg, dts; + test_time = nothing, + only_first_index = false) Solve `prob` with `alg` and fixed time steps taken from `dts`, and compute the errors at `test_time`. If`test_time` is not specified the error is computed From dbaa027df2d76d9e3833dbc3bb6aeb1c9853c70b Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Tue, 2 Jul 2024 10:00:07 +0200 Subject: [PATCH 32/48] format --- src/interpolation.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/interpolation.jl b/src/interpolation.jl index 1c35e763..e86cdf84 100644 --- a/src/interpolation.jl +++ b/src/interpolation.jl @@ -54,7 +54,7 @@ function interp_summary(::Type{cacheType}, end function _ode_interpolant(Θ, dt, u0, u1, k, - cache:: MPRKCaches, + cache::MPRKCaches, idxs, # Optionally specialize for ::Nothing and others T::Type{Val{0}}, differential_vars::Nothing) @@ -62,7 +62,7 @@ function _ode_interpolant(Θ, dt, u0, u1, k, end function _ode_interpolant!(out, Θ, dt, u0, u1, k, - cache:: MPRKCaches, + cache::MPRKCaches, idxs, # Optionally specialize for ::Nothing and others T::Type{Val{0}}, differential_vars::Nothing) @@ -70,7 +70,7 @@ function _ode_interpolant!(out, Θ, dt, u0, u1, k, end function _ode_interpolant(Θ, dt, u0, u1, k, - cache:: MPRKCaches, + cache::MPRKCaches, idxs, # Optionally specialize for ::Nothing and others T::Type{Val{1}}, differential_vars::Nothing) @@ -78,7 +78,7 @@ function _ode_interpolant(Θ, dt, u0, u1, k, end function _ode_interpolant!(out, Θ, dt, u0, u1, k, - cache:: MPRKCaches, + cache::MPRKCaches, idxs, # Optionally specialize for ::Nothing and others T::Type{Val{1}}, differential_vars::Nothing) From 3d699d7b7d0ca90249c77dc2a736f10d8269e3ca Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Tue, 2 Jul 2024 10:55:17 +0200 Subject: [PATCH 33/48] fixed typos in experimental_orders_of_convergence --- test/runtests.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index 5fc632b2..8065adec 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -13,7 +13,7 @@ using LinearSolve: RFLUFactorization, LUFactorization, KrylovJL_GMRES using Aqua: Aqua """ - experimental_order_of_convergence(prob, alg, dts; + experimental_orders_of_convergence(prob, alg, dts; test_time = nothing, only_first_index = false) @@ -25,8 +25,8 @@ Return the associated experimental orders of convergence. If `only_first_index == true`, only the first solution component is used to compute the error. """ -function experimental_order_of_convergence(prob, alg, dts; test_time = nothing, - only_first_index = false) +function experimental_orders_of_convergence(prob, alg, dts; test_time = nothing, + only_first_index = false) @assert length(dts) > 1 errors = zeros(eltype(dts), length(dts)) analytic = t -> prob.f.analytic(prob.u0, prob.p, t) @@ -57,11 +57,11 @@ function experimental_order_of_convergence(prob, alg, dts; test_time = nothing, end end - return experimental_order_of_convergence(errors, dts) + return experimental_orders_of_convergence(errors, dts) end """ - experimental_order_of_convergence(errors, dts) + experimental_orders_of_convergence(errors, dts) Compute the experimental orders of convergence for given `errors` and time step sizes `dts`. From d666c1a2ea7e4edd7581b43ac23022a6e9e4550d Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Tue, 2 Jul 2024 11:45:01 +0200 Subject: [PATCH 34/48] fixed typo function calls --- test/runtests.jl | 32 ++++++++++++++++++-------------- 1 file changed, 18 insertions(+), 14 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index 8065adec..6a566632 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -715,7 +715,7 @@ const prob_pds_linmod_nonconservative_inplace = PDSProblem(linmodP!, linmodD!, [ alg = MPRK22(1.0) for prob in problems prob = problems[1] - orders = experimental_order_of_convergence(prob, alg, dts) + orders = experimental_orders_of_convergence_convergence(prob, alg, dts) @test check_order(orders, PositiveIntegrators.alg_order(alg)) test_times = [ @@ -723,13 +723,15 @@ const prob_pds_linmod_nonconservative_inplace = PDSProblem(linmodP!, linmodD!, [ 1.23456789, 1 + 1 / pi, 1 + exp(-1), ] for test_time in test_times - orders = experimental_order_of_convergence(prob, alg, dts; - test_time) + orders = experimental_orders_of_convergence_convergence(prob, alg, + dts; + test_time) @test check_order(orders, PositiveIntegrators.alg_order(alg), atol = 0.2) - orders = experimental_order_of_convergence(prob, alg, dts; - test_time, - only_first_index = true) + orders = experimental_orders_of_convergence_convergence(prob, alg, + dts; + test_time, + only_first_index = true) @test check_order(orders, PositiveIntegrators.alg_order(alg), atol = 0.2) end @@ -747,7 +749,7 @@ const prob_pds_linmod_nonconservative_inplace = PDSProblem(linmodP!, linmodD!, [ @testset "$alg" for alg in algs alg = MPRK22(1.0) for prob in problems - orders = experimental_order_of_convergence(prob, alg, dts) + orders = experimental_orders_of_convergence_convergence(prob, alg, dts) @test check_order(orders, PositiveIntegrators.alg_order(alg)) test_times = [ @@ -755,13 +757,15 @@ const prob_pds_linmod_nonconservative_inplace = PDSProblem(linmodP!, linmodD!, [ 1.23456789, 1 + 1 / pi, 1 + exp(-1), ] for test_time in test_times - orders = experimental_order_of_convergence(prob, alg, dts; - test_time) + orders = experimental_orders_of_convergence_convergence(prob, alg, + dts; + test_time) @test check_order(orders, PositiveIntegrators.alg_order(alg), atol = 0.2) - orders = experimental_order_of_convergence(prob, alg, dts; - test_time, - only_first_index = true) + orders = experimental_orders_of_convergence_convergence(prob, alg, + dts; + test_time, + only_first_index = true) @test check_order(orders, PositiveIntegrators.alg_order(alg), atol = 0.2) end @@ -778,7 +782,7 @@ const prob_pds_linmod_nonconservative_inplace = PDSProblem(linmodP!, linmodD!, [ algs = (MPRK43I(1.0, 0.5), MPRK43I(0.5, 0.75), MPRK43II(0.5), MPRK43II(2.0 / 3.0), SSPMPRK43()) for alg in algs, prob in problems - orders = experimental_order_of_convergence(prob, alg, dts) + orders = experimental_orders_of_convergence_convergence(prob, alg, dts) @test check_order(orders, PositiveIntegrators.alg_order(alg), atol = 0.2) end end @@ -792,7 +796,7 @@ const prob_pds_linmod_nonconservative_inplace = PDSProblem(linmodP!, linmodD!, [ algs = (MPRK43I(1.0, 0.5), MPRK43I(0.5, 0.75), MPRK43II(0.5), MPRK43II(2.0 / 3.0), SSPMPRK43()) for alg in algs, prob in problems - orders = experimental_order_of_convergence(prob, alg, dts) + orders = experimental_orders_of_convergence_convergence(prob, alg, dts) @test check_order(orders, PositiveIntegrators.alg_order(alg), atol = 0.2) end end From 01dc1d9229921655a12376bf67ecfed4b287d1cc Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Tue, 2 Jul 2024 12:19:38 +0200 Subject: [PATCH 35/48] fixed another typo in function calls --- test/runtests.jl | 36 ++++++++++++++++++------------------ 1 file changed, 18 insertions(+), 18 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index 6a566632..b67366b8 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -715,7 +715,7 @@ const prob_pds_linmod_nonconservative_inplace = PDSProblem(linmodP!, linmodD!, [ alg = MPRK22(1.0) for prob in problems prob = problems[1] - orders = experimental_orders_of_convergence_convergence(prob, alg, dts) + orders = experimental_orders_of_convergence(prob, alg, dts) @test check_order(orders, PositiveIntegrators.alg_order(alg)) test_times = [ @@ -723,15 +723,15 @@ const prob_pds_linmod_nonconservative_inplace = PDSProblem(linmodP!, linmodD!, [ 1.23456789, 1 + 1 / pi, 1 + exp(-1), ] for test_time in test_times - orders = experimental_orders_of_convergence_convergence(prob, alg, - dts; - test_time) + orders = experimental_orders_of_convergence(prob, alg, + dts; + test_time) @test check_order(orders, PositiveIntegrators.alg_order(alg), atol = 0.2) - orders = experimental_orders_of_convergence_convergence(prob, alg, - dts; - test_time, - only_first_index = true) + orders = experimental_orders_of_convergence(prob, alg, + dts; + test_time, + only_first_index = true) @test check_order(orders, PositiveIntegrators.alg_order(alg), atol = 0.2) end @@ -749,7 +749,7 @@ const prob_pds_linmod_nonconservative_inplace = PDSProblem(linmodP!, linmodD!, [ @testset "$alg" for alg in algs alg = MPRK22(1.0) for prob in problems - orders = experimental_orders_of_convergence_convergence(prob, alg, dts) + orders = experimental_orders_of_convergence(prob, alg, dts) @test check_order(orders, PositiveIntegrators.alg_order(alg)) test_times = [ @@ -757,15 +757,15 @@ const prob_pds_linmod_nonconservative_inplace = PDSProblem(linmodP!, linmodD!, [ 1.23456789, 1 + 1 / pi, 1 + exp(-1), ] for test_time in test_times - orders = experimental_orders_of_convergence_convergence(prob, alg, - dts; - test_time) + orders = experimental_orders_of_convergence(prob, alg, + dts; + test_time) @test check_order(orders, PositiveIntegrators.alg_order(alg), atol = 0.2) - orders = experimental_orders_of_convergence_convergence(prob, alg, - dts; - test_time, - only_first_index = true) + orders = experimental_orders_of_convergence(prob, alg, + dts; + test_time, + only_first_index = true) @test check_order(orders, PositiveIntegrators.alg_order(alg), atol = 0.2) end @@ -782,7 +782,7 @@ const prob_pds_linmod_nonconservative_inplace = PDSProblem(linmodP!, linmodD!, [ algs = (MPRK43I(1.0, 0.5), MPRK43I(0.5, 0.75), MPRK43II(0.5), MPRK43II(2.0 / 3.0), SSPMPRK43()) for alg in algs, prob in problems - orders = experimental_orders_of_convergence_convergence(prob, alg, dts) + orders = experimental_orders_of_convergence(prob, alg, dts) @test check_order(orders, PositiveIntegrators.alg_order(alg), atol = 0.2) end end @@ -796,7 +796,7 @@ const prob_pds_linmod_nonconservative_inplace = PDSProblem(linmodP!, linmodD!, [ algs = (MPRK43I(1.0, 0.5), MPRK43I(0.5, 0.75), MPRK43II(0.5), MPRK43II(2.0 / 3.0), SSPMPRK43()) for alg in algs, prob in problems - orders = experimental_orders_of_convergence_convergence(prob, alg, dts) + orders = experimental_orders_of_convergence(prob, alg, dts) @test check_order(orders, PositiveIntegrators.alg_order(alg), atol = 0.2) end end From 5caccfce4a48993f927a72b33fa10c88e37f9500 Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Tue, 2 Jul 2024 14:00:57 +0200 Subject: [PATCH 36/48] Update src/mprk.jl Co-authored-by: Hendrik Ranocha --- src/mprk.jl | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/src/mprk.jl b/src/mprk.jl index d7e7f99d..888e52de 100644 --- a/src/mprk.jl +++ b/src/mprk.jl @@ -149,8 +149,15 @@ struct MPE{F, T} <: OrdinaryDiffEqAlgorithm small_constant::T end -function MPE(; linsolve = LUFactorization(), small_constant = floatmin()) - MPE(linsolve, small_constant) +function MPE(; linsolve = LUFactorization(), small_constant = nothing) + if isnothing(small_constant) + small_constant_function = floatmin + elseif small_constant isa Number + small_constant_function = Returns(small_constant) + else # assume small_constant isa Function + small_constant_function = small_constant + end + MPE(linsolve, small_constant_function) end alg_order(::MPE) = 1 From 5801417e149c713042fc0589a3650fecdde70475 Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Tue, 2 Jul 2024 14:06:47 +0200 Subject: [PATCH 37/48] MPE now uses small_constant_function internally to set the parameter small_constant --- src/mprk.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/mprk.jl b/src/mprk.jl index 888e52de..79c9f623 100644 --- a/src/mprk.jl +++ b/src/mprk.jl @@ -133,7 +133,7 @@ The modified Patankar-Euler method requires the special structure of a 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 +You can also choose the parameter `small_constant` which is added to all Patankar-weight denominators to avoid divisions by zero. ## References @@ -175,7 +175,7 @@ function alg_cache(alg::MPE, u, rate_prototype, ::Type{uEltypeNoUnits}, if !(f isa PDSFunction || f isa ConservativePDSFunction) throw(ArgumentError("MPE can only be applied to production-destruction systems")) end - MPEConstantCache(convert(uEltypeNoUnits, alg.small_constant)) + MPEConstantCache(alg.small_constant_function(uEltypeNoUnits)) end function initialize!(integrator, cache::MPEConstantCache) @@ -237,7 +237,7 @@ function alg_cache(alg::MPE, u, rate_prototype, ::Type{uEltypeNoUnits}, ::Val{true}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} P = p_prototype(u, f) σ = zero(u) - tab = MPEConstantCache(alg.small_constant) + tab = MPEConstantCache(alg.small_constant_function(uEltypeNoUnits)) if f isa ConservativePDSFunction # We use P to store the evaluation of the PDS From 5414efbd34fa5a011930a105d8bace2569377198 Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Tue, 2 Jul 2024 14:11:01 +0200 Subject: [PATCH 38/48] MPE now has a field called small_constant_function --- src/mprk.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/mprk.jl b/src/mprk.jl index 79c9f623..ecb2529b 100644 --- a/src/mprk.jl +++ b/src/mprk.jl @@ -146,7 +146,7 @@ to avoid divisions by zero. """ struct MPE{F, T} <: OrdinaryDiffEqAlgorithm linsolve::F - small_constant::T + small_constant_function::T end function MPE(; linsolve = LUFactorization(), small_constant = nothing) From 628d6842bfff2e2a6fb82322e6c960d9aeff73cd Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Tue, 2 Jul 2024 14:40:58 +0200 Subject: [PATCH 39/48] All MPRK schemes internally use small_constant_function --- src/mprk.jl | 58 +++++++++++++++++++++++++++++++++++------------------ 1 file changed, 38 insertions(+), 20 deletions(-) diff --git a/src/mprk.jl b/src/mprk.jl index ecb2529b..0c338752 100644 --- a/src/mprk.jl +++ b/src/mprk.jl @@ -368,15 +368,22 @@ to avoid divisions by zero. Applied Numerical Mathematics 182 (2022): 117-147. [DOI: 10.1016/j.apnum.2022.07.014](https://doi.org/10.1016/j.apnum.2022.07.014) """ -struct MPRK22{T, F} <: OrdinaryDiffEqAdaptiveAlgorithm +struct MPRK22{T, F, T2} <: OrdinaryDiffEqAdaptiveAlgorithm alpha::T linsolve::F - small_constant::T + small_constant_function::T2 end function MPRK22(alpha; linsolve = LUFactorization(), - small_constant = floatmin(typeof(alpha))) - MPRK22{typeof(alpha), typeof(linsolve)}(alpha, linsolve, small_constant) + small_constant = nothing) + if isnothing(small_constant) + small_constant_function = floatmin + elseif small_constant isa Number + small_constant_function = Returns(small_constant) + else # assume small_constant isa Function + small_constant_function = small_constant + end + MPRK22{typeof(alpha), typeof(linsolve)}(alpha, linsolve, small_constant_function) end alg_order(::MPRK22) = 2 @@ -415,8 +422,7 @@ function alg_cache(alg::MPRK22, u, rate_prototype, ::Type{uEltypeNoUnits}, end a21, b1, b2 = get_constant_parameters(alg) - small_constant = alg.small_constant - MPRK22ConstantCache(a21, b1, b2, small_constant) + MPRK22ConstantCache(a21, b1, b2, alg.small_constant_function(uEltypeNoUnits)) end function initialize!(integrator, cache::MPRK22ConstantCache) @@ -525,8 +531,7 @@ function alg_cache(alg::MPRK22, u, rate_prototype, ::Type{uEltypeNoUnits}, uprev, uprev2, f, t, dt, reltol, p, calck, ::Val{true}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} a21, b1, b2 = get_constant_parameters(alg) - small_constant = alg.small_constant - tab = MPRK22ConstantCache(a21, b1, b2, small_constant) + tab = MPRK22ConstantCache(a21, b1, b2, alg.small_constant_function(uEltypeNoUnits)) tmp = zero(u) P = p_prototype(u, f) # We use P2 to store the last evaluation of the PDS @@ -732,16 +737,23 @@ to avoid divisions by zero. BIT Numerical Mathematics 58 (2018): 691–728. [DOI: 10.1007/s10543-018-0705-1](https://doi.org/10.1007/s10543-018-0705-1) """ -struct MPRK43I{T, F} <: OrdinaryDiffEqAdaptiveAlgorithm +struct MPRK43I{T, F, T2} <: OrdinaryDiffEqAdaptiveAlgorithm alpha::T beta::T linsolve::F - small_constant::T + small_constant_function::T2 end function MPRK43I(alpha, beta; linsolve = LUFactorization(), - small_constant = floatmin(typeof(alpha))) - MPRK43I{typeof(alpha), typeof(linsolve)}(alpha, beta, linsolve, small_constant) + small_constant = nothing) + if isnothing(small_constant) + small_constant_function = floatmin + elseif small_constant isa Number + small_constant_function = Returns(small_constant) + else # assume small_constant isa Function + small_constant_function = small_constant + end + MPRK43I{typeof(alpha), typeof(linsolve)}(alpha, beta, linsolve, small_constant_function) end alg_order(::MPRK43I) = 3 @@ -822,14 +834,21 @@ to avoid divisions by zero. BIT Numerical Mathematics 58 (2018): 691–728. [DOI: 10.1007/s10543-018-0705-1](https://doi.org/10.1007/s10543-018-0705-1) """ -struct MPRK43II{T, F} <: OrdinaryDiffEqAdaptiveAlgorithm +struct MPRK43II{T, F, T2} <: OrdinaryDiffEqAdaptiveAlgorithm gamma::T linsolve::F - small_constant::T + small_constant_function::T2 end -function MPRK43II(gamma; linsolve = LUFactorization(), small_constant = floatmin()) - MPRK43II{typeof(gamma), typeof(linsolve)}(gamma, linsolve, small_constant) +function MPRK43II(gamma; linsolve = LUFactorization(), small_constant = nothing) + if isnothing(small_constant) + small_constant_function = floatmin + elseif small_constant isa Number + small_constant_function = Returns(small_constant) + else # assume small_constant isa Function + small_constant_function = small_constant + end + MPRK43II{typeof(gamma), typeof(linsolve)}(gamma, linsolve, small_constant_function) end alg_order(::MPRK43II) = 3 @@ -887,9 +906,8 @@ function alg_cache(alg::Union{MPRK43I, MPRK43II}, u, rate_prototype, ::Type{uElt throw(ArgumentError("MPRK43 can only be applied to production-destruction systems")) end a21, a31, a32, b1, b2, b3, c2, c3, beta1, beta2, q1, q2 = get_constant_parameters(alg) - small_constant = alg.small_constant MPRK43ConstantCache(a21, a31, a32, b1, b2, b3, c2, c3, - beta1, beta2, q1, q2, small_constant) + beta1, beta2, q1, q2, alg.small_constant_function(uEltypeNoUnits)) end function initialize!(integrator, cache::MPRK43ConstantCache) @@ -1046,9 +1064,9 @@ function alg_cache(alg::Union{MPRK43I, MPRK43II}, u, rate_prototype, ::Type{uElt uprev, uprev2, f, t, dt, reltol, p, calck, ::Val{true}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} a21, a31, a32, b1, b2, b3, c2, c3, beta1, beta2, q1, q2 = get_constant_parameters(alg) - small_constant = alg.small_constant tab = MPRK43ConstantCache(a21, a31, a32, b1, b2, b3, c2, c3, - beta1, beta2, q1, q2, small_constant) + beta1, beta2, q1, q2, + alg.small_constant_function(uEltypeNoUnits)) tmp = zero(u) tmp2 = zero(u) P = p_prototype(u, f) From d5c38a7fc3fb2d93ca4e5dc8bad7d1332e26c001 Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Tue, 2 Jul 2024 14:47:47 +0200 Subject: [PATCH 40/48] fixed bug in constructors of MPRK schemes --- src/mprk.jl | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/src/mprk.jl b/src/mprk.jl index 0c338752..aa49afa5 100644 --- a/src/mprk.jl +++ b/src/mprk.jl @@ -383,7 +383,9 @@ function MPRK22(alpha; linsolve = LUFactorization(), else # assume small_constant isa Function small_constant_function = small_constant end - MPRK22{typeof(alpha), typeof(linsolve)}(alpha, linsolve, small_constant_function) + MPRK22{typeof(alpha), typeof(linsolve), typeof(small_constant_function)}(alpha, + linsolve, + small_constant_function) end alg_order(::MPRK22) = 2 @@ -753,7 +755,9 @@ function MPRK43I(alpha, beta; linsolve = LUFactorization(), else # assume small_constant isa Function small_constant_function = small_constant end - MPRK43I{typeof(alpha), typeof(linsolve)}(alpha, beta, linsolve, small_constant_function) + MPRK43I{typeof(alpha), typeof(linsolve), typeof(small_constant_function)}(alpha, beta, + linsolve, + small_constant_function) end alg_order(::MPRK43I) = 3 @@ -848,7 +852,9 @@ function MPRK43II(gamma; linsolve = LUFactorization(), small_constant = nothing) else # assume small_constant isa Function small_constant_function = small_constant end - MPRK43II{typeof(gamma), typeof(linsolve)}(gamma, linsolve, small_constant_function) + MPRK43II{typeof(gamma), typeof(linsolve), typeof(small_constant_function)}(gamma, + linsolve, + small_constant_function) end alg_order(::MPRK43II) = 3 From b74059b71aa0414cb190fc21bd0ba179cc8a9928 Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Tue, 2 Jul 2024 14:58:35 +0200 Subject: [PATCH 41/48] All SSPMPRK schemes internally use small_constant_function --- src/sspmprk.jl | 48 ++++++++++++++++++++++++++++++++---------------- 1 file changed, 32 insertions(+), 16 deletions(-) diff --git a/src/sspmprk.jl b/src/sspmprk.jl index 219bb7fc..f698eb9a 100644 --- a/src/sspmprk.jl +++ b/src/sspmprk.jl @@ -1,6 +1,6 @@ ### SSPMPRK ##################################################################################### """ - SSPMPRK22(α, β; [linsolve = ...]) + SSPMPRK22(α, β; [linsolve = ..., small_constant = ...]) A family of second-order modified Patankar-Runge-Kutta algorithms for production-destruction systems. Each member of this family is a one-step, two-stage method which is @@ -36,16 +36,25 @@ to avoid divisions by zero. ESAIM: Mathematical Modelling and Numerical Analysis 57 (2023):1063–1086 [DOI: 10.1051/m2an/2023005](https://doi.org/10.1051/m2an/2023005) """ -struct SSPMPRK22{T, F} <: OrdinaryDiffEqAdaptiveAlgorithm +struct SSPMPRK22{T, F, T2} <: OrdinaryDiffEqAdaptiveAlgorithm alpha::T beta::T linsolve::F - small_constant::T + small_constant_function::T2 end function SSPMPRK22(alpha, beta; linsolve = LUFactorization(), - small_constant = floatmin(alpha)) - SSPMPRK22{typeof(alpha), typeof(linsolve)}(alpha, beta, linsolve, small_constant) + small_constant = nothing) + if isnothing(small_constant) + small_constant_function = floatmin + elseif small_constant isa Number + small_constant_function = Returns(small_constant) + else # assume small_constant isa Function + small_constant_function = small_constant + end + SSPMPRK22{typeof(alpha), typeof(linsolve), typeof(small_constant_function)}(alpha, beta, + linsolve, + small_constant_function) end alg_order(::SSPMPRK22) = 2 @@ -95,8 +104,8 @@ function alg_cache(alg::SSPMPRK22, u, rate_prototype, ::Type{uEltypeNoUnits}, end a21, a10, a20, b10, b20, b21, s = get_constant_parameters(alg) - small_constant = alg.small_constant - SSPMPRK22ConstantCache(a21, a10, a20, b10, b20, b21, s, small_constant) + SSPMPRK22ConstantCache(a21, a10, a20, b10, b20, b21, s, + alg.small_constant_function(uEltypeNoUnits)) end function initialize!(integrator, cache::SSPMPRK22ConstantCache) @@ -204,8 +213,8 @@ function alg_cache(alg::SSPMPRK22, u, rate_prototype, ::Type{uEltypeNoUnits}, uprev, uprev2, f, t, dt, reltol, p, calck, ::Val{true}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} a21, a10, a20, b10, b20, b21, s = get_constant_parameters(alg) - small_constant = alg.small_constant - tab = SSPMPRK22ConstantCache(a21, a10, a20, b10, b20, b21, s, small_constant) + tab = SSPMPRK22ConstantCache(a21, a10, a20, b10, b20, b21, s, + alg.small_constant_function(uEltypeNoUnits)) tmp = zero(u) P = p_prototype(u, f) # We use P2 to store the last evaluation of the PDS @@ -378,7 +387,7 @@ function perform_step!(integrator, cache::SSPMPRK22ConservativeCache, repeat_ste end """ - SSPMPRK43([linsolve = ...]) + SSPMPRK43([linsolve = ..., small_constant = ...]) A third-order modified Patankar-Runge-Kutta algorithm for production-destruction systems. This scheme is a one-step, two-stage method which is @@ -415,11 +424,19 @@ to avoid divisions by zero. """ struct SSPMPRK43{F, T} <: OrdinaryDiffEqAlgorithm linsolve::F - small_constant::T + small_constant_function::T end function SSPMPRK43(; linsolve = LUFactorization(), small_constant = 1e-50) - SSPMPRK43{typeof(linsolve), typeof(small_constant)}(linsolve, small_constant) + if isnothing(small_constant) + small_constant_function = floatmin + elseif small_constant isa Number + small_constant_function = Returns(small_constant) + else # assume small_constant isa Function + small_constant_function = small_constant + end + SSPMPRK43{typeof(linsolve), typeof(small_constant_function)}(linsolve, + small_constant_function) end alg_order(::SSPMPRK43) = 3 @@ -496,11 +513,10 @@ function alg_cache(alg::SSPMPRK43, u, rate_prototype, ::Type{uEltypeNoUnits}, 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) - small_constant = alg.small_constant SSPMPRK43ConstantCache(n1, n2, z, η1, η2, η3, η4, η5, η6, s, α10, α20, α21, α30, α31, α32, β10, β20, β21, β30, - β31, β32, c3, small_constant) + β31, β32, c3, alg.small_constant_function(uEltypeNoUnits)) end function initialize!(integrator, cache::SSPMPRK43ConstantCache) @@ -662,10 +678,10 @@ function alg_cache(alg::SSPMPRK43, u, rate_prototype, ::Type{uEltypeNoUnits}, uprev, uprev2, f, t, dt, reltol, p, calck, ::Val{true}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} 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) - small_constant = alg.small_constant tab = SSPMPRK43ConstantCache(n1, n2, z, η1, η2, η3, η4, η5, η6, s, α10, α20, α21, α30, α31, α32, - β10, β20, β21, β30, β31, β32, c3, small_constant) + β10, β20, β21, β30, β31, β32, c3, + alg.small_constant_function(uEltypeNoUnits)) tmp = zero(u) tmp2 = zero(u) P = p_prototype(u, f) From 73f64b4a98965c1f1b1a7d8ecd8d48af24132a81 Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Tue, 2 Jul 2024 15:38:11 +0200 Subject: [PATCH 42/48] Corrected comments regarding the error estimator in MPRK22 --- src/mprk.jl | 15 ++++++--------- 1 file changed, 6 insertions(+), 9 deletions(-) diff --git a/src/mprk.jl b/src/mprk.jl index aa49afa5..98e664b5 100644 --- a/src/mprk.jl +++ b/src/mprk.jl @@ -494,9 +494,8 @@ function perform_step!(integrator, cache::MPRK22ConstantCache, repeat_step = fal u = sol.u integrator.stats.nsolve += 1 - # If a21 = 1.0, then σ is the MPE approximation, i.e. a first order approximation - # of the solution, and can be used for error estimation. Moreover, MPE is suited for stiff problems. - # TODO: Find first order approximation if a21≠ 1.0. + # If a21 = 1, then σ is the MPE approximation, i.e. suited for stiff problems. + # Otherwise, this is not clear. tmp = u - σ atmp = calculate_residuals(tmp, uprev, u, integrator.opts.abstol, integrator.opts.reltol, integrator.opts.internalnorm, t) @@ -634,9 +633,8 @@ function perform_step!(integrator, cache::MPRK22Cache, repeat_step = false) integrator.stats.nsolve += 1 # Now σ stores the error estimate - # If a21 = 1.0, then σ is the MPE approximation, i.e. a first order approximation - # of the solution, and can be used for error estimation. Moreover, MPE is suited for stiff problems. - # TODO: Find first order approximation if a21≠ 1.0. + # If a21 = 1, then σ is the MPE approximation, i.e. suited for stiff problems. + # Otherwise, this is not clear. @.. broadcast=false σ=u - σ # Now tmp stores error residuals @@ -694,9 +692,8 @@ function perform_step!(integrator, cache::MPRK22ConservativeCache, repeat_step = integrator.stats.nsolve += 1 # Now σ stores the error estimate - # If a21 = 1.0, then σ is the MPE approximation, i.e. a first order approximation - # of the solution, and can be used for error estimation. Moreover, MPE is suited for stiff problems. - # TODO: Find first order approximation if a21≠ 1.0. + # If a21 = 1, then σ is the MPE approximation, i.e. suited for stiff problems. + # Otherwise, this is not clear. @.. broadcast=false σ=u - σ # Now tmp stores error residuals From f496304c550b6d69e45429b2b3da12fb3ca55ca5 Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Tue, 2 Jul 2024 15:40:52 +0200 Subject: [PATCH 43/48] corrected typos in docstrings --- src/mprk.jl | 2 +- src/sspmprk.jl | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/mprk.jl b/src/mprk.jl index 98e664b5..d0587a95 100644 --- a/src/mprk.jl +++ b/src/mprk.jl @@ -824,7 +824,7 @@ These modified Patankar-Runge-Kutta methods require the special structure of a 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 +You can also choose the parameter `small_constant` which is added to all Patankar-weight denominators to avoid divisions by zero. ## References diff --git a/src/sspmprk.jl b/src/sspmprk.jl index f698eb9a..9841befe 100644 --- a/src/sspmprk.jl +++ b/src/sspmprk.jl @@ -21,7 +21,7 @@ This modified Patankar-Runge-Kutta method requires the special structure of a 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 +You can also choose the parameter `small_constant` which is added to all Patankar-weight denominators to avoid divisions by zero. ## References @@ -407,7 +407,7 @@ This modified Patankar-Runge-Kutta method requires the special structure of a 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 +You can also choose the parameter `small_constant` which is added to all Patankar-weight denominators to avoid divisions by zero. ## References From 86a791b065f52e107c1ca033ba53b8292184e4dc Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Tue, 2 Jul 2024 16:23:56 +0200 Subject: [PATCH 44/48] 'Different matrix types' tests no use adaptive time stepping (if supported by the scheme) --- test/runtests.jl | 26 ++++++++++++-------------- 1 file changed, 12 insertions(+), 14 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index b67366b8..72de462c 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -552,14 +552,12 @@ const prob_pds_linmod_nonconservative_inplace = PDSProblem(linmodP!, linmodD!, [ prob_sparse_op = ConservativePDSProblem(prod, u0, tspan; p_prototype = P_sparse) - sol_tridiagonal_ip = solve(prob_tridiagonal_ip, alg; dt, - adaptive = false) - sol_tridiagonal_op = solve(prob_tridiagonal_op, alg; dt, - adaptive = false) - sol_dense_ip = solve(prob_dense_ip, alg; dt, adaptive = false) - sol_dense_op = solve(prob_dense_op, alg; dt, adaptive = false) - sol_sparse_ip = solve(prob_sparse_ip, alg; dt, adaptive = false) - sol_sparse_op = solve(prob_sparse_op, alg; dt, adaptive = false) + sol_tridiagonal_ip = solve(prob_tridiagonal_ip, alg; dt) + sol_tridiagonal_op = solve(prob_tridiagonal_op, alg; dt) + sol_dense_ip = solve(prob_dense_ip, alg; dt) + sol_dense_op = solve(prob_dense_op, alg; dt) + sol_sparse_ip = solve(prob_sparse_ip, alg; dt) + sol_sparse_op = solve(prob_sparse_op, alg; dt) @test sol_tridiagonal_ip.t ≈ sol_tridiagonal_op.t @test sol_dense_ip.t ≈ sol_dense_op.t @@ -676,17 +674,17 @@ const prob_pds_linmod_nonconservative_inplace = PDSProblem(linmodP!, linmodD!, [ p_prototype = P_sparse) sol_tridiagonal_ip = solve(prob_tridiagonal_ip, alg; - dt, adaptive = false) + dt) sol_tridiagonal_op = solve(prob_tridiagonal_op, alg; - dt, adaptive = false) + dt) sol_dense_ip = solve(prob_dense_ip, alg; - dt, adaptive = false) + dt) sol_dense_op = solve(prob_dense_op, alg; - dt, adaptive = false) + dt) sol_sparse_ip = solve(prob_sparse_ip, alg; - dt, adaptive = false) + dt) sol_sparse_op = solve(prob_sparse_op, alg; - dt, adaptive = false) + dt) @test sol_tridiagonal_ip.t ≈ sol_tridiagonal_op.t @test sol_dense_ip.t ≈ sol_dense_op.t From 10352b4e85f0599e1a25e09135b56195c0ecd14b Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Tue, 2 Jul 2024 22:14:01 +0200 Subject: [PATCH 45/48] We now test 'Different matrix types' with constant step size and adaptive time stepping. For adaptive time stepping isapprox is used with a larger relative tolerance. --- test/runtests.jl | 237 ++++++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 225 insertions(+), 12 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index 72de462c..d3f503b0 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -552,12 +552,14 @@ const prob_pds_linmod_nonconservative_inplace = PDSProblem(linmodP!, linmodD!, [ prob_sparse_op = ConservativePDSProblem(prod, u0, tspan; p_prototype = P_sparse) - sol_tridiagonal_ip = solve(prob_tridiagonal_ip, alg; dt) - sol_tridiagonal_op = solve(prob_tridiagonal_op, alg; dt) - sol_dense_ip = solve(prob_dense_ip, alg; dt) - sol_dense_op = solve(prob_dense_op, alg; dt) - sol_sparse_ip = solve(prob_sparse_ip, alg; dt) - sol_sparse_op = solve(prob_sparse_op, alg; dt) + sol_tridiagonal_ip = solve(prob_tridiagonal_ip, alg; dt, + adaptive = false) + sol_tridiagonal_op = solve(prob_tridiagonal_op, alg; dt, + adaptive = false) + sol_dense_ip = solve(prob_dense_ip, alg; dt, adaptive = false) + sol_dense_op = solve(prob_dense_op, alg; dt, adaptive = false) + sol_sparse_ip = solve(prob_sparse_ip, alg; dt, adaptive = false) + sol_sparse_op = solve(prob_sparse_op, alg; dt, adaptive = false) @test sol_tridiagonal_ip.t ≈ sol_tridiagonal_op.t @test sol_dense_ip.t ≈ sol_dense_op.t @@ -574,6 +576,89 @@ const prob_pds_linmod_nonconservative_inplace = PDSProblem(linmodP!, linmodD!, [ end 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 (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), SSPMPRK43()) + 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_tridiagonal_ip = ConservativePDSProblem(prod!, u0, tspan; + p_prototype = P_tridiagonal) + prob_tridiagonal_op = ConservativePDSProblem(prod, u0, tspan; + p_prototype = P_tridiagonal) + prob_dense_ip = ConservativePDSProblem(prod!, u0, tspan; + p_prototype = P_dense) + prob_dense_op = ConservativePDSProblem(prod, u0, tspan; + p_prototype = P_dense) + prob_sparse_ip = ConservativePDSProblem(prod!, u0, tspan; + p_prototype = P_sparse) + prob_sparse_op = ConservativePDSProblem(prod, u0, tspan; + p_prototype = P_sparse) + + sol_tridiagonal_ip = solve(prob_tridiagonal_ip, alg; dt) + sol_tridiagonal_op = solve(prob_tridiagonal_op, alg; dt) + sol_dense_ip = solve(prob_dense_ip, alg; dt) + sol_dense_op = solve(prob_dense_op, alg; dt) + sol_sparse_ip = solve(prob_sparse_ip, alg; dt) + sol_sparse_op = solve(prob_sparse_op, alg; dt) + + @test isapprox(sol_tridiagonal_ip.t, sol_tridiagonal_op.t; rtol) + @test isapprox(sol_dense_ip.t, sol_dense_op.t; rtol) + @test isapprox(sol_sparse_ip.t, sol_sparse_op.t; rtol) + @test isapprox(sol_tridiagonal_ip.t, sol_dense_ip.t; rtol) + @test isapprox(sol_tridiagonal_ip.t, sol_sparse_ip.t; rtol) + + @test isapprox(sol_tridiagonal_ip.u, sol_tridiagonal_op.u; rtol) + @test isapprox(sol_dense_ip.u, sol_dense_op.u; rtol) + @test isapprox(sol_sparse_ip.u, sol_sparse_op.u; rtol) + @test isapprox(sol_tridiagonal_ip.u, sol_dense_ip.u; rtol) + @test isapprox(sol_tridiagonal_ip.u, sol_sparse_ip.u; rtol) + end + end + end + # Here we check that in-place and out-of-place implementations # deliver the same results @testset "Different matrix types (nonconservative)" begin @@ -674,17 +759,17 @@ const prob_pds_linmod_nonconservative_inplace = PDSProblem(linmodP!, linmodD!, [ p_prototype = P_sparse) sol_tridiagonal_ip = solve(prob_tridiagonal_ip, alg; - dt) + dt, adaptive = false) sol_tridiagonal_op = solve(prob_tridiagonal_op, alg; - dt) + dt, adaptive = false) sol_dense_ip = solve(prob_dense_ip, alg; - dt) + dt, adaptive = false) sol_dense_op = solve(prob_dense_op, alg; - dt) + dt, adaptive = false) sol_sparse_ip = solve(prob_sparse_ip, alg; - dt) + dt, adaptive = false) sol_sparse_op = solve(prob_sparse_op, alg; - dt) + dt, adaptive = false) @test sol_tridiagonal_ip.t ≈ sol_tridiagonal_op.t @test sol_dense_ip.t ≈ sol_dense_op.t @@ -701,6 +786,134 @@ const prob_pds_linmod_nonconservative_inplace = PDSProblem(linmodP!, linmodD!, [ end end + @testset "Different matrix types (nonconservative, 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 + for i in 1:length(u) + P[i, i] = i * u[i] + end + return nothing + end + dest_1! = (D, u, p, t) -> begin + fill!(D, zero(eltype(D))) + for i in 1:length(u) + D[i] = (i + 1) * 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 + for i in 1:length(u) + P[i, i] = (i - 1) * u[i] + end + return nothing + end + dest_2! = (D, u, p, t) -> begin + fill!(D, zero(eltype(D))) + for i in 1:length(u) + D[i] = i * u[i] + 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 + for i in 1:length(u) + P[i, i] = (i + 1) * u[i] + end + return nothing + end + dest_3! = (D, u, p, t) -> begin + fill!(D, zero(eltype(D))) + for i in 1:length(u) + D[i] = (i - 1) * u[i] + 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] + D = u0 + tspan = (0.0, 1.0) + dt = 0.25 + + rtol = sqrt(eps(Float32)) + @testset "$alg" 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 (prod!, dest!) in zip((prod_1!, prod_2!, prod_3!), + (dest_1!, dest_2!, dest_3!)) + prod! = prod_3! + dest! = dest_3! + prod = (u, p, t) -> begin + P = similar(u, (length(u), length(u))) + prod!(P, u, p, t) + return P + end + dest = (u, p, t) -> begin + D = similar(u) + dest!(D, u, p, t) + return D + end + prob_tridiagonal_ip = PDSProblem(prod!, dest!, u0, tspan; + p_prototype = P_tridiagonal) + prob_tridiagonal_op = PDSProblem(prod, dest, u0, tspan; + p_prototype = P_tridiagonal) + prob_dense_ip = PDSProblem(prod!, dest!, u0, tspan; + p_prototype = P_dense) + prob_dense_op = PDSProblem(prod, dest, u0, tspan; + p_prototype = P_dense) + prob_sparse_ip = PDSProblem(prod!, dest!, u0, tspan; + p_prototype = P_sparse) + prob_sparse_op = PDSProblem(prod, dest, u0, tspan; + p_prototype = P_sparse) + + sol_tridiagonal_ip = solve(prob_tridiagonal_ip, alg; + dt) + sol_tridiagonal_op = solve(prob_tridiagonal_op, alg; + dt) + sol_dense_ip = solve(prob_dense_ip, alg; + dt) + sol_dense_op = solve(prob_dense_op, alg; + dt) + sol_sparse_ip = solve(prob_sparse_ip, alg; + dt) + sol_sparse_op = solve(prob_sparse_op, alg; + dt) + + @test isapprox(sol_tridiagonal_ip.t, sol_tridiagonal_op.t; rtol) + @test isapprox(sol_dense_ip.t, sol_dense_op.t; rtol) + @test isapprox(sol_sparse_ip.t, sol_sparse_op.t; rtol) + @test isapprox(sol_tridiagonal_ip.t, sol_dense_ip.t; rtol) + @test isapprox(sol_tridiagonal_ip.t, sol_sparse_ip.t; rtol) + + @test isapprox(sol_tridiagonal_ip.u, sol_tridiagonal_op.u; rtol) + @test isapprox(sol_dense_ip.u, sol_dense_op.u; rtol) + @test isapprox(sol_sparse_ip.u, sol_sparse_op.u; rtol) + @test isapprox(sol_tridiagonal_ip.u, sol_dense_ip.u; rtol) + @test isapprox(sol_tridiagonal_ip.u, sol_sparse_ip.u; rtol) + 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 019e5e05e52ff1dab09cb2f42a5969115b882d30 Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Tue, 2 Jul 2024 22:29:55 +0200 Subject: [PATCH 46/48] Additional explanation of small_constant in docstring --- src/mprk.jl | 12 ++++++++---- src/sspmprk.jl | 6 ++++-- 2 files changed, 12 insertions(+), 6 deletions(-) diff --git a/src/mprk.jl b/src/mprk.jl index d0587a95..02434d9e 100644 --- a/src/mprk.jl +++ b/src/mprk.jl @@ -134,7 +134,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. +to avoid divisions by zero. You can pass a value explicitly, otherwise `small_constant` is set to +`floatmin` of the floating point type used. ## References @@ -345,7 +346,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. +to avoid divisions by zero. You can pass a value explicitly, otherwise `small_constant` is set to +`floatmin` of the floating point type used. ## References @@ -726,7 +728,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. +to avoid divisions by zero. You can pass a value explicitly, otherwise `small_constant` is set to +`floatmin` of the floating point type used. ## References @@ -825,7 +828,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. +to avoid divisions by zero. You can pass a value explicitly, otherwise `small_constant` is set to +`floatmin` of the floating point type used. ## References diff --git a/src/sspmprk.jl b/src/sspmprk.jl index 9841befe..6380fd2b 100644 --- a/src/sspmprk.jl +++ b/src/sspmprk.jl @@ -22,7 +22,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. +to avoid divisions by zero. You can pass a value explicitly, otherwise `small_constant` is set to +`floatmin` of the floating point type used. ## References @@ -408,7 +409,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. +to avoid divisions by zero. You can pass a value explicitly, otherwise `small_constant` is set to +`floatmin` of the floating point type used. ## References From 0a88e06dd44eb94c23d65df2349bc5182ce8110c Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Wed, 3 Jul 2024 12:01:07 +0200 Subject: [PATCH 47/48] Use 1st order approximation for error estimates in SSPMPRK22 --- src/sspmprk.jl | 47 +++++++++++++++++++++++++++++++---------------- 1 file changed, 31 insertions(+), 16 deletions(-) diff --git a/src/sspmprk.jl b/src/sspmprk.jl index 6380fd2b..9f69fa2a 100644 --- a/src/sspmprk.jl +++ b/src/sspmprk.jl @@ -77,11 +77,13 @@ function get_constant_parameters(alg::SSPMPRK22) s = (b20 + b21 + a21 * b10^2) / (b10 * (b20 + b21)) + τ = 1 + (a21 * b10^2) / (b20 + b21) + # This should never happen if any(<(0), (a21, a10, a20, b10, b20, b21)) throw(ArgumentError("SSPMPRK22 requires nonnegative SSP coefficients.")) end - return a21, a10, a20, b10, b20, b21, s + return a21, a10, a20, b10, b20, b21, s, τ end struct SSPMPRK22ConstantCache{T} <: OrdinaryDiffEqConstantCache @@ -92,6 +94,7 @@ struct SSPMPRK22ConstantCache{T} <: OrdinaryDiffEqConstantCache b20::T b21::T s::T + τ::T small_constant::T end @@ -104,8 +107,8 @@ function alg_cache(alg::SSPMPRK22, u, rate_prototype, ::Type{uEltypeNoUnits}, throw(ArgumentError("SSPMPRK22 can only be applied to production-destruction systems")) end - a21, a10, a20, b10, b20, b21, s = get_constant_parameters(alg) - SSPMPRK22ConstantCache(a21, a10, a20, b10, b20, b21, s, + a21, a10, a20, b10, b20, b21, s, τ = get_constant_parameters(alg) + SSPMPRK22ConstantCache(a21, a10, a20, b10, b20, b21, s, τ, alg.small_constant_function(uEltypeNoUnits)) end @@ -114,7 +117,7 @@ end function perform_step!(integrator, cache::SSPMPRK22ConstantCache, repeat_step = false) @unpack alg, t, dt, uprev, f, p = integrator - @unpack a21, a10, a20, b10, b20, b21, s, small_constant = cache + @unpack a21, a10, a20, b10, b20, b21, s, τ, small_constant = cache f = integrator.f @@ -175,10 +178,17 @@ function perform_step!(integrator, cache::SSPMPRK22ConstantCache, repeat_step = u = sol.u integrator.stats.nsolve += 1 + # Unless τ = 1, σ is not a first order approximation, since + # σ = uprev + τ * dt *(P^n − D^n) + O(dt^2), see (2.7) in + # https://doi.org/10.1007/s10915-018-0852-1. + # But we can compute a 1st order approximation σ2, as follows. + # σ2 may become negative, but still can be used for error estimation. + σ2 = (σ - uprev) / τ + uprev + # If a21 = 0 or b10 = 0, then σ is a first order approximation of the solution and # can be used for error estimation. # TODO: Find first order approximation, if a21*b10 ≠ 0. - tmp = u - σ + tmp = u - σ2 atmp = calculate_residuals(tmp, uprev, u, integrator.opts.abstol, integrator.opts.reltol, integrator.opts.internalnorm, t) integrator.EEst = integrator.opts.internalnorm(atmp, t) @@ -213,8 +223,8 @@ function alg_cache(alg::SSPMPRK22, u, rate_prototype, ::Type{uEltypeNoUnits}, ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, dt, reltol, p, calck, ::Val{true}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} - a21, a10, a20, b10, b20, b21, s = get_constant_parameters(alg) - tab = SSPMPRK22ConstantCache(a21, a10, a20, b10, b20, b21, s, + a21, a10, a20, b10, b20, b21, s, τ = get_constant_parameters(alg) + tab = SSPMPRK22ConstantCache(a21, a10, a20, b10, b20, b21, s, τ, alg.small_constant_function(uEltypeNoUnits)) tmp = zero(u) P = p_prototype(u, f) @@ -253,7 +263,7 @@ end function perform_step!(integrator, cache::SSPMPRK22Cache, repeat_step = false) @unpack t, dt, uprev, u, f, p = integrator @unpack tmp, P, P2, D, D2, σ, linsolve = cache - @unpack a21, a10, a20, b10, b20, b21, s, small_constant = cache.tab + @unpack a21, a10, a20, b10, b20, b21, s, τ, small_constant = cache.tab # We use P2 to store the last evaluation of the PDS # as well as to store the system matrix of the linear system @@ -311,12 +321,14 @@ function perform_step!(integrator, cache::SSPMPRK22Cache, repeat_step = false) u .= linres integrator.stats.nsolve += 1 - # If a21 = 0 or b10 = 0, then σ is a first order approximation of the solution and - # can be used for error estimation. - # TODO: Find first order approximation, if a21*b10 ≠ 0. + # Unless τ = 1, σ is not a first order approximation, since + # σ = uprev + τ * dt *(P^n − D^n) + O(dt^2), see (2.7) in + # https://doi.org/10.1007/s10915-018-0852-1. + # But we can compute a 1st order approximation as σ2 = (σ - uprev) / τ + uprev. + # σ2 may become negative, but still can be used for error estimation. # Now σ stores the error estimate - @.. broadcast=false σ=u - σ + @.. broadcast=false σ=u - (σ - uprev) / τ - uprev # Now tmp stores error residuals calculate_residuals!(tmp, σ, uprev, u, integrator.opts.abstol, @@ -374,11 +386,14 @@ function perform_step!(integrator, cache::SSPMPRK22ConservativeCache, repeat_ste u .= linres integrator.stats.nsolve += 1 - # If a21 = 0 or b10 = 0, then σ is a first order approximation of the solution and - # can be used for error estimation. - # TODO: Find first order approximation, if a21*b10 ≠ 0. + # Unless τ = 1, σ is not a first order approximation, since + # σ = uprev + τ * dt *(P^n − D^n) + O(dt^2), see (2.7) in + # https://doi.org/10.1007/s10915-018-0852-1. + # But we can compute a 1st order approximation as σ2 = (σ - uprev) / τ + uprev. + # σ2 may become negative, but still can be used for error estimation. + # Now σ stores the error estimate - @.. broadcast=false σ=u - σ + @.. broadcast=false σ=u - (σ - uprev) / τ - uprev # Now tmp stores error residuals calculate_residuals!(tmp, σ, uprev, u, integrator.opts.abstol, From bd6cdca688a3d18a1cb7da3d595e616d490e9112 Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Wed, 3 Jul 2024 12:54:41 +0200 Subject: [PATCH 48/48] =?UTF-8?q?Fixed=20bug:=20undefined=20=CF=84?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/sspmprk.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/sspmprk.jl b/src/sspmprk.jl index 9f69fa2a..0a8ed96f 100644 --- a/src/sspmprk.jl +++ b/src/sspmprk.jl @@ -340,7 +340,7 @@ end function perform_step!(integrator, cache::SSPMPRK22ConservativeCache, repeat_step = false) @unpack t, dt, uprev, u, f, p = integrator @unpack tmp, P, P2, σ, linsolve = cache - @unpack a21, a10, a20, b10, b20, b21, s, small_constant = cache.tab + @unpack a21, a10, a20, b10, b20, b21, s, τ, small_constant = cache.tab # We use P2 to store the last evaluation of the PDS # as well as to store the system matrix of the linear system