Compose step-size knobs into BasicController; replace DummyController for BDF/JVODE

lib/ImplicitDiscreteSolve/Project.toml

@@ -29,7 +29,7 @@ NonlinearSolveBase = "2.0.0"
 OrdinaryDiffEqSDIRK = "2"
 SafeTestsets = "0.1.0"
 SciMLBase = "3"
-OrdinaryDiffEqCore = "4"
+OrdinaryDiffEqCore = "4.1"
 ConcreteStructs = "0.2.3"
 NonlinearSolveFirstOrder = "1.9.0, 2"
 SymbolicIndexingInterface = "0.3.38"

lib/ImplicitDiscreteSolve/src/controller.jl

@@ -25,33 +25,48 @@ The baseline algorithm has been derived in Peter Deuflhard's book "Newton Method
 Nonlinear Problems" in Section 5.1.3 (Adaptive pathfollowing algorithms). Please note
 that some implementation details deviate from the original algorithm.
 """
-Base.@kwdef struct KantorovichTypeController{T} <: AbstractController
+struct KantorovichTypeController{B <: OrdinaryDiffEqCore.CommonControllerOptions, T} <: AbstractController
+    basic::B
     Θmin::T
     p::Int64
-    Θreject::T = 0.95
-    Θbar::T = 0.5
-    γ::T = 0.95
-    qmin::T = 1 / 5
-    qmax::T = 5.0
-    strict::Bool = true
+    Θreject::T
+    Θbar::T
+    γ::T
+    strict::Bool
+end
+
+function KantorovichTypeController(;
+        Θmin, p, Θreject = 0.95, Θbar = 0.5, γ = 0.95,
+        qmin = 1 // 5, qmax = 5, strict = true,
+        kwargs...,
+    )
+    T = promote_type(typeof(Θmin), typeof(Θreject), typeof(Θbar), typeof(γ))
+    basic = OrdinaryDiffEqCore.CommonControllerOptions(; qmin, qmax, kwargs...)
+    return KantorovichTypeController{typeof(basic), T}(
+        basic, T(Θmin), Int64(p), T(Θreject), T(Θbar), T(γ), strict,
+    )
 end
 
 mutable struct KantorovichTypeControllerCache{T, E} <: AbstractControllerCache
-    controller::KantorovichTypeController{T}
+    controller::KantorovichTypeController{OrdinaryDiffEqCore.CommonControllerOptions{T}, T}
     EEst::E
 end
 
 function OrdinaryDiffEqCore.default_controller(
         QT, alg::IDSolve,
     )
-    return KantorovichTypeController{QT}(; Θmin = QT(1 // 8), p = 1)
+    return KantorovichTypeController(; Θmin = QT(1 // 8), p = 1)
 end
 
-function OrdinaryDiffEqCore.setup_controller_cache(alg, cache, controller::KantorovichTypeController{T}, ::Type{E}) where {T, E}
-    return KantorovichTypeControllerCache{T, E}(
-        controller,
-        oneunit(E),
+function OrdinaryDiffEqCore.setup_controller_cache(alg, cache, controller::KantorovichTypeController, ::Type{E}) where {E}
+    QT = OrdinaryDiffEqCore._resolved_QT(controller.basic)
+    basic = OrdinaryDiffEqCore.resolve_basic(controller.basic, alg, QT)
+    resolved = KantorovichTypeController{typeof(basic), QT}(
+        basic, QT(controller.Θmin), controller.p,
+        QT(controller.Θreject), QT(controller.Θbar), QT(controller.γ), controller.strict,
     )
+    T = QT
+    return KantorovichTypeControllerCache{T, E}(resolved, oneunit(E))
 end
 
 function OrdinaryDiffEqCore.stepsize_controller!(
@@ -62,7 +77,8 @@ function OrdinaryDiffEqCore.stepsize_controller!(
 
     # Adapt dt with a priori estimate (Eq. 5.24)
     (; Θks) = integrator.cache
-    (; Θbar, γ, Θmin, qmin, qmax, p) = controller
+    (; Θbar, γ, Θmin, p) = controller
+    (; qmin, qmax) = controller.basic
 
     Θ₀ = length(Θks) > 0 ? max(first(Θks), Θmin) : Θmin
     q = clamp(γ * (g(Θbar) / (g(Θ₀)))^(1 / p), qmin, qmax)
@@ -84,7 +100,8 @@ function OrdinaryDiffEqCore.step_reject_controller!(
 
     # Shorten dt according to (Eq. 5.24)
     (; Θks) = integrator.cache
-    (; Θbar, Θreject, γ, Θmin, qmin, qmax, p) = controller
+    (; Θbar, Θreject, γ, Θmin, p) = controller
+    (; qmin, qmax) = controller.basic
     for Θk in Θks
         if Θk > Θreject
             q = clamp(γ * (g(Θbar) / g(Θk))^(1 / p), qmin, qmax)

lib/OrdinaryDiffEqBDF/Project.toml

@@ -61,7 +61,7 @@ OrdinaryDiffEqDifferentiation = "3"
 OrdinaryDiffEqSDIRK = "2"
 TruncatedStacktraces = "1.4"
 SciMLBase = "3"
-OrdinaryDiffEqCore = "4"
+OrdinaryDiffEqCore = "4.1"
 ArrayInterface = "7.19"
 Enzyme = "0.13"
 Preferences = "1.4"

lib/OrdinaryDiffEqBDF/src/OrdinaryDiffEqBDF.jl

@@ -14,11 +14,16 @@ import OrdinaryDiffEqCore: alg_order, calculate_residuals!,
     constvalue, isadaptive, error_constant,
     has_special_newton_error,
     trivial_limiter!,
-    issplit, qsteady_min_default, qsteady_max_default,
+    issplit, qmin_default, qmax_default, gamma_default,
+    qsteady_min_default, qsteady_max_default,
     get_current_alg_order, get_current_adaptive_order,
     stepsize_controller!,
     step_accept_controller!,
     step_reject_controller!, post_newton_controller!,
+    accept_step_controller, get_EEst, set_EEst!,
+    setup_controller_cache, get_qmax, get_gamma, get_qsteady_min, get_qsteady_max,
+    get_failfactor, CommonControllerOptions, resolve_basic, _resolved_QT,
+    AbstractControllerCache,
     DAEAlgorithm, _unwrap_val, DummyController,
     get_fsalfirstlast, generic_solver_docstring, _ad_chunksize_int, _ad_fdtype, _fixup_ad,
     _ode_interpolant, _ode_interpolant!, has_stiff_interpolation,

lib/OrdinaryDiffEqBDF/src/controllers.jl

@@ -1,5 +1,72 @@
+"""
+    BDFController(; qmin, qmax, qsteady_min, qsteady_max, gamma, qmax_first_step,
+                  failfactor)
+
+Step-size controller for the variable-order BDF family (`QNDF`, `FBDF`,
+`DFBDF`). Composes the standard step-size knobs via [`CommonControllerOptions`](@ref);
+the adaptive logic is integrated into the algorithm itself, so the cache
+falls through to alg-level dispatch (the same way the legacy
+`DummyControllerCache` did) but exposes the knobs as a real, settable
+controller. Pass it explicitly with `solve(prob, alg; controller = BDFController(...))`,
+or rely on the default constructed by `default_controller(QT, alg)`.
+"""
+struct BDFController{B <: CommonControllerOptions} <: AbstractController
+    basic::B
+end
+
+BDFController(; kwargs...) = BDFController(CommonControllerOptions(; kwargs...))
+
+BDFController(alg; kwargs...) = BDFController(Float64, alg; kwargs...)
+BDFController(::Type{QT}, alg; kwargs...) where {QT} =
+    BDFController(resolve_basic(CommonControllerOptions(; kwargs...), alg, QT))
+
+mutable struct BDFControllerCache{T, E, C} <: AbstractControllerCache
+    controller::BDFController{CommonControllerOptions{T}}
+    cache::C
+    EEst::E
+end
+
+function setup_controller_cache(
+        alg::Union{QNDF, FBDF, DFBDF}, cache, controller::BDFController, ::Type{E},
+    ) where {E}
+    QT = _resolved_QT(controller.basic)
+    basic = resolve_basic(controller.basic, alg, QT)
+    resolved = BDFController(basic)
+    return BDFControllerCache{QT, E, typeof(cache)}(resolved, cache, oneunit(E))
+end
+
+# The BDF stepsize/accept/reject logic lives at the algorithm level — the
+# controller cache just delegates back, mirroring how DummyControllerCache did.
+@inline OrdinaryDiffEqCore.stepsize_controller!(integrator, ::BDFControllerCache, alg) =
+    stepsize_controller!(integrator, alg)
+@inline OrdinaryDiffEqCore.step_accept_controller!(integrator, ::BDFControllerCache, alg, q) =
+    step_accept_controller!(integrator, alg, q)
+@inline OrdinaryDiffEqCore.step_reject_controller!(integrator, ::BDFControllerCache, alg) =
+    step_reject_controller!(integrator, alg)
+@inline OrdinaryDiffEqCore.post_newton_controller!(integrator, ::BDFControllerCache, alg) =
+    post_newton_controller!(integrator, alg)
+@inline OrdinaryDiffEqCore.accept_step_controller(
+    integrator, cache::BDFControllerCache, alg,
+) = get_EEst(cache) <= 1
+
+# Per-algorithm defaults for `CommonControllerOptions` knobs that don't match the
+# generic IController defaults. These match the historical alg-struct kwargs
+# `QNDF(qmax=5//1, qsteady_min=9//10, qsteady_max=12//10)` etc. and the formerly
+# hard-coded `zₛ = 1.2` step-size safety factor in the BDF stepsize logic.
+qmax_default(::Union{QNDF, FBDF, DFBDF}) = 5 // 1
+qsteady_min_default(::Union{QNDF, QNDF1, QNDF2, FBDF, DFBDF}) = 9 // 10
+qsteady_max_default(::Union{QNDF, QNDF1, QNDF2, FBDF, DFBDF}) = 12 // 10
+gamma_default(::Union{QNDF, FBDF, DFBDF}) = 12 // 10
+
 function default_controller(QT, alg::Union{QNDF, FBDF, DFBDF})
-    return DummyController()
+    # Thread the alg-level kwargs through to the CommonControllerOptions so that
+    # `QNDF(qmax = 20)` keeps working (qmax = 20 ends up on the controller).
+    return BDFController(
+        QT, alg;
+        qmax = alg.qmax,
+        qsteady_min = alg.qsteady_min,
+        qsteady_max = alg.qsteady_max,
+    )
 end
 
 # QNDF
@@ -18,15 +85,15 @@ function step_accept_controller!(integrator, cache::Union{QNDFCache, QNDFConstan
     cache.consfailcnt = 0
     cache.nconsteps += 1
     if iszero(OrdinaryDiffEqCore.get_EEst(integrator))
-        return integrator.dt * get_current_qmax(integrator, alg.qmax)
+        return integrator.dt * get_current_qmax(integrator, get_qmax(integrator))
     else
         est = OrdinaryDiffEqCore.get_EEst(integrator)
         estₖ₋₁ = cache.EEst1
         estₖ₊₁ = cache.EEst2
         h = integrator.dt
         k = cache.order
         prefer_const_step = cache.nconsteps < cache.order + 2
-        zₛ = 1.2 # equivalent to integrator.opts.gamma
+        zₛ = get_gamma(integrator)
         zᵤ = 0.1
         Fᵤ = 10
         expo = 1 / (k + 1)
@@ -86,7 +153,7 @@ function step_accept_controller!(integrator, cache::Union{QNDFCache, QNDFConstan
             return integrator.dt
         end
     end
-    if q <= alg.qsteady_max && q >= alg.qsteady_min
+    if q <= get_qsteady_max(integrator) && q >= get_qsteady_min(integrator)
         return integrator.dt
     end
     return integrator.dt / q
@@ -100,7 +167,7 @@ function bdf_step_reject_controller!(integrator, cache, EEst1)
     if cache.consfailcnt > 1
         h = h / 2
     end
-    zₛ = 1.2  # equivalent to integrator.opts.gamma
+    zₛ = get_gamma(integrator)
     expo = 1 / (k + 1)
     z = zₛ * ((OrdinaryDiffEqCore.get_EEst(integrator))^expo)
     F = inv(z)
@@ -154,7 +221,7 @@ function post_newton_controller!(integrator, cache::Union{FBDFCache, FBDFConstan
     if cache.order > 1 && cache.nlsolver.nfails >= 3
         cache.order -= 1
     end
-    integrator.dt = integrator.dt / integrator.opts.failfactor
+    integrator.dt = integrator.dt / get_failfactor(integrator)
     cache.consfailcnt += 1
     cache.nconsteps = 0
     return nothing
@@ -299,7 +366,7 @@ function stepsize_controller!(
     end
 
     if iszero(terk)
-        q = inv(get_current_qmax(integrator, alg.qmax))
+        q = inv(get_current_qmax(integrator, get_qmax(integrator)))
     else
         # CVODE-style step size formula: eta = 1 / (BIAS2 * dsm)^(1/(k+1))
         # where dsm = terk / (alpha0 * (k+1)) and alpha0 is the BDF leading coefficient.
@@ -320,7 +387,7 @@ function step_accept_controller!(
         q
     ) where {max_order}
     cache.consfailcnt = 0
-    if q <= alg.qsteady_max && q >= alg.qsteady_min
+    if q <= get_qsteady_max(integrator) && q >= get_qsteady_min(integrator)
         q = one(q)
     end
     cache.nconsteps += 1
@@ -348,7 +415,7 @@ function post_newton_controller!(integrator, cache::Union{DFBDFCache, DFBDFConst
     if cache.order > 1 && cache.nlsolver.nfails >= 3
         cache.order -= 1
     end
-    integrator.dt = integrator.dt / integrator.opts.failfactor
+    integrator.dt = integrator.dt / get_failfactor(integrator)
     cache.consfailcnt += 1
     cache.nconsteps = 0
     return nothing
@@ -465,7 +532,7 @@ function stepsize_controller!(
         cache.order = k
     end
     if iszero(terk)
-        q = inv(get_current_qmax(integrator, alg.qmax))
+        q = inv(get_current_qmax(integrator, get_qmax(integrator)))
     else
         # CVODE-style step size formula matching FBDF change
         alpha0 = cache.bdf_coeffs[k, 1]
@@ -483,7 +550,7 @@ function step_accept_controller!(
         q
     ) where {max_order}
     cache.consfailcnt = 0
-    if q <= alg.qsteady_max && q >= alg.qsteady_min
+    if q <= get_qsteady_max(integrator) && q >= get_qsteady_min(integrator)
         q = one(q)
     end
     cache.nconsteps += 1

lib/OrdinaryDiffEqCore/Project.toml

@@ -1,7 +1,7 @@
 name = "OrdinaryDiffEqCore"
 uuid = "bbf590c4-e513-4bbe-9b18-05decba2e5d8"
 authors = ["ParamThakkar123 <[email protected]>"]
-version = "4.0.2"
+version = "4.1.0"
 
 [deps]
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

lib/OrdinaryDiffEqCore/src/alg_utils.jl

@@ -468,6 +468,16 @@ qsteady_max_default(alg) = 1
 qsteady_max_default(alg::OrdinaryDiffEqAdaptiveImplicitAlgorithm) = 6 // 5
 # But don't re-use Jacobian if not adaptive: too risky and cannot pull back
 qsteady_max_default(alg::OrdinaryDiffEqImplicitAlgorithm) = isadaptive(alg) ? 1 // 1 : 0
+
+# qmax_first_step is the upper bound on dt growth on the very first accepted
+# step — see https://github.com/SciML/DifferentialEquations.jl/issues/299.
+# 10000 mirrors the historical Sundials CVODE behavior.
+qmax_first_step_default(alg) = 10000
+
+# failfactor is the post-Newton-failure dt shrink factor used by
+# `post_newton_controller!`.  Default of 2 matches the historical
+# `integrator.opts.failfactor` default.
+failfactor_default(alg) = 2
 #TODO
 #SciMLBase.nlsolve_default(::QNDF, ::Val{κ}) = 1//2