Shifted methods to the correct subpackage

lib/OrdinaryDiffEqLowStorageRK/src/OrdinaryDiffEqLowStorageRK.jl

@@ -103,6 +103,6 @@ export ORK256, CarpenterKennedy2N54, SHLDDRK64, HSLDDRK64, DGLDDRK73_C, DGLDDRK8
        ParsaniKetchesonDeconinck3S53, ParsaniKetchesonDeconinck3S173,
        ParsaniKetchesonDeconinck3S94, ParsaniKetchesonDeconinck3S184,
        ParsaniKetchesonDeconinck3S105, ParsaniKetchesonDeconinck3S205,
-       RDPK3Sp35, RDPK3SpFSAL35, RDPK3Sp49, RDPK3SpFSAL49, RDPK3Sp510, RDPK3SpFSAL510,
-       KYK2014DGSSPRK_3S2, RK46NL
+       RDPK3Sp35, RDPK3SpFSAL35, RDPK3Sp49, RDPK3SpFSAL49, RDPK3Sp510, RDPK3SpFSAL510, 
+       RK46NL, SHLDDRK_2N, SHLDDRK52
 end

lib/OrdinaryDiffEqLowStorageRK/src/alg_utils.jl

@@ -1,4 +1,3 @@
-alg_order(alg::KYK2014DGSSPRK_3S2) = 2
 alg_order(alg::ORK256) = 2
 alg_order(alg::CarpenterKennedy2N54) = 4
 alg_order(alg::NDBLSRK124) = 4
@@ -41,6 +40,8 @@ alg_order(alg::CKLLSRK54_3N_4R) = 4
 alg_order(alg::CKLLSRK75_4M_5R) = 5
 alg_order(alg::CKLLSRK54_3M_4R) = 4
 alg_order(alg::CKLLSRK65_4M_4R) = 5
+alg_order(alg::SHLDDRK_2N) = 4
+alg_order(alg::SHLDDRK52) = 2
 
 isfsal(alg::ORK256) = false
 isfsal(alg::CarpenterKennedy2N54) = false
@@ -83,22 +84,6 @@ uses_uprev(alg::CKLLSRK65_4M_4R, adaptive::Bool) = adaptive
 uses_uprev(alg::CKLLSRK85_4FM_4R, adaptive::Bool) = adaptive
 uses_uprev(alg::CKLLSRK75_4M_5R, adaptive::Bool) = adaptive
 
-"""
-    ssp_coefficient(alg)
-
-Return the SSP coefficient of the ODE algorithm `alg`. If one time step of size
-`dt` with `alg` can be written as a convex combination of explicit Euler steps
-with step sizes `cᵢ * dt`, the SSP coefficient is the minimal value of `1/cᵢ`.
-
-# Examples
-
-```julia-repl
-julia> ssp_coefficient(SSPRK104())
-6
-```
-"""
-ssp_coefficient(alg::KYK2014DGSSPRK_3S2) = 0.8417
-
 function default_controller(alg::RDPK3Sp35, cache, qoldinit, args...)
     QT = typeof(qoldinit)
     return PIDController(map(Base.Fix1(convert, QT), (0.64, -0.31, 0.04))...)

lib/OrdinaryDiffEqLowStorageRK/src/algorithms.jl

@@ -51,6 +51,32 @@ function DGLDDRK73_C(stage_limiter!, step_limiter! = trivial_limiter!;
         williamson_condition)
 end
 
+@doc explicit_rk_docstring("TBD", "SHLDDRK_2N")
+Base.@kwdef struct SHLDDRK_2N{StageLimiter, StepLimiter, Thread} <: OrdinaryDiffEqAlgorithm
+    stage_limiter!::StageLimiter = trivial_limiter!
+    step_limiter!::StepLimiter = trivial_limiter!
+    thread::Thread = False()
+end
+# for backwards compatibility
+function SHLDDRK_2N(stage_limiter!, step_limiter! = trivial_limiter!)
+    SHLDDRK_2N(stage_limiter!,
+        step_limiter!,
+        False())
+end
+
+@doc explicit_rk_docstring("TBD", "SHLDDRK52")
+Base.@kwdef struct SHLDDRK52{StageLimiter, StepLimiter, Thread} <: OrdinaryDiffEqAlgorithm
+    stage_limiter!::StageLimiter = trivial_limiter!
+    step_limiter!::StepLimiter = trivial_limiter!
+    thread::Thread = False()
+end
+# for backwards compatibility
+function SHLDDRK52(stage_limiter!, step_limiter! = trivial_limiter!)
+    SHLDDRK52(stage_limiter!,
+        step_limiter!,
+        False())
+end
+
 @doc explicit_rk_docstring(
     "A fourth-order, five-stage explicit low-storage method of Carpenter and Kennedy
 (free 3rd order Hermite interpolant). Fixed timestep only. Designed for
@@ -888,19 +914,4 @@ function NDBLSRK134(stage_limiter!, step_limiter! = trivial_limiter!;
         williamson_condition)
 end
 
-#SSP Optimized Runge-Kutta Methods
-
-@doc explicit_rk_docstring("TBD",
-    "KYK2014DGSSPRK_3S2")
-Base.@kwdef struct KYK2014DGSSPRK_3S2{StageLimiter, StepLimiter, Thread} <:
-                   OrdinaryDiffEqAlgorithm
-    stage_limiter!::StageLimiter = trivial_limiter!
-    step_limiter!::StepLimiter = trivial_limiter!
-    thread::Thread = False()
-end
-# for backwards compatibility
-function KYK2014DGSSPRK_3S2(stage_limiter!, step_limiter! = trivial_limiter!)
-    KYK2014DGSSPRK_3S2(stage_limiter!,
-        step_limiter!,
-        False())
-end
+#SSP Optimized Runge-Kutta Methods

lib/OrdinaryDiffEqLowStorageRK/src/low_storage_rk_caches.jl

@@ -119,83 +119,6 @@ struct RK46NLConstantCache{T, T2} <: OrdinaryDiffEqConstantCache
     end
 end
 
-@cache struct KYK2014DGSSPRK_3S2_Cache{uType, rateType, TabType, StageLimiter, StepLimiter,
-    Thread} <:
-              LowStorageRKMutableCache
-    u::uType
-    uprev::uType
-    k::rateType
-    fsalfirst::rateType
-    tab::TabType
-    #temporary values for Shu-Osher
-    u_1::uType
-    u_2::uType
-    kk_1::rateType
-    kk_2::rateType
-    stage_limiter!::StageLimiter
-    step_limiter!::StepLimiter
-    thread::Thread
-end
-
-struct KYK2014DGSSPRK_3S2_ConstantCache{T, T2} <: OrdinaryDiffEqConstantCache
-    #These are not α and β for RK but for Shu-Osher
-    #see top of page 317 in
-    #Optimal Strong-Stability-Preserving Runge–Kutta Time Discretizations for
-    #Discontinuous Garlekin Methods, Kubatko, Yaeger, Ketcheson 2014
-    α_10::T
-    α_20::T
-    α_21::T
-    α_30::T
-    α_32::T
-    β_10::T
-    β_21::T
-    β_30::T
-    β_32::T
-    #Shu-Osher is normally stated for autonomous systems, the times
-    #are calculated by hand for this scheme
-    c_1::T
-    c_2::T
-
-    function KYK2014DGSSPRK_3S2_ConstantCache(T, T2)
-        α_10 = T(1.0)
-        α_20 = T(0.087353119859156)
-        α_21 = T(0.912646880140844)
-        α_30 = T(0.344956917166841)
-        α_32 = T(0.655043082833159)
-        β_10 = T(0.528005024856522)
-        β_21 = T(0.481882138633993)
-        β_30 = T(0.022826837460491)
-        β_32 = T(0.345866039233415)
-        c_1 = β_10
-        c_2 = α_21 * β_10 + β_21 # ==0.96376427726
-        new{T, T2}(α_10, α_20, α_21, α_30, α_32, β_10, β_21, β_30, β_32, c_1, c_2)
-    end
-end
-
-function alg_cache(alg::KYK2014DGSSPRK_3S2, u, rate_prototype, ::Type{uEltypeNoUnits},
-        ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t,
-        dt, reltol, p, calck,
-        ::Val{true}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
-    u_1 = zero(u)
-    u_2 = zero(u)
-    kk_1 = zero(rate_prototype)
-    kk_2 = zero(rate_prototype)
-    k = zero(rate_prototype)
-    fsalfirst = zero(rate_prototype)
-    tab = KYK2014DGSSPRK_3S2_ConstantCache(constvalue(uBottomEltypeNoUnits),
-        constvalue(tTypeNoUnits))
-    KYK2014DGSSPRK_3S2_Cache(u, uprev, k, fsalfirst, tab, u_1, u_2, kk_1, kk_2,
-        alg.stage_limiter!, alg.step_limiter!, alg.thread)
-end
-
-function alg_cache(alg::KYK2014DGSSPRK_3S2, u, rate_prototype, ::Type{uEltypeNoUnits},
-        ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t,
-        dt, reltol, p, calck,
-        ::Val{false}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
-    KYK2014DGSSPRK_3S2_ConstantCache(constvalue(uBottomEltypeNoUnits),
-        constvalue(tTypeNoUnits))
-end
-
 function alg_cache(alg::RK46NL, u, rate_prototype, ::Type{uEltypeNoUnits},
         ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t,
         dt, reltol, p, calck,
@@ -244,6 +167,164 @@ function CarpenterKennedy2N54ConstantCache(T, T2)
     LowStorageRK2NConstantCache{4, T, T2}(A2end, B1, B2end, c2end)
 end
 
+@cache mutable struct SHLDDRK_2NCache{uType, rateType, TabType, StageLimiter, StepLimiter,
+    Thread} <: SSPRKMutableCache
+    u::uType
+    uprev::uType
+    k::rateType
+    tmp::uType
+    fsalfirst::rateType
+    tab::TabType
+    step::Int
+    stage_limiter!::StageLimiter
+    step_limiter!::StepLimiter
+    thread::Thread
+end
+
+mutable struct SHLDDRK_2NConstantCache{T1, T2} <: SSPRKConstantCache
+    α21::T1
+    α31::T1
+    α41::T1
+    α51::T1
+    β11::T1
+    β21::T1
+    β31::T1
+    β41::T1
+    β51::T1
+    c21::T2
+    c31::T2
+    c41::T2
+    c51::T2
+
+    α22::T1
+    α32::T1
+    α42::T1
+    α52::T1
+    α62::T1
+    β12::T1
+    β22::T1
+    β32::T1
+    β42::T1
+    β52::T1
+    β62::T1
+    c22::T2
+    c32::T2
+    c42::T2
+    c52::T2
+    c62::T2
+
+    step::Int
+end
+
+function SHLDDRK_2NConstantCache(T1, T2)
+    α21 = T1(-0.6051226)
+    α31 = T1(-2.0437564)
+    α41 = T1(-0.7406999)
+    α51 = T1(-4.4231765)
+    β11 = T1(0.2687454)
+    β21 = T1(0.8014706)
+    β31 = T1(0.5051570)
+    β41 = T1(0.5623568)
+    β51 = T1(0.0590065)
+    c21 = T2(0.2687454)
+    c31 = T2(0.5852280)
+    c41 = T2(0.6827066)
+    c51 = T2(1.1646854)
+
+    α22 = T1(-0.4412737)
+    α32 = T1(-1.0739820)
+    α42 = T1(-1.7063570)
+    α52 = T1(-2.7979293)
+    α62 = T1(-4.0913537)
+    β12 = T1(0.1158488)
+    β22 = T1(0.3728769)
+    β32 = T1(0.7379536)
+    β42 = T1(0.5798110)
+    β52 = T1(1.0312849)
+    β62 = T1(0.15)
+    c22 = T2(0.1158485)
+    c32 = T2(0.3241850)
+    c42 = T2(0.6193208)
+    c52 = T2(0.8034472)
+    c62 = T2(0.9184166)
+    SHLDDRK_2NConstantCache(
+        α21, α31, α41, α51, β11, β21, β31, β41, β51, c21, c31, c41, c51,
+        α22, α32, α42, α52, α62, β12, β22, β32, β42, β52, β62, c22, c32,
+        c42, c52, c62, 1)
+end
+
+function alg_cache(alg::SHLDDRK_2N, u, rate_prototype, ::Type{uEltypeNoUnits},
+        ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t,
+        dt, reltol, p, calck,
+        ::Val{false}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
+    SHLDDRK_2NConstantCache(constvalue(uBottomEltypeNoUnits), constvalue(tTypeNoUnits))
+end
+
+function alg_cache(alg::SHLDDRK_2N, u, rate_prototype, ::Type{uEltypeNoUnits},
+        ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t,
+        dt, reltol, p, calck,
+        ::Val{true}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
+    tmp = zero(u)
+    k = zero(rate_prototype)
+    fsalfirst = zero(rate_prototype)
+    tab = SHLDDRK_2NConstantCache(constvalue(uBottomEltypeNoUnits),
+        constvalue(tTypeNoUnits))
+    SHLDDRK_2NCache(u, uprev, k, tmp, fsalfirst, tab, 1, alg.stage_limiter!,
+        alg.step_limiter!, alg.thread)
+end
+
+struct SHLDDRK52ConstantCache{T1, T2} <: SSPRKConstantCache
+    α2::T1
+    α3::T1
+    α4::T1
+    α5::T1
+    β1::T1
+    β2::T1
+    β3::T1
+    β4::T1
+    β5::T1
+    c2::T2
+    c3::T2
+    c4::T2
+    c5::T2
+end
+
+function SHLDDRK52ConstantCache(T1, T2)
+    α2 = T1(-0.6913065)
+    α3 = T1(-2.655155)
+    α4 = T1(-0.8147688)
+    α5 = T1(-0.6686587)
+    β1 = T1(0.1)
+    β2 = T1(0.75)
+    β3 = T1(0.7)
+    β4 = T1(0.479313)
+    β5 = T1(0.310392)
+    c2 = T2(0.1)
+    c3 = T2(0.3315201)
+    c4 = T2(0.4577796)
+    c5 = T2(0.8666528)
+    SHLDDRK52ConstantCache(α2, α3, α4, α5, β1, β2, β3, β4, β5, c2, c3, c4, c5)
+end
+
+function alg_cache(alg::SHLDDRK52, u, rate_prototype, ::Type{uEltypeNoUnits},
+        ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t,
+        dt, reltol, p, calck,
+        ::Val{false}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
+    SHLDDRK52ConstantCache(constvalue(uBottomEltypeNoUnits), constvalue(tTypeNoUnits))
+end
+
+function alg_cache(alg::SHLDDRK52, u, rate_prototype, ::Type{uEltypeNoUnits},
+        ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t,
+        dt, reltol, p, calck,
+        ::Val{true}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
+    tmp = zero(u)
+    k = zero(rate_prototype)
+    fsalfirst = zero(rate_prototype)
+    tab = SHLDDRK52ConstantCache(constvalue(uBottomEltypeNoUnits), constvalue(tTypeNoUnits))
+    SHLDDRK52Cache(u, uprev, k, tmp, fsalfirst, tab, alg.stage_limiter!, alg.step_limiter!,
+        alg.thread)
+end
+
 function alg_cache(alg::CarpenterKennedy2N54, u, rate_prototype, ::Type{uEltypeNoUnits},
         ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t,
         dt, reltol, p, calck,