Skip to content

remove redundant controller divisions #3745

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
oscardssmith wants to merge 3 commits into
base: master
Choose a base branch
from
Open

remove redundant controller divisions #3745

Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
7a07973
remove redundant controller divisions
oscardssmith Jun 11, 2026
7b6d392
fixes
oscardssmith Jun 11, 2026
13745cb
fix incorrect comment
oscardssmith Jun 11, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
85 changes: 42 additions & 43 deletions lib/OrdinaryDiffEqCore/src/integrators/controllers.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -554,25 +554,25 @@
EEst = DiffEqBase.value(get_EEst(integrator))

if iszero(EEst)
q = inv(qmax)
dt_factor = qmax
else
expo = 1 / (get_current_adaptive_order(alg, integrator.cache) + 1)
qtmp = fastpower(EEst, expo) / gamma
@fastmath q = DiffEqBase.value(max(inv(qmax), min(inv(qmin), qtmp)))
expo = inv(get_current_adaptive_order(alg, integrator.cache) + 1)
dt_factor = gamma / fastpower(EEst, expo)
@fastmath dt_factor = DiffEqBase.value(clamp(dt_factor, qmin, qmax))
# TODO: Shouldn't this be in `step_accept_controller!` as for the PI controller?
cache.dtreject = DiffEqBase.value(integrator.dt) / q
cache.dtreject = DiffEqBase.value(integrator.dt) * dt_factor
end
return q
return dt_factor
end

# TODO change signature to remove the q input
function step_accept_controller!(integrator, cache::IControllerCache, alg, q)
# TODO change signature to remove the dt_factor input
function step_accept_controller!(integrator, cache::IControllerCache, alg, dt_factor)
(; qsteady_min, qsteady_max) = cache.controller.basic

if qsteady_min <= q <= qsteady_max
q = one(q)
if qsteady_min <= inv(dt_factor) <= qsteady_max
dt_factor = one(dt_factor)
end
return integrator.dt / q # new dt
return integrator.dt * dt_factor # new dt
end

function step_reject_controller!(integrator, cache::IControllerCache, alg)
Expand Down Expand Up @@ -696,27 +696,27 @@
EEst = DiffEqBase.value(get_EEst(integrator))

if iszero(EEst)
q = inv(qmax)
dt_factor = qmax
else
q11 = fastpower(EEst, beta1)
q = q11 / fastpower(errold, beta2)
cache.q11 = q11
@fastmath q = clamp(q / gamma, inv(qmax), inv(qmin))
@fastmath dt_factor = clamp(gamma / q, qmin, qmax)
end
return q
return dt_factor
end

function step_accept_controller!(integrator, cache::PIControllerCache, alg, q)
function step_accept_controller!(integrator, cache::PIControllerCache, alg, dt_factor)
(; controller) = cache
(; qsteady_min, qsteady_max) = controller.basic
qoldinit = controller.qoldinit
EEst = DiffEqBase.value(get_EEst(integrator))

if qsteady_min <= q <= qsteady_max
q = one(q)
if qsteady_min <= inv(dt_factor) <= qsteady_max
dt_factor = one(dt_factor)
end
cache.errold = max(EEst, qoldinit)
return integrator.dt / q # new dt
return integrator.dt * dt_factor # new dt
end

function step_reject_controller!(integrator, cache::PIControllerCache, alg)
Expand All @@ -729,7 +729,7 @@
return integrator.dt
end
end
return integrator.dt /= min(inv(qmin), q11 / gamma)
return integrator.dt *= max(qmin, gamma / q11)
end

function reinit_controller!(integrator::SciMLBase.DEIntegrator, cache::PIControllerCache{T}) where {T}
Expand Down Expand Up @@ -924,14 +924,14 @@
err1, err2, err3 = cache.err

k = min(alg_order(alg), alg_adaptive_order(alg)) + 1
dt_factor = err1^(beta1 / k) * err2^(beta2 / k) * err3^(beta3 / k)
dt_factor = fastpower(err1, beta1 / k) * fastpower(err2, beta2 / k) * fastpower(err3, beta3 / k)
if isnan(dt_factor)
@warn "unlimited dt_factor" dt_factor err1 err2 err3 beta1 beta2 beta3 k
end
cache.dt_factor = controller.limiter(dt_factor)

# Note: No additional limiting of the form
# dt_factor = max(qmin, min(qmax, dt_factor))
# dt_factor = clamp(dt_factor, qmin, qmax)
# is necessary since the `limiter` should take care of that. The default limiter
# ensures
# 0.21 ≈ limiter(0) <= dt_factor <= limiter(Inf) ≈ 2.57
Expand Down Expand Up @@ -980,9 +980,9 @@
PredictiveController()

The Gustafsson acceleration algorithm accelerates changes so that way algorithms
can more swiftly change to handle quick transients. This algorithm is thus
well-suited for stiff solvers where this can be expected, and is the default
for algorithms like the (E)SDIRK methods.
can more swiftly change to handle quick transients.
This algorithm is gennerally less efficient than a PI controller, but is used in adaptive order methods
(e.g. AdaptiveRadau) where a simpler stepsize controller makes an order controller eaiser.

Check warning on line 985 in lib/OrdinaryDiffEqCore/src/integrators/controllers.jl

View workflow job for this annotation

GitHub Actions / Spell Check with Typos / Spell Check with Typos 

"eaiser" should be "easier".
Comment on lines +984 to +985

@ChrisRackauckas ChrisRackauckas Jun 12, 2026

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That's not true though? This isn't a Proportional Controller, it's a PredictiveController which is different.

@oscardssmith oscardssmith Jun 12, 2026

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

PI is our default controller though, it's just the 2nd part that is wrong, right?

@ChrisRackauckas ChrisRackauckas Jun 12, 2026

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It's not generally less efficient for all solvers, and it's not used with AdaptiveRadau

@oscardssmith oscardssmith Jun 12, 2026

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

what solvers is it more efficient for? I don't think we use it by default for any of them.

@ChrisRackauckas ChrisRackauckas Jun 12, 2026

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

we probably should use them on some sdirks

@oscardssmith oscardssmith Jun 12, 2026

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

which/why?

@ChrisRackauckas ChrisRackauckas Jun 12, 2026

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It attempts to predict Newton instabilities, so possibly all. It's just an underexplored topic.


```julia
(; qmin, qmax, gamma) = controller
Expand All @@ -992,8 +992,8 @@
(1 + 2 * integrator.cache.nlsolver.maxiters) * gamma /
(niters + 2 * integrator.cache.nlsolver.maxiters))
expo = 1 / (get_current_adaptive_order(alg, integrator.cache) + 1)
qtmp = fastpower(get_EEst(integrator), expo) / fac
@fastmath q = max(inv(qmax), min(inv(qmin), qtmp))
qtmp = fac / fastpower(get_EEst(integrator), expo)
@fastmath q = clamp(qtmp, qmin, qmax)
cache.qold = q
q
```
Expand All @@ -1011,17 +1011,17 @@
if integrator.success_iter > 0
expo = 1 / (get_current_adaptive_order(alg, integrator.cache) + 1)
qgus = (dtacc / integrator.dt) * fastpower((get_EEst(integrator)^2) / erracc, expo)
qgus = max(inv(qmax), min(inv(qmin), qgus / gamma))
qacc = max(q, qgus)
qgus = clamp(gamma / qgus, qmin, qmax)
qacc = min(q, qgus)
else
qacc = q
end
if qsteady_min <= qacc <= qsteady_max
if qsteady_min <= inv(qacc) <= qsteady_max
qacc = one(qacc)
end
cache.dtacc = integrator.dt
cache.erracc = max(1e-2, get_EEst(integrator))
integrator.dt / qacc
integrator.dt * qacc
```

When it rejects, it's the same as the [`IController`](@ref):
Expand All @@ -1033,7 +1033,7 @@
if integrator.success_iter == 0
integrator.dt *= 0.1
else
integrator.dt = integrator.dt / cache.qold
integrator.dt = integrator.dt * cache.qold
end
```
"""
Expand Down Expand Up @@ -1084,7 +1084,7 @@
qmax = get_current_qmax(integrator, qmax)
EEst = DiffEqBase.value(get_EEst(integrator))
if iszero(EEst)
q = inv(qmax)
dt_factor = qmax
else
if fac_default_gamma(alg)
fac = gamma
Expand All @@ -1098,11 +1098,11 @@
fac = min(gamma, (1 + 2 * maxiters) * gamma / (iter + 2 * maxiters))
end
expo = 1 / (get_current_adaptive_order(alg, integrator.cache) + 1)
qtmp = fastpower(EEst, expo) / fac
@fastmath q = DiffEqBase.value(max(inv(qmax), min(inv(qmin), qtmp)))
cache.qold = q
dt_factor = fac / fastpower(EEst, expo)
@fastmath dt_factor = DiffEqBase.value(clamp(dt_factor, qmin, qmax))
cache.qold = dt_factor
end
return q
return dt_factor
end

function step_accept_controller!(integrator, cache::PredictiveControllerCache, alg, q)
Expand All @@ -1114,20 +1114,19 @@

if integrator.success_iter > 0
expo = 1 / (get_current_adaptive_order(alg, integrator.cache) + 1)
qgus = (dtacc / integrator.dt) *
fastpower((EEst^2) / erracc, expo)
qgus = max(inv(qmax), min(inv(qmin), qgus / gamma))
qacc = max(q, qgus)
qgus = (dtacc / integrator.dt) * fastpower((EEst^2) / erracc, expo)
qgus = clamp(gamma / qgus, qmin, qmax)
qacc = min(q, qgus)
else
qacc = q
end
if qsteady_min <= qacc <= qsteady_max
if qsteady_min <= inv(qacc) <= qsteady_max
qacc = one(qacc)
end
cache.dtacc = DiffEqBase.value(integrator.dt)
cache.erracc = max(1.0e-2, EEst)

return integrator.dt / qacc
return integrator.dt * qacc
end

function step_reject_controller!(integrator, cache::PredictiveControllerCache, alg)
Expand All @@ -1141,7 +1140,7 @@
return integrator.dt
end
end
return integrator.dt = success_iter == 0 ? 0.1 * dt : dt / qold
return integrator.dt = success_iter == 0 ? 0.1 * dt : dt * qold
end

"""
Expand Down
Loading