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" diff --git a/src/PositiveIntegrators.jl b/src/PositiveIntegrators.jl index 8929959c..f5fd7c26 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, SSPMPRK43 export prob_pds_linmod, prob_pds_linmod_inplace, prob_pds_nonlinmod, prob_pds_robertson, prob_pds_brusselator, prob_pds_sir, @@ -61,6 +62,11 @@ 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") + +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..e86cdf84 --- /dev/null +++ b/src/interpolation.jl @@ -0,0 +1,86 @@ +# 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 +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 <: MPRKCaches} + "1st order linear" +end + +function _ode_interpolant(Θ, dt, u0, u1, k, + cache::MPRKCaches, + 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::MPRKCaches, + 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::MPRKCaches, + 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::MPRKCaches, + 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 f54cc418..e1524bc1 100644 --- a/src/mprk.jl +++ b/src/mprk.jl @@ -114,52 +114,9 @@ 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 = ...]) + 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 @@ -176,6 +133,9 @@ 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. You can pass a value explicitly, otherwise `small_constant` is set to +`floatmin` of the floating point type used. ## References @@ -185,12 +145,20 @@ 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_function::T end -function MPE(; linsolve = LUFactorization()) - MPE(linsolve) +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 @@ -208,7 +176,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(alg.small_constant_function(uEltypeNoUnits)) end function initialize!(integrator, cache::MPEConstantCache) @@ -270,7 +238,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_function(uEltypeNoUnits)) if f isa ConservativePDSFunction # We use P to store the evaluation of the PDS @@ -356,48 +324,9 @@ 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 = ...]) + MPRK22(α; [linsolve = ..., small_constant = ...]) 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 @@ -416,6 +345,9 @@ 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. You can pass a value explicitly, otherwise `small_constant` is set to +`floatmin` of the floating point type used. ## References @@ -438,13 +370,24 @@ as keyword argument `linsolve`. 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_function::T2 end -function MPRK22(alpha; linsolve = LUFactorization()) - MPRK22{typeof(alpha), typeof(linsolve)}(alpha, linsolve) +function MPRK22(alpha; 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 + MPRK22{typeof(alpha), typeof(linsolve), typeof(small_constant_function)}(alpha, + linsolve, + small_constant_function) end alg_order(::MPRK22) = 2 @@ -483,7 +426,7 @@ function alg_cache(alg::MPRK22, u, rate_prototype, ::Type{uEltypeNoUnits}, end a21, b1, b2 = get_constant_parameters(alg) - MPRK22ConstantCache(a21, b1, b2, floatmin(uEltypeNoUnits)) + MPRK22ConstantCache(a21, b1, b2, alg.small_constant_function(uEltypeNoUnits)) end function initialize!(integrator, cache::MPRK22ConstantCache) @@ -553,9 +496,8 @@ 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, 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) @@ -592,7 +534,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) - tab = MPRK22ConstantCache(a21, b1, b2, floatmin(uEltypeNoUnits)) + 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 @@ -693,6 +635,8 @@ function perform_step!(integrator, cache::MPRK22Cache, repeat_step = false) integrator.stats.nsolve += 1 # Now σ stores the error estimate + # 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 @@ -750,6 +694,8 @@ function perform_step!(integrator, cache::MPRK22ConservativeCache, repeat_step = integrator.stats.nsolve += 1 # Now σ stores the error estimate + # 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 @@ -761,7 +707,7 @@ end ### MPRK43 ##################################################################################### """ - MPRK43I(α, β; [linsolve = ...]) + MPRK43I(α, β; [linsolve = ..., small_constant = ...]) A family of third-order modified Patankar-Runge-Kutta schemes for production-destruction systems, which is based on the two-parameter family of third order explicit Runge--Kutta schemes. @@ -781,6 +727,9 @@ 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. You can pass a value explicitly, otherwise `small_constant` is set to +`floatmin` of the floating point type used. ## References @@ -790,14 +739,25 @@ as keyword argument `linsolve`. 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_function::T2 end -function MPRK43I(alpha, beta; linsolve = LUFactorization()) - MPRK43I{typeof(alpha), typeof(linsolve)}(alpha, beta, linsolve) +function MPRK43I(alpha, beta; 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 + MPRK43I{typeof(alpha), typeof(linsolve), typeof(small_constant_function)}(alpha, beta, + linsolve, + small_constant_function) end alg_order(::MPRK43I) = 3 @@ -846,7 +806,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 @@ -867,6 +827,9 @@ 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. You can pass a value explicitly, otherwise `small_constant` is set to +`floatmin` of the floating point type used. ## References @@ -876,13 +839,23 @@ as keyword argument `linsolve`. 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_function::T2 end -function MPRK43II(gamma; linsolve = LUFactorization()) - MPRK43II{typeof(gamma), typeof(linsolve)}(gamma, linsolve) +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), typeof(small_constant_function)}(gamma, + linsolve, + small_constant_function) end alg_order(::MPRK43II) = 3 @@ -941,7 +914,7 @@ function alg_cache(alg::Union{MPRK43I, MPRK43II}, u, rate_prototype, ::Type{uElt end a21, a31, a32, b1, b2, b3, c2, c3, beta1, beta2, q1, q2 = get_constant_parameters(alg) MPRK43ConstantCache(a21, a31, a32, b1, b2, b3, c2, c3, - beta1, beta2, q1, q2, floatmin(uEltypeNoUnits)) + beta1, beta2, q1, q2, alg.small_constant_function(uEltypeNoUnits)) end function initialize!(integrator, cache::MPRK43ConstantCache) @@ -1099,7 +1072,8 @@ function alg_cache(alg::Union{MPRK43I, MPRK43II}, u, rate_prototype, ::Type{uElt ::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)) + beta1, beta2, q1, q2, + alg.small_constant_function(uEltypeNoUnits)) tmp = zero(u) tmp2 = zero(u) P = p_prototype(u, f) @@ -1338,48 +1312,3 @@ function perform_step!(integrator, cache::MPRK43ConservativeCache, repeat_step = end ######################################################################################## - -# interpolation specializations -function interp_summary(::Type{cacheType}, - dense::Bool) where { - cacheType <: - Union{MPRK22ConstantCache, MPRK22Cache, - MPRK43ConstantCache, MPRK43Cache}} - "1st order linear" -end - -function _ode_interpolant(Θ, dt, u0, u1, k, - cache::Union{MPRK22ConstantCache, MPRK22Cache, - MPRK43ConstantCache, MPRK43Cache}, - 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{MPRK22ConstantCache, MPRK22Cache, - MPRK43ConstantCache, MPRK43Cache}, - 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{MPRK22ConstantCache, MPRK22Cache, - MPRK43ConstantCache, MPRK43Cache}, - 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{MPRK22ConstantCache, MPRK22Cache, - MPRK43ConstantCache, MPRK43Cache}, - 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/sspmprk.jl b/src/sspmprk.jl new file mode 100644 index 00000000..0a8ed96f --- /dev/null +++ b/src/sspmprk.jl @@ -0,0 +1,941 @@ +### SSPMPRK ##################################################################################### +""" + 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 +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`. +You can also choose the parameter `small_constant` which is added to all Patankar-weight denominators +to avoid divisions by zero. You can pass a value explicitly, otherwise `small_constant` is set to +`floatmin` of the floating point type used. + +## 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, T2} <: OrdinaryDiffEqAdaptiveAlgorithm + alpha::T + beta::T + linsolve::F + small_constant_function::T2 +end + +function SSPMPRK22(alpha, beta; 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 + SSPMPRK22{typeof(alpha), typeof(linsolve), typeof(small_constant_function)}(alpha, beta, + linsolve, + small_constant_function) +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)) + + τ = 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, τ +end + +struct SSPMPRK22ConstantCache{T} <: OrdinaryDiffEqConstantCache + a21::T + a10::T + a20::T + b10::T + b20::T + b21::T + s::T + τ::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, τ, + alg.small_constant_function(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 + + # 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 - σ2 + 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, τ, + alg.small_constant_function(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("SSPMPRK22 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 + + # 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 - (σ - uprev) / τ - uprev + + # 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 + + # 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 - (σ - uprev) / τ - uprev + + # 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 + +""" + 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 +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`. +You can also choose the parameter `small_constant` which is added to all Patankar-weight denominators +to avoid divisions by zero. You can pass a value explicitly, otherwise `small_constant` is set to +`floatmin` of the floating point type used. + +## 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, T} <: OrdinaryDiffEqAlgorithm + linsolve::F + small_constant_function::T +end + +function SSPMPRK43(; linsolve = LUFactorization(), small_constant = 1e-50) + 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 +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 + η1 = 3.777285888379173E-02 + η2 = 1.0 / 3.0 + η3 = 1.868649805549811E-01 + η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 + α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, η5, η6, 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 + η5::T + η6::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, η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, alg.small_constant_function(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, η5, η6, 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 + + # 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) + 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 ./ ρ + # 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 + β32 * 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 SSPMPRK43Cache{uType, PType, tabType, F} <: OrdinaryDiffEqMutableCache + tmp::uType + tmp2::uType + P::PType + P2::PType + P3::PType + D::uType + D2::uType + D3::uType + σ::uType + ρ::uType + tab::tabType + linsolve::F +end + +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::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, η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, + alg.small_constant_function(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) + ρ = zero(u) + + if f isa ConservativePDSFunction + linprob = LinearProblem(P3, _vec(tmp)) + linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true, + assumptions = LinearSolve.OperatorAssumptions(true)) + SSPMPRK43ConservativeCache(tmp, tmp2, P, P2, P3, σ, ρ, tab, linsolve) + elseif f isa PDSFunction + D = zero(u) + D2 = zero(u) + D3 = zero(u) + + linprob = LinearProblem(P3, _vec(tmp)) + linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true, + assumptions = LinearSolve.OperatorAssumptions(true)) + + SSPMPRK43Cache(tmp, tmp2, P, P2, P3, D, D2, D3, σ, ρ, tab, linsolve) + else + throw(ArgumentError("SSPMPRK43 can only be applied to production-destruction systems")) + end +end + +function initialize!(integrator, cache::Union{SSPMPRK43ConservativeCache, SSPMPRK43Cache}) +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, η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 + + f.p(P, uprev, p, t) # evaluate production terms + f.d(D, uprev, p, t) # evaluate nonconservative destruction terms + @.. 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 + @.. broadcast=false tmp=α10 * 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 + tmp2 .= u + integrator.stats.nsolve += 1 + + @.. broadcast=false ρ=n1 * u + n2 * u^2 / σ + @.. broadcast=false ρ=ρ + small_constant + + 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 + @.. 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) + + # Same as linres = P3 \ tmp + linsolve.A = P3 + linres = solve!(linsolve) + + u .= linres + integrator.stats.nsolve += 1 + + @.. broadcast=false σ=σ^(1 - s) * tmp2^s + @.. broadcast=false σ=σ + small_constant + + @.. broadcast=false P3=η3 * P + η4 * P2 + @.. broadcast=false D3=η3 * D + η4 * D2 + + # tmp holds the right hand side of the linear system + @.. broadcast=false tmp=η1 * uprev + η2 * tmp2 + @inbounds for i in eachindex(tmp) + # 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) + + # Same as linres = P3 \ tmp + linsolve.A = P3 + linres = solve!(linsolve) + 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=β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 + @.. broadcast=false tmp=α30 * uprev + α31 * tmp2 + α32 * u + @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 + + #TODO: Figure out if a second order approximation of the solution + # is hidden somewhere. + #= + # 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::SSPMPRK43ConservativeCache, repeat_step = false) + @unpack t, dt, uprev, u, f, p = integrator + @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=β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 + linsolve.A = P3 + linres = solve!(linsolve) + + u .= linres + tmp2 .= u + integrator.stats.nsolve += 1 + + @.. broadcast=false ρ=n1 * u + n2 * u^2 / σ + @.. broadcast=false ρ=ρ + small_constant + + f.p(P2, u, p, t + β10 * dt) # evaluate production terms + @.. broadcast=false P3=β20 * P + β21 * P2 + integrator.stats.nf += 1 + + @.. broadcast=false tmp=α20 * uprev + α21 * tmp2 + build_mprk_matrix!(P3, P3, ρ, dt) + + # Same as linres = P3 \ tmp + linsolve.A = P3 + linres = solve!(linsolve) + + u .= linres + integrator.stats.nsolve += 1 + + @.. broadcast=false σ=σ^(1 - s) * tmp2^s + @.. broadcast=false σ=σ + small_constant + + @.. broadcast=false P3=η3 * P + η4 * P2 + @.. broadcast=false tmp=η1 * uprev + η2 * tmp2 + + build_mprk_matrix!(P3, P3, σ, dt) + + # Same as linres = P3 \ tmp + linsolve.A = P3 + linres = solve!(linsolve) + 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=β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 + linsolve.A = P3 + linres = solve!(linsolve) + + 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 - σ + + # 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 diff --git a/test/runtests.jl b/test/runtests.jl index bd9bd703..d3f503b0 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -13,65 +13,60 @@ using LinearSolve: RFLUFactorization, LUFactorization, KrylovJL_GMRES using Aqua: Aqua """ - experimental_order_of_convergence(prob, alg, dts, test_times; + experimental_orders_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 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; - 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) - 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) + return experimental_orders_of_convergence(errors, dts) end """ - experimental_order_of_convergence(errors, dts) + experimental_orders_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) +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) @@ -80,7 +75,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, @@ -324,6 +343,27 @@ 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)) + @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 "SSPMPRK43 can only be applied to production-destruction systems" solve(prob_oop, + SSPMPRK43(), + dt = 0.1) + @test_throws "SSPMPRK43 can only be applied to production-destruction systems" solve(prob_ip, + SSPMPRK43(), + dt = 0.1) end # Here we check that MPE equals implicit Euler (IE) for a linear PDS @@ -380,6 +420,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 +468,9 @@ 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...), + (; kwargs...) -> SSPMPRK43(; kwargs...)) for alg in algs # Check different linear solvers @@ -478,7 +531,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)) + MPRK43II(2.0 / 3.0), MPRK43II(0.5), + SSPMPRK22(0.5, 1.0), SSPMPRK43()) for prod! in (prod_1!, prod_2!, prod_3!) prod = (u, p, t) -> begin P = similar(u, (length(u), length(u))) @@ -522,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 @@ -594,7 +731,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)) + 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 = (u, p, t) -> begin @@ -648,28 +786,166 @@ 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 - algs = (MPE(), MPRK22(0.5), MPRK22(1.0), MPRK22(2.0)) - dts = 0.5 .^ (6:11) + algs = (MPE(), MPRK22(0.5), MPRK22(1.0), MPRK22(2.0), SSPMPRK22(0.5, 1.0)) + 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_orders_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_orders_of_convergence(prob, alg, + dts; + test_time) + @test check_order(orders, PositiveIntegrators.alg_order(alg), + atol = 0.2) + 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 end end end @@ -677,24 +953,33 @@ 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)) - dts = 0.5 .^ (6:11) + algs = (MPE(), MPRK22(0.5), MPRK22(1.0), MPRK22(2.0), SSPMPRK22(0.5, 1.0)) + 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_orders_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_orders_of_convergence(prob, alg, + dts; + test_time) + @test check_order(orders, PositiveIntegrators.alg_order(alg), + atol = 0.2) + 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 end end end @@ -702,37 +987,35 @@ 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) - for alg in [ - MPRK43I(1.0, 0.5), - MPRK43I(0.5, 0.75), - MPRK43II(0.5), - MPRK43II(2.0 / 3.0), - ], prob in problems - eoc = experimental_order_of_convergence(prob, alg, dts) - @test isapprox(eoc, PositiveIntegrators.alg_order(alg); atol = 0.2) + 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(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)) + 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_orders_of_convergence(prob, alg, dts) + @test check_order(orders, PositiveIntegrators.alg_order(alg), atol = 0.2) end end @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) @@ -748,7 +1031,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)) + 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_nonconservative, prob_pds_linmod_nonconservative_inplace) @@ -795,7 +1079,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), SSPMPRK43()) for alg in algs sol = solve(prob_ip, alg; dt = dt, adaptive = false) @@ -843,7 +1127,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), SSPMPRK43()) dt = 1e-3 for alg in algs