Added FIRK solvers

lib/OrdinaryDiffEqExtrapolation/test/runtests.jl

@@ -1,3 +1,3 @@
 using SafeTestsets
 
-@time @safetestset "Extrapolation Tests" include("ode_extrapolation_tests.jl")
+@time @safetestset "Extrapolation Tests" include("ode_extrapolation_tests.jl")

lib/OrdinaryDiffEqFIRK/Project.toml

@@ -0,0 +1,24 @@
+name = "OrdinaryDiffEqFIRK"
+uuid = "5960d6e9-dd7a-4743-88e7-cf307b64f125"
+authors = ["ParamThakkar123 <[email protected]>"]
+version = "1.0.0"
+
+[deps]
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
+RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
+LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
+SciMLOperators = "c0aeaf25-5076-4817-a8d5-81caf7dfa961"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+
+[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/OrdinaryDiffEqFIRK/src/OrdinaryDiffEqFIRK.jl

@@ -0,0 +1,32 @@
+module OrdinaryDiffEqFIRK
+
+import OrdinaryDiffEq: alg_order, calculate_residuals!,
+                       initialize!, perform_step!, @unpack, unwrap_alg,
+                       calculate_residuals,
+                       OrdinaryDiffEqAlgorithm, OrdinaryDiffEqNewtonAdaptiveAlgorithm,
+                       OrdinaryDiffEqMutableCache, OrdinaryDiffEqConstantCache,
+                       OrdinaryDiffEqAdaptiveAlgorithm, CompiledFloats, uses_uprev,
+                       alg_cache, _vec, _reshape, @cache, isfsal, full_cache,
+                       constvalue, _unwrap_val, du_alias_or_new,
+                       explicit_rk_docstring, trivial_limiter!,
+                       _ode_interpolant!, _ode_addsteps!, AbstractController,
+                       NLStatus, qmax_default, alg_adaptive_order, DEFAULT_PRECS,
+                       UJacobianWrapper, build_J_W, build_jac_config, UDerivativeWrapper,
+                       Convergence, calc_J!, dolinsolve, FastConvergence, calc_J,
+                       mul!
+using MuladdMacro, DiffEqBase, RecursiveArrayTools
+using SciMLOperators: AbstractSciMLOperator
+using LinearAlgebra: I, UniformScaling
+import LinearSolve
+import FastBroadcast: @..
+include("algorithms.jl")
+include("alg_utils.jl")
+include("controllers.jl")
+include("firk_caches.jl")
+include("firk_tableaus.jl")
+include("firk_perform_step.jl")
+include("integrator_interface.jl")
+
+export RadauIIA3, RadauIIA5
+
+end

lib/OrdinaryDiffEqFIRK/src/alg_utils.jl

@@ -0,0 +1,9 @@
+qmax_default(alg::Union{RadauIIA3, RadauIIA5, RadauIIA7}) = 8
+
+alg_order(alg::RadauIIA3) = 3
+alg_order(alg::RadauIIA5) = 5
+alg_order(alg::RadauIIA7) = 7
+
+alg_adaptive_order(alg::RadauIIA3) = 1
+alg_adaptive_order(alg::RadauIIA5) = 3
+alg_adaptive_order(alg::RadauIIA7) = 5

lib/OrdinaryDiffEqFIRK/src/algorithms.jl

@@ -0,0 +1,155 @@
+# FIRK Methods
+
+"""
+@article{hairer1999stiff,
+title={Stiff differential equations solved by Radau methods},
+author={Hairer, Ernst and Wanner, Gerhard},
+journal={Journal of Computational and Applied Mathematics},
+volume={111},
+number={1-2},
+pages={93--111},
+year={1999},
+publisher={Elsevier}
+}
+
+RadauIIA3: Fully-Implicit Runge-Kutta Method
+An A-B-L stable fully implicit Runge-Kutta method with internal tableau complex basis transform for efficiency.
+"""
+struct RadauIIA3{CS, AD, F, P, FDT, ST, CJ, Tol, C1, C2, StepLimiter} <:
+       OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ}
+    linsolve::F
+    precs::P
+    extrapolant::Symbol
+    κ::Tol
+    maxiters::Int
+    fast_convergence_cutoff::C1
+    new_W_γdt_cutoff::C2
+    controller::Symbol
+    step_limiter!::StepLimiter
+end
+
+function RadauIIA3(; chunk_size = Val{0}(), autodiff = Val{true}(),
+        standardtag = Val{true}(), concrete_jac = nothing,
+        diff_type = Val{:forward},
+        linsolve = nothing, precs = DEFAULT_PRECS,
+        extrapolant = :dense, fast_convergence_cutoff = 1 // 5,
+        new_W_γdt_cutoff = 1 // 5,
+        controller = :Predictive, κ = nothing, maxiters = 10,
+        step_limiter! = trivial_limiter!)
+    RadauIIA3{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve),
+        typeof(precs), diff_type, _unwrap_val(standardtag), _unwrap_val(concrete_jac),
+        typeof(κ), typeof(fast_convergence_cutoff),
+        typeof(new_W_γdt_cutoff), typeof(step_limiter!)}(linsolve,
+        precs,
+        extrapolant,
+        κ,
+        maxiters,
+        fast_convergence_cutoff,
+        new_W_γdt_cutoff,
+        controller,
+        step_limiter!)
+end
+
+"""
+@article{hairer1999stiff,
+title={Stiff differential equations solved by Radau methods},
+author={Hairer, Ernst and Wanner, Gerhard},
+journal={Journal of Computational and Applied Mathematics},
+volume={111},
+number={1-2},
+pages={93--111},
+year={1999},
+publisher={Elsevier}
+}
+
+RadauIIA5: Fully-Implicit Runge-Kutta Method
+An A-B-L stable fully implicit Runge-Kutta method with internal tableau complex basis transform for efficiency.
+"""
+struct RadauIIA5{CS, AD, F, P, FDT, ST, CJ, Tol, C1, C2, StepLimiter} <:
+       OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ}
+    linsolve::F
+    precs::P
+    smooth_est::Bool
+    extrapolant::Symbol
+    κ::Tol
+    maxiters::Int
+    fast_convergence_cutoff::C1
+    new_W_γdt_cutoff::C2
+    controller::Symbol
+    step_limiter!::StepLimiter
+end
+
+function RadauIIA5(; chunk_size = Val{0}(), autodiff = Val{true}(),
+        standardtag = Val{true}(), concrete_jac = nothing,
+        diff_type = Val{:forward},
+        linsolve = nothing, precs = DEFAULT_PRECS,
+        extrapolant = :dense, fast_convergence_cutoff = 1 // 5,
+        new_W_γdt_cutoff = 1 // 5,
+        controller = :Predictive, κ = nothing, maxiters = 10, smooth_est = true,
+        step_limiter! = trivial_limiter!)
+    RadauIIA5{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve),
+        typeof(precs), diff_type, _unwrap_val(standardtag), _unwrap_val(concrete_jac),
+        typeof(κ), typeof(fast_convergence_cutoff),
+        typeof(new_W_γdt_cutoff), typeof(step_limiter!)}(linsolve,
+        precs,
+        smooth_est,
+        extrapolant,
+        κ,
+        maxiters,
+        fast_convergence_cutoff,
+        new_W_γdt_cutoff,
+        controller,
+        step_limiter!)
+end
+
+"""
+@article{hairer1999stiff,
+title={Stiff differential equations solved by Radau methods},
+author={Hairer, Ernst and Wanner, Gerhard},
+journal={Journal of Computational and Applied Mathematics},
+volume={111},
+number={1-2},
+pages={93--111},
+year={1999},
+publisher={Elsevier}
+}
+
+RadauIIA7: Fully-Implicit Runge-Kutta Method
+An A-B-L stable fully implicit Runge-Kutta method with internal tableau complex basis transform for efficiency.
+"""
+struct RadauIIA7{CS, AD, F, P, FDT, ST, CJ, Tol, C1, C2, StepLimiter} <:
+       OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ}
+    linsolve::F
+    precs::P
+    smooth_est::Bool
+    extrapolant::Symbol
+    κ::Tol
+    maxiters::Int
+    fast_convergence_cutoff::C1
+    new_W_γdt_cutoff::C2
+    controller::Symbol
+    step_limiter!::StepLimiter
+end
+
+function RadauIIA7(; chunk_size = Val{0}(), autodiff = Val{true}(),
+        standardtag = Val{true}(), concrete_jac = nothing,
+        diff_type = Val{:forward},
+        linsolve = nothing, precs = DEFAULT_PRECS,
+        extrapolant = :dense, fast_convergence_cutoff = 1 // 5,
+        new_W_γdt_cutoff = 1 // 5,
+        controller = :Predictive, κ = nothing, maxiters = 10, smooth_est = true,
+        step_limiter! = trivial_limiter!)
+    RadauIIA7{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve),
+        typeof(precs), diff_type, _unwrap_val(standardtag), _unwrap_val(concrete_jac),
+        typeof(κ), typeof(fast_convergence_cutoff),
+        typeof(new_W_γdt_cutoff), typeof(step_limiter!)}(linsolve,
+        precs,
+        smooth_est,
+        extrapolant,
+        κ,
+        maxiters,
+        fast_convergence_cutoff,
+        new_W_γdt_cutoff,
+        controller,
+        step_limiter!)
+end

lib/OrdinaryDiffEqFIRK/src/controllers.jl

@@ -0,0 +1,54 @@
+struct PredictiveController <: AbstractController
+end
+
+@inline function stepsize_controller!(integrator, controller::PredictiveController, alg)
+    @unpack qmin, qmax, gamma = integrator.opts
+    EEst = DiffEqBase.value(integrator.EEst)
+
+    if iszero(EEst)
+        q = inv(qmax)
+    else
+        if fac_default_gamma(alg)
+            fac = gamma
+        else
+            if alg isa Union{RadauIIA3, RadauIIA5, RadauIIA7}
+                @unpack iter = integrator.cache
+                @unpack maxiters = alg
+            else
+                @unpack iter, maxiters = integrator.cache.nlsolver
+            end
+            fac = min(gamma, (1 + 2 * maxiters) * gamma / (iter + 2 * maxiters))
+        end
+        expo = 1 / (get_current_adaptive_order(alg, integrator.cache) + 1)
+        qtmp = DiffEqBase.fastpow(EEst, expo) / fac
+        @fastmath q = DiffEqBase.value(max(inv(qmax), min(inv(qmin), qtmp)))
+        integrator.qold = q
+    end
+    q
+end
+
+function step_accept_controller!(integrator, controller::PredictiveController, alg, q)
+    @unpack qmin, qmax, gamma, qsteady_min, qsteady_max = integrator.opts
+    EEst = DiffEqBase.value(integrator.EEst)
+
+    if integrator.success_iter > 0
+        expo = 1 / (get_current_adaptive_order(alg, integrator.cache) + 1)
+        qgus = (integrator.dtacc / integrator.dt) *
+               DiffEqBase.fastpow((EEst^2) / integrator.erracc, expo)
+        qgus = max(inv(qmax), min(inv(qmin), qgus / gamma))
+        qacc = max(q, qgus)
+    else
+        qacc = q
+    end
+    if qsteady_min <= qacc <= qsteady_max
+        qacc = one(qacc)
+    end
+    integrator.dtacc = integrator.dt
+    integrator.erracc = max(1e-2, EEst)
+    return integrator.dt / qacc
+end
+
+function step_reject_controller!(integrator, controller::PredictiveController, alg)
+    @unpack dt, success_iter, qold = integrator
+    integrator.dt = success_iter == 0 ? 0.1 * dt : dt / qold
+end

src/caches/firk_caches.jl → lib/OrdinaryDiffEqFIRK/src/firk_caches.jl

@@ -200,7 +200,6 @@ mutable struct RadauIIA5Cache{uType, cuType, uNoUnitsType, rateType, JType, W1Ty
     status::NLStatus
     step_limiter!::StepLimiter
 end
-TruncatedStacktraces.@truncate_stacktrace RadauIIA5Cache 1
 
 function alg_cache(alg::RadauIIA5, u, rate_prototype, ::Type{uEltypeNoUnits},
         ::Type{uBottomEltypeNoUnits},
@@ -370,7 +369,6 @@ mutable struct RadauIIA7Cache{uType, cuType, uNoUnitsType, rateType, JType, W1Ty
     status::NLStatus
     step_limiter!::StepLimiter
 end
-TruncatedStacktraces.@truncate_stacktrace RadauIIA7Cache 1
 
 function alg_cache(alg::RadauIIA7, u, rate_prototype, ::Type{uEltypeNoUnits},
         ::Type{uBottomEltypeNoUnits},

lib/OrdinaryDiffEqFIRK/src/integrator_interface.jl

@@ -0,0 +1,4 @@
+@inline function DiffEqBase.get_tmp_cache(integrator, alg::Union{RadauIIA3, RadauIIA5, RadauIIA7},
+        cache::OrdinaryDiffEqMutableCache)
+    (cache.tmp, cache.atmp)
+end

test/algconvergence/ode_firk_tests.jl → lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl

@@ -57,4 +57,4 @@ for iip in (true, false)
         @test sol.stats.njacs < sol.stats.nw # W reuse
     end
     @test length(sol) < 5000 # the error estimate is not very good
-end
+end

lib/OrdinaryDiffEqFIRK/test/runtests.jl

@@ -0,0 +1,3 @@
+using SafeTestsets
+
+@time @safetestset "FIRK Tests" include("ode_firk_tests.jl")

lib/OrdinaryDiffEqFeagin/test/runtests.jl

@@ -1 +1,4 @@
 
+using SafeTestsets
+
+@time @safetestset "Feagin Tests" include("ode_feagin_tests.jl")

src/OrdinaryDiffEq.jl

@@ -148,7 +148,6 @@ include("generic_rosenbrock.jl")
 include("caches/basic_caches.jl")
 include("caches/low_order_rk_caches.jl")
 include("caches/high_order_rk_caches.jl")
-include("caches/firk_caches.jl")
 include("caches/linear_caches.jl")
 include("caches/linear_nonlinear_caches.jl")
 include("caches/rosenbrock_caches.jl")
@@ -162,7 +161,6 @@ 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/firk_tableaus.jl")
 include("tableaus/qprk_tableaus.jl")
 
 include("integrators/type.jl")
@@ -181,7 +179,6 @@ include("perform_step/exponential_rk_perform_step.jl")
 include("perform_step/explicit_rk_perform_step.jl")
 include("perform_step/low_order_rk_perform_step.jl")
 include("perform_step/high_order_rk_perform_step.jl")
-include("perform_step/firk_perform_step.jl")
 include("perform_step/rosenbrock_perform_step.jl")
 include("perform_step/composite_perform_step.jl")
 include("perform_step/adams_bashforth_moulton_perform_step.jl")
@@ -284,11 +281,16 @@ include("../lib/OrdinaryDiffEqDefault/src/OrdinaryDiffEqDefault.jl")
 using ..OrdinaryDiffEqDefault
 export DefaultODEAlgorithm
 
+include("../lib/OrdinaryDiffEqFIRK/src/OrdinaryDiffEqFIRK.jl")
+using ..OrdinaryDiffEqFIRK
+export RadauIIA3, RadauIIA5, RadauIIA7
+
 using ..OrdinaryDiffEqBDF: reinitFBDF!, error_constant, estimate_terk!,
                            calc_Lagrange_interp!,
                            calc_finite_difference_weights, estimate_terk,
                            calc_Lagrange_interp,
                            bdf_step_reject_controller!
+
 include("nlsolve/newton.jl")
 include("perform_step/dae_perform_step.jl")
 
@@ -428,8 +430,6 @@ export FunctionMap, Euler, Heun, Ralston, Midpoint, RK4, ExplicitRK, OwrenZen3,
        FRK65, PFRK87,
        RKM, MSRK5, MSRK6, Stepanov5, SIR54, QPRK98, PSRK4p7q6, PSRK3p6q5, PSRK3p5q4
 
-export RadauIIA3, RadauIIA5, RadauIIA7
-
 export MagnusMidpoint, LinearExponential, MagnusLeapfrog, LieEuler, CayleyEuler,
        MagnusGauss4, MagnusNC6, MagnusGL6, MagnusGL8, MagnusNC8, MagnusGL4,
        MagnusAdapt4, RKMK2, RKMK4, LieRK4, CG2, CG3, CG4a