Add FW15Stage10 explicit 15-stage Runge-Kutta method of order 10

lib/OrdinaryDiffEqHighOrderRK/src/OrdinaryDiffEqHighOrderRK.jl

@@ -32,6 +32,6 @@ include("interpolants.jl")
 include("high_order_rk_addsteps.jl")
 include("high_order_rk_perform_step.jl")
 
-export TanYam7, DP8, PFRK87, TsitPap8
+export TanYam7, DP8, PFRK87, TsitPap8, FW15Stage10
 
 end

lib/OrdinaryDiffEqHighOrderRK/src/alg_utils.jl

@@ -2,6 +2,7 @@ alg_order(alg::TanYam7) = 7
 alg_order(alg::DP8) = 8
 alg_order(alg::TsitPap8) = 8
 alg_order(alg::PFRK87) = 8
+alg_order(alg::FW15Stage10) = 10
 
 qmin_default(alg::DP8) = 1 // 3
 

lib/OrdinaryDiffEqHighOrderRK/src/algorithms.jl

@@ -76,3 +76,23 @@ end
 function PFRK87(stage_limiter!, step_limiter! = trivial_limiter!; omega = 0.0)
     PFRK87(stage_limiter!, step_limiter!, False(), omega)
 end
+
+@doc explicit_rk_docstring(
+    "15-stage explicit Runge-Kutta method of order 10.",
+    "FW15Stage10",
+    references = """@article{FW2025,
+    title={15-stage Runge-Kutta pairs of order 10(9)},
+    author={F. White and Others},
+    journal={arXiv preprint arXiv:2504.17329},
+    year={2025}
+    }""")
+Base.@kwdef struct FW15Stage10{StageLimiter, StepLimiter, Thread} <:
+                   OrdinaryDiffEqAdaptiveAlgorithm
+    stage_limiter!::StageLimiter = trivial_limiter!
+    step_limiter!::StepLimiter = trivial_limiter!
+    thread::Thread = False()
+end
+# for backwards compatibility
+function FW15Stage10(stage_limiter!, step_limiter! = trivial_limiter!)
+    FW15Stage10(stage_limiter!, step_limiter!, False())
+end

lib/OrdinaryDiffEqHighOrderRK/src/high_order_rk_caches.jl

@@ -266,3 +266,69 @@ function alg_cache(alg::PFRK87, u, rate_prototype, ::Type{uEltypeNoUnits},
         ::Val{false}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
     PFRK87ConstantCache(constvalue(uBottomEltypeNoUnits), constvalue(tTypeNoUnits))
 end
+
+@cache struct FW15Stage10Cache{uType, rateType, uNoUnitsType, TabType, StageLimiter, StepLimiter,
+    Thread} <: HighOrderRKMutableCache
+    u::uType
+    uprev::uType
+    k1::rateType
+    k2::rateType
+    k3::rateType
+    k4::rateType
+    k5::rateType
+    k6::rateType
+    k7::rateType
+    k8::rateType
+    k9::rateType
+    k10::rateType
+    k11::rateType
+    k12::rateType
+    k13::rateType
+    k14::rateType
+    k15::rateType
+    utilde::uType
+    tmp::uType
+    atmp::uNoUnitsType
+    fsallast::rateType
+    tab::TabType
+    stage_limiter!::StageLimiter
+    step_limiter!::StepLimiter
+    thread::Thread
+end
+get_fsalfirstlast(cache::FW15Stage10Cache, u) = (cache.k1, cache.fsallast)
+
+function alg_cache(alg::FW15Stage10, u, rate_prototype, ::Type{uEltypeNoUnits},
+        ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t,
+        dt, reltol, p, calck,
+        ::Val{true}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
+    tab = FW15Stage10ConstantCache(constvalue(uBottomEltypeNoUnits), constvalue(tTypeNoUnits))
+    k1 = zero(rate_prototype)
+    k2 = zero(rate_prototype)
+    k3 = zero(rate_prototype)
+    k4 = zero(rate_prototype)
+    k5 = zero(rate_prototype)
+    k6 = zero(rate_prototype)
+    k7 = zero(rate_prototype)
+    k8 = zero(rate_prototype)
+    k9 = zero(rate_prototype)
+    k10 = zero(rate_prototype)
+    k11 = zero(rate_prototype)
+    k12 = zero(rate_prototype)
+    k13 = zero(rate_prototype)
+    k14 = zero(rate_prototype)
+    k15 = zero(rate_prototype)
+    utilde = zero(u)
+    tmp = zero(u)
+    atmp = similar(u, uEltypeNoUnits)
+    recursivefill!(atmp, false)
+    fsallast = zero(rate_prototype)
+    FW15Stage10Cache(u, uprev, k1, k2, k3, k4, k5, k6, k7, k8, k9, k10, k11, k12, k13, k14, k15,
+        utilde, tmp, atmp, fsallast, tab, alg.stage_limiter!, alg.step_limiter!, alg.thread)
+end
+
+function alg_cache(alg::FW15Stage10, u, rate_prototype, ::Type{uEltypeNoUnits},
+        ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t,
+        dt, reltol, p, calck,
+        ::Val{false}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
+    FW15Stage10ConstantCache(constvalue(uBottomEltypeNoUnits), constvalue(tTypeNoUnits))
+end

lib/OrdinaryDiffEqHighOrderRK/src/high_order_rk_perform_step.jl

@@ -1143,3 +1143,112 @@ end
     OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1)
     return nothing
 end
+
+function initialize!(integrator, cache::FW15Stage10ConstantCache)
+    integrator.kshortsize = 2
+    integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
+    integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t) # Pre-start fsal
+    OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1)
+
+    # Avoid undefined entries if k is an array of arrays
+    integrator.fsallast = zero(integrator.fsalfirst)
+    integrator.k[1] = integrator.fsalfirst
+    integrator.k[2] = integrator.fsallast
+end
+
+function initialize!(integrator, cache::FW15Stage10Cache)
+    integrator.kshortsize = 2
+    resize!(integrator.k, integrator.kshortsize)
+    integrator.k[1] = integrator.fsalfirst
+    integrator.k[2] = integrator.fsallast
+    integrator.f(integrator.fsalfirst, integrator.uprev, integrator.p, integrator.t) # Pre-start fsal
+    OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1)
+end
+
+@muladd function perform_step!(integrator, cache::FW15Stage10ConstantCache, repeat_step = false)
+    @unpack t, dt, uprev, u, f, p = integrator
+    @unpack c1, c2, c3, c4, c5, c6, c7, c8, c9, c10, c11, c12, c13, c14, c15, 
+            a21, a32, a41, a43, a51, a53, a54, 
+            a61, a63, a64, a65, a71, a74, a75, a76, a81, a84, a85, a86, a87, a91, a94, a95, a96, a97, a98, a101, a104, a105, a106, a107, a108, a109, 
+            a111, a114, a115, a116, a117, a118, a119, a1110, a121, a124, a125, a126, a127, a128, a129, a1210, a1211, a131, a134, a135, a136, a137, a138, a139, a1310, a1311, a1312, a141, a144, a145, a146, a147, a148, a149, a1410, a1411, a1412, a1413, a151, a154, a155, a156, a157, a158, a159, a1510, a1511, a1512, a1513, a1514, 
+            b1, b9, b10, b11, b12, b13, b14, b15 = cache
+    k1 = integrator.fsalfirst
+    k2 = f(uprev + dt * (a21 * k1), p, t + c2 * dt)
+    k3 = f(uprev + dt * (a32 * k2), p, t + c3 * dt)
+    k4 = f(uprev + dt * (a41 * k1 + a43 * k3), p, t + c4 * dt)
+    k5 = f(uprev + dt * (a51 * k1 + a53 * k3 + a54 * k4), p, t + c5 * dt)
+    k6 = f(uprev + dt * (a61 * k1 + a63 * k3 + a64 * k4 + a65 * k5), p, t + c6 * dt)
+    k7 = f(uprev + dt * (a71 * k1 + a74 * k4 + a75 * k5 + a76 * k6), p, t + c7 * dt)
+    k8 = f(uprev + dt * (a81 * k1 + a84 * k4 + a85 * k5 + a86 * k6 + a87 * k7), p, t + c8 * dt)
+    k9 = f(uprev + dt * (a91 * k1 + a94 * k4 + a95 * k5 + a96 * k6 + a97 * k7 + a98 * k8), p, t + c9 * dt)
+    k10 = f(uprev + dt * (a101 * k1 + a104 * k4 + a105 * k5 + a106 * k6 + a107 * k7 + a108 * k8 + a109 * k9), p, t + c10 * dt)
+    k11 = f(uprev + dt * (a111 * k1 + a114 * k4 + a115 * k5 + a116 * k6 + a117 * k7 + a118 * k8 + a119 * k9 + a1110 * k10), p, t + c11 * dt)
+    k12 = f(uprev + dt * (a121 * k1 + a124 * k4 + a125 * k5 + a126 * k6 + a127 * k7 + a128 * k8 + a129 * k9 + a1210 * k10 + a1211 * k11), p, t + c12 * dt)
+    k13 = f(uprev + dt * (a131 * k1 + a134 * k4 + a135 * k5 + a136 * k6 + a137 * k7 + a138 * k8 + a139 * k9 + a1310 * k10 + a1311 * k11 + a1312 * k12), p, t + c13 * dt)
+    k14 = f(uprev + dt * (a141 * k1 + a144 * k4 + a145 * k5 + a146 * k6 + a147 * k7 + a148 * k8 + a149 * k9 + a1410 * k10 + a1411 * k11 + a1412 * k12 + a1413 * k13), p, t + c14 * dt)
+    k15 = f(uprev + dt * (a151 * k1 + a154 * k4 + a155 * k5 + a156 * k6 + a157 * k7 + a158 * k8 + a159 * k9 + a1510 * k10 + a1511 * k11 + a1512 * k12 + a1513 * k13 + a1514 * k14), p, t + c15 * dt)
+    u = uprev + dt * (b1 * k1 + b9 * k9 + b10 * k10 + b11 * k11 + b12 * k12 + b13 * k13 + b14 * k14 + b15 * k15)
+    integrator.fsallast = k15
+    integrator.k[1] = integrator.fsalfirst
+    integrator.k[2] = integrator.fsallast
+    integrator.u = u
+end
+
+@muladd function perform_step!(integrator, cache::FW15Stage10Cache, repeat_step = false)
+    @unpack t, dt, uprev, u, f, p = integrator
+    @unpack k1, k2, k3, k4, k5, k6, k7, k8, k9, k10, k11, k12, k13, k14, k15,
+            utilde, tmp, atmp, fsallast, stage_limiter!, step_limiter!, thread = cache
+    @unpack c1, c2, c3, c4, c5, c6, c7, c8, c9, c10, c11, c12, c13, c14, c15, 
+            a21, a32, a41, a43, a51, a53, a54, 
+            a61, a63, a64, a65, a71, a74, a75, a76, a81, a84, a85, a86, a87, a91, a94, a95, a96, a97, a98, a101, a104, a105, a106, a107, a108, a109, 
+            a111, a114, a115, a116, a117, a118, a119, a1110, a121, a124, a125, a126, a127, a128, a129, a1210, a1211, a131, a134, a135, a136, a137, a138, a139, a1310, a1311, a1312, a141, a144, a145, a146, a147, a148, a149, a1410, a1411, a1412, a1413, a151, a154, a155, a156, a157, a158, a159, a1510, a1511, a1512, a1513, a1514, 
+            b1, b9, b10, b11, b12, b13, b14, b15 = cache.tab
+
+    f(k1, uprev, p, t)
+    @.. broadcast=false thread=thread tmp = uprev + dt * (a21 * k1)
+    stage_limiter!(tmp, integrator, p, t + c2 * dt)
+    f(k2, tmp, p, t + c2 * dt)
+    @.. broadcast=false thread=thread tmp = uprev + dt * (a32 * k2)
+    stage_limiter!(tmp, integrator, p, t + c3 * dt)
+    f(k3, tmp, p, t + c3 * dt)
+    @.. broadcast=false thread=thread tmp = uprev + dt * (a41 * k1 + a43 * k3)
+    stage_limiter!(tmp, integrator, p, t + c4 * dt)
+    f(k4, tmp, p, t + c4 * dt)
+    @.. broadcast=false thread=thread tmp = uprev + dt * (a51 * k1 + a53 * k3 + a54 * k4)
+    stage_limiter!(tmp, integrator, p, t + c5 * dt)
+    f(k5, tmp, p, t + c5 * dt)
+    @.. broadcast=false thread=thread tmp = uprev + dt * (a61 * k1 + a63 * k3 + a64 * k4 + a65 * k5)
+    stage_limiter!(tmp, integrator, p, t + c6 * dt)
+    f(k6, tmp, p, t + c6 * dt)
+    @.. broadcast=false thread=thread tmp = uprev + dt * (a71 * k1 + a74 * k4 + a75 * k5 + a76 * k6)
+    stage_limiter!(tmp, integrator, p, t + c7 * dt)
+    f(k7, tmp, p, t + c7 * dt)
+    @.. broadcast=false thread=thread tmp = uprev + dt * (a81 * k1 + a84 * k4 + a85 * k5 + a86 * k6 + a87 * k7)
+    stage_limiter!(tmp, integrator, p, t + c8 * dt)
+    f(k8, tmp, p, t + c8 * dt)
+    @.. broadcast=false thread=thread tmp = uprev + dt * (a91 * k1 + a94 * k4 + a95 * k5 + a96 * k6 + a97 * k7 + a98 * k8)
+    stage_limiter!(tmp, integrator, p, t + c9 * dt)
+    f(k9, tmp, p, t + c9 * dt)
+    @.. broadcast=false thread=thread tmp = uprev + dt * (a101 * k1 + a104 * k4 + a105 * k5 + a106 * k6 + a107 * k7 + a108 * k8 + a109 * k9)
+    stage_limiter!(tmp, integrator, p, t + c10 * dt)
+    f(k10, tmp, p, t + c10 * dt)
+    @.. broadcast=false thread=thread tmp = uprev + dt * (a111 * k1 + a114 * k4 + a115 * k5 + a116 * k6 + a117 * k7 + a118 * k8 + a119 * k9 + a1110 * k10)
+    stage_limiter!(tmp, integrator, p, t + c11 * dt)
+    f(k11, tmp, p, t + c11 * dt)
+    @.. broadcast=false thread=thread tmp = uprev + dt * (a121 * k1 + a124 * k4 + a125 * k5 + a126 * k6 + a127 * k7 + a128 * k8 + a129 * k9 + a1210 * k10 + a1211 * k11)
+    stage_limiter!(tmp, integrator, p, t + c12 * dt)
+    f(k12, tmp, p, t + c12 * dt)
+    @.. broadcast=false thread=thread tmp = uprev + dt * (a131 * k1 + a134 * k4 + a135 * k5 + a136 * k6 + a137 * k7 + a138 * k8 + a139 * k9 + a1310 * k10 + a1311 * k11 + a1312 * k12)
+    stage_limiter!(tmp, integrator, p, t + c13 * dt)
+    f(k13, tmp, p, t + c13 * dt)
+    @.. broadcast=false thread=thread tmp = uprev + dt * (a141 * k1 + a144 * k4 + a145 * k5 + a146 * k6 + a147 * k7 + a148 * k8 + a149 * k9 + a1410 * k10 + a1411 * k11 + a1412 * k12 + a1413 * k13)
+    stage_limiter!(tmp, integrator, p, t + c14 * dt)
+    f(k14, tmp, p, t + c14 * dt)
+    @.. broadcast=false thread=thread tmp = uprev + dt * (a151 * k1 + a154 * k4 + a155 * k5 + a156 * k6 + a157 * k7 + a158 * k8 + a159 * k9 + a1510 * k10 + a1511 * k11 + a1512 * k12 + a1513 * k13 + a1514 * k14)
+    stage_limiter!(tmp, integrator, p, t + c15 * dt)
+    f(k15, tmp, p, t + c15 * dt)
+    @.. broadcast=false thread=thread u = uprev + dt * (b1 * k1 + b9 * k9 + b10 * k10 + b11 * k11 + b12 * k12 + b13 * k13 + b14 * k14 + b15 * k15)
+    stage_limiter!(u, integrator, p, t + dt)
+    step_limiter!(u, integrator, p, t + dt)
+    f(fsallast, u, p, t + dt)
+end