Added QPRK solvers

lib/OrdinaryDiffEqQPRK/Project.toml

@@ -0,0 +1,21 @@
+name = "OrdinaryDiffEqQPRK"
+uuid = "04162be5-8125-4266-98ed-640baecc6514"
+authors = ["ParamThakkar123 <paramthakkar864@gmail.com>"]
+version = "1.0.0"
+
+[deps]
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221"
+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/OrdinaryDiffEqQPRK/src/OrdinaryDiffEqQPRK.jl

@@ -0,0 +1,18 @@
+module OrdinaryDiffEqQPRK
+
+import OrdinaryDiffEq: OrdinaryDiffEqAdaptiveAlgorithm, OrdinaryDiffEqConstantCache,
+                       explicit_rk_docstring, @cache, @unpack, OrdinaryDiffEqMutableCache,
+                       OrdinaryDiffEqAdaptiveAlgorithm, @fold, @OnDemandTableauExtract,
+                       trivial_limiter!, alg_cache, alg_order, initialize!, perform_step!
+using Static: False
+using MuladdMacro, FastBroadcast
+
+include("algorithms.jl")
+include("alg_utils.jl")
+include("qprk_caches.jl")
+include("qprk_tableaus.jl")
+include("qprk_perform_step.jl")
+
+export QPRK98
+
+end

lib/OrdinaryDiffEqQPRK/src/alg_utils.jl

@@ -0,0 +1 @@
+alg_order(alg::QPRK98) = 9

lib/OrdinaryDiffEqQPRK/src/algorithms.jl

@@ -0,0 +1,15 @@
+@doc explicit_rk_docstring(
+    "Runge–Kutta pairs of orders 9(8) for use in quadruple precision computations", "QPRK98",
+    references = "Kovalnogov VN, Fedorov RV, Karpukhina TV, Simos TE, Tsitouras C. Runge–Kutta pairs 
+    of orders 9 (8) for use in quadruple precision computations. Numerical Algorithms, 2023. 
+    doi: https://doi.org/10.1007/s11075-023-01632-8")
+Base.@kwdef struct QPRK98{StageLimiter, StepLimiter, Thread} <:
+                   OrdinaryDiffEqAdaptiveAlgorithm
+    stage_limiter!::StageLimiter = trivial_limiter!
+    step_limiter!::StepLimiter = trivial_limiter!
+    thread::Thread = False()
+end
+# for backwards compatibility
+function QPRK98(stage_limiter!, step_limiter! = trivial_limiter!)
+    QPRK98(stage_limiter!, step_limiter!, False())
+end

src/caches/qprk_caches.jl → lib/OrdinaryDiffEqQPRK/src/qprk_caches.jl

@@ -65,4 +65,4 @@ function alg_cache(::QPRK98, u, rate_prototype, ::Type{uEltypeNoUnits},
         dt, reltol, p, calck,
         ::Val{false}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
     QPRK98ConstantCache()
-end
+end

src/perform_step/qprk_perform_step.jl → lib/OrdinaryDiffEqQPRK/src/qprk_perform_step.jl

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

src/tableaus/qprk_tableaus.jl → lib/OrdinaryDiffEqQPRK/src/qprk_tableaus.jl

@@ -232,4 +232,4 @@ end
         b14_13, b15_1, b15_6, b15_7, b15_8, b15_9, b15_10, b15_11, b15_12, b15_13,
         b15_14, b16_1, b16_6, b16_7, b16_8, b16_9, b16_10, b16_11, b16_12, b16_13,
         b16_14)
-end
+end

src/OrdinaryDiffEq.jl

@@ -159,12 +159,10 @@ include("caches/adams_bashforth_moulton_caches.jl")
 include("caches/nordsieck_caches.jl")
 include("caches/prk_caches.jl")
 include("caches/pdirk_caches.jl")
-include("caches/qprk_caches.jl")
 
 include("tableaus/low_order_rk_tableaus.jl")
 include("tableaus/high_order_rk_tableaus.jl")
 include("tableaus/rosenbrock_tableaus.jl")
-include("tableaus/qprk_tableaus.jl")
 
 include("integrators/type.jl")
 include("integrators/controllers.jl")
@@ -187,7 +185,6 @@ include("perform_step/adams_bashforth_moulton_perform_step.jl")
 include("perform_step/nordsieck_perform_step.jl")
 include("perform_step/prk_perform_step.jl")
 include("perform_step/pdirk_perform_step.jl")
-include("perform_step/qprk_perform_step.jl")
 
 include("dense/generic_dense.jl")
 include("dense/interpolants.jl")
@@ -288,6 +285,10 @@ include("../lib/OrdinaryDiffEqFIRK/src/OrdinaryDiffEqFIRK.jl")
 using ..OrdinaryDiffEqFIRK
 export RadauIIA3, RadauIIA5, RadauIIA9
 
+include("../lib/OrdinaryDiffEqQPRK/src/OrdinaryDiffEqQPRK.jl")
+using ..OrdinaryDiffEqQPRK
+export QPRK98
+
 PrecompileTools.@compile_workload begin
     function lorenz(du, u, p, t)
         du[1] = 10.0(u[2] - u[1])
@@ -420,7 +421,7 @@ export FunctionMap, Euler, Heun, Ralston, Midpoint, RK4, ExplicitRK, OwrenZen3,
        OwrenZen5,
        BS3, BS5, DP5, Tsit5, DP8, TanYam7, TsitPap8, CompositeAlgorithm, Anas5, RKO65,
        FRK65, PFRK87,
-       RKM, MSRK5, MSRK6, Stepanov5, SIR54, QPRK98, PSRK4p7q6, PSRK3p6q5, PSRK3p5q4
+       RKM, MSRK5, MSRK6, Stepanov5, SIR54, PSRK4p7q6, PSRK3p6q5, PSRK3p5q4
 
 export MagnusMidpoint, LinearExponential, MagnusLeapfrog, LieEuler, CayleyEuler,
        MagnusGauss4, MagnusNC6, MagnusGL6, MagnusGL8, MagnusNC8, MagnusGL4,

src/alg_utils.jl

@@ -507,8 +507,6 @@ alg_order(alg::Alshina2) = 2
 alg_order(alg::Alshina3) = 3
 alg_order(alg::Alshina6) = 6
 
-alg_order(alg::QPRK98) = 9
-
 alg_maximum_order(alg) = alg_order(alg)
 alg_maximum_order(alg::CompositeAlgorithm) = maximum(alg_order(x) for x in alg.algs)
 

src/algorithms/explicit_rk.jl

@@ -482,20 +482,4 @@ This requires that simultaneous calls to f are thread-safe.
 """
 Base.@kwdef struct KuttaPRK2p5{TO} <: OrdinaryDiffEqAlgorithm
     threading::TO = true
-end
-
-@doc explicit_rk_docstring(
-    "Runge–Kutta pairs of orders 9(8) for use in quadruple precision computations", "QPRK98",
-    references = "Kovalnogov VN, Fedorov RV, Karpukhina TV, Simos TE, Tsitouras C. Runge–Kutta pairs 
-    of orders 9 (8) for use in quadruple precision computations. Numerical Algorithms, 2023. 
-    doi: https://doi.org/10.1007/s11075-023-01632-8")
-Base.@kwdef struct QPRK98{StageLimiter, StepLimiter, Thread} <:
-                   OrdinaryDiffEqAdaptiveAlgorithm
-    stage_limiter!::StageLimiter = trivial_limiter!
-    step_limiter!::StepLimiter = trivial_limiter!
-    thread::Thread = False()
-end
-# for backwards compatibility
-function QPRK98(stage_limiter!, step_limiter! = trivial_limiter!)
-    QPRK98(stage_limiter!, step_limiter!, False())
-end
+end