Added High Order and Low Order RK methods

lib/OrdinaryDiffEqHighOrderRK/Project.toml

@@ -0,0 +1,22 @@
+name = "OrdinaryDiffEqHighOrderRK"
+uuid = "d28bc4f8-55e1-4f49-af69-84c1a99f0f58"
+authors = ["ParamThakkar123 <[email protected]>"]
+version = "1.0.0"
+
+[deps]
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221"
+RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
+FastBroadcast = "7034ab61-46d4-4ed7-9d0f-46aef9175898"
+
+[compat]
+julia = "1.10"
+
+[extras]
+DiffEqDevTools = "f3b72e0c-5b89-59e1-b016-84e28bfd966d"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[targets]
+test = ["DiffEqDevTools", "Random", "SafeTestsets", "Test"]

lib/OrdinaryDiffEqHighOrderRK/src/OrdinaryDiffEqHighOrderRK.jl

@@ -0,0 +1,22 @@
+module OrdinaryDiffEqHighOrderRK
+
+import OrdinaryDiffEq: alg_order, qmax_default, qmin_default, beta2_default, beta1_default,
+                       explicit_rk_docstring, OrdinaryDiffEqAdaptiveAlgorithm, trivial_limiter!,
+                       _ode_addsteps!, @unpack, @cache, OrdinaryDiffEqMutableCache, constvalue,
+                       alg_cache, uses_uprev, initialize!, perform_step!, OrdinaryDiffEqConstantCache,
+                       calculate_residuals!, calculate_residuals
+import Static: False
+import MuladdMacro: @muladd
+import FastBroadcast: @..
+import RecursiveArrayTools: recursivefill!
+
+include("algorithms.jl")
+include("alg_utils.jl")
+include("high_order_rk_caches.jl")
+include("high_order_rk_tableaus.jl")
+include("high_order_rk_addsteps.jl")
+include("high_order_rk_perform_step.jl")
+
+export TanYam7, DP8, PFRK87, TsitPap8
+
+end

lib/OrdinaryDiffEqHighOrderRK/src/alg_utils.jl

@@ -0,0 +1,12 @@
+alg_order(alg::TanYam7) = 7
+alg_order(alg::DP8) = 8
+alg_order(alg::TsitPap8) = 8
+alg_order(alg::PFRK87) = 8
+
+qmin_default(alg::DP8) = 1 // 3
+
+qmax_default(alg::DP8) = 6
+
+beta2_default(alg::DP8) = 0 // 1
+
+beta1_default(alg::DP8, beta2) = typeof(beta2)(1 // alg_order(alg)) - beta2 / 5

lib/OrdinaryDiffEqHighOrderRK/src/algorithms.jl

@@ -0,0 +1,61 @@
+@doc explicit_rk_docstring(
+    "Tanaka-Yamashita 7 Runge-Kutta method. (7th order interpolant).",
+    "TanYam7",
+    references = "Tanaka M., Muramatsu S., Yamashita S., (1992), On the Optimization of Some Nine-Stage
+    Seventh-order Runge-Kutta Method, Information Processing Society of Japan,
+    33 (12), pp. 1512-1526.")
+Base.@kwdef struct TanYam7{StageLimiter, StepLimiter, Thread} <:
+                   OrdinaryDiffEqAdaptiveAlgorithm
+    stage_limiter!::StageLimiter = trivial_limiter!
+    step_limiter!::StepLimiter = trivial_limiter!
+    thread::Thread = False()
+end
+# for backwards compatibility
+function TanYam7(stage_limiter!, step_limiter! = trivial_limiter!)
+    TanYam7(stage_limiter!, step_limiter!, False())
+end
+
+@doc explicit_rk_docstring(
+    "Hairer's 8/5/3 adaption of the Dormand-Prince Runge-Kutta method. (7th order interpolant).",
+    "DP8",
+    references = "E. Hairer, S.P. Norsett, G. Wanner, (1993) Solving Ordinary Differential Equations I.
+    Nonstiff Problems. 2nd Edition. Springer Series in Computational Mathematics,
+    Springer-Verlag.")
+Base.@kwdef struct DP8{StageLimiter, StepLimiter, Thread} <: OrdinaryDiffEqAdaptiveAlgorithm
+    stage_limiter!::StageLimiter = trivial_limiter!
+    step_limiter!::StepLimiter = trivial_limiter!
+    thread::Thread = False()
+end
+# for backwards compatibility
+function DP8(stage_limiter!, step_limiter! = trivial_limiter!)
+    DP8(stage_limiter!, step_limiter!, False())
+end
+
+@doc explicit_rk_docstring("Phase-fitted Runge-Kutta of 8th order.", "PFRK87",
+    extra_keyword_description = """- `omega`: a periodicity phase estimate,
+                                   when accurate this method results in zero numerical dissipation.
+                    """,
+    extra_keyword_default = "omega = 0.0")
+Base.@kwdef struct PFRK87{StageLimiter, StepLimiter, Thread, T} <:
+                   OrdinaryDiffEqAdaptiveAlgorithm
+    stage_limiter!::StageLimiter = trivial_limiter!
+    step_limiter!::StepLimiter = trivial_limiter!
+    thread::Thread = False()
+    omega::T = 0.0
+end
+# for backwards compatibility
+function PFRK87(stage_limiter!, step_limiter! = trivial_limiter!; omega = 0.0)
+    PFRK87(stage_limiter!, step_limiter!, False(), omega)
+end
+
+@doc explicit_rk_docstring("Tsitouras-Papakostas 8/7 Runge-Kutta method.", "TsitPap8")
+Base.@kwdef struct TsitPap8{StageLimiter, StepLimiter, Thread} <:
+                   OrdinaryDiffEqAdaptiveAlgorithm
+    stage_limiter!::StageLimiter = trivial_limiter!
+    step_limiter!::StepLimiter = trivial_limiter!
+    thread::Thread = False()
+end
+# for backwards compatibility
+function TsitPap8(stage_limiter!, step_limiter! = trivial_limiter!)
+    TsitPap8(stage_limiter!, step_limiter!, False())
+end

src/dense/high_order_rk_addsteps.jl → lib/OrdinaryDiffEqHighOrderRK/src/high_order_rk_addsteps.jl

@@ -256,4 +256,4 @@ end
     end
   end
 end
-=#
+=#

src/caches/high_order_rk_caches.jl → lib/OrdinaryDiffEqHighOrderRK/src/high_order_rk_caches.jl

@@ -262,4 +262,4 @@ function alg_cache(alg::PFRK87, u, rate_prototype, ::Type{uEltypeNoUnits},
         dt, reltol, p, calck,
         ::Val{false}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
     PFRK87ConstantCache(constvalue(uBottomEltypeNoUnits), constvalue(tTypeNoUnits))
-end
+end

src/perform_step/high_order_rk_perform_step.jl → lib/OrdinaryDiffEqHighOrderRK/src/high_order_rk_perform_step.jl

@@ -1150,4 +1150,4 @@ end
     f(k, u, p, t + dt)
     integrator.stats.nf += 1
     return nothing
-end
+end

src/tableaus/high_order_rk_tableaus.jl → lib/OrdinaryDiffEqHighOrderRK/src/high_order_rk_tableaus.jl

@@ -1403,4 +1403,4 @@ function PFRK87ConstantCache(T1::Type, T2::Type)
         α1211, α1311, β1, β6, β7, β8, β9, β10, β11, β12, β13, β1tilde,
         β6tilde, β7tilde, β8tilde, β9tilde, β10tilde, β11tilde, β12tilde,
         c2, c3, c4, c5, c6, c7, c8, c9, c10, c11, c12, c13)
-end
+end

lib/OrdinaryDiffEqLowOrderRK/Project.toml

@@ -0,0 +1,22 @@
+name = "OrdinaryDiffEqLowOrderRK"
+uuid = "1344f307-1e59-4825-a18e-ace9aa3fa4c6"
+authors = ["ParamThakkar123 <[email protected]>"]
+version = "1.0.0"
+
+[deps]
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221"
+RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
+FastBroadcast = "7034ab61-46d4-4ed7-9d0f-46aef9175898"
+
+[compat]
+julia = "1.10"
+
+[extras]
+DiffEqDevTools = "f3b72e0c-5b89-59e1-b016-84e28bfd966d"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[targets]
+test = ["DiffEqDevTools", "Random", "SafeTestsets", "Test"]

lib/OrdinaryDiffEqLowOrderRK/src/OrdinaryDiffEqLowOrderRK.jl

@@ -0,0 +1,30 @@
+module OrdinaryDiffEqLowOrderRK
+
+import OrdinaryDiffEq: alg_order, isfsal, beta2_default, beta1_default, alg_stability_size,
+                       ssp_coefficient, OrdinaryDiffEqAlgorithm, OrdinaryDiffEqExponentialAlgorithm,
+                       explicit_rk_docstring, trivial_limiter!, OrdinaryDiffEqAdaptiveAlgorithm,
+                       unwrap_alg, @unpack, initialize!, perform_step!, calculate_residuals,
+                       calculate_residuals!, _ode_addsteps!, @OnDemandTableauExtract, constvalue,
+                       OrdinaryDiffEqMutableCache, uses_uprev, OrdinaryDiffEqConstantCache
+import DiffEqBase: @tight_loop_macros
+import MuladdMacro: @muladd
+import FastBroadcast: @..
+import RecursiveArrayTools: recursivefill!
+
+include("algorithms.jl")
+include("alg_utils.jl")
+include("low_order_rk_caches.jl")
+include("low_order_rk_tableaus.jl")
+include("interp_func.jl")
+include("low_order_rk_perform_step.jl")
+include("low_order_rk_addsteps.jl")
+include("split_perform_step.jl")
+include("fixed_timestep_perform_step.jl")
+
+export Euler, SplitEuler, Heun, Ralston, Midpoint, RK4,
+       BS3, OwrenZen3, OwrenZen4, OwrenZen5, BS5, Tsit5,
+       DP5, Anas5, RKO65, FRK65, RKM, MSRK5, MSRK6,
+       PSRK4p7q6, PSRK3p5q4, PSRK3p6q5, Stepanov5, SIR54,
+       Alshina2, Alshina3, Alshina6
+
+end

lib/OrdinaryDiffEqLowOrderRK/src/alg_utils.jl

@@ -0,0 +1,49 @@
+alg_order(alg::Euler) = 1
+alg_order(alg::SplitEuler) = 1
+alg_order(alg::Heun) = 2
+alg_order(alg::Ralston) = 2
+alg_order(alg::Midpoint) = 2
+alg_order(alg::BS3) = 3
+alg_order(alg::OwrenZen3) = 3
+alg_order(alg::BS5) = 5
+alg_order(alg::OwrenZen4) = 4
+alg_order(alg::OwrenZen5) = 5
+alg_order(alg::Tsit5) = 5
+alg_order(alg::DP5) = 5
+alg_order(alg::Anas5) = 5
+alg_order(alg::RKO65) = 5
+alg_order(alg::FRK65) = 6
+alg_order(alg::RK4) = 4
+alg_order(alg::RKM) = 4
+alg_order(alg::MSRK5) = 5
+alg_order(alg::MSRK6) = 6
+alg_order(alg::PSRK4p7q6) = 4
+alg_order(alg::PSRK3p6q5) = 3
+alg_order(alg::PSRK3p5q4) = 3
+alg_order(alg::Stepanov5) = 5
+alg_order(alg::SIR54) = 5
+alg_order(alg::Alshina2) = 2
+alg_order(alg::Alshina3) = 3
+alg_order(alg::Alshina6) = 6
+
+isfsal(alg::FRK65) = true
+isfsal(alg::RKO65) = false
+isfsal(alg::PSRK3p5q4) = false
+isfsal(alg::PSRK3p6q5) = false
+isfsal(alg::PSRK4p7q6) = false
+
+beta2_default(alg::DP5) = 4 // 100
+
+beta1_default(alg::DP5, beta2) = typeof(beta2)(1 // alg_order(alg)) - 3beta2 / 4
+
+alg_stability_size(alg::Tsit5) = 3.5068
+alg_stability_size(alg::DP5) = 3.3066
+
+ssp_coefficient(alg::Euler) = 1
+
+function DiffEqBase.prepare_alg(
+        alg::SplitEuler,
+        u0::AbstractArray,
+        p, prob)
+    alg
+end