Multirate: add implicit MRI-GARK-ESDIRK46a (order 4) (#961)

lib/OrdinaryDiffEqMultirate/Project.toml

@@ -1,7 +1,7 @@
 name = "OrdinaryDiffEqMultirate"
 uuid = "d4b830b4-ac80-426b-8507-16693d424963"
 authors = ["singhharsh1708 <hs1663531@gmail.com>"]
-version = "2.6.1"
+version = "2.7.0"
 
 [deps]
 ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"

lib/OrdinaryDiffEqMultirate/src/OrdinaryDiffEqMultirate.jl

@@ -29,6 +29,6 @@ include("multirate_caches.jl")
 include("multirate_perform_step.jl")
 
 export MREEF, MRAB, MRIGARKERK22a, MRIGARKERK22b, MRIGARKERK33a, MRIGARKERK45a,
-    MRIGARKIRK21a, MRIGARKESDIRK34a, MIS
+    MRIGARKIRK21a, MRIGARKESDIRK34a, MRIGARKESDIRK46a, MIS
 
 end

lib/OrdinaryDiffEqMultirate/src/alg_utils.jl

@@ -29,16 +29,15 @@ end
 
 alg_order(::MRIGARKIRK21a) = 2
 alg_order(::MRIGARKESDIRK34a) = 3
-isfsal(::Union{MRIGARKIRK21a, MRIGARKESDIRK34a}) = false
+alg_order(::MRIGARKESDIRK46a) = 4
+isfsal(::MRIGARKImplicitAlg) = false
 
-function prepare_alg(
-        alg::Union{MRIGARKIRK21a, MRIGARKESDIRK34a}, u0::AbstractArray, p, prob
-    )
+function prepare_alg(alg::MRIGARKImplicitAlg, u0::AbstractArray, p, prob)
     alg.m >= 1 || throw(ArgumentError("$(nameof(typeof(alg))): `m` must be ≥ 1"))
     return alg
 end
 
-nlsolve_f(f, ::Union{MRIGARKIRK21a, MRIGARKESDIRK34a}) = f isa SplitFunction ? f.f2 : f
+nlsolve_f(f, ::MRIGARKImplicitAlg) = f isa SplitFunction ? f.f2 : f
 
 alg_order(::MIS) = 2
 isfsal(::MIS) = false

lib/OrdinaryDiffEqMultirate/src/algorithms.jl

@@ -267,6 +267,55 @@ function MRIGARKESDIRK34a(;
     return MRIGARKESDIRK34a(m, linsolve, nlsolve, autodiff, _unwrap_val(concrete_jac))
 end
 
+@doc generic_solver_docstring(
+    "Multirate Infinitesimal GARK — 4th-order solve-decoupled implicit (MRI-GARK-ESDIRK46a).
+
+Solves a SplitODE `du/dt = f1(u,t) + f2(u,t)` where `f1` is the fast component
+and `f2` is the slow component (SciML convention). 11-stage, 4th-order
+solve-decoupled implicit MRI-GARK method of Sandu 2019: the fast component is
+integrated over five sub-intervals with `m` explicit-RK4 inner micro-steps, each
+followed by a `Δc = 0` stage that applies an ESDIRK slow correction requiring a
+nonlinear solve in the slow rate `f2`. The five implicit stages share the common
+diagonal coefficient `1/4`, so a single Jacobian/`W` is reused. Unlike the lower
+order members of the family this method carries a 3rd-order embedded pair, giving
+a proper local error estimate for adaptive stepping.",
+    "MRIGARKESDIRK46a",
+    "Multirate infinitesimal GARK solve-decoupled implicit method.",
+    """@article{sandu2019class,
+    title={A class of multirate infinitesimal {GARK} methods},
+    author={Sandu, Adrian},
+    journal={SIAM Journal on Numerical Analysis},
+    volume={57},
+    number={5},
+    pages={2300--2327},
+    year={2019}}""",
+    """
+    - `m`: number of inner RK4 micro-steps per fast stage.
+    """,
+    """
+    m::Int,
+    """
+)
+struct MRIGARKESDIRK46a{AD, F, F2, CJ} <: OrdinaryDiffEqNewtonAdaptiveAlgorithm
+    m::Int
+    linsolve::F
+    nlsolve::F2
+    autodiff::AD
+    concrete_jac::CJ
+end
+
+function MRIGARKESDIRK46a(;
+        m::Int, autodiff = AutoForwardDiff(), concrete_jac = nothing,
+        linsolve = nothing, nlsolve = NLNewton()
+    )
+    autodiff = _fixup_ad(autodiff)
+    return MRIGARKESDIRK46a(m, linsolve, nlsolve, autodiff, _unwrap_val(concrete_jac))
+end
+
+const MRIGARKImplicitAlg = Union{
+    MRIGARKIRK21a, MRIGARKESDIRK34a, MRIGARKESDIRK46a,
+}
+
 @doc generic_solver_docstring(
     "Multirate Infinitesimal Step (MIS).
 

lib/OrdinaryDiffEqMultirate/src/multirate_caches.jl

@@ -168,6 +168,7 @@ end
     kk::Vector{uType}
     z::Vector{uType}
     fS::Vector{rateType}
+    zemb::uType
     fsalfirst::rateType
     k::rateType
     nlsolver::N
@@ -178,8 +179,6 @@ get_fsalfirstlast(cache::MRIGARKImplicitCache, u) = (cache.fsalfirst, cache.k)
 
 _mrigark_impl_γ(tab) = tab.γ0[findfirst(!iszero, tab.γ0)]
 
-const MRIGARKImplicitAlg = Union{MRIGARKIRK21a, MRIGARKESDIRK34a}
-
 function alg_cache(
         alg::MRIGARKImplicitAlg, u, rate_prototype,
         ::Type{uEltypeNoUnits}, ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits},
@@ -203,10 +202,12 @@ function alg_cache(
     kk = [zero(u) for _ in 1:(tab.q)]
     z = [zero(u) for _ in 1:(s + 1)]
     fS = [zero(rate_prototype) for _ in 1:s]
+    zemb = zero(u)
     fsalfirst = zero(rate_prototype)
     k = zero(rate_prototype)
     return MRIGARKImplicitCache(
-        u, uprev, tmp, atmp, v, vtmp, f1eval, kk, z, fS, fsalfirst, k, nlsolver, tab
+        u, uprev, tmp, atmp, v, vtmp, f1eval, kk, z, fS, zemb,
+        fsalfirst, k, nlsolver, tab
     )
 end
 

lib/OrdinaryDiffEqMultirate/src/multirate_perform_step.jl

@@ -547,8 +547,8 @@ end
 
 function perform_step!(integrator, cache::MRIGARKImplicitCache, repeat_step = false)
     (; t, dt, uprev, u, f, p) = integrator
-    (; tmp, atmp, z, fS, nlsolver, tab) = cache
-    (; Δc, W0, W1, γ0, q) = tab
+    (; tmp, atmp, z, fS, zemb, nlsolver, tab) = cache
+    (; Δc, W0, W1, Wemb0, Wemb1, γ0, q) = tab
     alg = unwrap_alg(integrator, true)
     m = alg.m
     s = length(Δc)
@@ -585,7 +585,16 @@ function perform_step!(integrator, cache::MRIGARKImplicitCache, repeat_step = fa
     @.. broadcast = false u = z[s + 1]
 
     return if integrator.opts.adaptive
-        @.. broadcast = false tmp = u - z[s]
+        if isempty(Wemb0)
+            @.. broadcast = false tmp = u - z[s]
+        else
+            w1emb = isempty(Wemb1) ? nothing : Wemb1
+            _mrigark_substage!(
+                zemb, z[s], cache, f, p, t, dt, Δc[s], cprev - Δc[s], m, q,
+                Wemb0, w1emb, fS, s, stats
+            )
+            @.. broadcast = false tmp = u - zemb
+        end
         calculate_residuals!(
             atmp, tmp, uprev, u,
             integrator.opts.abstol, integrator.opts.reltol,
@@ -600,7 +609,7 @@ end
     )
     (; t, dt, uprev, f, p) = integrator
     (; nlsolver, tab) = cache
-    (; Δc, W0, W1, γ0, q) = tab
+    (; Δc, W0, W1, Wemb0, Wemb1, γ0, q) = tab
     alg = unwrap_alg(integrator, true)
     m = alg.m
     s = length(Δc)
@@ -640,7 +649,16 @@ end
     integrator.u = z[s + 1]
 
     if integrator.opts.adaptive
-        utilde = @.. broadcast = false integrator.u - z[s]
+        utilde = if isempty(Wemb0)
+            @.. broadcast = false integrator.u - z[s]
+        else
+            w1emb = isempty(Wemb1) ? nothing : Wemb1
+            zemb = _mrigark_substage(
+                z[s], f, p, t, dt, Δc[s], cprev - Δc[s], m, q,
+                Wemb0, w1emb, fS, s, stats
+            )
+            @.. broadcast = false integrator.u - zemb
+        end
         atmp = calculate_residuals(
             utilde, uprev, integrator.u,
             integrator.opts.abstol, integrator.opts.reltol,

lib/OrdinaryDiffEqMultirate/src/multirate_tableaus.jl

@@ -109,12 +109,90 @@ function MRIGARKESDIRK34aTableau(::Type{T}) where {T}
     return MRIGARKTableau{T}(Δc, W0, W1, T[], T[], γ0, 3)
 end
 
+# Stage 11 carries an all-zero coupling row: the method is stiffly accurate, so the
+# final stage repeats stage 10. It is kept because the embedding row couples f2 at
+# that stage value.
+function MRIGARKESDIRK46aTableau(::Type{T}) where {T}
+    γ = T(1 // 4)
+    Δc = T[1 // 5, 0, 1 // 5, 0, 1 // 5, 0, 1 // 5, 0, 1 // 5, 0, 0]
+    W0 = zeros(T, 11, 11)
+    W0[1, 1] = 1 // 5
+    W0[2, 1] = -1 // 4
+    W0[3, 1] = 1771023115159 // 1929363690800
+    W0[3, 3] = -1385150376999 // 1929363690800
+    W0[4, 1] = 914009 // 345800
+    W0[4, 3] = -1000459 // 345800
+    W0[5, 1] = 18386293581909 // 36657910125200
+    W0[5, 3] = 5506531089 // 80566835440
+    W0[5, 5] = -178423463189 // 482340922700
+    W0[6, 1] = 36036097 // 8299200
+    W0[6, 3] = 4621 // 118560
+    W0[6, 5] = -38434367 // 8299200
+    W0[7, 1] = -247809665162987 // 146631640500800
+    W0[7, 3] = 10604946373579 // 14663164050080
+    W0[7, 5] = 10838126175385 // 5865265620032
+    W0[7, 7] = -24966656214317 // 36657910125200
+    W0[8, 1] = 38519701 // 11618880
+    W0[8, 3] = 10517363 // 9682400
+    W0[8, 5] = -23284701 // 19364800
+    W0[8, 7] = -10018609 // 2904720
+    W0[9, 1] = -52907807977903 // 33838070884800
+    W0[9, 3] = 74846944529257 // 73315820250400
+    W0[9, 5] = 365022522318171 // 146631640500800
+    W0[9, 7] = -20513210406809 // 109973730375600
+    W0[9, 9] = -2918009798 // 1870301537
+    W0[10, 1] = 19 // 100
+    W0[10, 3] = -73 // 300
+    W0[10, 5] = 127 // 300
+    W0[10, 7] = 127 // 300
+    W0[10, 9] = -313 // 300
+    W1 = zeros(T, 11, 11)
+    W1[3, 1] = -1674554930619 // 964681845400
+    W1[3, 3] = 1674554930619 // 964681845400
+    W1[4, 1] = -1007739 // 172900
+    W1[4, 3] = 1007739 // 172900
+    W1[5, 1] = -8450070574289 // 18328955062600
+    W1[5, 3] = -39429409169 // 40283417720
+    W1[5, 5] = 173621393067 // 120585230675
+    W1[6, 1] = -122894383 // 16598400
+    W1[6, 3] = 14501 // 237120
+    W1[6, 5] = 121879313 // 16598400
+    W1[7, 1] = 32410002731287 // 15434909526400
+    W1[7, 3] = -46499276605921 // 29326328100160
+    W1[7, 5] = -34914135774643 // 11730531240064
+    W1[7, 7] = 45128506783177 // 18328955062600
+    W1[8, 1] = -128357303 // 23237760
+    W1[8, 3] = -35433927 // 19364800
+    W1[8, 5] = 71038479 // 38729600
+    W1[8, 7] = 8015933 // 1452360
+    W1[9, 1] = 136721604296777 // 67676141769600
+    W1[9, 3] = -349632444539303 // 146631640500800
+    W1[9, 5] = -1292744859249609 // 293263281001600
+    W1[9, 7] = 8356250416309 // 54986865187800
+    W1[9, 9] = 17282943803 // 3740603074
+    W1[10, 1] = 3 // 25
+    W1[10, 3] = -29 // 300
+    W1[10, 5] = 71 // 300
+    W1[10, 7] = 71 // 300
+    W1[10, 9] = -149 // 300
+    Wemb0 = zeros(T, 11)
+    Wemb0[1] = -1 // 4
+    Wemb0[3] = 5595 // 8804
+    Wemb0[5] = -2445 // 8804
+    Wemb0[7] = -4225 // 8804
+    Wemb0[9] = 2205 // 4402
+    Wemb0[11] = -567 // 4402
+    γ0 = T[0, γ, 0, γ, 0, γ, 0, γ, 0, γ, 0]
+    return MRIGARKTableau{T}(Δc, W0, W1, Wemb0, T[], γ0, 4)
+end
+
 mri_gark_tableau(::MRIGARKERK22a, ::Type{T}) where {T} = MRIGARKERK22aTableau(T)
 mri_gark_tableau(::MRIGARKERK22b, ::Type{T}) where {T} = MRIGARKERK22bTableau(T)
 mri_gark_tableau(::MRIGARKERK33a, ::Type{T}) where {T} = MRIGARKERK33aTableau(T)
 mri_gark_tableau(::MRIGARKERK45a, ::Type{T}) where {T} = MRIGARKERK45aTableau(T)
 mri_gark_tableau(::MRIGARKIRK21a, ::Type{T}) where {T} = MRIGARKIRK21aTableau(T)
 mri_gark_tableau(::MRIGARKESDIRK34a, ::Type{T}) where {T} = MRIGARKESDIRK34aTableau(T)
+mri_gark_tableau(::MRIGARKESDIRK46a, ::Type{T}) where {T} = MRIGARKESDIRK46aTableau(T)
 
 # ── MIS: multirate infinitesimal step ─────────────────────────────────────────
 

lib/OrdinaryDiffEqMultirate/test/mri_gark_tests.jl

@@ -294,3 +294,73 @@ end
         @test sim.𝒪est[:l∞] ≈ 3 atol = 0.3
     end
 end
+
+@testset "MRIGARKESDIRK46a" begin
+    @testset "Construction" begin
+        @test MRIGARKESDIRK46a(m = 10).m == 10
+        @test MRIGARKESDIRK46a(m = 20).m == 20
+    end
+
+    @testset "Scalar out-of-place" begin
+        prob = SplitODEProblem(
+            (u, p, t) -> -0.9 * u, (u, p, t) -> -0.1 * u, 1.0, (0.0, 1.0)
+        )
+        sol = solve(prob, MRIGARKESDIRK46a(m = 10), dt = 0.05, adaptive = false)
+        @test abs(sol.u[end] - exp(-1.0)) < 1.0e-3
+
+        sol_a = solve(prob, MRIGARKESDIRK46a(m = 10), reltol = 1.0e-6, abstol = 1.0e-6)
+        @test abs(sol_a.u[end] - exp(-1.0)) < 1.0e-3
+    end
+
+    @testset "Vector in-place" begin
+        f1!(du, u, p, t) = (du .= -0.9 .* u)
+        f2!(du, u, p, t) = (du .= -0.1 .* u)
+        u0 = [1.0, 2.0, 3.0]
+        prob = SplitODEProblem(f1!, f2!, u0, (0.0, 1.0))
+
+        sol = solve(prob, MRIGARKESDIRK46a(m = 10), dt = 0.05, adaptive = false)
+        @test norm(sol.u[end] - u0 .* exp(-1.0)) < 1.0e-2
+
+        sol_a = solve(prob, MRIGARKESDIRK46a(m = 10), reltol = 1.0e-6, abstol = 1.0e-6)
+        @test norm(sol_a.u[end] - u0 .* exp(-1.0)) < 1.0e-3
+    end
+
+    @testset "Stiff slow component" begin
+        analytic(u0, p, t) = u0 * exp(-51t)
+        prob = SplitODEProblem(
+            SplitFunction(
+                (u, p, t) -> -u, (u, p, t) -> -50u; analytic = analytic
+            ),
+            1.0, (0.0, 1.0)
+        )
+        sol = solve(prob, MRIGARKESDIRK46a(m = 4), dt = 0.02, adaptive = false)
+        @test abs(sol.u[end] - exp(-51.0)) < 1.0e-3
+        @test solve(prob, MRIGARKESDIRK46a(m = 4), reltol = 1.0e-6, abstol = 1.0e-8).retcode ==
+            ReturnCode.Success
+    end
+
+    @testset "Embedded error estimate" begin
+        # The final coupling row is zero, so `u - z[s]` degenerates to 0 here: the
+        # step size can only respond to the tolerance via the tableau embedding.
+        f1!(du, u, p, t) = (du .= -2 .* u)
+        f2!(du, u, p, t) = (du .= -u)
+        prob = SplitODEProblem(f1!, f2!, [1.0], (0.0, 1.0))
+        loose = solve(prob, MRIGARKESDIRK46a(m = 8), reltol = 1.0e-4, abstol = 1.0e-4)
+        tight = solve(prob, MRIGARKESDIRK46a(m = 8), reltol = 1.0e-8, abstol = 1.0e-8)
+        @test tight.stats.naccept > loose.stats.naccept
+        @test abs(tight.u[end][1] - exp(-3.0)) < abs(loose.u[end][1] - exp(-3.0))
+    end
+
+    @testset "Convergence" begin
+        analytic(u0, p, t) = u0 * exp(-3t)
+        prob = SplitODEProblem(
+            SplitFunction(
+                (u, p, t) -> -2u, (u, p, t) -> -u; analytic = analytic
+            ),
+            1.0, (0.0, 1.0)
+        )
+        dts = 1 ./ 2 .^ (8:-1:4)
+        sim = test_convergence(dts, prob, MRIGARKESDIRK46a(m = 8))
+        @test sim.𝒪est[:l∞] ≈ 4 atol = 0.3
+    end
+end