fixes

Author: Your Name

lib/OrdinaryDiffEqLowStorageRK/src/alg_utils.jl

@@ -43,6 +43,7 @@ alg_order(alg::CKLLSRK65_4M_4R) = 5
 alg_order(alg::SHLDDRK_2N) = 4
 alg_order(alg::SHLDDRK52) = 2
 
+isfsal(alg::RK46NL) = false
 isfsal(alg::ORK256) = false
 isfsal(alg::CarpenterKennedy2N54) = false
 isfsal(alg::DGLDDRK84_F) = false

lib/OrdinaryDiffEqLowStorageRK/src/low_storage_rk_caches.jl

@@ -47,6 +47,23 @@ function ORK256ConstantCache(::Type{T}, ::Type{T2}) where {T, T2}
     return LowStorageRK2NConstantCache(A2end, B1, B2end, c2end)
 end
 
+function RK46NLConstantCache(::Type{T}, ::Type{T2}) where {T, T2}
+    A2end = (
+        T(-0.737101392796), T(-1.634740794343), T(-0.74473900378),
+        T(-1.469897351522), T(-2.813971388035),
+    )
+    B1 = T(0.032918605146)
+    B2end = (
+        T(0.8232569982), T(0.3815309489), T(0.200092213184),
+        T(1.718581042715), T(0.27),
+    )
+    c2end = (
+        T2(0.032918605146), T2(0.249351723343), T2(0.466911705055),
+        T2(0.582030414044), T2(0.847252983783),
+    )
+    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,
@@ -56,45 +73,6 @@ function alg_cache(
     return RK46NLConstantCache(constvalue(uBottomEltypeNoUnits), constvalue(tTypeNoUnits))
 end
 
-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
-
 function alg_cache(
         alg::RK46NL, u, rate_prototype, ::Type{uEltypeNoUnits},
         ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t,

lib/OrdinaryDiffEqLowStorageRK/src/low_storage_rk_perform_step.jl

@@ -922,81 +922,21 @@ end
 @muladd function perform_step!(integrator, cache::RK46NLCache, repeat_step = false)
     (; t, dt, uprev, u, f, p) = integrator
     (; k, fsalfirst, tmp, stage_limiter!, step_limiter!, thread) = cache
-    (; α2, α3, α4, α5, α6, β1, β2, β3, β4, β5, β6, c2, c3, c4, c5, c6) = cache.tab
+    (; A2end, B1, B2end, c2end) = 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
+    @.. broadcast = false thread = thread u = uprev + B1 * tmp
+    for i in eachindex(A2end)
+        stage_limiter!(u, integrator, p, t + c2end[i] * dt)
+        f(k, u, p, t + c2end[i] * dt)
+        OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1)
+        @.. broadcast = false thread = thread tmp = A2end[i] * tmp + dt * k
+        @.. broadcast = false thread = thread u = u + B2end[i] * tmp
+    end
     stage_limiter!(u, integrator, p, t + dt)
     step_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)