refactor(IMEX): migrate KenCarp/Kvaerno to IMEXTableau and remove per-method implementations

Author: Harsh Singh

Project.toml

@@ -10,6 +10,7 @@ ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
 CommonSolve = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2"
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
+DiffEqDevTools = "f3b72e0c-5b89-59e1-b016-84e28bfd966d"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 EnumX = "4e289a0a-7415-4d19-859d-a7e5c4648b56"
 ExponentialUtilities = "d4d017d3-3776-5f7e-afef-a10c40355c18"
@@ -27,6 +28,7 @@ Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
 MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221"
 NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
+ODEProblemLibrary = "fdc4e326-1af4-4b90-96e7-779fcce2daa5"
 OrdinaryDiffEqAdamsBashforthMoulton = "89bda076-bce5-4f1c-845f-551c83cdda9a"
 OrdinaryDiffEqBDF = "6ad6398a-0878-4a85-9266-38940aa047c8"
 OrdinaryDiffEqCore = "bbf590c4-e513-4bbe-9b18-05decba2e5d8"
@@ -73,38 +75,6 @@ StaticArrayInterface = "0d7ed370-da01-4f52-bd93-41d350b8b718"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 TruncatedStacktraces = "781d530d-4396-4725-bb49-402e4bee1e77"
 
-[sources]
-OrdinaryDiffEqAdamsBashforthMoulton = {path = "lib/OrdinaryDiffEqAdamsBashforthMoulton"}
-OrdinaryDiffEqBDF = {path = "lib/OrdinaryDiffEqBDF"}
-OrdinaryDiffEqCore = {path = "lib/OrdinaryDiffEqCore"}
-OrdinaryDiffEqDefault = {path = "lib/OrdinaryDiffEqDefault"}
-OrdinaryDiffEqDifferentiation = {path = "lib/OrdinaryDiffEqDifferentiation"}
-OrdinaryDiffEqExplicitRK = {path = "lib/OrdinaryDiffEqExplicitRK"}
-OrdinaryDiffEqExponentialRK = {path = "lib/OrdinaryDiffEqExponentialRK"}
-OrdinaryDiffEqExtrapolation = {path = "lib/OrdinaryDiffEqExtrapolation"}
-OrdinaryDiffEqFIRK = {path = "lib/OrdinaryDiffEqFIRK"}
-OrdinaryDiffEqFeagin = {path = "lib/OrdinaryDiffEqFeagin"}
-OrdinaryDiffEqFunctionMap = {path = "lib/OrdinaryDiffEqFunctionMap"}
-OrdinaryDiffEqHighOrderRK = {path = "lib/OrdinaryDiffEqHighOrderRK"}
-OrdinaryDiffEqIMEXMultistep = {path = "lib/OrdinaryDiffEqIMEXMultistep"}
-OrdinaryDiffEqLinear = {path = "lib/OrdinaryDiffEqLinear"}
-OrdinaryDiffEqLowOrderRK = {path = "lib/OrdinaryDiffEqLowOrderRK"}
-OrdinaryDiffEqLowStorageRK = {path = "lib/OrdinaryDiffEqLowStorageRK"}
-OrdinaryDiffEqNonlinearSolve = {path = "lib/OrdinaryDiffEqNonlinearSolve"}
-OrdinaryDiffEqNordsieck = {path = "lib/OrdinaryDiffEqNordsieck"}
-OrdinaryDiffEqPDIRK = {path = "lib/OrdinaryDiffEqPDIRK"}
-OrdinaryDiffEqPRK = {path = "lib/OrdinaryDiffEqPRK"}
-OrdinaryDiffEqQPRK = {path = "lib/OrdinaryDiffEqQPRK"}
-OrdinaryDiffEqRKN = {path = "lib/OrdinaryDiffEqRKN"}
-OrdinaryDiffEqRosenbrock = {path = "lib/OrdinaryDiffEqRosenbrock"}
-OrdinaryDiffEqSDIRK = {path = "lib/OrdinaryDiffEqSDIRK"}
-OrdinaryDiffEqSSPRK = {path = "lib/OrdinaryDiffEqSSPRK"}
-OrdinaryDiffEqStabilizedIRK = {path = "lib/OrdinaryDiffEqStabilizedIRK"}
-OrdinaryDiffEqStabilizedRK = {path = "lib/OrdinaryDiffEqStabilizedRK"}
-OrdinaryDiffEqSymplecticRK = {path = "lib/OrdinaryDiffEqSymplecticRK"}
-OrdinaryDiffEqTsit5 = {path = "lib/OrdinaryDiffEqTsit5"}
-OrdinaryDiffEqVerner = {path = "lib/OrdinaryDiffEqVerner"}
-
 [compat]
 ADTypes = "1.16"
 Adapt = "4.3"
@@ -172,10 +142,10 @@ SciMLBase = "2.143"
 SciMLOperators = "1.8"
 SciMLStructures = "1.7"
 SimpleNonlinearSolve = "2.7"
+StableRNGs = "1"
 Static = "1.2"
 StaticArrayInterface = "1.8"
 StaticArrays = "1.9.14"
-StableRNGs = "1"
 StructArrays = "0.6, 0.7"
 Symbolics = "6, 7"
 TruncatedStacktraces = "1.4"
@@ -198,14 +168,104 @@ Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 RecursiveFactorization = "f2c3362d-daeb-58d1-803e-2bc74f2840b4"
 SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
-StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 SparseConnectivityTracer = "9f842d2f-2579-4b1d-911e-f412cf18a3f5"
 SparseMatrixColorings = "0a514795-09f3-496d-8182-132a7b665d35"
+StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 StructArrays = "09ab397b-f2b6-538f-b94a-2f83cf4a842a"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
+[sources.OrdinaryDiffEqAdamsBashforthMoulton]
+path = "lib/OrdinaryDiffEqAdamsBashforthMoulton"
+
+[sources.OrdinaryDiffEqBDF]
+path = "lib/OrdinaryDiffEqBDF"
+
+[sources.OrdinaryDiffEqCore]
+path = "lib/OrdinaryDiffEqCore"
+
+[sources.OrdinaryDiffEqDefault]
+path = "lib/OrdinaryDiffEqDefault"
+
+[sources.OrdinaryDiffEqDifferentiation]
+path = "lib/OrdinaryDiffEqDifferentiation"
+
+[sources.OrdinaryDiffEqExplicitRK]
+path = "lib/OrdinaryDiffEqExplicitRK"
+
+[sources.OrdinaryDiffEqExponentialRK]
+path = "lib/OrdinaryDiffEqExponentialRK"
+
+[sources.OrdinaryDiffEqExtrapolation]
+path = "lib/OrdinaryDiffEqExtrapolation"
+
+[sources.OrdinaryDiffEqFIRK]
+path = "lib/OrdinaryDiffEqFIRK"
+
+[sources.OrdinaryDiffEqFeagin]
+path = "lib/OrdinaryDiffEqFeagin"
+
+[sources.OrdinaryDiffEqFunctionMap]
+path = "lib/OrdinaryDiffEqFunctionMap"
+
+[sources.OrdinaryDiffEqHighOrderRK]
+path = "lib/OrdinaryDiffEqHighOrderRK"
+
+[sources.OrdinaryDiffEqIMEXMultistep]
+path = "lib/OrdinaryDiffEqIMEXMultistep"
+
+[sources.OrdinaryDiffEqLinear]
+path = "lib/OrdinaryDiffEqLinear"
+
+[sources.OrdinaryDiffEqLowOrderRK]
+path = "lib/OrdinaryDiffEqLowOrderRK"
+
+[sources.OrdinaryDiffEqLowStorageRK]
+path = "lib/OrdinaryDiffEqLowStorageRK"
+
+[sources.OrdinaryDiffEqNonlinearSolve]
+path = "lib/OrdinaryDiffEqNonlinearSolve"
+
+[sources.OrdinaryDiffEqNordsieck]
+path = "lib/OrdinaryDiffEqNordsieck"
+
+[sources.OrdinaryDiffEqPDIRK]
+path = "lib/OrdinaryDiffEqPDIRK"
+
+[sources.OrdinaryDiffEqPRK]
+path = "lib/OrdinaryDiffEqPRK"
+
+[sources.OrdinaryDiffEqQPRK]
+path = "lib/OrdinaryDiffEqQPRK"
+
+[sources.OrdinaryDiffEqRKN]
+path = "lib/OrdinaryDiffEqRKN"
+
+[sources.OrdinaryDiffEqRosenbrock]
+path = "lib/OrdinaryDiffEqRosenbrock"
+
+[sources.OrdinaryDiffEqSDIRK]
+path = "lib/OrdinaryDiffEqSDIRK"
+
+[sources.OrdinaryDiffEqSSPRK]
+path = "lib/OrdinaryDiffEqSSPRK"
+
+[sources.OrdinaryDiffEqStabilizedIRK]
+path = "lib/OrdinaryDiffEqStabilizedIRK"
+
+[sources.OrdinaryDiffEqStabilizedRK]
+path = "lib/OrdinaryDiffEqStabilizedRK"
+
+[sources.OrdinaryDiffEqSymplecticRK]
+path = "lib/OrdinaryDiffEqSymplecticRK"
+
+[sources.OrdinaryDiffEqTsit5]
+path = "lib/OrdinaryDiffEqTsit5"
+
+[sources.OrdinaryDiffEqVerner]
+path = "lib/OrdinaryDiffEqVerner"
+
 [targets]
 test = ["ComponentArrays", "AlgebraicMultigrid", "IncompleteLU", "DiffEqCallbacks", "DifferentiationInterface", "DiffEqDevTools", "ExplicitImports", "ODEProblemLibrary", "ElasticArrays", "JLArrays", "Random", "SafeTestsets", "StableRNGs", "StructArrays", "Test", "Unitful", "Pkg", "NLsolve", "RecursiveFactorization", "SparseConnectivityTracer", "SparseMatrixColorings", "Statistics"]

lib/OrdinaryDiffEqSDIRK/src/algorithms.jl

@@ -79,6 +79,9 @@ function SDIRK_docstring(
     )
 end
 
+abstract type OrdinaryDiffEqNewtonAdaptiveIMEXAlgorithm{CS, AD, FDT, ST, CJ} <:
+OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ} end
+
 @doc SDIRK_docstring(
     "A 1st order implicit solver. A-B-L-stable. Adaptive timestepping through a divided differences estimate. Strong-stability preserving (SSP). Good for highly stiff equations.",
     "ImplicitEuler";
@@ -486,7 +489,7 @@ end
     """
 )
 struct Kvaerno3{CS, AD, F, F2, P, FDT, ST, CJ, StepLimiter} <:
-    OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ}
+    OrdinaryDiffEqNewtonAdaptiveIMEXAlgorithm{CS, AD, FDT, ST, CJ}
     linsolve::F
     nlsolve::F2
     precs::P
@@ -538,7 +541,7 @@ end
     """
 )
 struct KenCarp3{CS, AD, F, F2, P, FDT, ST, CJ, StepLimiter} <:
-    OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ}
+    OrdinaryDiffEqNewtonAdaptiveIMEXAlgorithm{CS, AD, FDT, ST, CJ}
     linsolve::F
     nlsolve::F2
     precs::P
@@ -1045,7 +1048,7 @@ end
     """
 )
 struct Kvaerno4{CS, AD, F, F2, P, FDT, ST, CJ, StepLimiter} <:
-    OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ}
+    OrdinaryDiffEqNewtonAdaptiveIMEXAlgorithm{CS, AD, FDT, ST, CJ}
     linsolve::F
     nlsolve::F2
     precs::P
@@ -1101,7 +1104,7 @@ end
     """
 )
 struct Kvaerno5{CS, AD, F, F2, P, FDT, ST, CJ, StepLimiter} <:
-    OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ}
+    OrdinaryDiffEqNewtonAdaptiveIMEXAlgorithm{CS, AD, FDT, ST, CJ}
     linsolve::F
     nlsolve::F2
     precs::P
@@ -1153,7 +1156,7 @@ end
     """
 )
 struct KenCarp4{CS, AD, F, F2, P, FDT, ST, CJ, StepLimiter} <:
-    OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ}
+    OrdinaryDiffEqNewtonAdaptiveIMEXAlgorithm{CS, AD, FDT, ST, CJ}
     linsolve::F
     nlsolve::F2
     precs::P
@@ -1207,14 +1210,15 @@ end
     controller = :PI,
     """
 )
-struct KenCarp47{CS, AD, F, F2, P, FDT, ST, CJ} <:
-    OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ}
+struct KenCarp47{CS, AD, F, F2, P, FDT, ST, CJ, StepLimiter} <:
+    OrdinaryDiffEqNewtonAdaptiveIMEXAlgorithm{CS, AD, FDT, ST, CJ}
     linsolve::F
     nlsolve::F2
     precs::P
     smooth_est::Bool
     extrapolant::Symbol
     controller::Symbol
+    step_limiter!::StepLimiter
     autodiff::AD
 end
 function KenCarp47(;
@@ -1223,17 +1227,17 @@ function KenCarp47(;
         diff_type = Val{:forward}(),
         linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(),
         smooth_est = true, extrapolant = :linear,
-        controller = :PI
+        controller = :PI, step_limiter! = trivial_limiter!
     )
     AD_choice, chunk_size, diff_type = _process_AD_choice(autodiff, chunk_size, diff_type)
 
     return KenCarp47{
         _unwrap_val(chunk_size), typeof(AD_choice), typeof(linsolve),
         typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag),
-        _unwrap_val(concrete_jac),
+        _unwrap_val(concrete_jac), typeof(step_limiter!),
     }(
         linsolve, nlsolve, precs, smooth_est, extrapolant,
-        controller, AD_choice
+        controller, step_limiter!, AD_choice
     )
 end
 
@@ -1259,7 +1263,7 @@ end
     """
 )
 struct KenCarp5{CS, AD, F, F2, P, FDT, ST, CJ, StepLimiter} <:
-    OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ}
+    OrdinaryDiffEqNewtonAdaptiveIMEXAlgorithm{CS, AD, FDT, ST, CJ}
     linsolve::F
     nlsolve::F2
     precs::P
@@ -1311,14 +1315,15 @@ end
     controller = :PI,
     """
 )
-struct KenCarp58{CS, AD, F, F2, P, FDT, ST, CJ} <:
-    OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ}
+struct KenCarp58{CS, AD, F, F2, P, FDT, ST, CJ, StepLimiter} <:
+    OrdinaryDiffEqNewtonAdaptiveIMEXAlgorithm{CS, AD, FDT, ST, CJ}
     linsolve::F
     nlsolve::F2
     precs::P
     smooth_est::Bool
     extrapolant::Symbol
     controller::Symbol
+    step_limiter!::StepLimiter
     autodiff::AD
 end
 function KenCarp58(;
@@ -1327,17 +1332,17 @@ function KenCarp58(;
         diff_type = Val{:forward}(),
         linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(),
         smooth_est = true, extrapolant = :linear,
-        controller = :PI
+        controller = :PI, step_limiter! = trivial_limiter!
     )
     AD_choice, chunk_size, diff_type = _process_AD_choice(autodiff, chunk_size, diff_type)
 
     return KenCarp58{
         _unwrap_val(chunk_size), typeof(AD_choice), typeof(linsolve),
         typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag),
-        _unwrap_val(concrete_jac),
+        _unwrap_val(concrete_jac), typeof(step_limiter!),
     }(
         linsolve, nlsolve, precs, smooth_est, extrapolant,
-        controller, AD_choice
+        controller, step_limiter!, AD_choice
     )
 end
 
@@ -1594,9 +1599,6 @@ function ESDIRK659L2SA(;
     )
 end
 
-abstract type OrdinaryDiffEqNewtonAdaptiveIMEXAlgorithm{CS, AD, FDT, ST, CJ} <:
-OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ} end
-
 @doc SDIRK_docstring(
     "3rd order L-stable IMEX ARK method. Uses a generic tableau-driven implementation that supports both split and non-split forms.",
     "ARS343";