Merge pull request #2274 from ParamThakkar123/Verner

Added Verner Methods

lib/OrdinaryDiffEqDefault/Project.toml

@@ -0,0 +1,21 @@
+name = "OrdinaryDiffEqDefault"
+uuid = "50262376-6c5a-4cf5-baba-aaf4f84d72d7"
+authors = ["ParamThakkar123 <paramthakkar864@gmail.com>"]
+version = "1.0.0"
+
+[deps]
+EnumX = "4e289a0a-7415-4d19-859d-a7e5c4648b56"
+LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+
+[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"
+
+[compat]
+julia = "1.10"
+
+[targets]
+test = ["DiffEqDevTools", "Random", "SafeTestsets", "Test"]

lib/OrdinaryDiffEqDefault/src/OrdinaryDiffEqDefault.jl

@@ -0,0 +1,15 @@
+module OrdinaryDiffEqDefault
+
+using OrdinaryDiffEq: Vern7, Vern8, Vern9, Vern6, Tsit5, Rosenbrock23, Rodas5P, FBDF,
+    alg_stability_size, beta2_default, beta1_default, AutoSwitchCache, ODEIntegrator,
+    CompositeAlgorithm, OrdinaryDiffEqAlgorithm, OrdinaryDiffEqMutableCache, AutoAlgSwitch
+import OrdinaryDiffEq: is_mass_matrix_alg, default_autoswitch
+import LinearSolve
+using LinearAlgebra: I, isdiag
+using EnumX
+
+include("default_alg.jl")
+
+export DefaultODEAlgorithm
+
+end # module OrdinaryDiffEqDefault

lib/OrdinaryDiffEqDefault/src/default_alg.jl

@@ -0,0 +1,121 @@
+EnumX.@enumx DefaultSolverChoice begin
+    Tsit5 = 1
+    Vern7 = 2
+    Rosenbrock23 = 3
+    Rodas5P = 4
+    FBDF = 5
+    KrylovFBDF = 6
+end
+
+const NUM_NONSTIFF = 2
+const NUM_STIFF = 4
+const LOW_TOL = 1e-6
+const MED_TOL = 1e-2
+const EXTREME_TOL = 1e-9
+const SMALLSIZE = 50
+const MEDIUMSIZE = 500
+const STABILITY_SIZES = (alg_stability_size(Tsit5()), alg_stability_size(Vern7()))
+const DEFAULTBETA2S = (
+    beta2_default(Tsit5()), beta2_default(Vern7()), beta2_default(Rosenbrock23()),
+    beta2_default(Rodas5P()), beta2_default(FBDF()), beta2_default(FBDF()))
+const DEFAULTBETA1S = (
+    beta1_default(Tsit5(), DEFAULTBETA2S[1]), beta1_default(Vern7(), DEFAULTBETA2S[2]),
+    beta1_default(Rosenbrock23(), DEFAULTBETA2S[3]), beta1_default(
+        Rodas5P(), DEFAULTBETA2S[4]),
+    beta1_default(FBDF(), DEFAULTBETA2S[5]), beta1_default(FBDF(), DEFAULTBETA2S[6]))
+
+callbacks_exists(integrator) = !isempty(integrator.opts.callbacks)
+current_nonstiff(current) = ifelse(current <= NUM_NONSTIFF, current, current - NUM_STIFF)
+
+function DefaultODEAlgorithm(; lazy = true, stiffalgfirst = false, kwargs...)
+    nonstiff = (Tsit5(), Vern7(lazy = lazy))
+    stiff = (Rosenbrock23(; kwargs...), Rodas5P(; kwargs...), FBDF(; kwargs...),
+        FBDF(; linsolve = LinearSolve.KrylovJL_GMRES(), kwargs...))
+    AutoAlgSwitch(nonstiff, stiff; stiffalgfirst)
+end
+
+function is_stiff(integrator, alg, ntol, stol, is_stiffalg, current)
+    eigen_est, dt = integrator.eigen_est, integrator.dt
+    stiffness = abs(eigen_est * dt /
+                    STABILITY_SIZES[nonstiffchoice(integrator.opts.reltol)]) # `abs` here is just for safety
+    tol = is_stiffalg ? stol : ntol
+    os = oneunit(stiffness)
+    bool = stiffness > os * tol
+
+    if !bool
+        integrator.alg.choice_function.successive_switches += 1
+    else
+        integrator.alg.choice_function.successive_switches = 0
+    end
+
+    integrator.do_error_check = (integrator.alg.choice_function.successive_switches >
+                                 integrator.alg.choice_function.switch_max || !bool) ||
+                                is_stiffalg
+    bool
+end
+
+function nonstiffchoice(reltol)
+    x = if reltol < LOW_TOL
+        DefaultSolverChoice.Vern7
+    else
+        DefaultSolverChoice.Tsit5
+    end
+    Int(x)
+end
+
+function stiffchoice(reltol, len, mass_matrix)
+    x = if len > MEDIUMSIZE
+        DefaultSolverChoice.KrylovFBDF
+    elseif len > SMALLSIZE
+        DefaultSolverChoice.FBDF
+    else
+        if reltol < LOW_TOL || !isdiag(mass_matrix)
+            DefaultSolverChoice.Rodas5P
+        else
+            DefaultSolverChoice.Rosenbrock23
+        end
+    end
+    Int(x)
+end
+
+function default_autoswitch(AS::AutoSwitchCache, integrator)
+    len = length(integrator.u)
+    reltol = integrator.opts.reltol
+
+    # Choose the starting method
+    if AS.current == 0
+        choice = if AS.stiffalgfirst || integrator.f.mass_matrix != I
+            stiffchoice(reltol, len, integrator.f.mass_matrix)
+        else
+            nonstiffchoice(reltol)
+        end
+        AS.current = choice
+        return AS.current
+    end
+
+    dt = integrator.dt
+    # Successive stiffness test positives are counted by a positive integer,
+    # and successive stiffness test negatives are counted by a negative integer
+    AS.count = is_stiff(integrator, AS.nonstiffalg, AS.nonstifftol, AS.stifftol,
+        AS.is_stiffalg, AS.current) ?
+               AS.count < 0 ? 1 : AS.count + 1 :
+               AS.count > 0 ? -1 : AS.count - 1
+    if integrator.f.mass_matrix != I
+        #don't change anything
+    elseif (!AS.is_stiffalg && AS.count > AS.maxstiffstep)
+        integrator.dt = dt * AS.dtfac
+        AS.is_stiffalg = true
+        AS.current = stiffchoice(reltol, len, integrator.f.mass_matrix)
+    elseif (AS.is_stiffalg && AS.count < -AS.maxnonstiffstep)
+        integrator.dt = dt / AS.dtfac
+        AS.is_stiffalg = false
+        AS.current = nonstiffchoice(reltol)
+    end
+    return AS.current
+end
+
+# hack for the default alg
+function is_mass_matrix_alg(alg::CompositeAlgorithm{
+    <:Any, <:Tuple{Tsit5, Vern7, Rosenbrock23, Rodas5P, FBDF, FBDF}})
+    true
+end

lib/OrdinaryDiffEqVerner/Project.toml

@@ -0,0 +1,21 @@
+name = "OrdinaryDiffEqVerner"
+uuid = "79d7bb75-1356-48c1-b8c0-6832512096c2"
+authors = ["ParamThakkar123 <paramthakkar864@gmail.com>"]
+version = "1.0.0"
+
+[deps]
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+TruncatedStacktraces = "781d530d-4396-4725-bb49-402e4bee1e77"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+
+[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"
+
+[compat]
+julia = "1.10"
+
+[targets]
+test = ["DiffEqDevTools", "Random", "SafeTestsets", "Test"]

lib/OrdinaryDiffEqVerner/src/OrdinaryDiffEqVerner.jl

@@ -0,0 +1,35 @@
+module OrdinaryDiffEqVerner
+
+import OrdinaryDiffEq: alg_order, calculate_residuals!,
+                       initialize!, perform_step!, @unpack, unwrap_alg,
+                       calculate_residuals, alg_stability_size,
+                       OrdinaryDiffEqAlgorithm,
+                       CompositeAlgorithm, AbstractController, PIDController,
+                       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_interpolant!, _ode_addsteps!, @fold,
+                       @OnDemandTableauExtract, AutoAlgSwitch,
+                       DerivativeOrderNotPossibleError
+using DiffEqBase, FastBroadcast, Polyester, MuladdMacro, RecursiveArrayTools
+using DiffEqBase: @def, @tight_loop_macros
+using Static: False
+using TruncatedStacktraces
+using LinearAlgebra: norm
+
+include("algorithms.jl")
+include("alg_utils.jl")
+include("verner_tableaus.jl")
+include("verner_caches.jl")
+include("verner_addsteps.jl")
+include("interp_func.jl")
+include("interpolants.jl")
+include("controllers.jl")
+include("verner_rk_perform_step.jl")
+
+export Vern6, Vern7, Vern8, Vern9
+export AutoVern6, AutoVern7, AutoVern8, AutoVern9
+
+end

lib/OrdinaryDiffEqVerner/src/alg_utils.jl

@@ -0,0 +1,13 @@
+isfsal(alg::Vern7) = false
+isfsal(alg::Vern8) = false
+isfsal(alg::Vern9) = false
+
+alg_order(alg::Vern6) = 6
+alg_order(alg::Vern7) = 7
+alg_order(alg::Vern8) = 8
+alg_order(alg::Vern9) = 9
+
+alg_stability_size(alg::Vern6) = 4.8553
+alg_stability_size(alg::Vern7) = 4.6400
+alg_stability_size(alg::Vern8) = 5.8641
+alg_stability_size(alg::Vern9) = 4.4762

lib/OrdinaryDiffEqVerner/src/algorithms.jl

@@ -0,0 +1,119 @@
+@doc explicit_rk_docstring(
+    "Verner's “Most Efficient” 6/5 Runge-Kutta method. (lazy 6th order interpolant).",
+    "Vern6",
+    references = "@article{verner2010numerically,
+    title={Numerically optimal Runge--Kutta pairs with interpolants},
+    author={Verner, James H},
+    journal={Numerical Algorithms},
+    volume={53},
+    number={2-3},
+    pages={383--396},
+    year={2010},
+    publisher={Springer}
+    }",
+    extra_keyword_description = """- `lazy`: determines if the lazy interpolant is used.
+                    """,
+    extra_keyword_default = "lazy = true")
+Base.@kwdef struct Vern6{StageLimiter, StepLimiter, Thread} <:
+                   OrdinaryDiffEqAdaptiveAlgorithm
+    stage_limiter!::StageLimiter = trivial_limiter!
+    step_limiter!::StepLimiter = trivial_limiter!
+    thread::Thread = False()
+    lazy::Bool = true
+end
+TruncatedStacktraces.@truncate_stacktrace Vern6 3
+# for backwards compatibility
+function Vern6(stage_limiter!, step_limiter! = trivial_limiter!; lazy = true)
+    Vern6(stage_limiter!, step_limiter!, False(), lazy)
+end
+
+@doc explicit_rk_docstring(
+    "Verner's “Most Efficient” 7/6 Runge-Kutta method. (lazy 7th order interpolant).",
+    "Vern7",
+    references = "@article{verner2010numerically,
+    title={Numerically optimal Runge--Kutta pairs with interpolants},
+    author={Verner, James H},
+    journal={Numerical Algorithms},
+    volume={53},
+    number={2-3},
+    pages={383--396},
+    year={2010},
+    publisher={Springer}
+    }",
+    extra_keyword_description = """- `lazy`: determines if the lazy interpolant is used.
+                    """,
+    extra_keyword_default = "lazy = true")
+Base.@kwdef struct Vern7{StageLimiter, StepLimiter, Thread} <:
+                   OrdinaryDiffEqAdaptiveAlgorithm
+    stage_limiter!::StageLimiter = trivial_limiter!
+    step_limiter!::StepLimiter = trivial_limiter!
+    thread::Thread = False()
+    lazy::Bool = true
+end
+TruncatedStacktraces.@truncate_stacktrace Vern7 3
+# for backwards compatibility
+function Vern7(stage_limiter!, step_limiter! = trivial_limiter!; lazy = true)
+    Vern7(stage_limiter!, step_limiter!, False(), lazy)
+end
+
+@doc explicit_rk_docstring(
+    "Verner's “Most Efficient” 8/7 Runge-Kutta method. (lazy 8th order interpolant).",
+    "Vern8",
+    references = "@article{verner2010numerically,
+    title={Numerically optimal Runge--Kutta pairs with interpolants},
+    author={Verner, James H},
+    journal={Numerical Algorithms},
+    volume={53},
+    number={2-3},
+    pages={383--396},
+    year={2010},
+    publisher={Springer}
+    }",
+    extra_keyword_description = """- `lazy`: determines if the lazy interpolant is used.
+                    """,
+    extra_keyword_default = "lazy = true")
+Base.@kwdef struct Vern8{StageLimiter, StepLimiter, Thread} <:
+                   OrdinaryDiffEqAdaptiveAlgorithm
+    stage_limiter!::StageLimiter = trivial_limiter!
+    step_limiter!::StepLimiter = trivial_limiter!
+    thread::Thread = False()
+    lazy::Bool = true
+end
+TruncatedStacktraces.@truncate_stacktrace Vern8 3
+# for backwards compatibility
+function Vern8(stage_limiter!, step_limiter! = trivial_limiter!; lazy = true)
+    Vern8(stage_limiter!, step_limiter!, False(), lazy)
+end
+
+@doc explicit_rk_docstring(
+    "Verner's “Most Efficient” 9/8 Runge-Kutta method. (lazy9th order interpolant).",
+    "Vern9",
+    references = "@article{verner2010numerically,
+    title={Numerically optimal Runge--Kutta pairs with interpolants},
+    author={Verner, James H},
+    journal={Numerical Algorithms},
+    volume={53},
+    number={2-3},
+    pages={383--396},
+    year={2010},
+    publisher={Springer}
+    }",
+    extra_keyword_description = """- `lazy`: determines if the lazy interpolant is used.
+                    """, extra_keyword_default = "lazy = true")
+Base.@kwdef struct Vern9{StageLimiter, StepLimiter, Thread} <:
+                   OrdinaryDiffEqAdaptiveAlgorithm
+    stage_limiter!::StageLimiter = trivial_limiter!
+    step_limiter!::StepLimiter = trivial_limiter!
+    thread::Thread = False()
+    lazy::Bool = true
+end
+TruncatedStacktraces.@truncate_stacktrace Vern9 3
+# for backwards compatibility
+function Vern9(stage_limiter!, step_limiter! = trivial_limiter!; lazy = true)
+    Vern9(stage_limiter!, step_limiter!, False(), lazy)
+end
+
+AutoVern6(alg; lazy = true, kwargs...) = AutoAlgSwitch(Vern6(lazy = lazy), alg; kwargs...)
+AutoVern7(alg; lazy = true, kwargs...) = AutoAlgSwitch(Vern7(lazy = lazy), alg; kwargs...)
+AutoVern8(alg; lazy = true, kwargs...) = AutoAlgSwitch(Vern8(lazy = lazy), alg; kwargs...)
+AutoVern9(alg; lazy = true, kwargs...) = AutoAlgSwitch(Vern9(lazy = lazy), alg; kwargs...)

lib/OrdinaryDiffEqVerner/src/controllers.jl

@@ -0,0 +1,8 @@
+# checks whether the controller should accept a step based on the error estimate
+@inline function accept_step_controller(integrator, controller::AbstractController)
+    return integrator.EEst <= 1
+end
+
+@inline function accept_step_controller(integrator, controller::PIDController)
+    return integrator.qold >= controller.accept_safety
+end

lib/OrdinaryDiffEqVerner/src/interp_func.jl

@@ -0,0 +1,31 @@
+function DiffEqBase.interp_summary(::Type{cacheType},
+        dense::Bool) where {
+        cacheType <:
+        Union{Vern6Cache, Vern6ConstantCache
+}}
+    dense ? "specialized 6th order lazy interpolation" : "1st order linear"
+end
+
+function DiffEqBase.interp_summary(::Type{cacheType},
+        dense::Bool) where {
+        cacheType <:
+        Union{Vern7Cache, Vern7ConstantCache
+}}
+    dense ? "specialized 7th order lazy interpolation" : "1st order linear"
+end
+
+function DiffEqBase.interp_summary(::Type{cacheType},
+        dense::Bool) where {
+        cacheType <:
+        Union{Vern8Cache, Vern8ConstantCache
+}}
+    dense ? "specialized 8th order lazy interpolation" : "1st order linear"
+end
+
+function DiffEqBase.interp_summary(::Type{cacheType},
+        dense::Bool) where {
+        cacheType <:
+        Union{Vern9Cache, Vern9ConstantCache
+}}
+    dense ? "specialized 9th order lazy interpolation" : "1st order linear"
+end