Skip to content

Commit 2b6b3d7

Browse files
Add get_qmin to the controller interface
`qmin` (alongside `qmax`, `gamma`, `beta1/beta2`, `qsteady_*`, `qoldinit`) moved off `DEOptions` and onto the controller object in v7. 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. Surface this through the controller interface rather than a one-off `hasfield` walk: a new `get_qmin(integrator)` (with implementations dispatched on each concrete controller cache) returns the minimum step-size shrinkage factor used by the integrator's controller. `handle_step_rejection!` calls it instead of `integrator.opts.qmin`. Per-cache implementations: - `IControllerCache`, `PIControllerCache`, `PredictiveControllerCache` return their controller's `qmin` field. - `PIDControllerCache` has no `qmin` (it limits dt via `limiter` + `accept_safety` instead) and returns the historical default `1 // 5` so the out-of-domain shrink path still has something to multiply by. - `DummyControllerCache` (BDF / Nordsieck — algorithms that own the step-size logic) reads `integrator.alg.qmin` if present, else falls back to the historical default. - `CompositeControllerCache` delegates to the currently active sub-cache, mirroring how `stepsize_controller!` and friends dispatch. Refresh the `PredictiveController` docstring (which still showed the old `integrator.opts.qmin/qmax/qsteady_*/gamma` interface) to match the actual v7 implementation that destructures from `cache.controller`, and update two stale `# equivalent to integrator.opts.gamma` comments in `lib/OrdinaryDiffEqBDF/src/controllers.jl`. Reported in NumericalMathematics/PositiveIntegrators.jl#194 (comment) Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
1 parent 3a72e5e commit 2b6b3d7

3 files changed

Lines changed: 77 additions & 25 deletions

File tree

lib/OrdinaryDiffEqBDF/src/controllers.jl

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -26,7 +26,7 @@ function step_accept_controller!(integrator, cache::Union{QNDFCache, QNDFConstan
2626
h = integrator.dt
2727
k = cache.order
2828
prefer_const_step = cache.nconsteps < cache.order + 2
29-
zₛ = 1.2 # equivalent to integrator.opts.gamma
29+
zₛ = 1.2 # equivalent to controller `gamma`
3030
zᵤ = 0.1
3131
Fᵤ = 10
3232
expo = 1 / (k + 1)
@@ -100,7 +100,7 @@ function bdf_step_reject_controller!(integrator, cache, EEst1)
100100
if cache.consfailcnt > 1
101101
h = h / 2
102102
end
103-
zₛ = 1.2 # equivalent to integrator.opts.gamma
103+
zₛ = 1.2 # equivalent to controller `gamma`
104104
expo = 1 / (k + 1)
105105
z = zₛ * ((OrdinaryDiffEqCore.get_EEst(integrator))^expo)
106106
F = inv(z)

lib/OrdinaryDiffEqCore/src/integrators/controllers.jl

Lines changed: 74 additions & 22 deletions
Original file line numberDiff line numberDiff line change
@@ -157,6 +157,35 @@ See also: https://github.com/SciML/DifferentialEquations.jl/issues/299
157157
return qmax
158158
end
159159

160+
"""
161+
get_qmin(integrator)
162+
get_qmin(integrator, controller_cache)
163+
164+
Return the minimum step-size shrinkage factor `qmin` for the integrator's
165+
controller. Used by the out-of-domain rejection path in
166+
`handle_step_rejection!` to retry with a smaller step after the
167+
`isoutofdomain` predicate fires.
168+
169+
Each concrete `AbstractControllerCache` implements its own
170+
`get_qmin(integrator, cache)` method:
171+
172+
- `IController`, `PIController`, `PredictiveController` return their
173+
`qmin` field.
174+
- `PIDController` has no `qmin` (it uses `limiter` + `accept_safety`
175+
instead) and returns the historical default `1 // 5` so the
176+
out-of-domain path still has something to shrink by.
177+
- `DummyController` (BDF / Nordsieck — algorithms that own the
178+
step-size logic) reads `integrator.alg.qmin` if present, else falls
179+
back to the historical default.
180+
- `CompositeController` delegates to the currently active sub-cache.
181+
"""
182+
@inline get_qmin(integrator::SciMLBase.DEIntegrator) =
183+
get_qmin(integrator, integrator.controller_cache)
184+
185+
# Default fallback for any AbstractControllerCache that hasn't overridden
186+
# this — used by the historical step-size shrink behavior.
187+
@inline get_qmin(integrator, ::AbstractControllerCache) = 1 // 5
188+
160189
reset_alg_dependent_opts!(controller::AbstractControllerCache, alg1, alg2) = nothing
161190

162191
reinit_controller!(integrator::SciMLBase.DEIntegrator, controller::AbstractControllerCache) = nothing
@@ -210,6 +239,10 @@ end
210239
post_newton_controller!(integrator, alg)
211240
@inline accept_step_controller(integrator, cache::DummyControllerCache, alg) =
212241
get_EEst(cache) <= 1
242+
# BDF / Nordsieck-style algorithms with integrated controllers store the
243+
# step-size knobs on the algorithm itself.
244+
@inline get_qmin(integrator, ::DummyControllerCache) =
245+
hasfield(typeof(integrator.alg), :qmin) ? integrator.alg.qmin : 1 // 5
213246

214247

215248
# Standard integral (I) step size controller
@@ -323,6 +356,8 @@ function step_reject_controller!(integrator, cache::IControllerCache, alg)
323356
return integrator.dt = cache.dtreject
324357
end
325358

359+
@inline get_qmin(integrator, cache::IControllerCache) = cache.controller.qmin
360+
326361
reinit_controller!(integrator::SciMLBase.DEIntegrator, cache::IControllerCache) = nothing
327362

328363
function sync_controllers!(cache1::IControllerCache, cache2::IControllerCache)
@@ -468,6 +503,8 @@ function step_reject_controller!(integrator, cache::PIControllerCache, alg)
468503
return integrator.dt /= min(inv(qmin), q11 / gamma)
469504
end
470505

506+
@inline get_qmin(integrator, cache::PIControllerCache) = cache.controller.qmin
507+
471508
function reinit_controller!(integrator::SciMLBase.DEIntegrator, cache::PIControllerCache{T}) where {T}
472509
cache.q11 = one(T)
473510
return cache.errold = T(cache.controller.qoldinit)
@@ -686,6 +723,10 @@ function step_reject_controller!(integrator, cache::PIDControllerCache, alg)
686723
return integrator.dt *= cache.dt_factor
687724
end
688725

726+
# PIDController has no `qmin`; it limits dt via `limiter` + `accept_safety`.
727+
# Fall back to the historical default for the out-of-domain rejection path.
728+
@inline get_qmin(integrator, ::PIDControllerCache) = 1 // 5
729+
689730
function sync_controllers!(cache1::PIDControllerCache, cache2::PIDControllerCache)
690731
cache1.err = cache2.err
691732
cache1.dt_factor = cache2.dt_factor
@@ -702,38 +743,42 @@ well-suited for stiff solvers where this can be expected, and is the default
702743
for algorithms like the (E)SDIRK methods.
703744
704745
```julia
705-
gamma = integrator.opts.gamma
706-
niters = integrator.cache.newton_iters
746+
(; qmin, qmax, gamma) = controller
747+
qmax = get_current_qmax(integrator, qmax)
748+
niters = integrator.cache.nlsolver.iter
707749
fac = min(gamma,
708-
(1 + 2 * integrator.alg.max_newton_iter) * gamma /
709-
(niters + 2 * integrator.alg.max_newton_iter))
710-
expo = 1 / (alg_order(integrator.alg) + 1)
711-
qtmp = (get_EEst(integrator)^expo) / fac
712-
@fastmath q = max(inv(integrator.opts.qmax), min(inv(integrator.opts.qmin), qtmp))
713-
if q <= integrator.opts.qsteady_max && q >= integrator.opts.qsteady_min
714-
q = one(q)
715-
end
716-
integrator.controller_cache.qold = q
750+
(1 + 2 * integrator.cache.nlsolver.maxiters) * gamma /
751+
(niters + 2 * integrator.cache.nlsolver.maxiters))
752+
expo = 1 / (get_current_adaptive_order(alg, integrator.cache) + 1)
753+
qtmp = fastpower(get_EEst(integrator), expo) / fac
754+
@fastmath q = max(inv(qmax), min(inv(qmin), qtmp))
755+
cache.qold = q
717756
q
718757
```
719758
720-
In this case, `niters` is the number of Newton iterations which was required
721-
in the most recent step of the algorithm. Note that these values are used
722-
differently depending on acceptance and rejectance. When the step is accepted,
723-
the following logic is applied:
759+
where `controller` and `cache` are the controller and its cache stored in
760+
`integrator.controller_cache`. `niters` is the number of Newton iterations
761+
which was required in the most recent step of the algorithm. Note that
762+
these values are used differently depending on acceptance and rejectance.
763+
When the step is accepted, the following logic is applied:
724764
725765
```julia
766+
(; dtacc, erracc) = cache
767+
(; qmin, qmax, gamma, qsteady_min, qsteady_max) = controller
768+
qmax = get_current_qmax(integrator, qmax)
726769
if integrator.success_iter > 0
727-
expo = 1 / (alg_adaptive_order(integrator.alg) + 1)
728-
qgus = (integrator.dtacc / integrator.dt) * (((get_EEst(integrator)^2) / integrator.erracc)^expo)
729-
qgus = max(inv(integrator.opts.qmax),
730-
min(inv(integrator.opts.qmin), qgus / integrator.opts.gamma))
770+
expo = 1 / (get_current_adaptive_order(alg, integrator.cache) + 1)
771+
qgus = (dtacc / integrator.dt) * fastpower((get_EEst(integrator)^2) / erracc, expo)
772+
qgus = max(inv(qmax), min(inv(qmin), qgus / gamma))
731773
qacc = max(q, qgus)
732774
else
733775
qacc = q
734776
end
735-
integrator.dtacc = integrator.dt
736-
integrator.erracc = max(1e-2, get_EEst(integrator))
777+
if qsteady_min <= qacc <= qsteady_max
778+
qacc = one(qacc)
779+
end
780+
cache.dtacc = integrator.dt
781+
cache.erracc = max(1e-2, get_EEst(integrator))
737782
integrator.dt / qacc
738783
```
739784
@@ -743,7 +788,7 @@ When it rejects, it's the same as the [`IController`](@ref):
743788
if integrator.success_iter == 0
744789
integrator.dt *= 0.1
745790
else
746-
integrator.dt = integrator.dt / integrator.controller_cache.qold
791+
integrator.dt = integrator.dt / cache.qold
747792
end
748793
```
749794
"""
@@ -871,6 +916,8 @@ function step_reject_controller!(integrator, cache::PredictiveControllerCache, a
871916
return integrator.dt = success_iter == 0 ? 0.1 * dt : dt / qold
872917
end
873918

919+
@inline get_qmin(integrator, cache::PredictiveControllerCache) = cache.controller.qmin
920+
874921
# For a composite algorithm the default strategy is to switch forth and back between the controllers of the individual algorithms
875922
struct CompositeController{T} <: AbstractController
876923
controllers::T
@@ -931,6 +978,11 @@ end
931978
return post_newton_controller!(integrator, @inbounds(cache.caches[current_idx]), alg)
932979
end
933980

981+
@inline function get_qmin(integrator, cache::CompositeControllerCache)
982+
current_idx = integrator.cache.current
983+
return get_qmin(integrator, @inbounds(cache.caches[current_idx]))
984+
end
985+
934986
function setup_controller_cache(alg::CompositeAlgorithm, caches::DefaultCache, controller::CompositeController, ::Type{E}) where {E}
935987
sub = (
936988
setup_controller_cache(alg.algs[1], caches, controller.controllers[1], E),

lib/OrdinaryDiffEqCore/src/integrators/integrator_utils.jl

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -68,7 +68,7 @@ function handle_step_rejection!(integrator)
6868
integrator.opts.verbose, :step_rejected
6969
)
7070
if integrator.isout
71-
integrator.dt = integrator.dt * integrator.opts.qmin
71+
integrator.dt = integrator.dt * get_qmin(integrator)
7272
elseif !integrator.force_stepfail
7373
step_reject_controller!(integrator, integrator.alg)
7474
end

0 commit comments

Comments
 (0)