Fold RK46NL into the generic LowStorageRK2N dispatch

RK46NL was the last plain-2N method still carrying a bespoke constant
cache struct, mutable cache struct, and four hand-unrolled
initialize!/perform_step! methods. Its recurrence is the same Williamson
form the shared LowStorageRK2N kernels already run, so the tableau moves
into a RK46NLConstantCache function returning a
LowStorageRK2NConstantCache with unchanged coefficients, RK46NL joins
the LowStorageRK2NAlgorithm union, and the bespoke caches and kernels
are deleted.

The struct gains the williamson_condition keyword (default false) that
the shared alg_cache reads, and isfsal/uses_uprev are declared false to
match the rest of the family. The generic in-place kernel re-evaluates f
at the start of each step instead of using FSAL storage, keeping the
per-step f-call count at six. The new field is fourth, so positional
construction now takes four arguments, the same break the nine sibling
2N methods already took when they gained the keyword.

Fixed-dt end states are bit-identical to the removed implementation for
in-place problems and for out-of-place problems whose state is an array.
Out-of-place problems with a scalar state differ in the last bits: 3 ulp
after 128 steps of u' = 1.01u, 1 ulp after 256 steps of a nonlinear
scalar problem, roughly four thousand times below the method's own
discretization error there. Both kernels expand to the same
muladd(dt, k, A*tmp), but the straight-line code folds the other product
into the fused multiply-add, and only scalar states are affected because
muladd over arrays never fuses.

RK46NL gains the williamson_condition testset the nine siblings already
have. Version goes to 3.4.0 for the new public keyword.

Author: singhharsh1708

lib/OrdinaryDiffEqLowStorageRK/Project.toml

@@ -1,7 +1,7 @@
 name = "OrdinaryDiffEqLowStorageRK"
 uuid = "b0944070-b475-4768-8dec-fb6eb410534d"
 authors = ["ParamThakkar123 <paramthakkar864@gmail.com>"]
-version = "3.3.0"
+version = "3.4.0"
 
 [deps]
 CommonSolve = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2"

lib/OrdinaryDiffEqLowStorageRK/src/alg_utils.jl

@@ -60,6 +60,7 @@ alg_order(alg::SHLDDRK_2N) = 4
 alg_order(alg::SHLDDRK52) = 2
 
 isfsal(alg::ORK256) = false
+isfsal(alg::RK46NL) = false
 isfsal(alg::CarpenterKennedy2N54) = false
 isfsal(alg::DGLDDRK84_F) = false
 isfsal(alg::DGLDDRK73_C) = false
@@ -73,6 +74,7 @@ isfsal(alg::NDBLSRK124) = false
 isfsal(alg::NDBLSRK144) = false
 
 uses_uprev(alg::ORK256, adaptive::Bool) = false
+uses_uprev(alg::RK46NL, adaptive::Bool) = false
 uses_uprev(alg::SHLDDRK64, adaptive::Bool) = false
 uses_uprev(alg::CarpenterKennedy2N54, adaptive::Bool) = false
 uses_uprev(alg::NDBLSRK124, adaptive::Bool) = false

lib/OrdinaryDiffEqLowStorageRK/src/algorithms.jl

@@ -195,12 +195,16 @@ end
 @doc explicit_rk_docstring(
     "6-stage, fourth order low-stage, low-dissipation, low-dispersion scheme.
 Fixed timestep only.", "RK46NL",
-    references = "Julien Berland, Christophe Bogey, Christophe Bailly. Low-Dissipation and Low-Dispersion Fourth-Order Runge-Kutta Algorithm. Computers & Fluids, 35(10), pp 1459-1463, 2006. doi: https://doi.org/10.1016/j.compfluid.2005.04.003"
+    references = "Julien Berland, Christophe Bogey, Christophe Bailly. Low-Dissipation and Low-Dispersion Fourth-Order Runge-Kutta Algorithm. Computers & Fluids, 35(10), pp 1459-1463, 2006. doi: https://doi.org/10.1016/j.compfluid.2005.04.003",
+    extra_keyword_description = """- `williamson_condition`: allows for an optimization that allows fusing broadcast expressions with the function call `f`. However, it only works for `Array` types.
+    """,
+    extra_keyword_default = "williamson_condition = false"
 )
 Base.@kwdef struct RK46NL{StageLimiter, StepLimiter, Thread} <: OrdinaryDiffEqAlgorithm
     stage_limiter!::StageLimiter = trivial_limiter!
     step_limiter!::StepLimiter = trivial_limiter!
     thread::Thread = Serial()
+    williamson_condition::Bool = false
 end
 
 @doc explicit_rk_docstring(

lib/OrdinaryDiffEqLowStorageRK/src/low_storage_rk_caches.jl

@@ -47,81 +47,30 @@ function ORK256ConstantCache(::Type{T}, ::Type{T2}) where {T, T2}
     return LowStorageRK2NConstantCache(A2end, B1, B2end, c2end)
 end
 
-function alg_cache(
-        alg::RK46NL, u, rate_prototype, ::Type{uEltypeNoUnits},
-        ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t,
-        dt, reltol, p, calck,
-        ::Val{false}, verbose
-    ) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
-    return RK46NLConstantCache(constvalue(uBottomEltypeNoUnits), constvalue(tTypeNoUnits))
-end
+function RK46NLConstantCache(::Type{T}, ::Type{T2}) where {T, T2}
+    A2 = convert(T, -0.737101392796)
+    A3 = convert(T, -1.634740794343)
+    A4 = convert(T, -0.74473900378)
+    A5 = convert(T, -1.469897351522)
+    A6 = convert(T, -2.813971388035)
+    A2end = (A2, A3, A4, A5, A6)
 
-struct RK46NLConstantCache{T, T2} <: OrdinaryDiffEqConstantCache
-    α2::T
-    α3::T
-    α4::T
-    α5::T
-    α6::T
-    β1::T
-    β2::T
-    β3::T
-    β4::T
-    β5::T
-    β6::T
-    c2::T2
-    c3::T2
-    c4::T2
-    c5::T2
-    c6::T2
-
-    function RK46NLConstantCache(::Type{T}, ::Type{T2}) where {T, T2}
-        α2 = T(-0.737101392796)
-        α3 = T(-1.634740794343)
-        α4 = T(-0.74473900378)
-        α5 = T(-1.469897351522)
-        α6 = T(-2.813971388035)
-        β1 = T(0.032918605146)
-        β2 = T(0.8232569982)
-        β3 = T(0.3815309489)
-        β4 = T(0.200092213184)
-        β5 = T(1.718581042715)
-        β6 = T(0.27)
-        c2 = T2(0.032918605146)
-        c3 = T2(0.249351723343)
-        c4 = T2(0.466911705055)
-        c5 = T2(0.582030414044)
-        c6 = T2(0.847252983783)
-        return new{T, T2}(α2, α3, α4, α5, α6, β1, β2, β3, β4, β5, β6, c2, c3, c4, c5, c6)
-    end
-end
+    B1 = convert(T, 0.032918605146)
+    B2 = convert(T, 0.8232569982)
+    B3 = convert(T, 0.3815309489)
+    B4 = convert(T, 0.200092213184)
+    B5 = convert(T, 1.718581042715)
+    B6 = convert(T, 0.27)
+    B2end = (B2, B3, B4, B5, B6)
 
-function alg_cache(
-        alg::RK46NL, u, rate_prototype, ::Type{uEltypeNoUnits},
-        ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t,
-        dt, reltol, p, calck,
-        ::Val{true}, verbose
-    ) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
-    tmp = zero(u)
-    k = zero(rate_prototype)
-    fsalfirst = zero(rate_prototype)
-    tab = RK46NLConstantCache(constvalue(uBottomEltypeNoUnits), constvalue(tTypeNoUnits))
-    return RK46NLCache(
-        u, uprev, k, tmp, fsalfirst, tab, alg.stage_limiter!, alg.step_limiter!,
-        alg.thread
-    )
-end
+    c2 = convert(T2, 0.032918605146)
+    c3 = convert(T2, 0.249351723343)
+    c4 = convert(T2, 0.466911705055)
+    c5 = convert(T2, 0.582030414044)
+    c6 = convert(T2, 0.847252983783)
+    c2end = (c2, c3, c4, c5, c6)
 
-@cache struct RK46NLCache{uType, rateType, TabType, StageLimiter, StepLimiter, Thread} <:
-    LowStorageRKMutableCache
-    u::uType
-    uprev::uType
-    k::rateType
-    tmp::uType
-    fsalfirst::rateType
-    tab::TabType
-    stage_limiter!::StageLimiter
-    step_limiter!::StepLimiter
-    thread::Thread
+    return LowStorageRK2NConstantCache(A2end, B1, B2end, c2end)
 end
 
 function CarpenterKennedy2N54ConstantCache(::Type{T}, ::Type{T2}) where {T, T2}
@@ -4357,11 +4306,12 @@ end
 
 # --- LowStorageRK2N ---
 const LowStorageRK2NAlgorithm = Union{
-    ORK256, CarpenterKennedy2N54, SHLDDRK64, DGLDDRK73_C, DGLDDRK84_C, DGLDDRK84_F,
+    ORK256, RK46NL, CarpenterKennedy2N54, SHLDDRK64, DGLDDRK73_C, DGLDDRK84_C, DGLDDRK84_F,
     NDBLSRK124, NDBLSRK134, NDBLSRK144,
 }
 
 _lowstorage_2n_tableau(::ORK256, ::Type{T}, ::Type{T2}) where {T, T2} = ORK256ConstantCache(T, T2)
+_lowstorage_2n_tableau(::RK46NL, ::Type{T}, ::Type{T2}) where {T, T2} = RK46NLConstantCache(T, T2)
 _lowstorage_2n_tableau(::CarpenterKennedy2N54, ::Type{T}, ::Type{T2}) where {T, T2} = CarpenterKennedy2N54ConstantCache(T, T2)
 _lowstorage_2n_tableau(::SHLDDRK64, ::Type{T}, ::Type{T2}) where {T, T2} = SHLDDRK64ConstantCache(T, T2)
 _lowstorage_2n_tableau(::DGLDDRK73_C, ::Type{T}, ::Type{T2}) where {T, T2} = DGLDDRK73_CConstantCache(T, T2)

lib/OrdinaryDiffEqLowStorageRK/src/low_storage_rk_perform_step.jl

@@ -903,96 +903,6 @@ end
     OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1)
 end
 
-function initialize!(integrator, cache::RK46NLCache)
-    (; k, fsalfirst) = cache
-
-    integrator.kshortsize = 1
-    resize!(integrator.k, integrator.kshortsize)
-    integrator.k[1] = integrator.fsalfirst
-    integrator.f(integrator.fsalfirst, integrator.uprev, integrator.p, integrator.t) # FSAL for interpolation
-    return OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1)
-end
-
-@muladd function perform_step!(integrator, cache::RK46NLCache, repeat_step = false)
-    (; t, dt, uprev, u, f, p) = integrator
-    (; k, fsalfirst, tmp, thread) = cache
-    stage_limiter! = integrator.opts.stage_limiter!
-    (; α2, α3, α4, α5, α6, β1, β2, β3, β4, β5, β6, c2, c3, c4, c5, c6) = cache.tab
-
-    # u1
-    @.. broadcast = false thread = thread tmp = dt * fsalfirst
-    @.. broadcast = false thread = thread u = uprev + β1 * tmp
-    stage_limiter!(u, integrator, p, t + c2 * dt)
-    # u2
-    f(k, u, p, t + c2 * dt)
-    @.. broadcast = false thread = thread tmp = α2 * tmp + dt * k
-    @.. broadcast = false thread = thread u = u + β2 * tmp
-    stage_limiter!(u, integrator, p, t + c3 * dt)
-    # u3
-    f(k, u, p, t + c3 * dt)
-    @.. broadcast = false thread = thread tmp = α3 * tmp + dt * k
-    @.. broadcast = false thread = thread u = u + β3 * tmp
-    stage_limiter!(u, integrator, p, t + c4 * dt)
-    # u4
-    f(k, u, p, t + c4 * dt)
-    @.. broadcast = false thread = thread tmp = α4 * tmp + dt * k
-    @.. broadcast = false thread = thread u = u + β4 * tmp
-    stage_limiter!(u, integrator, p, t + c5 * dt)
-    # u5 = u
-    f(k, u, p, t + c5 * dt)
-    @.. broadcast = false thread = thread tmp = α5 * tmp + dt * k
-    @.. broadcast = false thread = thread u = u + β5 * tmp
-    stage_limiter!(u, integrator, p, t + c6 * dt)
-
-    f(k, u, p, t + c6 * dt)
-    @.. broadcast = false thread = thread tmp = α6 * tmp + dt * k
-    @.. broadcast = false thread = thread u = u + β6 * tmp
-    stage_limiter!(u, integrator, p, t + dt)
-
-    f(k, u, p, t + dt)
-    OrdinaryDiffEqCore.increment_nf!(integrator.stats, 6)
-end
-
-function initialize!(integrator, cache::RK46NLConstantCache)
-    integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t) # Pre-start fsal
-    OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1)
-    integrator.kshortsize = 1
-    integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
-
-    # Avoid undefined entries if k is an array of arrays
-    integrator.fsallast = zero(integrator.fsalfirst)
-    return integrator.k[1] = integrator.fsalfirst
-end
-
-@muladd function perform_step!(integrator, cache::RK46NLConstantCache, repeat_step = false)
-    (; t, dt, uprev, u, f, p) = integrator
-    (; α2, α3, α4, α5, α6, β1, β2, β3, β4, β5, β6, c2, c3, c4, c5, c6) = cache
-
-    # u1
-    tmp = dt * integrator.fsalfirst
-    u = uprev + β1 * tmp
-    # u2
-    tmp = α2 * tmp + dt * f(u, p, t + c2 * dt)
-    u = u + β2 * tmp
-    # u3
-    tmp = α3 * tmp + dt * f(u, p, t + c3 * dt)
-    u = u + β3 * tmp
-    # u4
-    tmp = α4 * tmp + dt * f(u, p, t + c4 * dt)
-    u = u + β4 * tmp
-    # u5 = u
-    tmp = α5 * tmp + dt * f(u, p, t + c5 * dt)
-    u = u + β5 * tmp
-    # u6
-    tmp = α6 * tmp + dt * f(u, p, t + c6 * dt)
-    u = u + β6 * tmp
-
-    integrator.fsallast = f(u, p, t + dt) # For interpolation, then FSAL'd
-    OrdinaryDiffEqCore.increment_nf!(integrator.stats, 6)
-    integrator.k[1] = integrator.fsalfirst
-    integrator.u = u
-end
-
 function initialize!(integrator, cache::SHLDDRK52ConstantCache)
     integrator.kshortsize = 2
     integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)

lib/OrdinaryDiffEqLowStorageRK/test/ode_low_storage_rk_tests.jl

@@ -143,6 +143,64 @@ end
     @test sol_old.u[end] ≈ sol_new.u[end]
 end
 
+@testset "RK46NL" begin
+    alg = RK46NL(; williamson_condition = true)
+    alg2 = RK46NL(; williamson_condition = false)
+    dts = 1 ./ 2 .^ (7:-1:3)
+    for prob in test_problems_only_time
+        sim = test_convergence(dts, prob, alg)
+        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
+        sim = test_convergence(dts, prob, alg2)
+        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
+    end
+    for prob in test_problems_linear
+        sim = test_convergence(dts, prob, alg)
+        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
+        sim = test_convergence(dts, prob, alg2)
+        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
+    end
+    for prob in test_problems_nonlinear
+        sim = test_convergence(dts, prob, alg)
+        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
+        sim = test_convergence(dts, prob, alg2)
+        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
+    end
+    integ = init(
+        prob_ode_large, alg, dt = 1.0e-2, save_start = false, save_end = false,
+        save_everystep = false
+    )
+    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 3
+    integ = init(
+        prob_ode_large, alg, dt = 1.0e-2, save_start = false, save_end = false,
+        save_everystep = false, alias = ODEAliasSpecifier(alias_u0 = true)
+    )
+    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 2
+    integ = init(
+        prob_ode_large, alg2, dt = 1.0e-2, save_start = false, save_end = false,
+        save_everystep = false
+    )
+    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 4
+    integ = init(
+        prob_ode_large, alg2, dt = 1.0e-2, save_start = false, save_end = false,
+        save_everystep = false, alias = ODEAliasSpecifier(alias_u0 = true)
+    )
+    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 3
+    # test whether aliasing u0 is bad
+    new_prob_ode_nonlinear_inplace = ODEProblem(
+        prob_ode_nonlinear_inplace.f, [1.0],
+        (0.0, 0.5)
+    )
+    sol_old = solve(
+        prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
+        save_start = false
+    )
+    sol_new = solve(
+        new_prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
+        save_start = false, alias = ODEAliasSpecifier(alias_u0 = true)
+    )
+    @test sol_old.u[end] ≈ sol_new.u[end]
+end
+
 @testset "CarpenterKennedy2N54" begin
     alg = CarpenterKennedy2N54(; williamson_condition = true)
     alg2 = CarpenterKennedy2N54(; williamson_condition = false)