diff --git a/NEWS.md b/NEWS.md index 00c80934..d0b4dc63 100644 --- a/NEWS.md +++ b/NEWS.md @@ -12,6 +12,9 @@ for human readability. - The minimum required Julia version was updated to v1.10 +### Added + +- MPDeC methods of Öffner and Torlo ## Changes when updating to v0.2 from v0.1.x diff --git a/docs/src/api_reference.md b/docs/src/api_reference.md index ac1a6f9f..c9f26925 100644 --- a/docs/src/api_reference.md +++ b/docs/src/api_reference.md @@ -35,6 +35,7 @@ SSPMPRK22 MPRK43I MPRK43II SSPMPRK43 +MPDeC ``` ## Auxiliary functions diff --git a/docs/src/stratospheric_reaction_benchmark.md b/docs/src/stratospheric_reaction_benchmark.md index cb6c3a60..4b746a91 100644 --- a/docs/src/stratospheric_reaction_benchmark.md +++ b/docs/src/stratospheric_reaction_benchmark.md @@ -23,7 +23,7 @@ function stratreac_plot(sols, labels = fill("", length(sols)), sol_ref = nothing labels = [labels] end - tspan = prob.tspan + tspan = prob_pds_stratreac.tspan layout = (3, 2) linewidth = 2 xticks = (range(first(tspan), last(tspan), 4), range(12.0, 84.0, 4)) diff --git a/src/PositiveIntegrators.jl b/src/PositiveIntegrators.jl index 51306dec..737fac19 100644 --- a/src/PositiveIntegrators.jl +++ b/src/PositiveIntegrators.jl @@ -6,7 +6,7 @@ using Statistics: median using SparseArrays: SparseArrays, AbstractSparseMatrix, issparse, nonzeros, nzrange, rowvals, spdiagm -using StaticArrays: SVector, SMatrix, StaticArray, @SVector, @SMatrix +using StaticArrays: SVector, SMatrix, StaticArray, @SVector, @SMatrix, MMatrix using FastBroadcast: @.. using MuladdMacro: @muladd @@ -48,6 +48,7 @@ export ConservativePDSFunction, ConservativePDSProblem export MPE, MPRK22, MPRK43I, MPRK43II export SSPMPRK22, SSPMPRK43 +export MPDeC export prob_pds_linmod, prob_pds_linmod_inplace, prob_pds_nonlinmod, prob_pds_robertson, prob_pds_brusselator, prob_pds_sir, @@ -69,6 +70,9 @@ include("mprk.jl") # modified Patankar-Runge-Kutta based on the SSP formulation of RK methods (SSPMPRK) include("sspmprk.jl") +# MPDeC methods +include("mpdec.jl") + # interpolation for dense output include("interpolation.jl") diff --git a/src/interpolation.jl b/src/interpolation.jl index e86cdf84..d257009d 100644 --- a/src/interpolation.jl +++ b/src/interpolation.jl @@ -46,7 +46,8 @@ const MPRKCaches = Union{MPEConstantCache, MPECache, MPEConservativeCache, MPRK22ConstantCache, MPRK22Cache, MPRK22ConservativeCache, MPRK43ConstantCache, MPRK43Cache, MPRK43ConservativeCache, SSPMPRK22ConstantCache, SSPMPRK22Cache, SSPMPRK22ConservativeCache, - SSPMPRK43ConstantCache, SSPMPRK43Cache, SSPMPRK43ConservativeCache} + SSPMPRK43ConstantCache, SSPMPRK43Cache, SSPMPRK43ConservativeCache, + MPDeCConstantCache, MPDeCCache, MPDeCConservativeCache} function interp_summary(::Type{cacheType}, dense::Bool) where {cacheType <: MPRKCaches} diff --git a/src/mpdec.jl b/src/mpdec.jl new file mode 100644 index 00000000..c31d4ee0 --- /dev/null +++ b/src/mpdec.jl @@ -0,0 +1,766 @@ +""" + MPDeC(K; [nodes = :gausslobatto, linsolve = ..., small_constant = ...]) + +A family of arbitrary order modified Patankar-Runge-Kutta algorithms for +production-destruction systems. Each member of this family is an adaptive, one-step method which is +Kth order accurate, unconditionally positivity-preserving, and linearly +implicit. The integer K must be chosen to satisfy 2 ≤ K ≤ 10. +Available node choices are Lagrange or Gauss-Lobatto nodes, with the latter being the default. +These methods support adaptive time stepping, using the numerical solution obtained with one correction step less as a lower-order approximation to estimate the error. +The MPDeC schemes were introduced by Torlo and Öffner (2020) for autonomous conservative production-destruction systems and +further investigated in Torlo, Öffner and Ranocha (2022). + +For nonconservative production–destruction systems we use a straight forward extension +analogous to [`MPE`](@ref). +A general discussion of DeC schemes applied to non-autonomous differential equations +and using general integration nodes is given by Ong and Spiteri (2020). + +The MPDeC methods require 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 +`1e-300` in double precision computations or `floatmin` of the floating point type used. + +## References + +- Davide Torlo and Philipp Öffner. + "Arbitrary high-order, conservative and positivity preserving Patankar-type deferred correction schemes." + Applied Numerical Mathematics 153 (2020): 15-34. + [DOI: 10.1016/j.apnum.2020.01.025](https://doi.org/10.1016/j.apnum.2020.01.025) +- Davide Torlo, Philipp Öffner, and Hendrik Ranocha. + "Issues with positivity-preserving Patankar-type schemes." + 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) + +- Benjamin W. Ong and Raymond J. Spiteri. + "Deferred Correction Methods for Ordinary Differential Equations." + Journal of Scientific Computing 83 (2020): Article 60. + [DOI: 10.1007/s10915-020-01235-8](https://doi.org/10.1007/s10915-020-01235-8) +""" +struct MPDeC{N, F, T2} <: OrdinaryDiffEqAdaptiveAlgorithm + K::Int + M::Int + nodes::N + linsolve::F + small_constant_function::T2 +end + +function small_constant_function_MPDeC(type) + if type == Float64 + # small_constant is chosen such that + # the testset "Zero initial values" passes. + small_constant = 1e-300 + else + small_constant = floatmin(type) + end + return small_constant +end + +function MPDeC(K::Integer; + nodes = :gausslobatto, + linsolve = LUFactorization(), + small_constant = small_constant_function_MPDeC) + if small_constant isa Number + small_constant_function = Returns(small_constant) + else # assume small_constant isa Function + small_constant_function = small_constant + end + + if nodes === :lagrange + M = K - 1 + elseif nodes === :gausslobatto + M = cld(K, 2) + else + throw(ArgumentError("Unknown node choice for MPDeC. Possible options are :gauslobatto or :lagrange")) + end + + MPDeC{typeof(nodes), typeof(linsolve), typeof(small_constant_function)}(K, M, + nodes, + linsolve, + small_constant_function) +end + +alg_order(alg::MPDeC) = alg.K +isfsal(::MPDeC) = false + +function get_constant_parameters(alg::MPDeC, type) + oneType = one(type) + if alg.nodes === :lagrange + nodes = collect(i * oneType / alg.M for i in 0:(alg.M)) + # M is one less than the methods order + if alg.M == 1 + theta = [0 oneType/2; 0 oneType/2] + elseif alg.M == 2 + theta = [0 5oneType/24 oneType/6; + 0 oneType/3 2oneType/3; + 0 -oneType/24 oneType/6] + elseif alg.M == 3 + theta = [0 1oneType/8 1oneType/9 1oneType/8; + 0 19oneType/72 4oneType/9 3oneType/8; + 0 -5oneType/72 1oneType/9 3oneType/8; + 0 1oneType/72 0 1oneType/8] + elseif alg.M == 4 + theta = th2 = [0 251oneType/2880 29oneType/360 27oneType/320 7oneType/90; + 0 323oneType/1440 31oneType/90 51oneType/160 16oneType/45; + 0 -11oneType/120 1oneType/15 9oneType/40 2oneType/15; + 0 53oneType/1440 1oneType/90 21oneType/160 16oneType/45; + 0 -19oneType/2880 -1oneType/360 -3oneType/320 7oneType/90] + elseif alg.M == 5 + theta = [0 19oneType/288 14oneType/225 51oneType/800 14oneType/225 19oneType/288; + 0 1427oneType/7200 43oneType/150 219oneType/800 64oneType/225 25oneType/96; + 0 -133oneType/1200 7oneType/225 57oneType/400 8oneType/75 25oneType/144; + 0 241oneType/3600 7oneType/225 57oneType/400 64oneType/225 25oneType/144; + 0 -173oneType/7200 -1oneType/75 -21oneType/800 14oneType/225 25oneType/96; + 0 3oneType/800 1oneType/450 3oneType/800 0 19oneType/288] + elseif alg.M == 6 + theta = [0 19087oneType/362880 1139oneType/22680 137oneType/2688 143oneType/2835 3715oneType/72576 41oneType/840; + 0 2713oneType/15120 47oneType/189 27oneType/112 232oneType/945 725oneType/3024 9oneType/35; + 0 -15487oneType/120960 11oneType/7560 387oneType/4480 64oneType/945 2125oneType/24192 9oneType/280; + 0 293oneType/2835 166oneType/2835 17oneType/105 752oneType/2835 125oneType/567 34oneType/105; + 0 -6737oneType/120960 -269oneType/7560 -243oneType/4480 29oneType/945 3875oneType/24192 9oneType/280; + 0 263oneType/15120 11oneType/945 9oneType/560 8oneType/945 235oneType/3024 9oneType/35; + 0 -863oneType/362880 -37oneType/22680 -29oneType/13440 -4oneType/2835 -275oneType/72576 41oneType/840] + elseif alg.M == 7 + theta = [0 751oneType/17280 41oneType/980 265oneType/6272 278oneType/6615 265oneType/6272 41oneType/980 751oneType/17280; + 0 139849oneType/846720 1466oneType/6615 1359oneType/6272 1448oneType/6615 36725oneType/169344 54oneType/245 3577oneType/17280; + 0 -4511oneType/31360 -71oneType/2940 1377oneType/31360 8oneType/245 775oneType/18816 27oneType/980 49oneType/640; + 0 123133oneType/846720 68oneType/735 5927oneType/31360 1784oneType/6615 4625oneType/18816 68oneType/245 2989oneType/17280; + 0 -88547oneType/846720 -1927oneType/26460 -3033oneType/31360 -106oneType/6615 13625oneType/169344 27oneType/980 2989oneType/17280; + 0 1537oneType/31360 26oneType/735 1377oneType/31360 8oneType/245 1895oneType/18816 54oneType/245 49oneType/640; + 0 -11351oneType/846720 -29oneType/2940 -373oneType/31360 -64oneType/6615 -275oneType/18816 41oneType/980 3577oneType/17280; + 0 275oneType/169344 8oneType/6615 9oneType/6272 8oneType/6615 275oneType/169344 0 751oneType/17280] + elseif alg.M == 8 + theta = [0 1070017oneType/29030400 32377oneType/907200 12881oneType/358400 4063oneType/113400 41705oneType/1161216 401oneType/11200 149527oneType/4147200 989oneType/28350; + 0 2233547oneType/14515200 22823oneType/113400 35451oneType/179200 2822oneType/14175 115075oneType/580608 279oneType/1400 408317oneType/2073600 2944oneType/14175; + 0 -2302297oneType/14515200 -21247oneType/453600 1719oneType/179200 61oneType/28350 3775oneType/580608 9oneType/5600 24353oneType/2073600 -464oneType/14175; + 0 2797679oneType/14515200 15011oneType/113400 39967oneType/179200 4094oneType/14175 159175oneType/580608 403oneType/1400 542969oneType/2073600 5248oneType/14175; + 0 -31457oneType/181440 -2903oneType/22680 -351oneType/2240 -227oneType/2835 -125oneType/36288 -9oneType/280 343oneType/25920 -454oneType/2835; + 0 1573169oneType/14515200 9341oneType/113400 17217oneType/179200 1154oneType/14175 85465oneType/580608 333oneType/1400 368039oneType/2073600 5248oneType/14175; + 0 -645607oneType/14515200 -15577oneType/453600 -7031oneType/179200 -989oneType/28350 -24575oneType/580608 79oneType/5600 261023oneType/2073600 -464oneType/14175; + 0 156437oneType/14515200 953oneType/113400 243oneType/25600 122oneType/14175 5725oneType/580608 9oneType/1400 111587oneType/2073600 2944oneType/14175; + 0 -33953oneType/29030400 -119oneType/129600 -369oneType/358400 -107oneType/113400 -175oneType/165888 -9oneType/11200 -8183oneType/4147200 989oneType/28350] + + elseif alg.M == 9 + theta = [0 2857oneType/89600 3956oneType/127575 399oneType/12800 7oneType/225 81385oneType/2612736 7oneType/225 399oneType/12800 3956oneType/127575 2857oneType/89600; + 0 9449717oneType/65318400 37829oneType/204120 16381oneType/89600 3346oneType/18225 478525oneType/2612736 257oneType/1400 341383oneType/1866240 23552oneType/127575 15741oneType/89600; + 0 -1408913oneType/8164800 -34369oneType/510300 -31oneType/1600 -628oneType/25515 -7225oneType/326592 -17oneType/700 -24647oneType/1166400 -3712oneType/127575 27oneType/2240; + 0 200029oneType/816480 45331oneType/255150 13273oneType/50400 40648oneType/127575 50425oneType/163296 199oneType/630 178703oneType/583200 41984oneType/127575 1209oneType/5600; + 0 -8641823oneType/32659200 -103987oneType/510300 -2123oneType/8960 -20924oneType/127575 -131575oneType/1306368 -83oneType/700 -460649oneType/4665600 -3632oneType/25515 2889oneType/44800; + 0 6755041oneType/32659200 83291oneType/510300 8201oneType/44800 21076oneType/127575 298505oneType/1306368 211oneType/700 1251607oneType/4665600 41984oneType/127575 2889oneType/44800; + 0 -462127oneType/4082400 -9239oneType/102060 -5039oneType/50400 -11852oneType/127575 -16775oneType/163296 -299oneType/6300 4459oneType/116640 -3712oneType/127575 1209oneType/5600; + 0 335983oneType/8164800 8467oneType/255150 407oneType/11200 872oneType/25515 11975oneType/326592 11oneType/350 92617oneType/1166400 23552oneType/127575 27oneType/2240; + 0 -116687oneType/13063680 -3697oneType/510300 -101oneType/12800 -953oneType/127575 -20675oneType/2612736 -1oneType/140 -90013oneType/9331200 3956oneType/127575 15741oneType/89600; + 0 8183oneType/9331200 1oneType/1400 25oneType/32256 94oneType/127575 25oneType/32256 1oneType/1400 8183oneType/9331200 0 2857oneType/89600] + + else + error("MPDeC requires 2 ≤ K ≤ 10.") + end + else # alg.nodes === :gausslobatto + if alg.M == 1 + nodes = [0, oneType] + theta = [0 oneType/2; 0 oneType/2] + elseif alg.M == 2 + nodes = [0, oneType / 2, oneType] + theta = [0 5oneType/24 oneType/6; + 0 oneType/3 2oneType/3; + 0 -oneType/24 oneType/6] + elseif alg.M == 3 + nodes = [ + 0, + (1 - sqrt(oneType / 5)) / 2, + (1 + sqrt(oneType / 5)) / 2, + oneType + ] + theta = [0 (11 + sqrt(5oneType))/120 (11 - sqrt(5oneType))/120 oneType/12; + 0 (25 - sqrt(5oneType))/120 (13 * sqrt(5oneType) + 25)/120 5oneType/12; + 0 (25 - 13 * sqrt(5oneType))/120 (sqrt(5oneType) + 25)/120 5oneType/12; + 0 (sqrt(5oneType) - 1)/120 (-sqrt(5oneType) - 1)/120 oneType/12] + elseif alg.M == 4 + nodes = [ + 0, + (1 - sqrt(3oneType / 7)) / 2, + oneType / 2, + (1 + sqrt(3oneType / 7)) / 2, + oneType + ] + + theta = [0 (3sqrt(21oneType) + 119)/1960 13oneType/320 (119 - 3sqrt(21oneType))/1960 oneType/20; + 0 49oneType / 360-sqrt(21oneType) / 280 (7sqrt(21oneType)) / 192+49oneType / 360 (23sqrt(21oneType)) / 840+49oneType / 360 49oneType/180; + 0 8oneType / 45-(32sqrt(21oneType)) / 735 8oneType/45 (32sqrt(21oneType)) / 735+8oneType / 45 16oneType/45; + 0 49oneType / 360-(23sqrt(21oneType)) / 840 49oneType / 360-(7sqrt(21oneType)) / 192 sqrt(21oneType) / 280+49oneType / 360 49oneType/180; + 0 ((3sqrt(21oneType)) / 1960-3oneType / 280) 3oneType/320 (-(3sqrt(21oneType)) / 1960-3oneType / 280) 1oneType/20] + + elseif alg.M == 5 + nodes = [ + 0, + oneType / 2 - sqrt((2sqrt(7oneType)) / 21 + oneType / 3) / 2, + oneType / 2 - sqrt(oneType / 3 - (2sqrt(7oneType)) / 21) / 2, + oneType / 2 + sqrt(oneType / 3 - (2sqrt(7oneType)) / 21) / 2, + oneType / 2 + sqrt((2sqrt(7oneType)) / 21 + oneType / 3) / 2, + oneType + ] + + theta = [0 (sqrt(7oneType) / 756 + sqrt(6sqrt(7oneType) + 21) / 540 - sqrt(42sqrt(7oneType) + 147) / 3780+19oneType / 540) (19oneType / 540 - sqrt(21 - 6sqrt(7oneType)) / 540 - sqrt(147 - 42sqrt(7oneType)) / 3780-sqrt(7oneType) / 756) (sqrt(21 - 6sqrt(7oneType)) / 540-sqrt(7oneType) / 756+sqrt(147 - 42sqrt(7oneType))/3780+19oneType/540) (sqrt(7oneType) / 756-sqrt(6sqrt(7oneType) + 21) / 540+sqrt(42sqrt(7oneType) + 147)/3780+19oneType/540) oneType/30; + 0 (7oneType / 60 - sqrt(42sqrt(7oneType) + 147) / 1260-sqrt(7oneType) / 120) (-((sqrt(3oneType) - 2sqrt(21oneType)) * (63sqrt(21oneType) + 35sqrt(2sqrt(7oneType) + 7) + 70sqrt(6sqrt(7oneType) + 21) + 20sqrt(14sqrt(7oneType) + 49) - 23sqrt(42sqrt(7oneType) + 147)))/22680) (-((sqrt(3oneType) - 2sqrt(21oneType)) * (63sqrt(21oneType) + 35sqrt(2sqrt(7oneType) + 7) - 70sqrt(6sqrt(7oneType) + 21) + 20sqrt(14sqrt(7oneType) + 49) + 23sqrt(42sqrt(7oneType) + 147)))/22680) (sqrt(42sqrt(7oneType) + 147) / 60 - sqrt(6sqrt(7oneType) + 21) / 36 - sqrt(7oneType) / 120+7oneType / 60) (7oneType / 30-sqrt(7oneType) / 60); + 0 (-((sqrt(3oneType) + 2sqrt(21oneType)) * (21sqrt(2sqrt(7oneType) + 7) - 63sqrt(21oneType) + 70sqrt(6sqrt(7oneType) + 21) + 24 * sqrt(14sqrt(7oneType) + 49) - 25sqrt(42sqrt(7oneType) + 147)))/22680) (sqrt(7oneType) / 120 - sqrt(147 - 42sqrt(7oneType)) / 1260+7oneType / 60) (sqrt(7oneType)/120+sqrt(21 - 6sqrt(7oneType))/36+sqrt(147 - 42sqrt(7oneType))/60+7oneType/60) (((sqrt(3oneType) + 2sqrt(21oneType)) * (63sqrt(21oneType) + 21sqrt(2sqrt(7oneType) + 7) - 70sqrt(6sqrt(7oneType) + 21) + 24 * sqrt(14sqrt(7oneType) + 49) + 25sqrt(42sqrt(7oneType) + 147)))/22680) (sqrt(7oneType) / 60+7oneType / 30); + 0 (-((sqrt(3oneType) + 2sqrt(21oneType)) * (21sqrt(2sqrt(7oneType) + 7) - 63sqrt(21oneType) - 70sqrt(6sqrt(7oneType) + 21) + 24 * sqrt(14sqrt(7oneType) + 49) + 25sqrt(42sqrt(7oneType) + 147)))/22680) (sqrt(7oneType) / 120 - sqrt(21 - 6sqrt(7oneType)) / 36 - sqrt(147 - 42sqrt(7oneType)) / 60+7oneType / 60) (((sqrt(7oneType) + 14) * (32 * sqrt(2sqrt(7oneType) + 7) - 10sqrt(14sqrt(7oneType) + 49) + 567))/68040) (((sqrt(3oneType) + 2sqrt(21oneType)) * (63sqrt(21oneType) + 21sqrt(2sqrt(7oneType) + 7) + 70sqrt(6sqrt(7oneType) + 21) + 24 * sqrt(14sqrt(7oneType) + 49) - 25sqrt(42sqrt(7oneType) + 147)))/22680) (sqrt(7oneType) / 60+7oneType / 30); + 0 (sqrt(6sqrt(7oneType) + 21) / 36 - sqrt(7oneType) / 120 - sqrt(42sqrt(7oneType) + 147) / 60+7oneType / 60) (((sqrt(3oneType) - 2sqrt(21oneType)) * (35sqrt(2sqrt(7oneType) + 7) - 63sqrt(21oneType) - 70sqrt(6sqrt(7oneType) + 21) + 20sqrt(14sqrt(7oneType) + 49) + 23sqrt(42sqrt(7oneType) + 147)))/22680) (((sqrt(3oneType) - 2sqrt(21oneType)) * (35sqrt(2sqrt(7oneType) + 7) - 63sqrt(21oneType) + 70sqrt(6sqrt(7oneType) + 21) + 20sqrt(14sqrt(7oneType) + 49) - 23sqrt(42sqrt(7oneType) + 147)))/22680) (sqrt(42sqrt(7oneType) + 147) / 1260 - sqrt(7oneType) / 120+7oneType / 60) (7oneType / 30-sqrt(7oneType) / 60); + 0 (sqrt(6sqrt(7oneType) + 21) / 540 - sqrt(7oneType) / 756 - sqrt(42sqrt(7oneType) + 147) / 3780-oneType / 540) (sqrt(7oneType) / 756 - sqrt(21 - 6sqrt(7oneType)) / 540 - sqrt(147 - 42sqrt(7oneType)) / 3780-oneType / 540) (sqrt(7oneType) / 756 + sqrt(21 - 6sqrt(7oneType)) / 540 + sqrt(147 - 42sqrt(7oneType)) / 3780-oneType / 540) (sqrt(42sqrt(7oneType) + 147) / 3780 - sqrt(6sqrt(7oneType) + 21) / 540 - sqrt(7oneType) / 756-oneType / 540) oneType/30] + else + error("MPDeC requires 2 ≤ K ≤ 10.") + end + end + return nodes, theta +end + +struct MPDeCConstantCache{NType, T, T2} <: OrdinaryDiffEqConstantCache + K::Int + M::Int + nodes::NType + theta::T2 + small_constant::T +end + +# Out-of-place +function alg_cache(alg::MPDeC, 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("MPDeC can only be applied to production-destruction systems")) + end + + nodes, theta = get_constant_parameters(alg, uEltypeNoUnits) + MPDeCConstantCache(alg.K, alg.M, nodes, theta, + alg.small_constant_function(uEltypeNoUnits)) +end + +function initialize!(integrator, cache::MPDeCConstantCache) +end + +# out-of-place +function build_mpdec_matrix_and_rhs_oop(uprev, m, f, C, p, t, dt, nodes, theta, + small_constant) + N = length(uprev) + if f isa PDSFunction + # Additional destruction terms + Mmat, rhs = _build_mpdec_matrix_and_rhs_oop(uprev, m, f.p, C, p, t, dt, nodes, + theta, + small_constant, f.d) + else + # No additional destruction terms + Mmat, rhs = _build_mpdec_matrix_and_rhs_oop(uprev, m, f.p, C, p, t, dt, nodes, + theta, + small_constant) + end + + if uprev isa StaticArray + return SMatrix{N, N}(Mmat), SVector{N}(rhs) + else + return Mmat, rhs + end +end + +# out-of-place for dense arrays +@muladd function _build_mpdec_matrix_and_rhs_oop(uprev, m, prod, C, p, t, dt, nodes, theta, + small_constant, dest = nothing) + N, M_plus_1 = size(C) + + # Create linear system matrix and rhs + if uprev isa StaticArray + Mmat = MMatrix{N, N}(zeros(eltype(uprev), N, N)) + else + Mmat = zeros(eltype(uprev), N, N) + end + rhs = similar(uprev) + + # Initialize + oneMmat = one(eltype(Mmat)) + @inbounds for i in 1:N + Mmat[i, i] = oneMmat + end + rhs .= uprev + + σ = add_small_constant(C[:, m], small_constant) + + @fastmath @inbounds @simd for r in 1:M_plus_1 + th = theta[r, m] + dt_th = dt * th + P = prod(C[:, r], p, t + nodes[r] * dt) + if !isnothing(dest) + d = dest(C[:, r], p, t + nodes[r] * dt) + else + d = nothing + end + _build_mpdec_matrix_and_rhs!(Mmat, rhs, P, dt_th, σ, d) + end + + return Mmat, rhs +end + +# in-place for dense arrays +@muladd function build_mpdec_matrix_and_rhs_ip!(Mmat, rhs, m, prod, P, C, p, t, dt, σ, tmp, + nodes, theta, small_constant, + dest = nothing, d = nothing) + N, M_plus_1 = size(C) + + oneMmat = one(eltype(Mmat)) + zeroMmat = zero(eltype(Mmat)) + + # Initialize Mmat as identity matrix + if Mmat isa Tridiagonal + Mmat.d .= oneMmat + Mmat.du .= zeroMmat + Mmat.dl .= zeroMmat + elseif issparse(Mmat) + # Fill sparse matrix with zeros without changing the sparsity pattern, see + # https://github.com/JuliaSparse/SparseArrays.jl/issues/190#issuecomment-1186690008. + fill!(Mmat.nzval, false) + + M_rows = rowvals(Mmat) + M_vals = nonzeros(Mmat) + for j in 1:N + for idx_M in nzrange(Mmat, j) + i = M_rows[idx_M] + if i == j + M_vals[idx_M] = oneMmat + end + end + end + else + fill!(Mmat, zeroMmat) + @inbounds for i in 1:N + Mmat[i, i] = oneMmat + end + end + + σ .= C[:, m] .+ small_constant + + @fastmath @inbounds @simd for r in 1:M_plus_1 + th = theta[r, m] + dt_th = dt * th + + prod(P, C[:, r], p, t + nodes[r] * dt) + if !isnothing(dest) + dest(d, C[:, r], p, t + nodes[r] * dt) + end + + if issparse(Mmat) + _build_mpdec_matrix_and_rhs!(Mmat, rhs, P, dt_th, σ, d, tmp) + else + _build_mpdec_matrix_and_rhs!(Mmat, rhs, P, dt_th, σ, d) + end + end +end + +function _build_mpdec_matrix_and_rhs!(M, rhs, P, dt_th, σ, d = nothing) + Base.require_one_based_indexing(M, P, σ) + @assert size(M, 1) == size(M, 2) == size(P, 1) == size(P, 2) == length(σ) + + if dt_th ≥ 0 + @fastmath @inbounds @simd for I in CartesianIndices(P) + if !iszero(P[I]) + dt_th_P = dt_th * P[I] + if I[1] != I[2] + M[I] -= dt_th_P / σ[I[2]] + M[I[2], I[2]] += dt_th_P / σ[I[2]] + else # diagonal elements + rhs[I[1]] += dt_th_P + end + end + end + + if !isnothing(d) + @fastmath @inbounds @simd for i in eachindex(d) + if !iszero(d[i]) + M[i, i] += dt_th * d[i] / σ[i] + end + end + end + else # dt_th ≤ 0 + @fastmath @inbounds @simd for I in CartesianIndices(P) + if !iszero(P[I]) + dt_th_P = dt_th * P[I] + if I[1] != I[2] + M[I[2], I[1]] += dt_th_P / σ[I[1]] + M[I[1], I[1]] -= dt_th_P / σ[I[1]] + else # diagonal elements + M[I] -= dt_th_P / σ[I[1]] + end + end + end + + if !isnothing(d) + @fastmath @inbounds @simd for i in eachindex(d) + if !iszero(d[i]) + rhs[i] -= dt_th * d[i] + end + end + end + end +end + +# optimized version for Tridiagonal matrices +function _build_mpdec_matrix_and_rhs!(M::Tridiagonal, rhs, P::Tridiagonal, dt_th, σ, + d = nothing) + Base.require_one_based_indexing(M.dl, M.d, M.du, P.dl, P.d, P.du, σ) + @assert length(M.dl) + 1 == length(M.d) == length(M.du) + 1 == + length(P.dl) + 1 == length(P.d) == length(P.du) + 1 == length(σ) + + if dt_th ≥ 0 + @fastmath @inbounds @simd for i in eachindex(P.d, rhs) + rhs[i] += dt_th * P.d[i] + end + + for i in eachindex(M.dl, P.dl) + dt_th_P = dt_th * P.dl[i] + M.dl[i] -= dt_th_P / σ[i] + M.d[i] += dt_th_P / σ[i] + end + + for i in eachindex(M.du, P.du) + dt_th_P = dt_th * P.du[i] + M.du[i] -= dt_th_P / σ[i + 1] + M.d[i + 1] += dt_th_P / σ[i + 1] + end + + if !isnothing(d) + @fastmath @inbounds @simd for i in eachindex(M.d, σ, d) + M.d[i] += dt_th * d[i] / σ[i] + end + end + else # dt_th ≤ 0 + @fastmath @inbounds @simd for i in eachindex(M.d, P.d, σ) + M.d[i] -= dt_th * P.d[i] / σ[i] + end + + for i in eachindex(M.dl, P.dl) + dt_th_P = dt_th * P.dl[i] + M.du[i] += dt_th_P / σ[i + 1] + M.d[i + 1] -= dt_th_P / σ[i + 1] + end + + for i in eachindex(M.du, P.du) + dt_th_P = dt_th * P.du[i] + M.dl[i] += dt_th_P / σ[i] + M.d[i] -= dt_th_P / σ[i] + end + + if !isnothing(d) + @fastmath @inbounds @simd for i in eachindex(rhs, d) + rhs[i] -= dt_th * d[i] + end + end + end +end + +# optimized version for sparse matrices +function _build_mpdec_matrix_and_rhs!(M::AbstractSparseMatrix, rhs, P::AbstractSparseMatrix, + dt_th, σ, + d = nothing, tmp = nothing) + Base.require_one_based_indexing(M, P, σ) + @assert size(M, 1) == size(M, 2) == size(P, 1) == size(P, 2) == length(σ) + if !isnothing(d) + Base.require_one_based_indexing(d) + @assert length(σ) == length(d) + end + + M_rows = rowvals(M) + M_vals = nonzeros(M) + P_rows = rowvals(P) + P_vals = nonzeros(P) + n = size(M, 2) + + # We use tmp as a buffer for the diagonal elements of M. + fill!(tmp, zero(eltype(tmp))) + + if dt_th ≥ 0 + for j in 1:n # run through columns of P + for idx_P in nzrange(P, j) # run through rows of P + i = P_rows[idx_P] + dt_th_P = dt_th * P_vals[idx_P] + if i != j + for idx_M in nzrange(M, j) + if M_rows[idx_M] == i + M_vals[idx_M] -= dt_th_P / σ[j] # M_ij <- P_ij + break + end + end + tmp[j] += dt_th_P / σ[j] # M_jj <- P_ij = D_ji + else + rhs[i] += dt_th_P # rhs_i <- P_ii + end + end + end + + if !isnothing(d) + for i in eachindex(d) + tmp[i] += dt_th * d[i] / σ[i] # M_ii <- D_i + end + end + + for j in 1:n + for idx_M in nzrange(M, j) + i = M_rows[idx_M] + if i == j + M_vals[idx_M] += tmp[j] + break + end + end + end + else # dt ≤ 0 + for j in 1:n # j is column index + for idx_P in nzrange(P, j) + i = P_rows[idx_P] # i is row index + dt_th_P = dt_th * P_vals[idx_P] + if i != j + for idx_M in nzrange(M, i) + if M_rows[idx_M] == j + M_vals[idx_M] += dt_th_P / σ[i] # M_ji <- P_ij + break + end + end + tmp[i] -= dt_th_P / σ[i] + else + for idx_M in nzrange(M, j) + if i == M_rows[idx_M] + M_vals[idx_M] -= dt_th_P / σ[i] # M_ij <- P_ij + break + end + end + end + end + end + + for j in 1:n + for idx_M in nzrange(M, j) + i = M_rows[idx_M] + if i == j + M_vals[idx_M] += tmp[j] + break + end + end + end + + if !isnothing(d) + @.. broadcast=false rhs-=dt_th * d + end + end +end + +@muladd function perform_step!(integrator, cache::MPDeCConstantCache, repeat_step = false) + @unpack alg, t, dt, uprev, f, p = integrator + @unpack K, M, nodes, theta, small_constant = cache + + N = length(uprev) + + if uprev isa StaticArray + C = MMatrix{N, M + 1}(zeros(eltype(uprev), N, M + 1)) + C2 = MMatrix{N, M + 1}(zeros(eltype(uprev), N, M + 1)) + else + C = zeros(eltype(uprev), N, M + 1) + C2 = zeros(eltype(uprev), N, M + 1) + end + + for i in 1:(M + 1) + C2[:, i] = uprev + end + + for _ in 1:K + C .= C2 + for m in 2:(M + 1) + Mmat, rhs = build_mpdec_matrix_and_rhs_oop(uprev, m, f, C, p, t, dt, nodes, + theta, + small_constant) + # solve linear system + linprob = LinearProblem(Mmat, rhs) + sol = solve(linprob, alg.linsolve) + C2[:, m] = sol.u + integrator.stats.nsolve += 1 + end + end + u = C2[:, M + 1] + u1 = C[:, M + 1] # one order less accurate + + tmp = u - u1 + 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 MPDeCCache{uType, PType, CType, tabType, F} <: MPRKMutableCache + tmp::uType + P::PType + P2::PType + d::uType + σ::uType + C::CType + C2::CType + tab::tabType + linsolve_rhs::uType + linsolve::F +end + +struct MPDeCConservativeCache{uType, PType, CType, tabType, F} <: MPRKMutableCache + tmp::uType + P::PType + P2::PType + σ::uType + C::CType + C2::CType + tab::tabType + linsolve_rhs::uType + linsolve::F +end + +get_tmp_cache(integrator, ::MPDeC, cache::OrdinaryDiffEqMutableCache) = (cache.σ,) + +# In-place +function alg_cache(alg::MPDeC, u, rate_prototype, ::Type{uEltypeNoUnits}, + ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, + uprev, uprev2, f, t, dt, reltol, p, calck, + ::Val{true}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} + nodes, theta = get_constant_parameters(alg, uEltypeNoUnits) + tab = MPDeCConstantCache(alg.K, alg.M, nodes, theta, + alg.small_constant_function(uEltypeNoUnits)) + + tmp = zero(u) + P = p_prototype(u, f) # stores evaluation of the production matrix + P2 = p_prototype(u, f) # stores the linear system matrix + if issparse(P2) + # We need to ensure that evaluating the production function does + # not alter the sparsity pattern given by the production matrix prototype + f.p(P2, uprev, p, t) + @assert P.rowval == P2.rowval&&P.colptr == P2.colptr "Evaluation of the production terms must not alter the sparsity pattern given by the prototype." + + if alg.K > 2 + # Negative weights of MPDeC(K) , K > 2 require + # a symmetric sparsity pattern of the linear system matrix P2 + P2 = P + P' + end + end + d = zero(u) + σ = zero(u) + C = zeros(eltype(u), length(u), alg.M + 1) + C2 = zeros(eltype(u), length(u), alg.M + 1) + linsolve_rhs = zero(u) + + if f isa ConservativePDSFunction + # The right hand side of the linear system is always uprev. But using + # linsolve_rhs 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(P2, _vec(linsolve_rhs)) + linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true, + assumptions = LinearSolve.OperatorAssumptions(true)) + + MPDeCConservativeCache(tmp, P, P2, σ, C, C2, + tab, #MPDeCConstantCache + linsolve_rhs, + linsolve) + elseif f isa PDSFunction + linprob = LinearProblem(P2, _vec(linsolve_rhs)) + linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true, + assumptions = LinearSolve.OperatorAssumptions(true)) + + MPDeCCache(tmp, P, P2, d, σ, C, C2, + tab, #MPDeCConstantCache + linsolve_rhs, + linsolve) + else + throw(ArgumentError("MPDeC can only be applied to production-destruction systems")) + end +end + +function initialize!(integrator, cache::Union{MPDeCCache, MPDeCConservativeCache}) +end + +@muladd function perform_step!(integrator, cache::MPDeCCache, repeat_step = false) + @unpack t, dt, uprev, u, f, p = integrator + @unpack tmp, P, P2, d, σ, C, C2, linsolve_rhs, linsolve = cache + @unpack K, M, nodes, theta, small_constant = cache.tab + + # Initialize C matrices + for i in 1:(M + 1) + C2[:, i] .= uprev + end + + for _ in 1:K + C .= C2 + for m in 2:(M + 1) + linsolve_rhs .= uprev + build_mpdec_matrix_and_rhs_ip!(P2, linsolve_rhs, m, f.p, P, C, p, t, dt, σ, tmp, + nodes, + theta, + small_constant, f.d, d) + + # Same as linres = P2 \ linsolve_rhs + linsolve.A = P2 + linres = solve!(linsolve) + C2[:, m] .= linres + integrator.stats.nsolve += 1 + end + end + + u .= C2[:, M + 1] + σ .= C[:, M + 1] # one order less accurate + + # 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 + +@muladd function perform_step!(integrator, cache::MPDeCConservativeCache, + repeat_step = false) + @unpack t, dt, uprev, u, f, p = integrator + @unpack tmp, P, P2, σ, C, C2, linsolve_rhs, linsolve = cache + @unpack K, M, nodes, theta, small_constant = cache.tab + + # Initialize right hand side of linear system + linsolve_rhs .= uprev + + # Initialize C matrices + for i in 1:(M + 1) + C2[:, i] .= uprev + end + + for _ in 1:K + C .= C2 + for m in 2:(M + 1) + linsolve_rhs .= uprev + build_mpdec_matrix_and_rhs_ip!(P2, linsolve_rhs, m, f.p, P, C, p, t, dt, σ, tmp, + nodes, + theta, + small_constant) + + # Same as linres = P2 \ linsolve_rhs + linsolve.A = P2 + linres = solve!(linsolve) + C2[:, m] .= linres + integrator.stats.nsolve += 1 + end + end + + u .= C2[:, M + 1] + σ .= C[:, M + 1] # one order less accurate + + # 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 diff --git a/test/Project.toml b/test/Project.toml index 65a43072..7a9ef3ff 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -1,6 +1,7 @@ [deps] ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b" Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" +DoubleFloats = "497a8b3b-efae-58df-a0af-a86822472b78" ExplicitImports = "7d51a73a-1435-4ff3-83d9-f097790105c7" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae" @@ -21,6 +22,7 @@ Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" [compat] ADTypes = "1.12" Aqua = "0.8" +DoubleFloats = "1.4.3" ExplicitImports = "1.0.1" LinearAlgebra = "1" LinearSolve = "3.7.1" diff --git a/test/runtests.jl b/test/runtests.jl index 48be033e..64b6c2fa 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -3,7 +3,8 @@ using LinearAlgebra using SparseArrays using Statistics: mean, median -using StaticArrays: MVector, @SVector, SA +using DoubleFloats: Double64 +using StaticArrays: SMatrix, MVector, @SVector, SVector, SA using Unitful: @u_str, ustrip @@ -40,7 +41,7 @@ solution is computed using `ref_alg`. """ function experimental_orders_of_convergence(prob, alg, dts; test_time = nothing, only_first_index = false, - ref_alg = TRBDF2(autodiff = AutoFiniteDiff())) + ref_alg = Vern7()) @assert length(dts) > 1 errors = zeros(eltype(dts), length(dts)) @@ -354,22 +355,22 @@ end 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] + P[i, i + 1] = i * u[i + 1] 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] + P[i + 1, 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] + P[i, i + 1] = i * u[i + 1] + P[i + 1, i] = i * u[i] end return nothing end @@ -1103,6 +1104,21 @@ end @test_throws "SSPMPRK43 can only be applied to production-destruction systems" solve(prob_ip, SSPMPRK43(), dt = 0.1) + @test_throws "MPDeC can only be applied to production-destruction systems" solve(prob_ip, + MPDeC(2), + dt = 0.1) + @test_throws "MPDeC can only be applied to production-destruction systems" solve(prob_oop, + MPDeC(2), + dt = 0.1) + @test_throws "MPDeC requires 2 ≤ K ≤ 10." solve(prob_pds_linmod, MPDeC(123)) + @test_throws "MPDeC requires 2 ≤ K ≤ 10." solve(prob_pds_linmod, + MPDeC(123, nodes = :lagrange)) + P = spdiagm(1 => [1.0]) + function prod!(P, u, p, t) + P[2, 1] = one(eltype(P)) + end + prob = ConservativePDSProblem(prod!, ones(2, 1), (0.0, 1.0); p_prototype = P) + @test_throws "must not alter the sparsity pattern" solve(prob, MPDeC(2)) end # Here we check that algorithms which accept input parameters return constants @@ -1132,7 +1148,8 @@ end algs = (MPRK22(0.5f0), MPRK22(1.0f0), MPRK22(2.0f0), MPRK43I(1.0f0, 0.5f0), MPRK43I(0.5f0, 0.75f0), MPRK43II(0.5f0), MPRK43II(2.0f0 / 3.0f0), - SSPMPRK22(0.5f0, 1.0f0), SSPMPRK43()) + SSPMPRK22(0.5f0, 1.0f0), SSPMPRK43(), MPDeC(2), + MPDeC(2, nodes = :lagrange)) for alg in algs sol = solve(prob, alg; dt = 0.1f0, save_everystep = false) @test sol.t isa Vector{Float32} @@ -1213,6 +1230,28 @@ end end end + # Here we check that MPRK22(1.0) = MPDeC(2) + @testset "MPRK22(1.0) = MPDeC(2)" begin + # conservative PDS + sol1 = solve(prob_pds_linmod, MPRK22(1.0)) + sol2 = solve(prob_pds_linmod, MPDeC(2)) + sol3 = solve(prob_pds_linmod, MPDeC(2, nodes = :lagrange)) + sol4 = solve(prob_pds_linmod_inplace, MPRK22(1.0)) + sol5 = solve(prob_pds_linmod_inplace, MPDeC(2)) + sol6 = solve(prob_pds_linmod_inplace, MPDeC(2, nodes = :lagrange)) + @test sol1.u ≈ sol2.u ≈ sol3.u ≈ sol4.u ≈ sol5.u ≈ sol6.u + + # nonconservative PDS + sol1 = solve(prob_pds_linmod_nonconservative, MPRK22(1.0)) + sol2 = solve(prob_pds_linmod_nonconservative, MPDeC(2)) + sol3 = solve(prob_pds_linmod_nonconservative, MPDeC(2, nodes = :lagrange)) + sol4 = solve(prob_pds_linmod_nonconservative_inplace, MPRK22(1.0)) + sol5 = solve(prob_pds_linmod_nonconservative_inplace, MPDeC(2)) + sol6 = solve(prob_pds_linmod_nonconservative_inplace, + MPDeC(2, nodes = :lagrange)) + @test sol1.u ≈ sol2.u ≈ sol3.u ≈ sol4.u ≈ sol5.u ≈ sol6.u + end + # Here we check that different linear solvers can be used @testset "Different linear solvers" begin # problem data @@ -1244,15 +1283,19 @@ end prob_ip_2 = ConservativePDSProblem(linmodP!, u0, tspan, p; analytic = f_analytic) - algs = (MPE, (; kwargs...) -> MPRK22(1.0; kwargs...), - (; kwargs...) -> MPRK22(0.5; kwargs...), - (; kwargs...) -> MPRK22(2.0; kwargs...), - (; 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...) -> SSPMPRK22(0.5, 1.0; kwargs...), - (; kwargs...) -> SSPMPRK43(; kwargs...)) + algs = [MPE, (; kwargs...) -> MPRK22(1.0; kwargs...), + (; kwargs...) -> MPRK22(0.5; kwargs...), + (; kwargs...) -> MPRK22(2.0; kwargs...), + (; 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...) -> SSPMPRK22(0.5, 1.0; kwargs...), + (; kwargs...) -> SSPMPRK43(; kwargs...)] + for k in 2:10 + push!(algs, (; kwargs...) -> MPDeC(k; kwargs...), + (; kwargs...) -> MPDeC(k; nodes = :lagrange, kwargs...)) + end for alg in algs # Check different linear solvers @@ -1278,7 +1321,7 @@ end 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] + P[i, i + 1] = i * u[i + 1] end return nothing end @@ -1286,7 +1329,7 @@ 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] + P[i + 1, i] = i * u[i] end return nothing end @@ -1294,8 +1337,8 @@ 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] + P[i, i + 1] = i * u[i + 1] + P[i + 1, i] = i * u[i] end return nothing end @@ -1310,11 +1353,21 @@ end tspan = (0.0, 1.0) dt = 0.25 - @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()) + algs = [ + 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 k in 2:10 + push!(algs, MPDeC(k), MPDeC(k; nodes = :lagrange)) + end + @testset "$alg" for alg in algs for prod! in (prod_1!, prod_2!, prod_3!) prod = (u, p, t) -> begin P = similar(u, (length(u), length(u))) @@ -1362,7 +1415,7 @@ end 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] + P[i, i + 1] = i * u[i + 1] end return nothing end @@ -1370,7 +1423,7 @@ 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] + P[i + 1, i] = i * u[i] end return nothing end @@ -1378,8 +1431,8 @@ 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] + P[i, i + 1] = i * u[i + 1] + P[i + 1, i] = i * u[i] end return nothing end @@ -1396,11 +1449,14 @@ end 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()) + algs = [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)] + for k in 2:10 + push!(algs, MPDeC(k), MPDeC(k; nodes = :lagrange)) + end + @testset "$alg" for alg in algs for prod! in (prod_1!, prod_2!, prod_3!) prod = (u, p, t) -> begin P = similar(u, (length(u), length(u))) @@ -1420,12 +1476,12 @@ end 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) + sol_tridiagonal_op = solve(prob_tridiagonal_op, alg) + sol_dense_ip = solve(prob_dense_ip, alg) + sol_dense_op = solve(prob_dense_op, alg) + sol_sparse_ip = solve(prob_sparse_ip, alg) + sol_sparse_op = solve(prob_sparse_op, alg) @test isapprox(sol_tridiagonal_ip.t, sol_tridiagonal_op.t; rtol) @test isapprox(sol_dense_ip.t, sol_dense_op.t; rtol) @@ -1448,7 +1504,7 @@ end 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] + P[i, i + 1] = i * u[i + 1] end for i in 1:length(u) P[i, i] = i * u[i] @@ -1466,7 +1522,7 @@ 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] + P[i + 1, i] = i * u[i] end for i in 1:length(u) P[i, i] = (i - 1) * u[i] @@ -1484,8 +1540,8 @@ 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] + P[i, i + 1] = i * u[i + 1] + P[i + 1, i] = i * u[i] end for i in 1:length(u) P[i, i] = (i + 1) * u[i] @@ -1511,11 +1567,16 @@ end tspan = (0.0, 1.0) dt = 0.25 - @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()) + algs = [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 k in 2:10 + push!(algs, MPDeC(k), MPDeC(k; nodes = :lagrange)) + end + + @testset "$alg" for alg in algs for (prod!, dest!) in zip((prod_1!, prod_2!, prod_3!), (dest_1!, dest_2!, dest_3!)) prod = (u, p, t) -> begin @@ -1573,7 +1634,7 @@ end 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] + P[i, i + 1] = i * u[i + 1] end for i in 1:length(u) P[i, i] = i * u[i] @@ -1591,7 +1652,7 @@ 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] + P[i + 1, i] = i * u[i] end for i in 1:length(u) P[i, i] = (i - 1) * u[i] @@ -1609,8 +1670,8 @@ 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] + P[i, i + 1] = i * u[i + 1] + P[i + 1, i] = i * u[i] end for i in 1:length(u) P[i, i] = (i + 1) * u[i] @@ -1636,12 +1697,16 @@ end tspan = (0.0, 1.0) dt = 0.25 + algs = [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)] + for k in 2:10 + push!(algs, MPDeC(k), MPDeC(k; nodes = :lagrange)) + end + 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()) + @testset "$alg" for alg in algs for (prod!, dest!) in zip((prod_1!, prod_2!, prod_3!), (dest_1!, dest_2!, dest_3!)) prod! = prod_3! @@ -1669,12 +1734,12 @@ end 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) + sol_tridiagonal_ip = solve(prob_tridiagonal_ip, alg) + sol_tridiagonal_op = solve(prob_tridiagonal_op, alg) + sol_dense_ip = solve(prob_dense_ip, alg) + sol_dense_op = solve(prob_dense_op, alg) + sol_sparse_ip = solve(prob_sparse_ip, alg) + sol_sparse_op = solve(prob_sparse_op, alg) @test isapprox(sol_tridiagonal_ip.t, sol_tridiagonal_op.t; rtol) @test isapprox(sol_dense_ip.t, sol_dense_op.t; rtol) @@ -1699,7 +1764,7 @@ end prod_inner! = (P, u, p, t) -> begin fill!(P, zero(eltype(P))) for i in 1:(length(u) - 1) - P[i, i + 1] = i * u[i] + P[i, i + 1] = i * u[i + 1] end return nothing end @@ -1747,11 +1812,17 @@ end p_prototype = P_dense) prob_sparse2 = PDSProblem(prod_sparse!, dest!, u0, tspan; p_prototype = P_sparse) - # solve and test - for alg in (MPE(), MPRK22(0.5), MPRK22(1.0), MPRK43I(1.0, 0.5), - MPRK43I(0.5, 0.75), - MPRK43II(2.0 / 3.0), MPRK43II(0.5), SSPMPRK22(0.5, 1.0), - SSPMPRK43()) + + algs = [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 k in 2:10 + push!(algs, MPDeC(k), MPDeC(k; nodes = :lagrange)) + end + #solve and test + for alg in algs for prob in (prob_default, prob_tridiagonal, prob_dense, prob_sparse, prob_default2, prob_tridiagonal2, prob_dense2, prob_sparse2) @@ -1775,7 +1846,8 @@ 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), SSPMPRK22(0.5, 1.0)) + algs = (MPE(), MPRK22(0.5), MPRK22(1.0), MPRK22(2.0), SSPMPRK22(0.5, 1.0), + MPDeC(2), MPDeC(2, nodes = :lagrange)) dts = 0.5 .^ (4:15) problems = (prob_pds_linmod, prob_pds_linmod_array, prob_pds_linmod_mvector, prob_pds_linmod_inplace) @@ -1811,12 +1883,12 @@ end # 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), + MPDeC(2), MPDeC(2; nodes = :lagrange)) dts = 0.5 .^ (4:15) problems = (prob_pds_linmod_nonconservative, prob_pds_linmod_nonconservative_inplace) @testset "$alg" for alg in algs - alg = MPRK22(1.0) for prob in problems orders = experimental_orders_of_convergence(prob, alg, dts) @test check_order(orders, PositiveIntegrators.alg_order(alg)) @@ -1830,13 +1902,13 @@ end dts; test_time) @test check_order(orders, PositiveIntegrators.alg_order(alg), - atol = 0.2) + atol = 0.3) 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) + N = 2, atol = 0.3) end end end @@ -1849,7 +1921,8 @@ end 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()) + MPRK43II(0.5), MPRK43II(2.0 / 3.0), SSPMPRK43(), MPDeC(3), + MPDeC(3, nodes = :lagrange)) 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) @@ -1863,7 +1936,9 @@ end 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()) + MPRK43II(0.5), MPRK43II(2.0 / 3.0), SSPMPRK43(), + MPDeC(3), MPDeC(3; nodes = :lagrange), + MPDeC(4), MPDeC(4; nodes = :lagrange)) 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) @@ -1873,7 +1948,13 @@ 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), - SSPMPRK22(0.5, 1.0), SSPMPRK43()) + SSPMPRK22(0.5, 1.0), SSPMPRK43(), + MPDeC(2), MPDeC(2, nodes = :gausslobatto), + MPDeC(3), MPDeC(3, nodes = :gausslobatto), + MPDeC(4), MPDeC(4, nodes = :gausslobatto), + MPDeC(5), MPDeC(5, nodes = :gausslobatto), + MPDeC(6), MPDeC(6, nodes = :gausslobatto), + MPDeC(7), MPDeC(7, nodes = :gausslobatto)) dt = 0.5^6 problems = (prob_pds_linmod, prob_pds_linmod_array, prob_pds_linmod_mvector, prob_pds_linmod_inplace) @@ -1892,7 +1973,7 @@ end prod! = (P, u, p, t) -> begin fill!(P, zero(eltype(P))) P[1, 2] = sin(t)^2 * u[2] - P[2, 1] = cos(2 * t)^2 * u[1] + P[2, 1] = cos(t)^2 * u[1] return nothing end prod = (u, p, t) -> begin @@ -1900,14 +1981,16 @@ end prod!(P, u, p, t) return P end - u0 = [1.0; 0.0] + u0 = [1.1; 0.9] # values close to zero may decrease the order tspan = (0.0, 1.0) prob_oop = ConservativePDSProblem(prod, u0, tspan) #out-of-place prob_ip = ConservativePDSProblem(prod!, u0, tspan) #in-place dts = 0.5 .^ (4:15) algs = (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()) + MPRK43II(2.0 / 3.0), MPRK43II(0.5), SSPMPRK22(0.5, 1.0), SSPMPRK43(), + MPDeC(2), MPDeC(2, nodes = :lagrange), MPDeC(3), + MPDeC(3, nodes = :lagrange)) @testset "$alg" for alg in algs orders = experimental_orders_of_convergence(prob_oop, alg, dts) @test check_order(orders, PositiveIntegrators.alg_order(alg), atol = 0.2) @@ -1940,14 +2023,18 @@ end dest!(d, u, p, t) return d end - u0 = [1.0; 0.0] + u0 = [1.1; 0.9] # stay away from zero tspan = (0.0, 1.0) prob_oop = PDSProblem(prod, dest, u0, tspan) #out-of-place prob_ip = PDSProblem(prod!, dest!, u0, tspan) #in-place - dts = 0.5 .^ (4:15) + dts = 0.5 .^ (4:10) algs = (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()) + MPRK43II(2.0 / 3.0), MPRK43II(0.5), SSPMPRK22(0.5, 1.0), SSPMPRK43(), + MPDeC(2), MPDeC(2; nodes = :lagrange), MPDeC(3), + MPDeC(3; nodes = :lagrange), + MPDeC(4), MPDeC(4; nodes = :lagrange), MPDeC(5), + MPDeC(5; nodes = :lagrange)) @testset "$alg" for alg in algs orders = experimental_orders_of_convergence(prob_oop, alg, dts) @test check_order(orders, PositiveIntegrators.alg_order(alg), atol = 0.2) @@ -1956,6 +2043,179 @@ end end end + @testset "Convergence of higher order schemes (autonomous)" begin + function f_analytic(u0, p, t) + u₁⁰, u₂⁰ = u0 + a, b = p + c = a + b + return ((u₁⁰ + u₂⁰) * [b; a] + + exp(-c * t) * (a * u₁⁰ - b * u₂⁰) * [1; -1]) / c + end + + # linear model problem - conservative + function linmodP!(P, u, p, t) + P[1, 2] = u[2] + P[2, 1] = 5u[1] + return nothing + end + function linmodP(u, p, t) + P = zeros(eltype(u), 2, 2) + linmodP!(P, u, p, t) + return SMatrix{2, 2}(P) + end + + # linear model problem - nonconservative - in-place + function linmodP_noncons!(P, u, p, t) + P[1, 1] = u[2] / 2 + P[1, 2] = u[2] / 2 + P[2, 1] = 3u[1] + P[2, 2] = 2u[1] + return nothing + end + function linmodD_noncons!(D, u, p, t) + D[1] = 2u[1] + D[2] = u[2] / 2 + return nothing + end + function linmodP_noncons(u, p, t) + P = zeros(eltype(u), 2, 2) + linmodP_noncons!(P, u, p, t) + return SMatrix{2, 2}(P) + end + function linmodD_noncons(u, p, t) + D = zeros(eltype(u), 2, 1) + linmodD_noncons!(D, u, p, t) + return SVector{2}(D) + end + + u0 = [Double64(9) / 10; Double64(1) / 10] + p = [Double64(5); Double64(1)] + tspan = (Double64(0), Double64(1)) + + prob_ip = ConservativePDSProblem(linmodP!, u0, tspan, p; analytic = f_analytic) + prob_op = ConservativePDSProblem(linmodP, SVector{2}(u0), tspan, SVector{2}(p); + analytic = f_analytic) + prob_ip_noncons = PDSProblem(linmodP_noncons!, linmodD_noncons!, u0, tspan, p; + analytic = f_analytic) + prob_op_noncons = PDSProblem(linmodP_noncons, linmodD_noncons, SVector{2}(u0), + tspan, + SVector{2}(p); analytic = f_analytic) + + dts = 0.5 .^ (8:13) + + algs = MPDeC[] + for K in 4:10 + push!(algs, MPDeC(K)) + push!(algs, MPDeC(K, nodes = :lagrange)) + end + + atol = 0.3 + @testset "$alg" for alg in algs + orders = experimental_orders_of_convergence(prob_ip, alg, dts) + @test check_order(orders, PositiveIntegrators.alg_order(alg); atol) + + orders = experimental_orders_of_convergence(prob_op, alg, dts) + @test check_order(orders, PositiveIntegrators.alg_order(alg); atol) + + orders = experimental_orders_of_convergence(prob_ip_noncons, alg, dts) + @test check_order(orders, PositiveIntegrators.alg_order(alg); atol) + + orders = experimental_orders_of_convergence(prob_op_noncons, alg, dts) + @test check_order(orders, PositiveIntegrators.alg_order(alg); atol) + end + end + + @testset "Convergence of higher order schemes (nonautonomous)" begin + prod! = (P, u, p, t) -> begin + fill!(P, zero(eltype(P))) + P[1, 2] = sin(t)^2 * u[2] + P[2, 1] = cos(t)^2 * u[1] + return nothing + end + prod = (u, p, t) -> begin + P = similar(u, (length(u), length(u))) + prod!(P, u, p, t) + return SMatrix{2, 2}(P) + end + function fct!(du, u, p, t) + du[1] = sin(t)^2 * u[2] - cos(t)^2 * u[1] + du[2] = -sin(t)^2 * u[2] + cos(t)^2 * u[1] + return nothing + end + fct = (u, p, t) -> begin + f = similar(u) + fct!(f, u, p, t) + return SVector{2}(f) + end + + prod_noncons! = (P, u, p, t) -> begin + fill!(P, zero(eltype(P))) + P[1, 2] = sin(t)^2 * u[2] + P[2, 1] = cos(2 * t)^2 * u[1] + P[2, 2] = cos(t)^2 * u[2] + return nothing + end + dest_noncons! = (d, u, p, t) -> begin + fill!(d, zero(eltype(d))) + d[1] = sin(2 * t)^2 * u[1] + d[2] = sin(0.5 * t)^2 * u[2] + return nothing + end + prod_noncons = (u, p, t) -> begin + P = similar(u, (length(u), length(u))) + prod_noncons!(P, u, p, t) + return SMatrix{2, 2}(P) + end + dest_noncons = (u, p, t) -> begin + d = similar(u, (length(u),)) + dest_noncons!(d, u, p, t) + return SVector{2}(d) + end + fct_noncons! = (du, u, p, t) -> begin + du[1] = sin(t)^2 * u[2] - cos(2 * t)^2 * u[1] - sin(2 * t)^2 * u[1] + du[2] = -sin(t)^2 * u[2] + cos(2 * t)^2 * u[1] + cos(t)^2 * u[2] - + sin(0.5 * t)^2 * u[2] + return nothing + end + fct_noncons = (u, p, t) -> begin + f = similar(u) + fct_noncons!(f, u, p, t) + return SVector{2}(f) + end + + u0 = [Double64(11) / 10; Double64(9) / 10] # values close to zero may decrease the order + tspan = (Double64(0), Double64(1)) + prob_op = ConservativePDSProblem(prod, SVector{2}(u0), tspan; std_rhs = fct) #out-of-place + prob_ip = ConservativePDSProblem(prod!, u0, tspan; std_rhs = fct!) #in-place + prob_op_noncons = PDSProblem(prod_noncons, dest_noncons, SVector{2}(u0), tspan; + std_rhs = fct_noncons) #out-of-place + prob_ip_noncons = PDSProblem(prod_noncons!, dest_noncons!, u0, tspan; + std_rhs = fct_noncons!) #in-place + + dts = 0.5 .^ (5:10) + + algs = [] + for K in 4:10 + push!(algs, MPDeC(K)) + push!(algs, MPDeC(K, nodes = :lagrange)) + end + + atol = 0.3 + @testset "$alg" for alg in algs + orders = experimental_orders_of_convergence(prob_ip, alg, dts) + @test check_order(orders, PositiveIntegrators.alg_order(alg); atol) + + orders = experimental_orders_of_convergence(prob_op, alg, dts) + @test check_order(orders, PositiveIntegrators.alg_order(alg); atol) + + orders = experimental_orders_of_convergence(prob_ip_noncons, alg, dts) + @test check_order(orders, PositiveIntegrators.alg_order(alg); atol) + + orders = experimental_orders_of_convergence(prob_op_noncons, alg, dts) + @test check_order(orders, PositiveIntegrators.alg_order(alg); atol) + end + end + @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), @@ -2004,9 +2264,12 @@ end prob_oop = ConservativePDSProblem(prod, u0, tspan, p) prob_oop_2 = PDSProblem(prod, dest, u0, tspan, p) - 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), SSPMPRK43()) + 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), SSPMPRK43()] + for k in 2:10 + push!(algs, MPDeC(k), MPDeC(k; nodes = :lagrange)) + end for alg in algs sol = solve(prob_ip, alg; dt = dt, adaptive = false) @@ -2052,9 +2315,12 @@ end prob_oop = ConservativePDSProblem(prod, u0, tspan, p) prob_oop_2 = PDSProblem(prod, dest, u0, tspan, p) - 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), SSPMPRK43()) + 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), SSPMPRK43()] + for k in 2:10 + push!(algs, MPDeC(k), MPDeC(k; nodes = :lagrange)) + end dt = 1e-3 for alg in algs @@ -2073,8 +2339,11 @@ end # Here we check that the implemented schemes can solve the predefined PDS # (at least for specific parameters) @testset "PDS problem library (adaptive schemes)" begin - algs = (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)) + algs = [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)] + for k in 2:10 + push!(algs, MPDeC(k), MPDeC(k; nodes = :lagrange)) + end probs = (prob_pds_linmod, prob_pds_linmod_inplace, prob_pds_nonlinmod, prob_pds_robertson, prob_pds_bertolazzi, prob_pds_brusselator, prob_pds_npzd, @@ -2088,6 +2357,9 @@ end elseif prob == prob_pds_stratreac && alg == MPRK43I(0.5, 0.75) # Not successful on Julia 1.9 break + elseif prob == prob_pds_stratreac && alg == MPDeC(9; nodes = :lagrange) + # unstable + break end # later versions of OrdinaryDiffEq.jl use dtmin = 0 by default, # see https://github.com/SciML/OrdinaryDiffEq.jl/pull/2098 @@ -2149,11 +2421,15 @@ end prob_pds_npzd, prob_pds_sir, prob_pds_stratreac) - algs = (MPE, (; kwargs...) -> MPRK22(1.0; kwargs...), - (; kwargs...) -> MPRK43I(1.0, 0.5; kwargs...), - (; kwargs...) -> MPRK43II(0.5; kwargs...), - (; kwargs...) -> SSPMPRK22(0.5, 1.0; kwargs...), - (; kwargs...) -> SSPMPRK43(; kwargs...)) + algs = [MPE, (; kwargs...) -> MPRK22(1.0; kwargs...), + (; kwargs...) -> MPRK43I(1.0, 0.5; kwargs...), + (; kwargs...) -> MPRK43II(0.5; kwargs...), + (; kwargs...) -> SSPMPRK22(0.5, 1.0; kwargs...), + (; kwargs...) -> SSPMPRK43(; kwargs...)] + for k in 2:10 + push!(algs, (; kwargs...) -> MPDeC(k; kwargs...), + (; kwargs...) -> MPDeC(k; nodes = :lagrange, kwargs...)) + end for alg in algs for prob in probs sol1 = solve(prob_pds_linmod, alg(), dt = 0.1) @@ -2166,12 +2442,20 @@ end end end - #Here we check if the RK methods on which the MPRK schemes are based integrate + # Here we check if the RK methods on which the MPRK schemes are based integrate # u'(t) = q * t^(q-1) exactly for q from 1 to the order of the method. + # This is also true for MPDeC as long as the theta matrix is nonnegative, i.e. K = 2. + # Nevertheless, the results of most MPDeC schemes are good enough to pass this test @testset "Exact solutions (RK)" 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), SSPMPRK43()) + 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), SSPMPRK43()] + for k in 2:10 + push!(algs, MPDeC(k)) + if k != 9 + push!(algs, MPDeC(k, nodes = :lagrange)) + end + end @testset "$alg, $q" for alg in algs, q in 1:PositiveIntegrators.alg_order(alg) f(t) = q * t^(q - 1) function prod!(P, u, p, t)