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

In v7, `qmin` (alongside `qmax`, `gamma`, `beta1/beta2`, `qsteady_*`,
`qoldinit`) moved off `DEOptions` and onto the controller object. The
out-of-domain rejection path in `handle_step_rejection!` was still
reaching for the old `integrator.opts.qmin`, which throws on the v7
`DEOptions` struct — only the legacy DelayDiffEq constructor still
mentions it. The same was true of `integrator.opts.failfactor` in
`post_newton_controller!`.

Reported in
NumericalMathematics/PositiveIntegrators.jl#194 (comment)

Rather than papering over with a one-off `hasfield` walk, this lifts
the standard step-size knobs into a composable building block and
retires the `DummyController` workaround that BDF / Nordsieck were
using to keep the knobs on their algorithm structs.

A new `CommonControllerOptions{T}` struct holds `qmin`, `qmax`,
`qmax_first_step`, `gamma`, `qsteady_min`, `qsteady_max`, `failfactor`
— the seven scalars the integrator-level paths actually read. A single
type parameter `T` keeps the type signatures simple even if more knobs
are added later. All fields default to `nothing`; algorithm-specific
defaults are filled in by `resolve_basic` at `setup_controller_cache`
time. Concrete controllers (`IController`, `PIController`,
`PIDController`, `PredictiveController`, `ExtrapolationController`,
`KantorovichTypeController`, plus the new `BDFController` and
`JVODEController`) all embed a `CommonControllerOptions` as
`controller.basic`.

Seven generic accessors — `get_qmin`, `get_qmax`, `get_qmax_first_step`,
`get_gamma`, `get_qsteady_min`, `get_qsteady_max`, `get_failfactor` —
dispatch on `cache::AbstractControllerCache` and read through
`cache.controller.basic`. `CompositeControllerCache` overrides each one
to delegate to the active sub-cache. `DummyControllerCache` keeps an
alg-field fallback for any SDE algorithm still using it.

`integrator.opts.qmin` → `get_qmin(integrator)`,
`integrator.opts.failfactor` → `get_failfactor(integrator)`. Same in
the BDF post-Newton paths.

QNDF/FBDF/DFBDF used to keep `qmax`, `qsteady_min`, `qsteady_max` as
fields on the algorithm struct itself, with a `DummyController`
hard-wired into `default_controller`. The stepsize logic read
`alg.qmax` / `alg.qsteady_min` / `alg.qsteady_max` directly, plus a
hard-coded `zₛ = 1.2` magic-number gamma, so the controller surface was
unsettable.

`BDFController` embeds `CommonControllerOptions` and has a cache that
delegates back to alg-level dispatch (the existing BDF order-selection
logic is left intact). The hard-coded `zₛ = 1.2` is now
`get_gamma(integrator)`. `default_controller(QT, alg::Union{QNDF, FBDF, DFBDF})`
threads `alg.qmax` / `alg.qsteady_min` / `alg.qsteady_max` through to
the controller, so existing usage like `QNDF(qmax = 20)` keeps working.
Users can now also pass `controller = BDFController(qmin = …, gamma = …)`
to set knobs that previously had no surface (incl. `qmin` and `gamma`).
BDF-tuned per-algorithm defaults (`qmax = 5//1`, `qsteady_min = 9//10`,
`qsteady_max = 12//10`, `gamma = 12//10`) are encoded as
`qmax_default(::QNDF)` / `gamma_default(::QNDF)` etc.

Same pattern. `setη!` / `chooseη!` / `step_accept_controller!(::JVODE)`
now read `get_qmin(integrator)` / `get_qmax(integrator)` /
`get_qsteady_*(integrator)` instead of `alg.qmin` etc.

- `IController` / `PIController` / `PredictiveController` /
  `PIDController` shed their flat `qmin/qmax/...` fields and embed
  `CommonControllerOptions` instead. PI-specific knobs (`beta1`,
  `beta2`, `qoldinit`) and PID-specific knobs (`beta`, `accept_safety`,
  `limiter`) stay on the controller alongside `basic`.
- `ExtrapolationController` and `KantorovichTypeController` likewise
  embed `CommonControllerOptions`. Their stepsize logic reads
  `get_qmax(integrator)` / `get_qmin(integrator)` rather than direct
  field access.

`test/runtests.jl` walks transitive `[sources]` dependencies and
`Pkg.develop`s them. Pre-seed the `developed` set with the active
project so a `[sources]` entry that points back to it (e.g. via the
umbrella `OrdinaryDiffEq`'s transitive sources) is skipped — `Pkg.develop`
cannot develop the active project itself, and that error was the
"package X has the same name or UUID as the active project" failure
across the sublibrary CI matrix.

Reproducer (`isoutofdomain` predicate that fires once on the first
proposed step) plus a smoke test of every controller path (default
`solve`, user-supplied `BDFController`, `CommonControllerOptions`
construction, controller-composition invariants) — 21/21 pass on Julia
1.12.

- On master (without this fix): all algorithms error out — accessing
  `integrator.opts.qmin` throws because the v7 `DEOptions` struct
  doesn't have the field.
- With this fix: `Tsit5` / `Vern7` / `Rosenbrock23` / `FBDF` / `QNDF`
  all complete the isout-rejection problem successfully, and
  `BDFController(qmax = 3)` is honored end-to-end.

Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>

Author: ChrisRackauckas

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 <paramthakkar864@gmail.com>"]
-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