Merge branch 'master' of https://github.com/Shreyas-Ekanathan/Ordinar…

…yDiffEq.jl

Project.toml

@@ -1,7 +1,7 @@
 name = "OrdinaryDiffEq"
 uuid = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 authors = ["Chris Rackauckas <[email protected]>", "Yingbo Ma <[email protected]>"]
-version = "6.86.0"
+version = "6.87.0"
 
 [deps]
 ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
@@ -101,6 +101,7 @@ InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56"
 ODEProblemLibrary = "fdc4e326-1af4-4b90-96e7-779fcce2daa5"
+ParameterizedFunctions = "65888b18-ceab-5e60-b2b9-181511a3b968"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 PoissonRandom = "e409e4f3-bfea-5376-8464-e040bb5c01ab"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
@@ -114,4 +115,4 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [targets]
-test = ["Calculus", "ComponentArrays", "Symbolics", "AlgebraicMultigrid", "IncompleteLU", "DiffEqCallbacks", "DiffEqDevTools", "ODEProblemLibrary", "ElasticArrays", "InteractiveUtils", "PoissonRandom", "Printf", "Random", "ReverseDiff", "SafeTestsets", "SparseArrays", "Statistics", "Test", "Unitful", "ModelingToolkit", "Pkg", "NLsolve"]
+test = ["Calculus", "ComponentArrays", "Symbolics", "AlgebraicMultigrid", "IncompleteLU", "DiffEqCallbacks", "DiffEqDevTools", "ODEProblemLibrary", "ElasticArrays", "InteractiveUtils", "ParameterizedFunctions", "PoissonRandom", "Printf", "Random", "ReverseDiff", "SafeTestsets", "SparseArrays", "Statistics", "Test", "Unitful", "ModelingToolkit", "Pkg", "NLsolve"]

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,35 @@
+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,
+                       stepsize_controller!, islinearfunction, step_accept_controller!, step_reject_controller!,
+                       PredictiveController, alg_can_repeat_jac, NewtonAlgorithm, fac_default_gamma,
+                       get_new_W_γdt_cutoff, get_current_adaptive_order, VerySlowConvergence, 
+                       Divergence, isfirk
+using MuladdMacro, DiffEqBase, RecursiveArrayTools
+using SciMLOperators: AbstractSciMLOperator
+using LinearAlgebra: I, UniformScaling, mul!, lu
+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, RadauIIA9
+
+end

lib/OrdinaryDiffEqFIRK/src/alg_utils.jl

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

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}
+}
+
+RadauIIA9: Fully-Implicit Runge-Kutta Method
+An A-B-L stable fully implicit Runge-Kutta method with internal tableau complex basis transform for efficiency.
+"""
+struct RadauIIA9{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 RadauIIA9(; 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!)
+    RadauIIA9{_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,51 @@
+@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 isfirk(alg)
+                @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 RadauIIA9Cache{uType, cuType, uNoUnitsType, rateType, JType, W1Ty
     status::NLStatus
     step_limiter!::StepLimiter
 end
-TruncatedStacktraces.@truncate_stacktrace RadauIIA9Cache 1
 
 function alg_cache(alg::RadauIIA9, 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, RadauIIA9},
+        cache::OrdinaryDiffEqMutableCache)
+    (cache.tmp, cache.atmp)
+end

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

@@ -8,10 +8,10 @@ for prob in [prob_ode_linear, prob_ode_2Dlinear]
     @test sim21.𝒪est[:final]≈5 atol=testTol
 end
 
-sim21 = test_convergence(1 ./ 2 .^ (2.75:-0.5:0.25), prob_ode_linear, RadauIIA9())
+sim21 = test_convergence(1 ./ 2 .^ (2.5:-1:0.5), prob_ode_linear, RadauIIA9())
 @test sim21.𝒪est[:final]≈8 atol=testTol
 
-sim21 = test_convergence(1 ./ 2 .^ (2.75:-0.5:0.25), prob_ode_2Dlinear, RadauIIA0())
+sim21 = test_convergence(1 ./ 2 .^ (2.5:-1:0.5), prob_ode_2Dlinear, RadauIIA9())
 @test sim21.𝒪est[:final]≈8 atol=testTol
 
 # test adaptivity

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")