From 6739daa0ab7e88353f701677532f899c691f051d Mon Sep 17 00:00:00 2001 From: termi-official Date: Fri, 24 Oct 2025 23:56:31 +0200 Subject: [PATCH 01/21] [skip ci] First steps --- .../src/OrdinaryDiffEqBDF.jl | 2 +- lib/OrdinaryDiffEqBDF/src/controllers.jl | 2 +- lib/OrdinaryDiffEqCore/src/alg_utils.jl | 16 +++++-- .../src/integrators/controllers.jl | 10 ++-- .../src/integrators/type.jl | 10 ++-- lib/OrdinaryDiffEqCore/src/solve.jl | 48 ++++++++++++++----- .../src/OrdinaryDiffEqExtrapolation.jl | 2 +- .../src/alg_utils.jl | 2 +- .../src/OrdinaryDiffEqLowStorageRK.jl | 2 +- .../src/alg_utils.jl | 12 ++--- lib/OrdinaryDiffEqNordsieck/src/alg_utils.jl | 2 +- 11 files changed, 74 insertions(+), 34 deletions(-) diff --git a/lib/OrdinaryDiffEqBDF/src/OrdinaryDiffEqBDF.jl b/lib/OrdinaryDiffEqBDF/src/OrdinaryDiffEqBDF.jl index fb7ab3d305..9e99524d43 100644 --- a/lib/OrdinaryDiffEqBDF/src/OrdinaryDiffEqBDF.jl +++ b/lib/OrdinaryDiffEqBDF/src/OrdinaryDiffEqBDF.jl @@ -16,7 +16,7 @@ import OrdinaryDiffEqCore: alg_order, calculate_residuals!, trivial_limiter!, issplit, qsteady_min_default, qsteady_max_default, get_current_alg_order, get_current_adaptive_order, - default_controller, stepsize_controller!, + legacy_default_controller, stepsize_controller!, step_accept_controller!, step_reject_controller!, post_newton_controller!, u_modified!, DAEAlgorithm, _unwrap_val, DummyController, diff --git a/lib/OrdinaryDiffEqBDF/src/controllers.jl b/lib/OrdinaryDiffEqBDF/src/controllers.jl index cd65919803..57be463133 100644 --- a/lib/OrdinaryDiffEqBDF/src/controllers.jl +++ b/lib/OrdinaryDiffEqBDF/src/controllers.jl @@ -1,4 +1,4 @@ -function default_controller(alg::Union{QNDF, FBDF}, args...) +function legacy_default_controller(alg::Union{QNDF, FBDF}, args...) DummyController() end diff --git a/lib/OrdinaryDiffEqCore/src/alg_utils.jl b/lib/OrdinaryDiffEqCore/src/alg_utils.jl index 4539082cbb..3fd0caeb93 100644 --- a/lib/OrdinaryDiffEqCore/src/alg_utils.jl +++ b/lib/OrdinaryDiffEqCore/src/alg_utils.jl @@ -308,10 +308,20 @@ alg_adaptive_order(alg::Union{OrdinaryDiffEqAlgorithm, DAEAlgorithm}) = alg_orde # this is actually incorrect and is purposefully decreased as this tends # to track the real error much better -# this is actually incorrect and is purposefully decreased as this tends -# to track the real error much better -function default_controller(alg, cache, qoldinit, _beta1 = nothing, _beta2 = nothing) +function default_controller(alg) + if ispredictive(alg) + return PredictiveController() + elseif isstandard(alg) + return IController() + else + beta2 = beta2_default(alg) + beta1 = beta1_default(alg, beta2) + return PIController(beta1, beta2) + end +end + +function legacy_default_controller(alg, cache, qoldinit, _beta1 = nothing, _beta2 = nothing) if ispredictive(alg) return PredictiveController() elseif isstandard(alg) diff --git a/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl b/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl index 9c58c48432..e57d92332d 100644 --- a/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl +++ b/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl @@ -1,4 +1,6 @@ abstract type AbstractController end +abstract type AbstractLegacyController <: AbstractController end + using OrdinaryDiffEqCore @inline function stepsize_controller!(integrator, alg) @@ -54,7 +56,7 @@ the predicted step size. Solving Ordinary Differential Equations I Nonstiff Problems [DOI: 10.1007/978-3-540-78862-1](https://doi.org/10.1007/978-3-540-78862-1) """ -struct IController <: AbstractController +struct IController <: AbstractLegacyController end struct DummyController <: AbstractController @@ -127,7 +129,7 @@ the predicted step size. Solving Ordinary Differential Equations I Nonstiff Problems [DOI: 10.1007/978-3-540-78862-1](https://doi.org/10.1007/978-3-540-78862-1) """ -mutable struct PIController{QT} <: AbstractController +mutable struct PIController{QT} <: AbstractLegacyController beta1::QT beta2::QT end @@ -241,7 +243,7 @@ Some standard controller parameters suggested in the literature are Compressible Computational Fluid Dynamics # is bigger than this parameter [arXiv:2104.06836](https://arxiv.org/abs/2104.06836) # limiter of the dt factor (before clipping) """ -struct PIDController{QT, Limiter} <: AbstractController +struct PIDController{QT, Limiter} <: AbstractLegacyController beta::MVector{3, QT} # controller coefficients err::MVector{3, QT} # history of the error estimates accept_safety::QT # accept a step if the predicted change of the step size @@ -392,7 +394,7 @@ else end ``` """ -struct PredictiveController <: AbstractController +struct PredictiveController <: AbstractLegacyController end function post_newton_controller!(integrator, alg) diff --git a/lib/OrdinaryDiffEqCore/src/integrators/type.jl b/lib/OrdinaryDiffEqCore/src/integrators/type.jl index 7dc9e1a9c6..d9919042bd 100644 --- a/lib/OrdinaryDiffEqCore/src/integrators/type.jl +++ b/lib/OrdinaryDiffEqCore/src/integrators/type.jl @@ -6,12 +6,14 @@ mutable struct DEOptions{absType, relType, QT, tType, Controller, F1, F2, F3, F4 adaptive::Bool abstol::absType reltol::relType - gamma::QT - qmax::QT - qmin::QT + # TODO vvv remove this block as these are controller and not integrator parameters vvv + # gamma::QT + # qmax::QT + # qmin::QT qsteady_max::QT qsteady_min::QT - qoldinit::QT + # qoldinit::QT + # TODO ^^^ remove this block as these are controller and not integrator parameters ^^^ failfactor::QT dtmax::tType dtmin::tType diff --git a/lib/OrdinaryDiffEqCore/src/solve.jl b/lib/OrdinaryDiffEqCore/src/solve.jl index 6ddd866f87..9c5ac8a2b7 100644 --- a/lib/OrdinaryDiffEqCore/src/solve.jl +++ b/lib/OrdinaryDiffEqCore/src/solve.jl @@ -39,14 +39,14 @@ function SciMLBase.__init( gamma = gamma_default(alg), abstol = nothing, reltol = nothing, - qmin = qmin_default(alg), - qmax = qmax_default(alg), - qsteady_min = qsteady_min_default(alg), - qsteady_max = qsteady_max_default(alg), + qmin = nothing, + qmax = nothing, + qsteady_min = nothing, + qsteady_max = nothing, beta1 = nothing, beta2 = nothing, - qoldinit = anyadaptive(alg) ? 1 // 10^4 : 0, - controller = nothing, + qoldinit = nothing, + controller = any((qmin, qmax, qsteady_min, qsteady_max, beta1, beta2, qoldinit) .!== nothing) ? nothing : default_controller(alg), # We have to reconstruct the old controller before breaking release., fullnormalize = true, failfactor = 2, maxiters = anyadaptive(alg) ? 1000000 : typemax(Int), @@ -417,14 +417,37 @@ function SciMLBase.__init( throw(ArgumentError("Setting both the legacy PID parameters `beta1, beta2 = $((beta1, beta2))` and the `controller = $controller` is not allowed.")) end + # Deprecation warnings for users to break down which parameters they accidentally set. if (beta1 !== nothing || beta2 !== nothing) message = "Providing the legacy PID parameters `beta1, beta2` is deprecated. Use the keyword argument `controller` instead." Base.depwarn(message, :init) Base.depwarn(message, :solve) end + if (qmin !== nothing || qmax !== nothing) + message = "Providing the legacy PID parameters `qmin, qmax` is deprecated. Use the keyword argument `controller` instead." + Base.depwarn(message, :init) + Base.depwarn(message, :solve) + end + if (qsteady_min !== nothing || qsteady_max !== nothing) + message = "Providing the legacy PID parameters `qsteady_min, qsteady_max` is deprecated. Use the keyword argument `controller` instead." + Base.depwarn(message, :init) + Base.depwarn(message, :solve) + end + if (qoldinit !== nothing) + message = "Providing the legacy PID parameters `qoldinit` is deprecated. Use the keyword argument `controller` instead." + Base.depwarn(message, :init) + Base.depwarn(message, :solve) + end - if controller === nothing - controller = default_controller(_alg, cache, qoldinit, beta1, beta2) + # The folowing code provides an upgrade path for users by preserving the old behavior. + legacy_controller_parameters = (qmin, qmax, qsteady_min, qsteady_max, beta1, beta2, qoldinit) + if any(legacy_controller_parameters .== nothing) && controller === nothing # We have to reconstruct the old controller before breaking release. + qmin === nothing ? qmin_default(alg) : qmin + qmax === nothing ? qmax_default(alg) : qmax + qsteady_min === nothing ? qsteady_min_default(alg) : qsteady_min + qsteady_max === nothing ? qsteady_max_default(alg) : qsteady_max + qoldinit === nothing ? (anyadaptive(alg) ? 1 // 10^4 : 0) : qoldinit + controller = legacy_default_controller(_alg, cache, qoldinit, beta1, beta2) end save_end_user = save_end @@ -446,11 +469,14 @@ function SciMLBase.__init( typeof(saveat), typeof(d_discontinuities)}(maxiters, save_everystep, adaptive, abstol_internal, reltol_internal, - QT(gamma), QT(qmax), - QT(qmin), + # TODO vvv remove this block as these are controller and not integrator parameters vvv + # QT(gamma), + # QT(qmax), + # QT(qmin), QT(qsteady_max), QT(qsteady_min), - QT(qoldinit), + # QT(qoldinit), + # TODO ^^^remove this block as these are controller and not integrator parameters ^^^ QT(failfactor), tType(dtmax), tType(dtmin), controller, diff --git a/lib/OrdinaryDiffEqExtrapolation/src/OrdinaryDiffEqExtrapolation.jl b/lib/OrdinaryDiffEqExtrapolation/src/OrdinaryDiffEqExtrapolation.jl index 194b30e166..befa7152fc 100644 --- a/lib/OrdinaryDiffEqExtrapolation/src/OrdinaryDiffEqExtrapolation.jl +++ b/lib/OrdinaryDiffEqExtrapolation/src/OrdinaryDiffEqExtrapolation.jl @@ -3,7 +3,7 @@ module OrdinaryDiffEqExtrapolation import OrdinaryDiffEqCore: alg_order, alg_maximum_order, get_current_adaptive_order, get_current_alg_order, calculate_residuals!, accept_step_controller, - default_controller, beta2_default, beta1_default, gamma_default, + legacy_default_controller, beta2_default, beta1_default, gamma_default, initialize!, perform_step!, @unpack, @cache, unwrap_alg, isthreaded, step_accept_controller!, calculate_residuals, diff --git a/lib/OrdinaryDiffEqExtrapolation/src/alg_utils.jl b/lib/OrdinaryDiffEqExtrapolation/src/alg_utils.jl index ec2215e511..730eb11bc9 100644 --- a/lib/OrdinaryDiffEqExtrapolation/src/alg_utils.jl +++ b/lib/OrdinaryDiffEqExtrapolation/src/alg_utils.jl @@ -32,7 +32,7 @@ alg_maximum_order(alg::ImplicitHairerWannerExtrapolation) = 2(alg.max_order + 1) alg_maximum_order(alg::ImplicitEulerExtrapolation) = 2(alg.max_order + 1) alg_maximum_order(alg::ImplicitEulerBarycentricExtrapolation) = alg.max_order -function default_controller( +function legacy_default_controller( alg::Union{ExtrapolationMidpointDeuflhard, ImplicitDeuflhardExtrapolation, ExtrapolationMidpointHairerWanner, diff --git a/lib/OrdinaryDiffEqLowStorageRK/src/OrdinaryDiffEqLowStorageRK.jl b/lib/OrdinaryDiffEqLowStorageRK/src/OrdinaryDiffEqLowStorageRK.jl index 22596ff12b..fa12e2ce1d 100644 --- a/lib/OrdinaryDiffEqLowStorageRK/src/OrdinaryDiffEqLowStorageRK.jl +++ b/lib/OrdinaryDiffEqLowStorageRK/src/OrdinaryDiffEqLowStorageRK.jl @@ -7,7 +7,7 @@ import OrdinaryDiffEqCore: alg_order, alg_adaptive_order, calculate_residuals!, OrdinaryDiffEqAlgorithm, ispredictive, OrdinaryDiffEqMutableCache, OrdinaryDiffEqConstantCache, OrdinaryDiffEqAdaptiveAlgorithm, uses_uprev, - default_controller, PIDController, + legacy_default_controller, PIDController, alg_cache, _vec, _reshape, @cache, isfsal, full_cache, constvalue, _unwrap_val, trivial_limiter!, perform_step!, initialize!, diff --git a/lib/OrdinaryDiffEqLowStorageRK/src/alg_utils.jl b/lib/OrdinaryDiffEqLowStorageRK/src/alg_utils.jl index 1dc1a9e81a..9ce2293d54 100644 --- a/lib/OrdinaryDiffEqLowStorageRK/src/alg_utils.jl +++ b/lib/OrdinaryDiffEqLowStorageRK/src/alg_utils.jl @@ -84,32 +84,32 @@ uses_uprev(alg::CKLLSRK65_4M_4R, adaptive::Bool) = adaptive uses_uprev(alg::CKLLSRK85_4FM_4R, adaptive::Bool) = adaptive uses_uprev(alg::CKLLSRK75_4M_5R, adaptive::Bool) = adaptive -function default_controller(alg::RDPK3Sp35, cache, qoldinit, args...) +function legacy_default_controller(alg::RDPK3Sp35, cache, qoldinit, args...) QT = typeof(qoldinit) return PIDController(map(Base.Fix1(convert, QT), (0.64, -0.31, 0.04))...) end -function default_controller(alg::RDPK3SpFSAL35, cache, qoldinit, args...) +function legacy_default_controller(alg::RDPK3SpFSAL35, cache, qoldinit, args...) QT = typeof(qoldinit) return PIDController(map(Base.Fix1(convert, QT), (0.70, -0.23, 0.00))...) end -function default_controller(alg::RDPK3Sp49, cache, qoldinit, args...) +function legacy_default_controller(alg::RDPK3Sp49, cache, qoldinit, args...) QT = typeof(qoldinit) return PIDController(map(Base.Fix1(convert, QT), (0.25, -0.12, 0.00))...) end -function default_controller(alg::RDPK3SpFSAL49, cache, qoldinit, args...) +function legacy_default_controller(alg::RDPK3SpFSAL49, cache, qoldinit, args...) QT = typeof(qoldinit) return PIDController(map(Base.Fix1(convert, QT), (0.38, -0.18, 0.01))...) end -function default_controller(alg::RDPK3Sp510, cache, qoldinit, args...) +function legacy_default_controller(alg::RDPK3Sp510, cache, qoldinit, args...) QT = typeof(qoldinit) return PIDController(map(Base.Fix1(convert, QT), (0.47, -0.20, 0.06))...) end -function default_controller(alg::RDPK3SpFSAL510, cache, qoldinit, args...) +function legacy_default_controller(alg::RDPK3SpFSAL510, cache, qoldinit, args...) QT = typeof(qoldinit) return PIDController(map(Base.Fix1(convert, QT), (0.45, -0.13, 0.00))...) end diff --git a/lib/OrdinaryDiffEqNordsieck/src/alg_utils.jl b/lib/OrdinaryDiffEqNordsieck/src/alg_utils.jl index 7083340c99..067f482539 100644 --- a/lib/OrdinaryDiffEqNordsieck/src/alg_utils.jl +++ b/lib/OrdinaryDiffEqNordsieck/src/alg_utils.jl @@ -8,6 +8,6 @@ qsteady_max_default(alg::JVODE) = 3 // 2 get_current_alg_order(alg::JVODE, cache) = get_current_adaptive_order(alg, cache) -function default_controller(alg::Union{JVODE}, args...) +function legacy_default_controller(alg::Union{JVODE}, args...) DummyController() end From c20bbf17a072cb5576a9c661686f7d3dcbc1bd1e Mon Sep 17 00:00:00 2001 From: termi-official <9196588+termi-official@users.noreply.github.com> Date: Sun, 26 Oct 2025 19:04:53 +0100 Subject: [PATCH 02/21] Introduce caches and new controller types. --- lib/OrdinaryDiffEqCore/src/alg_utils.jl | 19 +- .../src/integrators/controllers.jl | 466 ++++++++++++++++-- .../src/integrators/type.jl | 13 +- lib/OrdinaryDiffEqCore/src/solve.jl | 50 +- .../src/alg_utils.jl | 2 +- .../src/controllers.jl | 4 +- .../src/OrdinaryDiffEqFIRK.jl | 1 + lib/OrdinaryDiffEqFIRK/src/controllers.jl | 71 ++- .../src/OrdinaryDiffEqLowStorageRK.jl | 2 +- .../src/alg_utils.jl | 12 +- .../src/OrdinaryDiffEqVerner.jl | 3 +- lib/OrdinaryDiffEqVerner/src/controllers.jl | 8 - 12 files changed, 558 insertions(+), 93 deletions(-) delete mode 100644 lib/OrdinaryDiffEqVerner/src/controllers.jl diff --git a/lib/OrdinaryDiffEqCore/src/alg_utils.jl b/lib/OrdinaryDiffEqCore/src/alg_utils.jl index 3fd0caeb93..8f918da0e1 100644 --- a/lib/OrdinaryDiffEqCore/src/alg_utils.jl +++ b/lib/OrdinaryDiffEqCore/src/alg_utils.jl @@ -309,30 +309,31 @@ alg_adaptive_order(alg::Union{OrdinaryDiffEqAlgorithm, DAEAlgorithm}) = alg_orde # this is actually incorrect and is purposefully decreased as this tends # to track the real error much better -function default_controller(alg) +function default_controller_v7(alg) if ispredictive(alg) - return PredictiveController() + return PredictiveController(alg) elseif isstandard(alg) - return IController() + return IController(alg) else - beta2 = beta2_default(alg) - beta1 = beta1_default(alg, beta2) - return PIController(beta1, beta2) + return PIController(alg) end end function legacy_default_controller(alg, cache, qoldinit, _beta1 = nothing, _beta2 = nothing) if ispredictive(alg) - return PredictiveController() + return LegacyPredictiveController() elseif isstandard(alg) - return IController() + return LegacyIController() else # Default is PI-controller QT = typeof(qoldinit) beta1, beta2 = _digest_beta1_beta2(alg, cache, Val(QT), _beta1, _beta2) - return PIController(beta1, beta2) + return LegacyPIController(beta1, beta2) end end +# TODO remove this when done +default_controller = legacy_default_controller + function _digest_beta1_beta2(alg, cache, ::Val{QT}, _beta1, _beta2) where {QT} if alg isa OrdinaryDiffEqCompositeAlgorithm beta2 = _beta2 === nothing ? diff --git a/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl b/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl index e57d92332d..6cdce66042 100644 --- a/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl +++ b/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl @@ -1,18 +1,26 @@ abstract type AbstractController end abstract type AbstractLegacyController <: AbstractController end -using OrdinaryDiffEqCore +abstract type AbstractControllerCache end + +# The legacy controllers do not have this concept. +setup_controller_cache(alg, cache, controller::AbstractLegacyController) = controller + +# checks whether the controller should accept a step based on the error estimate +@inline function accept_step_controller(integrator, controller::Union{<:AbstractLegacyController, <:AbstractControllerCache}) + return integrator.EEst <= 1 +end @inline function stepsize_controller!(integrator, alg) - stepsize_controller!(integrator, integrator.opts.controller, alg) + stepsize_controller!(integrator, integrator.controller_cache, alg) end @inline function step_accept_controller!(integrator, alg, q) - step_accept_controller!(integrator, integrator.opts.controller, alg, q) + step_accept_controller!(integrator, integrator.controller_cache, alg, q) end @inline function step_reject_controller!(integrator, alg) - step_reject_controller!(integrator, integrator.opts.controller, alg) + step_reject_controller!(integrator, integrator.controller_cache, alg) cache = integrator.cache if hasfield(typeof(cache), :nlsolve) nlsolve = cache.nlsolve @@ -25,9 +33,14 @@ reset_alg_dependent_opts!(controller::AbstractController, alg1, alg2) = nothing SciMLBase.reinit!(integrator::ODEIntegrator, controller::AbstractController) = nothing +function post_newton_controller!(integrator, alg) + integrator.dt = integrator.dt / integrator.opts.failfactor + nothing +end + # Standard integral (I) step size controller """ - IController() + LegacyIController() The standard (integral) controller is the most basic step size controller. This controller is usually the first one introduced in numerical analysis classes @@ -56,13 +69,13 @@ the predicted step size. Solving Ordinary Differential Equations I Nonstiff Problems [DOI: 10.1007/978-3-540-78862-1](https://doi.org/10.1007/978-3-540-78862-1) """ -struct IController <: AbstractLegacyController +struct LegacyIController <: AbstractLegacyController end struct DummyController <: AbstractController end -@inline function stepsize_controller!(integrator, controller::IController, alg) +@inline function stepsize_controller!(integrator, controller::LegacyIController, alg) @unpack qmin, qmax, gamma = integrator.opts EEst = DiffEqBase.value(integrator.EEst) @@ -78,7 +91,7 @@ end q end -function step_accept_controller!(integrator, controller::IController, alg, q) +function step_accept_controller!(integrator, controller::LegacyIController, alg, q) @unpack qsteady_min, qsteady_max = integrator.opts if qsteady_min <= q <= qsteady_max @@ -87,14 +100,81 @@ function step_accept_controller!(integrator, controller::IController, alg, q) integrator.dt / q # new dt end -function step_reject_controller!(integrator, controller::IController, alg) +function step_reject_controller!(integrator, controller::LegacyIController, alg) @unpack qold = integrator integrator.dt = qold end +struct IController{T} <: AbstractController + qmin::T + qmax::T + gamma::T + qsteady_min::T + qsteady_max::T +end + +function IController(alg; qmin = nothing, qmax = nothing, gamma = nothing, qsteady_min = nothing, qsteady_max = nothing) + IController( + qmin === nothing ? qmin_default(alg) : qmin, + qmax === nothing ? qmax_default(alg) : qmax, + gamma === nothing ? gamma_default(alg) : gamma, + qsteady_min === nothing ? qsteady_min_default(alg) : qsteady_min, + qsteady_max === nothing ? qsteady_max_default(alg) : qsteady_max, + ) +end + +mutable struct IControllerCache{C, T} <: AbstractControllerCache + controller::C + q::T + dtreject::T + # I believe this should go here or in the algorithm cache, but not in the integrator itself. + # EEst::T +end + +function setup_controller_cache(alg, cache, controller::IController) + IControllerCache( + controller, + 1 // 1, + 1 // 10^4,# TODO which value? + ) +end + +@inline function stepsize_controller!(integrator, cache::IControllerCache, alg) + @unpack qmin, qmax, gamma = cache.controller + EEst = DiffEqBase.value(integrator.EEst) + + if iszero(EEst) + cache.q = inv(qmax) + else + expo = 1 / (get_current_adaptive_order(alg, integrator.cache) + 1) + qtmp = fastpower(EEst, expo) / gamma + @fastmath cache.q = DiffEqBase.value(max(inv(qmax), min(inv(qmin), qtmp))) + # TODO: Shouldn't this be in `step_accept_controller!` as for the PI controller? + cache.dtreject = integrator.qold = DiffEqBase.value(integrator.dt) / cache.q + end + cache.q +end + +# TODO change signature to remove the q input +function step_accept_controller!(integrator, cache::IControllerCache, alg, q) + @unpack qsteady_min, qsteady_max = cache.controller + @assert q ≈ cache.q "Controller cache went out of sync with time stepping logic." + + if qsteady_min <= q <= qsteady_max + cache.q = q = one(q) + end + integrator.dt / q # new dt +end + +function step_reject_controller!(integrator, cache::IControllerCache, alg) + @assert cache.dtreject ≈ integrator.qold "Controller cache went out of sync with time stepping logic." + integrator.dt = cache.dtreject # TODO this does not look right. +end + + # PI step size controller """ - PIController(beta1, beta2) + LegacyPIController(beta1, beta2) The proportional-integral (PI) controller is a widespread step size controller with improved stability properties compared to the [`IController`](@ref). @@ -117,7 +197,7 @@ the predicted step size. !!! note The coefficients `beta1, beta2` are not scaled by the order of the method, - in contrast to the [`PIDController`](@ref). For the `PIController`, this + in contrast to the [`LegacyPIDController`](@ref). For the `LegacyPIController`, this scaling by the order must be done when the controller is constructed. ## References @@ -129,12 +209,12 @@ the predicted step size. Solving Ordinary Differential Equations I Nonstiff Problems [DOI: 10.1007/978-3-540-78862-1](https://doi.org/10.1007/978-3-540-78862-1) """ -mutable struct PIController{QT} <: AbstractLegacyController +mutable struct LegacyPIController{QT} <: AbstractLegacyController beta1::QT beta2::QT end -@inline function stepsize_controller!(integrator, controller::PIController, alg) +@inline function stepsize_controller!(integrator, controller::LegacyPIController, alg) @unpack qold = integrator @unpack qmin, qmax, gamma = integrator.opts @unpack beta1, beta2 = controller @@ -151,7 +231,7 @@ end q end -function step_accept_controller!(integrator, controller::PIController, alg, q) +function step_accept_controller!(integrator, controller::LegacyPIController, alg, q) @unpack qsteady_min, qsteady_max, qoldinit = integrator.opts EEst = DiffEqBase.value(integrator.EEst) @@ -162,13 +242,13 @@ function step_accept_controller!(integrator, controller::PIController, alg, q) return integrator.dt / q # new dt end -function step_reject_controller!(integrator, controller::PIController, alg) +function step_reject_controller!(integrator, controller::LegacyPIController, alg) @unpack q11 = integrator @unpack qmin, gamma = integrator.opts integrator.dt /= min(inv(qmin), q11 / gamma) end -function reset_alg_dependent_opts!(controller::PIController, alg1, alg2) +function reset_alg_dependent_opts!(controller::LegacyPIController, alg1, alg2) if controller.beta2 == beta2_default(alg1) controller.beta2 = beta2_default(alg2) end @@ -177,14 +257,107 @@ function reset_alg_dependent_opts!(controller::PIController, alg1, alg2) end end +struct PIController{T} <: AbstractController + beta1::T + beta2::T + qmin::T + qmax::T + gamma::T + qsteady_min::T + qsteady_max::T + qoldinit::T +end + +function PIController(alg; beta1 = nothing, beta2 = nothing, qmin = nothing, qmax = nothing, gamma = nothing, qsteady_min = nothing, qsteady_max = nothing, qoldinit = nothing) + beta2 = beta2 === nothing ? beta2_default(alg) : beta2 + beta1 = beta1 === nothing ? beta1_default(alg, beta2) : beta1 + qoldinit = qoldinit === nothing ? 1 // 10^4 : qoldinit + PIController{promote_type(typeof(qoldinit), typeof(beta1))}( + beta1, + beta2, + qmin === nothing ? qmin_default(alg) : qmin, + qmax === nothing ? qmax_default(alg) : qmax, + gamma === nothing ? gamma_default(alg) : gamma, + qsteady_min === nothing ? qsteady_min_default(alg) : qsteady_min, + qsteady_max === nothing ? qsteady_max_default(alg) : qsteady_max, + qoldinit, + ) +end + +mutable struct PIControllerCache{T} <: AbstractControllerCache + controller::PIController{T} + q::T + q11::T + qold::T +end + + +function setup_controller_cache(alg, cache, controller::PIController) + PIControllerCache( + controller, + 1 // 1, + 1 // 10^4,# TODO which value? + 1 // 10^4, + ) +end + +@inline function stepsize_controller!(integrator, cache::PIControllerCache, alg) + @unpack controller = cache + @unpack qold = integrator + @unpack qmin, qmax, gamma = controller + @unpack beta1, beta2 = controller + EEst = DiffEqBase.value(integrator.EEst) + + if iszero(EEst) + q = inv(qmax) + else + q11 = fastpower(EEst, convert(typeof(EEst), beta1)) + q = q11 / fastpower(qold, convert(typeof(EEst), beta2)) + integrator.q11 = q11 + @fastmath q = max(inv(qmax), min(inv(qmin), q / gamma)) + end + cache.q = q + q +end + +function step_accept_controller!(integrator, cache::PIControllerCache, alg, q) + @assert q ≈ cache.q "Controller cache went out of sync with time stepping logic (q=$q | cache.q=$(cache.q))." + @unpack controller = cache + @unpack qsteady_min, qsteady_max, qoldinit = controller + EEst = DiffEqBase.value(integrator.EEst) + + if qsteady_min <= q <= qsteady_max + q = one(q) + end + cache.qold = max(EEst, qoldinit) + return integrator.dt / q # new dt +end + +function step_reject_controller!(integrator, cache::PIControllerCache, alg) + @unpack controller = cache + @unpack q11 = cache + @unpack qmin, gamma = controller + integrator.dt /= min(inv(qmin), q11 / gamma) +end + +# FIXME multi-controller? +# function reset_alg_dependent_opts!(controller::PIControllerCache, alg1, alg2) +# if controller.beta2 == beta2_default(alg1) +# controller.beta2 = beta2_default(alg2) +# end +# if controller.beta1 == beta1_default(alg1, controller.beta2) +# controller.beta1 = beta1_default(alg2, controller.beta2) +# end +# end + # PID step size controller """ - PIDController(beta1, beta2, beta3=zero(beta1); + LegacyPIDController(beta1, beta2, beta3=zero(beta1); limiter=default_dt_factor_limiter, accept_safety=0.81) The proportional-integral-derivative (PID) controller is a generalization of the -[`PIController`](@ref) and can have improved stability and efficiency properties. +[`LegacyPIController`](@ref) and can have improved stability and efficiency properties. Construct a PID step size controller adapting the time step based on the formula @@ -217,17 +390,17 @@ Some standard controller parameters suggested in the literature are !!! note - In contrast to the [`PIController`](@ref), the coefficients `beta1, beta2, beta3` + In contrast to the [`LegacyPIController`](@ref), the coefficients `beta1, beta2, beta3` are scaled by the order of the method. Thus, standard controllers such as PI42 can use the same coefficients `beta1, beta2, beta3` for different algorithms. !!! note - In contrast to other controllers, the `PIDController` does not use the keyword + In contrast to other controllers, the `LegacyPIDController` does not use the keyword arguments `qmin, qmax` to limit the step size change or the safety factor `gamma`. These common keyword arguments are replaced by the `limiter` and `accept_safety` to guarantee a smooth behavior (Söderlind and Wang, 2006). - Because of this, a `PIDController` behaves different from a [`PIController`](@ref), + Because of this, a `LegacyPIDController` behaves different from a [`LegacyPIController`](@ref), even if `beta1, beta2` are adapted accordingly and `iszero(beta3)`. ## References @@ -243,7 +416,7 @@ Some standard controller parameters suggested in the literature are Compressible Computational Fluid Dynamics # is bigger than this parameter [arXiv:2104.06836](https://arxiv.org/abs/2104.06836) # limiter of the dt factor (before clipping) """ -struct PIDController{QT, Limiter} <: AbstractLegacyController +struct LegacyPIDController{QT, Limiter} <: AbstractLegacyController beta::MVector{3, QT} # controller coefficients err::MVector{3, QT} # history of the error estimates accept_safety::QT # accept a step if the predicted change of the step size @@ -251,17 +424,17 @@ struct PIDController{QT, Limiter} <: AbstractLegacyController limiter::Limiter # limiter of the dt factor (before clipping) end -function PIDController(beta1, beta2, beta3 = zero(beta1); +function LegacyPIDController(beta1, beta2, beta3 = zero(beta1); limiter = default_dt_factor_limiter, accept_safety = 0.81) beta = MVector(map(float, promote(beta1, beta2, beta3))...) QT = eltype(beta) err = MVector{3, QT}(true, true, true) - return PIDController(beta, err, convert(QT, accept_safety), limiter) + return LegacyPIDController(beta, err, convert(QT, accept_safety), limiter) end -function Base.show(io::IO, controller::PIDController) - print(io, "PIDController(beta=", controller.beta, +function Base.show(io::IO, controller::LegacyPIDController) + print(io, "LegacyPIDController(beta=", controller.beta, ", accept_safety=", controller.accept_safety, ", limiter=", controller.limiter, ")") @@ -269,7 +442,7 @@ end @inline default_dt_factor_limiter(x) = one(x) + atan(x - one(x)) -@inline function stepsize_controller!(integrator, controller::PIDController, alg) +@inline function stepsize_controller!(integrator, controller::LegacyPIDController, alg) @unpack qmax = integrator.opts beta1, beta2, beta3 = controller.beta @@ -313,16 +486,11 @@ end return dt_factor end -# 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) +@inline function accept_step_controller(integrator, controller::LegacyPIDController) return integrator.qold >= controller.accept_safety end -function step_accept_controller!(integrator, controller::PIDController, alg, dt_factor) +function step_accept_controller!(integrator, controller::LegacyPIDController, alg, dt_factor) @unpack qsteady_min, qsteady_max = integrator.opts if qsteady_min <= inv(dt_factor) <= qsteady_max @@ -335,10 +503,135 @@ function step_accept_controller!(integrator, controller::PIDController, alg, dt_ return integrator.dt * dt_factor # new dt end -function step_reject_controller!(integrator, controller::PIDController, alg) +function step_reject_controller!(integrator, controller::LegacyPIDController, alg) integrator.dt *= integrator.qold end + + +struct PIDController{T, Limiter} <: AbstractController + beta::SVector{3, T} # controller coefficients + accept_safety::T # accept a step if the predicted change of the step size + # is bigger than this parameter + limiter::Limiter # limiter of the dt factor (before clipping) + qmin::T + qmax::T + gamma::T + qsteady_min::T + qsteady_max::T + qoldinit::T +end + +function PIDController(alg; beta1 = nothing, beta2 = nothing, beta3 = nothing, accept_safety = 0.81, limiter = default_dt_factor_limiter, qmin = nothing, qmax = nothing, gamma = nothing, qsteady_min = nothing, qsteady_max = nothing, qoldinit = nothing) + beta2 = beta2 === nothing ? beta2_default(alg) : beta2 + beta1 = beta1 === nothing ? beta1_default(alg, beta2) : beta1 + beta3 = beta3 === nothing ? zero(beta1) : beta3 + beta = MVector(map(float, promote(beta1, beta2, beta3))...) + PIDController( + beta, + err, + accept_safety, + limiter, + qmin === nothing ? qmin_default(alg) : qmin, + qmax === nothing ? qmax_default(alg) : qmax, + gamma === nothing ? gamma_default(alg) : gamma, + qsteady_min === nothing ? qsteady_min_default(alg) : qsteady_min, + qsteady_max === nothing ? qsteady_max_default(alg) : qsteady_max, + qoldinit === nothing ? 1 // 10^4 : qoldinit, + ) +end + +function Base.show(io::IO, controller::PIDController) + print(io, "PIDController(beta=", controller.beta, + ", accept_safety=", controller.accept_safety, + ", limiter=", controller.limiter, + ")") +end + +struct PIDControllerCache{T, Limiter} <: AbstractControllerCache + controller::PIDController{T, Limiter} + err::MVector{3, T} # history of the error estimates + dt_factor::T +end + +function setup_controller_cache(alg, cache, controller::PIDController{QT}) where QT + err = MVector{3, QT}(true, true, true) + PIControllerCache( + controller, + err, + QT(1 // 10^4), + ) +end + +@inline function stepsize_controller!(integrator, cache::PIDControllerCache, alg) + @unpack controller = cache + beta1, beta2, beta3 = controller.beta + @unpack qmax = controller + + EEst = DiffEqBase.value(integrator.EEst) + + # If the error estimate is zero, we can increase the step size as much as + # desired. This additional check fixes problems of the code below when the + # error estimates become zero + # -> err1, err2, err3 become Inf + # -> err1^positive_number * err2^negative_number becomes NaN + # -> dt becomes NaN + # + # `EEst_min` is smaller than PETSC_SMALL used in the equivalent logic in PETSc. + # For example, `eps(Float64) ≈ 2.2e-16` but `PETSC_SMALL ≈ 1.0e-10` for `double`. + EEst_min = eps(typeof(EEst)) + # The code below is a bit more robust than + # ``` + # if iszero(EEst) + # EEst = eps(typeof(EEst)) + # end + # ``` + EEst = ifelse(EEst > EEst_min, EEst, EEst_min) + + cache.err[1] = inv(EEst) + 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) + if isnan(dt_factor) + @warn "unlimited dt_factor" dt_factor err1 err2 err3 beta1 beta2 beta3 k + end + dt_factor = controller.limiter(dt_factor) + + # Note: No additional limiting of the form + # dt_factor = max(qmin, min(qmax, dt_factor)) + # is necessary since the `limiter` should take care of that. The default limiter + # ensures + # 0.21 ≈ limiter(0) <= dt_factor <= limiter(Inf) ≈ 2.57 + # See Söderlind, Wang (2006), Section 6. + cache.qold = dt_factor + return dt_factor +end + +@inline function accept_step_controller(integrator, cache::PIDControllerCache) + return cache.qold >= cache.controller.accept_safety +end + +function step_accept_controller!(integrator, cache::PIDControllerCache, alg, dt_factor) + @assert dt_factor ≈ cache.q "Controller cache went out of sync with time stepping logic." + @unpack controller = cache + @unpack qsteady_min, qsteady_max = controller + + if qsteady_min <= inv(dt_factor) <= qsteady_max + dt_factor = one(dt_factor) + # cache.q = ...? + end + @inbounds begin + controller.err[3] = controller.err[2] + controller.err[2] = controller.err[1] + end + return integrator.dt * dt_factor # new dt +end + +function step_reject_controller!(integrator, cache::PIDControllerCache, alg) + integrator.dt *= cache.qold +end + # Gustafsson predictive step size controller """ PredictiveController() @@ -394,15 +687,10 @@ else end ``` """ -struct PredictiveController <: AbstractLegacyController -end - -function post_newton_controller!(integrator, alg) - integrator.dt = integrator.dt / integrator.opts.failfactor - nothing +struct LegacyPredictiveController <: AbstractController end -@inline function stepsize_controller!(integrator, controller::PredictiveController, alg) +@inline function stepsize_controller!(integrator, controller::LegacyPredictiveController, alg) @unpack qmin, qmax, gamma = integrator.opts EEst = DiffEqBase.value(integrator.EEst) if iszero(EEst) @@ -427,7 +715,7 @@ end q end -function step_accept_controller!(integrator, controller::PredictiveController, alg, q) +function step_accept_controller!(integrator, controller::LegacyPredictiveController, alg, q) @unpack qmin, qmax, gamma, qsteady_min, qsteady_max = integrator.opts EEst = DiffEqBase.value(integrator.EEst) @@ -450,7 +738,97 @@ function step_accept_controller!(integrator, controller::PredictiveController, a return integrator.dt / qacc end -function step_reject_controller!(integrator, controller::PredictiveController, alg) +function step_reject_controller!(integrator, controller::LegacyPredictiveController, alg) @unpack dt, success_iter, qold = integrator integrator.dt = success_iter == 0 ? 0.1 * dt : dt / qold end + + +struct PredictiveController{T} <: AbstractController + qmin::T + qmax::T + gamma::T + qsteady_min::T + qsteady_max::T +end + +function PredictiveController(alg; qmin = nothing, qmax = nothing, gamma = nothing, qsteady_min = nothing, qsteady_max = nothing) + PredictiveController( + qmin === nothing ? qmin_default(alg) : qmin, + qmax === nothing ? qmax_default(alg) : qmax, + gamma === nothing ? gamma_default(alg) : gamma, + qsteady_min === nothing ? qsteady_min_default(alg) : qsteady_min, + qsteady_max === nothing ? qsteady_max_default(alg) : qsteady_max, + ) +end + +mutable struct PredictiveControllerCache{T} <: AbstractControllerCache + controller::PredictiveController{T} + dtacc::T + erracc::T +end + +function setup_controller_cache(alg, cache, controller::PredictiveController) + PredictiveControllerCache( + controller, + 1 // 1, + 1 // 1, + ) +end + +@inline function stepsize_controller!(integrator, cache::PredictiveControllerCache, alg) + @unpack qmin, qmax, gamma = integrator.opts + EEst = DiffEqBase.value(integrator.EEst) + if iszero(EEst) + q = inv(qmax) + else + if fac_default_gamma(alg) + fac = gamma + else + if isfirk(alg) + @unpack iter = integrator.cache + @unpack maxiters = alg + else + @unpack iter, maxiters = integrator.cache.nlsolver + end + 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 + end + cache.q = q + q +end + +function step_accept_controller!(integrator, cache::PredictiveControllerCache, alg, q) + @assert q ≈ cache.q "Controller cache went out of sync with time stepping logic." + @unpack dtacc, erracc, controller = cache + @unpack qmin, qmax, gamma, qsteady_min, qsteady_max = controller + + EEst = DiffEqBase.value(integrator.EEst) + + 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) + else + qacc = q + end + if qsteady_min <= qacc <= qsteady_max + qacc = one(qacc) + end + cache.dtacc = integrator.dt + cache.erracc = max(1e-2, EEst) + + return integrator.dt / qacc +end + +function step_reject_controller!(integrator, cache::PredictiveControllerCache, alg) + @unpack dt, success_iter = integrator + @unpack qold = cache + integrator.dt = success_iter == 0 ? 0.1 * dt : dt / qold +end diff --git a/lib/OrdinaryDiffEqCore/src/integrators/type.jl b/lib/OrdinaryDiffEqCore/src/integrators/type.jl index d9919042bd..804fd275f3 100644 --- a/lib/OrdinaryDiffEqCore/src/integrators/type.jl +++ b/lib/OrdinaryDiffEqCore/src/integrators/type.jl @@ -7,12 +7,12 @@ mutable struct DEOptions{absType, relType, QT, tType, Controller, F1, F2, F3, F4 abstol::absType reltol::relType # TODO vvv remove this block as these are controller and not integrator parameters vvv - # gamma::QT - # qmax::QT - # qmin::QT + gamma::QT + qmax::QT + qmin::QT qsteady_max::QT qsteady_min::QT - # qoldinit::QT + qoldinit::QT # TODO ^^^ remove this block as these are controller and not integrator parameters ^^^ failfactor::QT dtmax::tType @@ -86,7 +86,7 @@ For more info see the linked documentation page. mutable struct ODEIntegrator{algType <: Union{OrdinaryDiffEqAlgorithm, DAEAlgorithm}, IIP, uType, duType, tType, pType, eigenType, EEstT, QT, tdirType, ksEltype, SolType, F, CacheType, O, FSALType, EventErrorType, - CallbackCacheType, IA, DV} <: + CallbackCacheType, IA, DV, CC} <: SciMLBase.AbstractODEIntegrator{algType, IIP, uType, tType} sol::SolType u::uType @@ -107,10 +107,13 @@ mutable struct ODEIntegrator{algType <: Union{OrdinaryDiffEqAlgorithm, DAEAlgori tdir::tdirType eigen_est::eigenType EEst::EEstT + # TODO vvv remove these qold::QT q11::QT erracc::QT dtacc::tType + # TODO ^^^ remove these + controller_cache::CC success_iter::Int iter::Int saveiter::Int diff --git a/lib/OrdinaryDiffEqCore/src/solve.jl b/lib/OrdinaryDiffEqCore/src/solve.jl index 9c5ac8a2b7..a06abb05fd 100644 --- a/lib/OrdinaryDiffEqCore/src/solve.jl +++ b/lib/OrdinaryDiffEqCore/src/solve.jl @@ -36,9 +36,9 @@ function SciMLBase.__init( dtmax = eltype(prob.tspan)((prob.tspan[end] - prob.tspan[1])), force_dtmin = false, adaptive = anyadaptive(alg), - gamma = gamma_default(alg), abstol = nothing, reltol = nothing, + gamma = nothing, qmin = nothing, qmax = nothing, qsteady_min = nothing, @@ -46,7 +46,7 @@ function SciMLBase.__init( beta1 = nothing, beta2 = nothing, qoldinit = nothing, - controller = any((qmin, qmax, qsteady_min, qsteady_max, beta1, beta2, qoldinit) .!== nothing) ? nothing : default_controller(alg), # We have to reconstruct the old controller before breaking release., + controller = any((gamma, qmin, qmax, qsteady_min, qsteady_max, beta1, beta2, qoldinit) .!== nothing) ? nothing : default_controller_v7(alg), # We have to reconstruct the old controller before breaking release., fullnormalize = true, failfactor = 2, maxiters = anyadaptive(alg) ? 1000000 : typemax(Int), @@ -440,23 +440,39 @@ function SciMLBase.__init( end # The folowing code provides an upgrade path for users by preserving the old behavior. - legacy_controller_parameters = (qmin, qmax, qsteady_min, qsteady_max, beta1, beta2, qoldinit) - if any(legacy_controller_parameters .== nothing) && controller === nothing # We have to reconstruct the old controller before breaking release. + legacy_controller_parameters = (gamma, qmin, qmax, qsteady_min, qsteady_max, beta1, beta2, qoldinit) + if controller === nothing # We have to reconstruct the old controller before breaking release. + if any(legacy_controller_parameters .== nothing) + gamma === nothing ? gamma_default(alg) : gamma + qmin === nothing ? qmin_default(alg) : qmin + qmax === nothing ? qmax_default(alg) : qmax + qsteady_min === nothing ? qsteady_min_default(alg) : qsteady_min + qsteady_max === nothing ? qsteady_max_default(alg) : qsteady_max + qoldinit === nothing ? (anyadaptive(alg) ? 1 // 10^4 : 0) : qoldinit + controller = legacy_default_controller(_alg, cache, qoldinit, beta1, beta2) + end + elseif controller isa AbstractLegacyController # Legacy controller has been passed + gamma === nothing ? gamma_default(alg) : gamma qmin === nothing ? qmin_default(alg) : qmin qmax === nothing ? qmax_default(alg) : qmax qsteady_min === nothing ? qsteady_min_default(alg) : qsteady_min qsteady_max === nothing ? qsteady_max_default(alg) : qsteady_max qoldinit === nothing ? (anyadaptive(alg) ? 1 // 10^4 : 0) : qoldinit - controller = legacy_default_controller(_alg, cache, qoldinit, beta1, beta2) + else + @unpack gamma, qmin, qmax, qsteady_min, qsteady_max, qoldinit = controller end + controller_cache = setup_controller_cache(_alg, cache, controller) + save_end_user = save_end save_end = save_end === nothing ? save_everystep || isempty(saveat) || saveat isa Number || prob.tspan[2] in saveat : save_end opts = DEOptions{typeof(abstol_internal), typeof(reltol_internal), - QT, tType, typeof(controller), + QT, tType, + # typeof(controller), + typeof(controller_cache), typeof(internalnorm), typeof(internalopnorm), typeof(save_end_user), typeof(callbacks_internal), @@ -470,16 +486,19 @@ function SciMLBase.__init( adaptive, abstol_internal, reltol_internal, # TODO vvv remove this block as these are controller and not integrator parameters vvv - # QT(gamma), - # QT(qmax), - # QT(qmin), + QT(gamma), + QT(qmax), + QT(qmin), QT(qsteady_max), QT(qsteady_min), - # QT(qoldinit), + QT(qoldinit), # TODO ^^^remove this block as these are controller and not integrator parameters ^^^ QT(failfactor), tType(dtmax), tType(dtmin), - controller, + # TODO vvv remove this vvv + # controller, + controller_cache, + # TODO ^^^ remove this ^^^ internalnorm, internalopnorm, save_idxs, tstops_internal, @@ -556,13 +575,18 @@ function SciMLBase.__init( FType, cacheType, typeof(opts), typeof(fsalfirst), typeof(last_event_error), typeof(callback_cache), - typeof(initializealg), typeof(differential_vars)}( + typeof(initializealg), typeof(differential_vars), + typeof(controller_cache)}( sol, u, du, k, t, tType(dt), f, p, uprev, uprev2, duprev, tprev, _alg, dtcache, dtchangeable, dtpropose, tdir, eigen_est, EEst, + # TODO vvv remove these QT(qoldinit), q11, - erracc, dtacc, success_iter, + erracc, dtacc, + # TODO ^^^ remove these + controller_cache, + success_iter, iter, saveiter, saveiter_dense, cache, callback_cache, kshortsize, force_stepfail, diff --git a/lib/OrdinaryDiffEqExtrapolation/src/alg_utils.jl b/lib/OrdinaryDiffEqExtrapolation/src/alg_utils.jl index 730eb11bc9..c607dc6dee 100644 --- a/lib/OrdinaryDiffEqExtrapolation/src/alg_utils.jl +++ b/lib/OrdinaryDiffEqExtrapolation/src/alg_utils.jl @@ -43,7 +43,7 @@ function legacy_default_controller( qoldinit, _beta1 = nothing, _beta2 = nothing) QT = typeof(qoldinit) beta1, beta2 = _digest_beta1_beta2(alg, cache, Val(QT), _beta1, _beta2) - return ExtrapolationController(beta1) + return LegacyExtrapolationController(beta1) end beta2_default(alg::ExtrapolationMidpointDeuflhard) = 0 // 1 diff --git a/lib/OrdinaryDiffEqExtrapolation/src/controllers.jl b/lib/OrdinaryDiffEqExtrapolation/src/controllers.jl index a62d61cb77..e47139b474 100644 --- a/lib/OrdinaryDiffEqExtrapolation/src/controllers.jl +++ b/lib/OrdinaryDiffEqExtrapolation/src/controllers.jl @@ -1,9 +1,9 @@ # Extrapolation methods -mutable struct ExtrapolationController{QT} <: AbstractController +mutable struct LegacyExtrapolationController{QT} <: AbstractLegacyController beta1::QT end -function reset_alg_dependent_opts!(controller::ExtrapolationController, alg1, alg2) +function reset_alg_dependent_opts!(controller::LegacyExtrapolationController, alg1, alg2) if controller.beta1 == beta1_default(alg1, beta2_default(alg1)) controller.beta1 = beta1_default(alg2, beta2_default(alg2)) end diff --git a/lib/OrdinaryDiffEqFIRK/src/OrdinaryDiffEqFIRK.jl b/lib/OrdinaryDiffEqFIRK/src/OrdinaryDiffEqFIRK.jl index 7ab3cffccf..5c7860651e 100644 --- a/lib/OrdinaryDiffEqFIRK/src/OrdinaryDiffEqFIRK.jl +++ b/lib/OrdinaryDiffEqFIRK/src/OrdinaryDiffEqFIRK.jl @@ -14,6 +14,7 @@ import OrdinaryDiffEqCore: alg_order, calculate_residuals!, qmax_default, alg_adaptive_order, DEFAULT_PRECS, stepsize_controller!, step_accept_controller!, step_reject_controller!, + LegacyPredictiveController, PredictiveControllerCache, PredictiveController, alg_can_repeat_jac, NewtonAlgorithm, fac_default_gamma, get_current_adaptive_order, get_fsalfirstlast, diff --git a/lib/OrdinaryDiffEqFIRK/src/controllers.jl b/lib/OrdinaryDiffEqFIRK/src/controllers.jl index f491a9e7d5..480503010f 100644 --- a/lib/OrdinaryDiffEqFIRK/src/controllers.jl +++ b/lib/OrdinaryDiffEqFIRK/src/controllers.jl @@ -1,5 +1,5 @@ function step_accept_controller!( - integrator, controller::PredictiveController, alg::AdaptiveRadau, q) + integrator, controller::LegacyPredictiveController, alg::AdaptiveRadau, q) @unpack qmin, qmax, gamma, qsteady_min, qsteady_max = integrator.opts @unpack cache = integrator @unpack num_stages, step, iter, hist_iter, index = cache @@ -43,7 +43,7 @@ function step_accept_controller!( end function step_reject_controller!( - integrator, controller::PredictiveController, alg::AdaptiveRadau) + integrator, controller::LegacyPredictiveController, alg::AdaptiveRadau) @unpack dt, success_iter, qold = integrator @unpack cache = integrator @unpack num_stages, step, iter, hist_iter = cache @@ -62,3 +62,70 @@ function step_reject_controller!( end end end + +function step_accept_controller!( + integrator, ccache::PredictiveControllerCache, alg::AdaptiveRadau, q) + @unpack controller = ccache + @unpack qmin, qmax, gamma, qsteady_min, qsteady_max = controller + @unpack cache = integrator + @unpack num_stages, step, iter, hist_iter, index = cache + + EEst = DiffEqBase.value(integrator.EEst) + + if integrator.success_iter > 0 + expo = 1 / (get_current_adaptive_order(alg, cache) + 1) + qgus = (ccache.dtacc / integrator.dt) * + fastpow((EEst^2) / ccache.erracc, expo) + qgus = max(inv(qmax), min(inv(qmin), qgus / gamma)) + qacc = max(q, qgus) + else + qacc = q + end + if qsteady_min <= qacc <= qsteady_max + qacc = one(qacc) + end + ccache.dtacc = integrator.dt + ccache.erracc = max(1e-2, EEst) + cache.step = step + 1 + hist_iter = hist_iter * 0.8 + iter * 0.2 + cache.hist_iter = hist_iter + max_stages = (alg.max_order - 1) ÷ 4 * 2 + 1 + min_stages = (alg.min_order - 1) ÷ 4 * 2 + 1 + if (step > 10) + if (hist_iter < 2.6 && num_stages < max_stages) + cache.num_stages += 2 + cache.index += 1 + cache.step = 1 + cache.hist_iter = iter + elseif ((hist_iter > 8 || cache.status == VerySlowConvergence || + cache.status == Divergence) && num_stages > min_stages) + cache.num_stages -= 2 + cache.index -= 1 + cache.step = 1 + cache.hist_iter = iter + end + end + return integrator.dt / qacc +end + +function step_reject_controller!( + integrator, ccache::PredictiveControllerCache, alg::AdaptiveRadau) + @unpack controller = ccache + @unpack dt, success_iter = integrator + @unpack cache = integrator + @unpack num_stages, step, iter, hist_iter = cache + integrator.dt = success_iter == 0 ? 0.1 * dt : dt / ccache.qold + cache.step = step + 1 + hist_iter = hist_iter * 0.8 + iter * 0.2 + cache.hist_iter = hist_iter + min_stages = (alg.min_order - 1) ÷ 4 * 2 + 1 + if (step > 10) + if ((hist_iter > 8 || cache.status == VerySlowConvergence || + cache.status == Divergence) && num_stages > min_stages) + cache.num_stages -= 2 + cache.index -= 1 + cache.step = 1 + cache.hist_iter = iter + end + end +end diff --git a/lib/OrdinaryDiffEqLowStorageRK/src/OrdinaryDiffEqLowStorageRK.jl b/lib/OrdinaryDiffEqLowStorageRK/src/OrdinaryDiffEqLowStorageRK.jl index fa12e2ce1d..315938b601 100644 --- a/lib/OrdinaryDiffEqLowStorageRK/src/OrdinaryDiffEqLowStorageRK.jl +++ b/lib/OrdinaryDiffEqLowStorageRK/src/OrdinaryDiffEqLowStorageRK.jl @@ -7,7 +7,7 @@ import OrdinaryDiffEqCore: alg_order, alg_adaptive_order, calculate_residuals!, OrdinaryDiffEqAlgorithm, ispredictive, OrdinaryDiffEqMutableCache, OrdinaryDiffEqConstantCache, OrdinaryDiffEqAdaptiveAlgorithm, uses_uprev, - legacy_default_controller, PIDController, + legacy_default_controller, LegacyPIDController, alg_cache, _vec, _reshape, @cache, isfsal, full_cache, constvalue, _unwrap_val, trivial_limiter!, perform_step!, initialize!, diff --git a/lib/OrdinaryDiffEqLowStorageRK/src/alg_utils.jl b/lib/OrdinaryDiffEqLowStorageRK/src/alg_utils.jl index 9ce2293d54..064403627f 100644 --- a/lib/OrdinaryDiffEqLowStorageRK/src/alg_utils.jl +++ b/lib/OrdinaryDiffEqLowStorageRK/src/alg_utils.jl @@ -86,30 +86,30 @@ uses_uprev(alg::CKLLSRK75_4M_5R, adaptive::Bool) = adaptive function legacy_default_controller(alg::RDPK3Sp35, cache, qoldinit, args...) QT = typeof(qoldinit) - return PIDController(map(Base.Fix1(convert, QT), (0.64, -0.31, 0.04))...) + return LegacyPIDController(map(Base.Fix1(convert, QT), (0.64, -0.31, 0.04))...) end function legacy_default_controller(alg::RDPK3SpFSAL35, cache, qoldinit, args...) QT = typeof(qoldinit) - return PIDController(map(Base.Fix1(convert, QT), (0.70, -0.23, 0.00))...) + return LegacyPIDController(map(Base.Fix1(convert, QT), (0.70, -0.23, 0.00))...) end function legacy_default_controller(alg::RDPK3Sp49, cache, qoldinit, args...) QT = typeof(qoldinit) - return PIDController(map(Base.Fix1(convert, QT), (0.25, -0.12, 0.00))...) + return LegacyPIDController(map(Base.Fix1(convert, QT), (0.25, -0.12, 0.00))...) end function legacy_default_controller(alg::RDPK3SpFSAL49, cache, qoldinit, args...) QT = typeof(qoldinit) - return PIDController(map(Base.Fix1(convert, QT), (0.38, -0.18, 0.01))...) + return LegacyPIDController(map(Base.Fix1(convert, QT), (0.38, -0.18, 0.01))...) end function legacy_default_controller(alg::RDPK3Sp510, cache, qoldinit, args...) QT = typeof(qoldinit) - return PIDController(map(Base.Fix1(convert, QT), (0.47, -0.20, 0.06))...) + return LegacyPIDController(map(Base.Fix1(convert, QT), (0.47, -0.20, 0.06))...) end function legacy_default_controller(alg::RDPK3SpFSAL510, cache, qoldinit, args...) QT = typeof(qoldinit) - return PIDController(map(Base.Fix1(convert, QT), (0.45, -0.13, 0.00))...) + return LegacyPIDController(map(Base.Fix1(convert, QT), (0.45, -0.13, 0.00))...) end diff --git a/lib/OrdinaryDiffEqVerner/src/OrdinaryDiffEqVerner.jl b/lib/OrdinaryDiffEqVerner/src/OrdinaryDiffEqVerner.jl index 4fd7a60f84..ee48edce7e 100644 --- a/lib/OrdinaryDiffEqVerner/src/OrdinaryDiffEqVerner.jl +++ b/lib/OrdinaryDiffEqVerner/src/OrdinaryDiffEqVerner.jl @@ -4,7 +4,7 @@ import OrdinaryDiffEqCore: alg_order, calculate_residuals!, initialize!, perform_step!, @unpack, unwrap_alg, calculate_residuals, alg_stability_size, OrdinaryDiffEqAlgorithm, - CompositeAlgorithm, AbstractController, PIDController, + CompositeAlgorithm, OrdinaryDiffEqMutableCache, OrdinaryDiffEqConstantCache, OrdinaryDiffEqAdaptiveAlgorithm, CompiledFloats, uses_uprev, alg_cache, _vec, _reshape, @cache, isfsal, full_cache, @@ -30,7 +30,6 @@ include("verner_caches.jl") include("verner_addsteps.jl") include("interp_func.jl") include("interpolants.jl") -include("controllers.jl") include("verner_rk_perform_step.jl") import PrecompileTools diff --git a/lib/OrdinaryDiffEqVerner/src/controllers.jl b/lib/OrdinaryDiffEqVerner/src/controllers.jl deleted file mode 100644 index e61ba34b6c..0000000000 --- a/lib/OrdinaryDiffEqVerner/src/controllers.jl +++ /dev/null @@ -1,8 +0,0 @@ -# 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 From d9dd6bb860f9c5a0abf354d92d68df794ebc4e80 Mon Sep 17 00:00:00 2001 From: termi-official <9196588+termi-official@users.noreply.github.com> Date: Sun, 26 Oct 2025 19:24:39 +0100 Subject: [PATCH 03/21] Fix precompilation errors --- lib/OrdinaryDiffEqCore/src/alg_utils.jl | 8 +-- .../src/integrators/controllers.jl | 58 ++++++++++++------- lib/OrdinaryDiffEqCore/src/solve.jl | 32 +++++----- .../src/OrdinaryDiffEqExtrapolation.jl | 2 +- 4 files changed, 57 insertions(+), 43 deletions(-) diff --git a/lib/OrdinaryDiffEqCore/src/alg_utils.jl b/lib/OrdinaryDiffEqCore/src/alg_utils.jl index 8f918da0e1..d773fc1345 100644 --- a/lib/OrdinaryDiffEqCore/src/alg_utils.jl +++ b/lib/OrdinaryDiffEqCore/src/alg_utils.jl @@ -309,13 +309,13 @@ alg_adaptive_order(alg::Union{OrdinaryDiffEqAlgorithm, DAEAlgorithm}) = alg_orde # this is actually incorrect and is purposefully decreased as this tends # to track the real error much better -function default_controller_v7(alg) +function default_controller_v7(QT, alg) if ispredictive(alg) - return PredictiveController(alg) + return PredictiveController(QT, alg) elseif isstandard(alg) - return IController(alg) + return IController(QT, alg) else - return PIController(alg) + return PIController(QT, alg) end end diff --git a/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl b/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl index 6cdce66042..c4063db42c 100644 --- a/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl +++ b/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl @@ -4,7 +4,7 @@ abstract type AbstractLegacyController <: AbstractController end abstract type AbstractControllerCache end # The legacy controllers do not have this concept. -setup_controller_cache(alg, cache, controller::AbstractLegacyController) = controller +setup_controller_cache(alg_cache, controller::AbstractLegacyController) = controller # checks whether the controller should accept a step based on the error estimate @inline function accept_step_controller(integrator, controller::Union{<:AbstractLegacyController, <:AbstractControllerCache}) @@ -113,8 +113,12 @@ struct IController{T} <: AbstractController qsteady_max::T end -function IController(alg; qmin = nothing, qmax = nothing, gamma = nothing, qsteady_min = nothing, qsteady_max = nothing) - IController( +function IController(alg; kwargs...) + IController(Float64, alg; kwargs...) +end + +function IController(QT, alg; qmin = nothing, qmax = nothing, gamma = nothing, qsteady_min = nothing, qsteady_max = nothing) + IController{QT}( qmin === nothing ? qmin_default(alg) : qmin, qmax === nothing ? qmax_default(alg) : qmax, gamma === nothing ? gamma_default(alg) : gamma, @@ -131,11 +135,11 @@ mutable struct IControllerCache{C, T} <: AbstractControllerCache # EEst::T end -function setup_controller_cache(alg, cache, controller::IController) +function setup_controller_cache(alg_cache, controller::IController{T}) where T IControllerCache( controller, - 1 // 1, - 1 // 10^4,# TODO which value? + T(1), + T(1 // 10^4), # TODO which value? ) end @@ -268,11 +272,15 @@ struct PIController{T} <: AbstractController qoldinit::T end -function PIController(alg; beta1 = nothing, beta2 = nothing, qmin = nothing, qmax = nothing, gamma = nothing, qsteady_min = nothing, qsteady_max = nothing, qoldinit = nothing) +function PIController(alg; kwargs...) + PIController(Float64, alg; kwargs...) +end + +function PIController(QT, alg; beta1 = nothing, beta2 = nothing, qmin = nothing, qmax = nothing, gamma = nothing, qsteady_min = nothing, qsteady_max = nothing, qoldinit = nothing) beta2 = beta2 === nothing ? beta2_default(alg) : beta2 beta1 = beta1 === nothing ? beta1_default(alg, beta2) : beta1 qoldinit = qoldinit === nothing ? 1 // 10^4 : qoldinit - PIController{promote_type(typeof(qoldinit), typeof(beta1))}( + PIController{QT}( beta1, beta2, qmin === nothing ? qmin_default(alg) : qmin, @@ -292,12 +300,12 @@ mutable struct PIControllerCache{T} <: AbstractControllerCache end -function setup_controller_cache(alg, cache, controller::PIController) +function setup_controller_cache(alg_cache, controller::PIController{T}) where T PIControllerCache( controller, - 1 // 1, - 1 // 10^4,# TODO which value? - 1 // 10^4, + T(1), + T(1 // 10^4),# TODO which value? + T(1 // 10^4), ) end @@ -522,12 +530,16 @@ struct PIDController{T, Limiter} <: AbstractController qoldinit::T end -function PIDController(alg; beta1 = nothing, beta2 = nothing, beta3 = nothing, accept_safety = 0.81, limiter = default_dt_factor_limiter, qmin = nothing, qmax = nothing, gamma = nothing, qsteady_min = nothing, qsteady_max = nothing, qoldinit = nothing) +function PIDController(alg; kwargs...) + PIDController(Float64, alg; kwargs...) +end + +function PIDController(QT, alg; beta1 = nothing, beta2 = nothing, beta3 = nothing, accept_safety = 0.81, limiter = default_dt_factor_limiter, qmin = nothing, qmax = nothing, gamma = nothing, qsteady_min = nothing, qsteady_max = nothing, qoldinit = nothing) beta2 = beta2 === nothing ? beta2_default(alg) : beta2 beta1 = beta1 === nothing ? beta1_default(alg, beta2) : beta1 beta3 = beta3 === nothing ? zero(beta1) : beta3 beta = MVector(map(float, promote(beta1, beta2, beta3))...) - PIDController( + PIDController{QT}( beta, err, accept_safety, @@ -554,7 +566,7 @@ struct PIDControllerCache{T, Limiter} <: AbstractControllerCache dt_factor::T end -function setup_controller_cache(alg, cache, controller::PIDController{QT}) where QT +function setup_controller_cache(alg_cache, controller::PIDController{QT}) where QT err = MVector{3, QT}(true, true, true) PIControllerCache( controller, @@ -752,8 +764,12 @@ struct PredictiveController{T} <: AbstractController qsteady_max::T end -function PredictiveController(alg; qmin = nothing, qmax = nothing, gamma = nothing, qsteady_min = nothing, qsteady_max = nothing) - PredictiveController( +function PredictiveController(alg; kwargs...) + PIController(Float64, alg; kwargs...) +end + +function PredictiveController(QT, alg; qmin = nothing, qmax = nothing, gamma = nothing, qsteady_min = nothing, qsteady_max = nothing) + PredictiveController{QT}( qmin === nothing ? qmin_default(alg) : qmin, qmax === nothing ? qmax_default(alg) : qmax, gamma === nothing ? gamma_default(alg) : gamma, @@ -768,11 +784,11 @@ mutable struct PredictiveControllerCache{T} <: AbstractControllerCache erracc::T end -function setup_controller_cache(alg, cache, controller::PredictiveController) - PredictiveControllerCache( +function setup_controller_cache(alg_cache, controller::PredictiveController{T}) where T + PredictiveControllerCache{T}( controller, - 1 // 1, - 1 // 1, + T(1), + T(1), ) end diff --git a/lib/OrdinaryDiffEqCore/src/solve.jl b/lib/OrdinaryDiffEqCore/src/solve.jl index a06abb05fd..582822527c 100644 --- a/lib/OrdinaryDiffEqCore/src/solve.jl +++ b/lib/OrdinaryDiffEqCore/src/solve.jl @@ -46,7 +46,7 @@ function SciMLBase.__init( beta1 = nothing, beta2 = nothing, qoldinit = nothing, - controller = any((gamma, qmin, qmax, qsteady_min, qsteady_max, beta1, beta2, qoldinit) .!== nothing) ? nothing : default_controller_v7(alg), # We have to reconstruct the old controller before breaking release., + controller = any((gamma, qmin, qmax, qsteady_min, qsteady_max, beta1, beta2, qoldinit) .!== nothing) ? nothing : default_controller_v7(typeof(one(eltype(prob.tspan))), alg), # We have to reconstruct the old controller before breaking release., fullnormalize = true, failfactor = 2, maxiters = anyadaptive(alg) ? 1000000 : typemax(Int), @@ -443,26 +443,24 @@ function SciMLBase.__init( legacy_controller_parameters = (gamma, qmin, qmax, qsteady_min, qsteady_max, beta1, beta2, qoldinit) if controller === nothing # We have to reconstruct the old controller before breaking release. if any(legacy_controller_parameters .== nothing) - gamma === nothing ? gamma_default(alg) : gamma - qmin === nothing ? qmin_default(alg) : qmin - qmax === nothing ? qmax_default(alg) : qmax - qsteady_min === nothing ? qsteady_min_default(alg) : qsteady_min - qsteady_max === nothing ? qsteady_max_default(alg) : qsteady_max - qoldinit === nothing ? (anyadaptive(alg) ? 1 // 10^4 : 0) : qoldinit + gamma = gamma === nothing ? gamma_default(alg) : gamma + qmin = qmin === nothing ? qmin_default(alg) : qmin + qmax = qmax === nothing ? qmax_default(alg) : qmax + qsteady_min = qsteady_min === nothing ? qsteady_min_default(alg) : qsteady_min + qsteady_max = qsteady_max === nothing ? qsteady_max_default(alg) : qsteady_max + qoldinit = qoldinit === nothing ? (anyadaptive(alg) ? 1 // 10^4 : 0) : qoldinit controller = legacy_default_controller(_alg, cache, qoldinit, beta1, beta2) end - elseif controller isa AbstractLegacyController # Legacy controller has been passed - gamma === nothing ? gamma_default(alg) : gamma - qmin === nothing ? qmin_default(alg) : qmin - qmax === nothing ? qmax_default(alg) : qmax - qsteady_min === nothing ? qsteady_min_default(alg) : qsteady_min - qsteady_max === nothing ? qsteady_max_default(alg) : qsteady_max - qoldinit === nothing ? (anyadaptive(alg) ? 1 // 10^4 : 0) : qoldinit - else - @unpack gamma, qmin, qmax, qsteady_min, qsteady_max, qoldinit = controller + else # Controller has been passed + gamma = hasfield(typeof(controller), :gamma) ? controller.gamma : gamma_default(alg) + qmin = hasfield(typeof(controller), :qmin) ? controller.qmin : qmin_default(alg) + qmax = hasfield(typeof(controller), :qmax) ? controller.qmax : qmax_default(alg) + qsteady_min = hasfield(typeof(controller), :qsteady_min) ? controller.qsteady_min : qsteady_min_default(alg) + qsteady_max = hasfield(typeof(controller), :qsteady_max) ? controller.qsteady_max : qsteady_max_default(alg) + qoldinit = hasfield(typeof(controller), :qoldinit) ? controller.qoldinit : (anyadaptive(alg) ? 1 // 10^4 : 0) end - controller_cache = setup_controller_cache(_alg, cache, controller) + controller_cache = setup_controller_cache(cache, controller) save_end_user = save_end save_end = save_end === nothing ? diff --git a/lib/OrdinaryDiffEqExtrapolation/src/OrdinaryDiffEqExtrapolation.jl b/lib/OrdinaryDiffEqExtrapolation/src/OrdinaryDiffEqExtrapolation.jl index befa7152fc..0d3c18a3a9 100644 --- a/lib/OrdinaryDiffEqExtrapolation/src/OrdinaryDiffEqExtrapolation.jl +++ b/lib/OrdinaryDiffEqExtrapolation/src/OrdinaryDiffEqExtrapolation.jl @@ -8,7 +8,7 @@ import OrdinaryDiffEqCore: alg_order, alg_maximum_order, get_current_adaptive_or isthreaded, step_accept_controller!, calculate_residuals, OrdinaryDiffEqMutableCache, OrdinaryDiffEqConstantCache, - reset_alg_dependent_opts!, AbstractController, + reset_alg_dependent_opts!, AbstractController, AbstractLegacyController, step_accept_controller!, step_reject_controller!, OrdinaryDiffEqAdaptiveAlgorithm, OrdinaryDiffEqAdaptiveImplicitAlgorithm, From 169c09abf8bac2530afa9a264dc42513efc5be7d Mon Sep 17 00:00:00 2001 From: termi-official <9196588+termi-official@users.noreply.github.com> Date: Sun, 26 Oct 2025 19:31:42 +0100 Subject: [PATCH 04/21] [skip ci] typo --- lib/OrdinaryDiffEqCore/src/solve.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/OrdinaryDiffEqCore/src/solve.jl b/lib/OrdinaryDiffEqCore/src/solve.jl index 582822527c..dffb8fddd0 100644 --- a/lib/OrdinaryDiffEqCore/src/solve.jl +++ b/lib/OrdinaryDiffEqCore/src/solve.jl @@ -439,7 +439,7 @@ function SciMLBase.__init( Base.depwarn(message, :solve) end - # The folowing code provides an upgrade path for users by preserving the old behavior. + # The following code provides an upgrade path for users by preserving the old behavior. legacy_controller_parameters = (gamma, qmin, qmax, qsteady_min, qsteady_max, beta1, beta2, qoldinit) if controller === nothing # We have to reconstruct the old controller before breaking release. if any(legacy_controller_parameters .== nothing) From f17a7a3f04c4bb1ee9ff96583f035ea733a02de8 Mon Sep 17 00:00:00 2001 From: termi-official <9196588+termi-official@users.noreply.github.com> Date: Mon, 27 Oct 2025 11:02:30 +0100 Subject: [PATCH 05/21] Fix some issues with recirpocals in the new implementation --- .../src/integrators/controllers.jl | 39 ++++++++++--------- 1 file changed, 21 insertions(+), 18 deletions(-) diff --git a/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl b/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl index c4063db42c..9313c8999f 100644 --- a/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl +++ b/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl @@ -294,9 +294,12 @@ end mutable struct PIControllerCache{T} <: AbstractControllerCache controller::PIController{T} + # Propsoed scaling factor for the time step length q::T + # Cached εₙ₊₁^β₁ q11::T - qold::T + # Previous EEst + errold::T end @@ -304,25 +307,25 @@ function setup_controller_cache(alg_cache, controller::PIController{T}) where T PIControllerCache( controller, T(1), - T(1 // 10^4),# TODO which value? + T(1), T(1 // 10^4), ) end @inline function stepsize_controller!(integrator, cache::PIControllerCache, alg) - @unpack controller = cache - @unpack qold = integrator + @unpack errold, controller = cache @unpack qmin, qmax, gamma = controller @unpack beta1, beta2 = controller EEst = DiffEqBase.value(integrator.EEst) if iszero(EEst) - q = inv(qmax) + q = qmax else - q11 = fastpower(EEst, convert(typeof(EEst), beta1)) - q = q11 / fastpower(qold, convert(typeof(EEst), beta2)) - integrator.q11 = q11 - @fastmath q = max(inv(qmax), min(inv(qmin), q / gamma)) + q11 = fastpower(EEst, -beta1) + q = q11 * fastpower(errold, -beta2) + cache.q11 = q11 + integrator.q11 = q11 # TODO remove + @fastmath q = clamp(q * gamma, qmin, qmax) end cache.q = q q @@ -337,15 +340,15 @@ function step_accept_controller!(integrator, cache::PIControllerCache, alg, q) if qsteady_min <= q <= qsteady_max q = one(q) end - cache.qold = max(EEst, qoldinit) - return integrator.dt / q # new dt + cache.errold = max(EEst, qoldinit) + integrator.qold = cache.errold # TODO remove + return integrator.dt * q # new dt end function step_reject_controller!(integrator, cache::PIControllerCache, alg) - @unpack controller = cache - @unpack q11 = cache + @unpack controller, q11 = cache @unpack qmin, gamma = controller - integrator.dt /= min(inv(qmin), q11 / gamma) + integrator.dt *= max(qmin, q11 * gamma) end # FIXME multi-controller? @@ -568,7 +571,7 @@ end function setup_controller_cache(alg_cache, controller::PIDController{QT}) where QT err = MVector{3, QT}(true, true, true) - PIControllerCache( + PIDControllerCache( controller, err, QT(1 // 10^4), @@ -616,16 +619,16 @@ end # ensures # 0.21 ≈ limiter(0) <= dt_factor <= limiter(Inf) ≈ 2.57 # See Söderlind, Wang (2006), Section 6. - cache.qold = dt_factor + cache.dt_factor = dt_factor return dt_factor end @inline function accept_step_controller(integrator, cache::PIDControllerCache) - return cache.qold >= cache.controller.accept_safety + return cache.dt_factor >= cache.controller.accept_safety end function step_accept_controller!(integrator, cache::PIDControllerCache, alg, dt_factor) - @assert dt_factor ≈ cache.q "Controller cache went out of sync with time stepping logic." + @assert dt_factor ≈ cache.dt_factor "Controller cache went out of sync with time stepping logic." @unpack controller = cache @unpack qsteady_min, qsteady_max = controller From dd3d7bf06b3c02eb1cb8c67a0d4c21dfe10f5785 Mon Sep 17 00:00:00 2001 From: termi-official <9196588+termi-official@users.noreply.github.com> Date: Mon, 27 Oct 2025 12:29:22 +0100 Subject: [PATCH 06/21] ifelse -> max --- lib/OrdinaryDiffEqCore/src/integrators/controllers.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl b/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl index 9313c8999f..cd000b9f7c 100644 --- a/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl +++ b/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl @@ -475,7 +475,7 @@ end # EEst = eps(typeof(EEst)) # end # ``` - EEst = ifelse(EEst > EEst_min, EEst, EEst_min) + EEst = max(EEst, EEst_min) controller.err[1] = inv(EEst) err1, err2, err3 = controller.err @@ -601,7 +601,7 @@ end # EEst = eps(typeof(EEst)) # end # ``` - EEst = ifelse(EEst > EEst_min, EEst, EEst_min) + EEst = max(EEst, EEst_min) cache.err[1] = inv(EEst) err1, err2, err3 = cache.err From 8e0e9793c7cd7fe06c82353836400ba2d39358ec Mon Sep 17 00:00:00 2001 From: termi-official Date: Mon, 27 Oct 2025 13:20:09 +0100 Subject: [PATCH 07/21] Reverse removing reciprocals --- lib/OrdinaryDiffEqCore/src/integrators/controllers.jl | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl b/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl index cd000b9f7c..95e26b8329 100644 --- a/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl +++ b/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl @@ -321,11 +321,12 @@ end if iszero(EEst) q = qmax else - q11 = fastpower(EEst, -beta1) - q = q11 * fastpower(errold, -beta2) + # Legacy code + q11 = fastpower(EEst, beta1) + q = q11 / fastpower(errold, beta2) cache.q11 = q11 integrator.q11 = q11 # TODO remove - @fastmath q = clamp(q * gamma, qmin, qmax) + @fastmath q = clamp(q / gamma, inv(qmax), inv(qmin)) end cache.q = q q @@ -342,13 +343,13 @@ function step_accept_controller!(integrator, cache::PIControllerCache, alg, q) end cache.errold = max(EEst, qoldinit) integrator.qold = cache.errold # TODO remove - return integrator.dt * q # new dt + return integrator.dt / q # new dt end function step_reject_controller!(integrator, cache::PIControllerCache, alg) @unpack controller, q11 = cache @unpack qmin, gamma = controller - integrator.dt *= max(qmin, q11 * gamma) + integrator.dt /= min(inv(qmin), q11 / gamma) end # FIXME multi-controller? From 9bc55b00d5887fabf0a734bde2109f62cc6cfea0 Mon Sep 17 00:00:00 2001 From: termi-official Date: Mon, 27 Oct 2025 13:34:36 +0100 Subject: [PATCH 08/21] Add new PIDController to los storage methods --- .../src/OrdinaryDiffEqLowStorageRK.jl | 3 ++- .../src/alg_utils.jl | 24 +++++++++++++++++++ 2 files changed, 26 insertions(+), 1 deletion(-) diff --git a/lib/OrdinaryDiffEqLowStorageRK/src/OrdinaryDiffEqLowStorageRK.jl b/lib/OrdinaryDiffEqLowStorageRK/src/OrdinaryDiffEqLowStorageRK.jl index 315938b601..e1f9ebc192 100644 --- a/lib/OrdinaryDiffEqLowStorageRK/src/OrdinaryDiffEqLowStorageRK.jl +++ b/lib/OrdinaryDiffEqLowStorageRK/src/OrdinaryDiffEqLowStorageRK.jl @@ -7,7 +7,8 @@ import OrdinaryDiffEqCore: alg_order, alg_adaptive_order, calculate_residuals!, OrdinaryDiffEqAlgorithm, ispredictive, OrdinaryDiffEqMutableCache, OrdinaryDiffEqConstantCache, OrdinaryDiffEqAdaptiveAlgorithm, uses_uprev, - legacy_default_controller, LegacyPIDController, + default_controller_v7, + legacy_default_controller, PIDController, alg_cache, _vec, _reshape, @cache, isfsal, full_cache, constvalue, _unwrap_val, trivial_limiter!, perform_step!, initialize!, diff --git a/lib/OrdinaryDiffEqLowStorageRK/src/alg_utils.jl b/lib/OrdinaryDiffEqLowStorageRK/src/alg_utils.jl index 064403627f..1a42d81c4a 100644 --- a/lib/OrdinaryDiffEqLowStorageRK/src/alg_utils.jl +++ b/lib/OrdinaryDiffEqLowStorageRK/src/alg_utils.jl @@ -113,3 +113,27 @@ function legacy_default_controller(alg::RDPK3SpFSAL510, cache, qoldinit, args... QT = typeof(qoldinit) return LegacyPIDController(map(Base.Fix1(convert, QT), (0.45, -0.13, 0.00))...) end + +function default_controller_v7(QT, alg::RDPK3Sp35) + return PIDController(map(Base.Fix1(convert, QT), (0.64, -0.31, 0.04))...) +end + +function default_controller_v7(QT, alg::RDPK3SpFSAL35) + return PIDController(map(Base.Fix1(convert, QT), (0.70, -0.23, 0.00))...) +end + +function default_controller_v7(QT, alg::RDPK3Sp49) + return PIDController(map(Base.Fix1(convert, QT), (0.25, -0.12, 0.00))...) +end + +function default_controller_v7(QT, alg::RDPK3SpFSAL49) + return PIDController(map(Base.Fix1(convert, QT), (0.38, -0.18, 0.01))...) +end + +function default_controller_v7(QT, alg::RDPK3Sp510) + return PIDController(map(Base.Fix1(convert, QT), (0.47, -0.20, 0.06))...) +end + +function default_controller_v7(QT, alg::RDPK3SpFSAL510) + return PIDController(map(Base.Fix1(convert, QT), (0.45, -0.13, 0.00))...) +end From 64ae15bed0c3d6640ea0ad894a496a9b46ab6f99 Mon Sep 17 00:00:00 2001 From: termi-official Date: Mon, 27 Oct 2025 13:52:43 +0100 Subject: [PATCH 09/21] Update BDF module --- lib/OrdinaryDiffEqBDF/src/OrdinaryDiffEqBDF.jl | 1 + lib/OrdinaryDiffEqBDF/src/controllers.jl | 4 ++++ lib/OrdinaryDiffEqCore/src/integrators/controllers.jl | 9 ++++++--- 3 files changed, 11 insertions(+), 3 deletions(-) diff --git a/lib/OrdinaryDiffEqBDF/src/OrdinaryDiffEqBDF.jl b/lib/OrdinaryDiffEqBDF/src/OrdinaryDiffEqBDF.jl index 9e99524d43..3910a77c9f 100644 --- a/lib/OrdinaryDiffEqBDF/src/OrdinaryDiffEqBDF.jl +++ b/lib/OrdinaryDiffEqBDF/src/OrdinaryDiffEqBDF.jl @@ -16,6 +16,7 @@ import OrdinaryDiffEqCore: alg_order, calculate_residuals!, trivial_limiter!, issplit, qsteady_min_default, qsteady_max_default, get_current_alg_order, get_current_adaptive_order, + default_controller_v7, legacy_default_controller, stepsize_controller!, step_accept_controller!, step_reject_controller!, post_newton_controller!, diff --git a/lib/OrdinaryDiffEqBDF/src/controllers.jl b/lib/OrdinaryDiffEqBDF/src/controllers.jl index 57be463133..fd46c2b02e 100644 --- a/lib/OrdinaryDiffEqBDF/src/controllers.jl +++ b/lib/OrdinaryDiffEqBDF/src/controllers.jl @@ -2,6 +2,10 @@ function legacy_default_controller(alg::Union{QNDF, FBDF}, args...) DummyController() end +function default_controller_v7(QT, alg::Union{QNDF, FBDF}, args...) + DummyController() +end + # QNBDF stepsize_controller!(integrator, alg::QNDF) = nothing diff --git a/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl b/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl index 95e26b8329..2ab453bdbb 100644 --- a/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl +++ b/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl @@ -38,6 +38,12 @@ function post_newton_controller!(integrator, alg) nothing end +# This is a helper for Nordsieck and BDF methods, which come with an integrated controller. +struct DummyController <: AbstractController +end + +setup_controller_cache(alg_cache, controller::DummyController) = controller + # Standard integral (I) step size controller """ LegacyIController() @@ -72,9 +78,6 @@ the predicted step size. struct LegacyIController <: AbstractLegacyController end -struct DummyController <: AbstractController -end - @inline function stepsize_controller!(integrator, controller::LegacyIController, alg) @unpack qmin, qmax, gamma = integrator.opts EEst = DiffEqBase.value(integrator.EEst) From 07e65d4831f74b4a022534ffa760622a8572f7ae Mon Sep 17 00:00:00 2001 From: termi-official Date: Mon, 27 Oct 2025 14:08:03 +0100 Subject: [PATCH 10/21] Update Nordsieck methods --- .../src/OrdinaryDiffEqNordsieck.jl | 5 +++-- lib/OrdinaryDiffEqNordsieck/src/alg_utils.jl | 6 +++++- lib/OrdinaryDiffEqNordsieck/src/algorithms.jl | 9 ++++++--- lib/OrdinaryDiffEqNordsieck/src/controllers.jl | 11 ++++------- lib/OrdinaryDiffEqNordsieck/src/nordsieck_caches.jl | 6 ++++-- 5 files changed, 22 insertions(+), 15 deletions(-) diff --git a/lib/OrdinaryDiffEqNordsieck/src/OrdinaryDiffEqNordsieck.jl b/lib/OrdinaryDiffEqNordsieck/src/OrdinaryDiffEqNordsieck.jl index 5db6344b12..79c871d576 100644 --- a/lib/OrdinaryDiffEqNordsieck/src/OrdinaryDiffEqNordsieck.jl +++ b/lib/OrdinaryDiffEqNordsieck/src/OrdinaryDiffEqNordsieck.jl @@ -1,7 +1,7 @@ module OrdinaryDiffEqNordsieck import OrdinaryDiffEqCore: alg_order, alg_adaptive_order, qsteady_max_default, - get_current_alg_order, + get_current_alg_order, DummyController, AbstractController, OrdinaryDiffEqAdaptiveAlgorithm, OrdinaryDiffEqAdamsVarOrderVarStepAlgorithm, alg_cache, OrdinaryDiffEqMutableCache, @@ -11,7 +11,8 @@ import OrdinaryDiffEqCore: alg_order, alg_adaptive_order, qsteady_max_default, calculate_residuals, calculate_residuals!, get_current_adaptive_order, get_fsalfirstlast, ode_interpolant, ode_interpolant!, trivial_limiter!, - generic_solver_docstring + generic_solver_docstring, + default_controller_v7, legacy_default_controller using MuladdMacro, FastBroadcast, RecursiveArrayTools import LinearAlgebra: rmul! import Static: False diff --git a/lib/OrdinaryDiffEqNordsieck/src/alg_utils.jl b/lib/OrdinaryDiffEqNordsieck/src/alg_utils.jl index 067f482539..f09630bf0c 100644 --- a/lib/OrdinaryDiffEqNordsieck/src/alg_utils.jl +++ b/lib/OrdinaryDiffEqNordsieck/src/alg_utils.jl @@ -8,6 +8,10 @@ qsteady_max_default(alg::JVODE) = 3 // 2 get_current_alg_order(alg::JVODE, cache) = get_current_adaptive_order(alg, cache) -function legacy_default_controller(alg::Union{JVODE}, args...) +function legacy_default_controller(alg::JVODE, args...) + DummyController() +end + +function default_controller_v7(QT, alg::JVODE, args...) DummyController() end diff --git a/lib/OrdinaryDiffEqNordsieck/src/algorithms.jl b/lib/OrdinaryDiffEqNordsieck/src/algorithms.jl index c8e49d97a7..11fca50bfc 100644 --- a/lib/OrdinaryDiffEqNordsieck/src/algorithms.jl +++ b/lib/OrdinaryDiffEqNordsieck/src/algorithms.jl @@ -15,17 +15,20 @@ struct AN5 <: OrdinaryDiffEqAdaptiveAlgorithm end `JVODE` is experimental, the solver `VCABM` is generally preferred. """ -struct JVODE{bType, aType} <: OrdinaryDiffEqAdamsVarOrderVarStepAlgorithm +struct JVODE{bType, aType, qType} <: OrdinaryDiffEqAdamsVarOrderVarStepAlgorithm algorithm::Symbol bias1::bType bias2::bType bias3::bType addon::aType + qmax::qType + qsteady_min::qType + qsteady_max::qType end function JVODE(algorithm = :Adams; bias1 = 6, bias2 = 6, bias3 = 10, - addon = 1 // 10^6) - JVODE(algorithm, bias1, bias2, bias3, addon) + addon = 1 // 10^6, qmax = float(10), qsteady_min = float(1), qsteady_max = float(3 // 2)) + JVODE(algorithm, bias1, bias2, bias3, addon, qmax, qsteady_min, qsteady_max) end """ !!! warning "Experimental" diff --git a/lib/OrdinaryDiffEqNordsieck/src/controllers.jl b/lib/OrdinaryDiffEqNordsieck/src/controllers.jl index 216eb36168..cd6745b67e 100644 --- a/lib/OrdinaryDiffEqNordsieck/src/controllers.jl +++ b/lib/OrdinaryDiffEqNordsieck/src/controllers.jl @@ -1,25 +1,22 @@ -struct DummyController <: AbstractController -end - # JVODE function stepsize_controller!(integrator, alg::JVODE) if iszero(integrator.EEst) - η = integrator.opts.qmax + η = alg.qmax else η = integrator.cache.η - integrator.qold = η + integrator.cache.ηold = η end η end function step_accept_controller!(integrator, alg::JVODE, η) q = inv(η) - if q <= integrator.opts.qsteady_max && q >= integrator.opts.qsteady_min + if q <= alg.qsteady_max && q >= alg.qsteady_min q = one(q) end return integrator.dt / q # dtnew end function step_reject_controller!(integrator, alg::JVODE) - integrator.dt *= integrator.qold + integrator.dt *= integrator.cache.η end diff --git a/lib/OrdinaryDiffEqNordsieck/src/nordsieck_caches.jl b/lib/OrdinaryDiffEqNordsieck/src/nordsieck_caches.jl index a9ecb45215..4c72b0f126 100644 --- a/lib/OrdinaryDiffEqNordsieck/src/nordsieck_caches.jl +++ b/lib/OrdinaryDiffEqNordsieck/src/nordsieck_caches.jl @@ -133,6 +133,7 @@ mutable struct JVODEConstantCache{zType, lType, dtsType, dType, tsit5Type, etaTy n_wait::Int # `η` is `dtₙ₊₁/dtₙ` η::etaType + ηold::etaType ηq::etaType η₊₁::etaType η₋₁::etaType @@ -154,7 +155,7 @@ function alg_cache(alg::JVODE, u, rate_prototype, ::Type{uEltypeNoUnits}, η = zero(dt / dt) JVODEConstantCache(z, l, m, c_LTE₊₁, c_LTE, c_LTE₋₁, c_conv, c_𝒟, prev_𝒟, - dts, Δ, tsit5tab, 2, 1, 1, 2, η, η, η, η, η) + dts, Δ, tsit5tab, 2, 1, 1, 2, η, η, η, η, η, η) end mutable struct JVODECache{ @@ -205,6 +206,7 @@ mutable struct JVODECache{ n_wait::Int # `η` is `dtₙ₊₁/dtₙ` η::etaType + ηold::etaType ηq::etaType η₊₁::etaType η₋₁::etaType @@ -262,7 +264,7 @@ function alg_cache(alg::JVODE, u, rate_prototype, ::Type{uEltypeNoUnits}, JVODECache(u, uprev, tmp, fsalfirst, ratetmp, z, l, m, c_LTE₊₁, c_LTE, c_LTE₋₁, c_conv, c_𝒟, prev_𝒟, - dts, Δ, atmp, tsit5cache, 2, 1, 1, 2, η, η, η, η, η) + dts, Δ, atmp, tsit5cache, 2, 1, 1, 2, η, η, η, η, η, η) end function get_fsalfirstlast(cache::Union{JVODECache, AN5Cache}, u) From f77aef9eca85da3f7d4ee472fa43c9a83e904336 Mon Sep 17 00:00:00 2001 From: termi-official Date: Mon, 27 Oct 2025 16:33:34 +0100 Subject: [PATCH 11/21] Fix low storage sublib --- .../src/integrators/controllers.jl | 41 ++++++++++--------- .../src/OrdinaryDiffEqLowStorageRK.jl | 2 +- .../src/alg_utils.jl | 12 +++--- 3 files changed, 29 insertions(+), 26 deletions(-) diff --git a/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl b/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl index 2ab453bdbb..b5d9ace551 100644 --- a/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl +++ b/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl @@ -44,6 +44,10 @@ end setup_controller_cache(alg_cache, controller::DummyController) = controller +@inline function accept_step_controller(integrator, controller::DummyController) + return integrator.EEst <= 1 +end + # Standard integral (I) step size controller """ LegacyIController() @@ -531,32 +535,31 @@ struct PIDController{T, Limiter} <: AbstractController limiter::Limiter # limiter of the dt factor (before clipping) qmin::T qmax::T - gamma::T qsteady_min::T qsteady_max::T - qoldinit::T end function PIDController(alg; kwargs...) PIDController(Float64, alg; kwargs...) end -function PIDController(QT, alg; beta1 = nothing, beta2 = nothing, beta3 = nothing, accept_safety = 0.81, limiter = default_dt_factor_limiter, qmin = nothing, qmax = nothing, gamma = nothing, qsteady_min = nothing, qsteady_max = nothing, qoldinit = nothing) - beta2 = beta2 === nothing ? beta2_default(alg) : beta2 - beta1 = beta1 === nothing ? beta1_default(alg, beta2) : beta1 - beta3 = beta3 === nothing ? zero(beta1) : beta3 - beta = MVector(map(float, promote(beta1, beta2, beta3))...) - PIDController{QT}( +function PIDController(QT, alg; beta = nothing, accept_safety = 0.81, limiter = default_dt_factor_limiter, qmin = nothing, qmax = nothing, qsteady_min = nothing, qsteady_max = nothing) + if beta === nothing + beta2 = QT(beta2_default(alg)) + beta1 = QT(beta1_default(alg, beta2)) + beta3 = QT(zero(beta1)) + else + beta1, beta2, beta3 = beta + end + beta = SVector(map(float, promote(beta1, beta2, beta3))...) + PIDController{QT, typeof(limiter)}( beta, - err, - accept_safety, + QT(accept_safety), limiter, - qmin === nothing ? qmin_default(alg) : qmin, - qmax === nothing ? qmax_default(alg) : qmax, - gamma === nothing ? gamma_default(alg) : gamma, - qsteady_min === nothing ? qsteady_min_default(alg) : qsteady_min, - qsteady_max === nothing ? qsteady_max_default(alg) : qsteady_max, - qoldinit === nothing ? 1 // 10^4 : qoldinit, + QT(qmin === nothing ? qmin_default(alg) : qmin), + QT(qmax === nothing ? qmax_default(alg) : qmax), + QT(qsteady_min === nothing ? qsteady_min_default(alg) : qsteady_min), + QT(qsteady_max === nothing ? qsteady_max_default(alg) : qsteady_max), ) end @@ -567,7 +570,7 @@ function Base.show(io::IO, controller::PIDController) ")") end -struct PIDControllerCache{T, Limiter} <: AbstractControllerCache +mutable struct PIDControllerCache{T, Limiter} <: AbstractControllerCache controller::PIDController{T, Limiter} err::MVector{3, T} # history of the error estimates dt_factor::T @@ -641,8 +644,8 @@ function step_accept_controller!(integrator, cache::PIDControllerCache, alg, dt_ # cache.q = ...? end @inbounds begin - controller.err[3] = controller.err[2] - controller.err[2] = controller.err[1] + cache.err[3] = cache.err[2] + cache.err[2] = cache.err[1] end return integrator.dt * dt_factor # new dt end diff --git a/lib/OrdinaryDiffEqLowStorageRK/src/OrdinaryDiffEqLowStorageRK.jl b/lib/OrdinaryDiffEqLowStorageRK/src/OrdinaryDiffEqLowStorageRK.jl index e1f9ebc192..416f52efd3 100644 --- a/lib/OrdinaryDiffEqLowStorageRK/src/OrdinaryDiffEqLowStorageRK.jl +++ b/lib/OrdinaryDiffEqLowStorageRK/src/OrdinaryDiffEqLowStorageRK.jl @@ -8,7 +8,7 @@ import OrdinaryDiffEqCore: alg_order, alg_adaptive_order, calculate_residuals!, OrdinaryDiffEqMutableCache, OrdinaryDiffEqConstantCache, OrdinaryDiffEqAdaptiveAlgorithm, uses_uprev, default_controller_v7, - legacy_default_controller, PIDController, + legacy_default_controller, PIDController, LegacyPIDController, alg_cache, _vec, _reshape, @cache, isfsal, full_cache, constvalue, _unwrap_val, trivial_limiter!, perform_step!, initialize!, diff --git a/lib/OrdinaryDiffEqLowStorageRK/src/alg_utils.jl b/lib/OrdinaryDiffEqLowStorageRK/src/alg_utils.jl index 1a42d81c4a..e849eb6469 100644 --- a/lib/OrdinaryDiffEqLowStorageRK/src/alg_utils.jl +++ b/lib/OrdinaryDiffEqLowStorageRK/src/alg_utils.jl @@ -115,25 +115,25 @@ function legacy_default_controller(alg::RDPK3SpFSAL510, cache, qoldinit, args... end function default_controller_v7(QT, alg::RDPK3Sp35) - return PIDController(map(Base.Fix1(convert, QT), (0.64, -0.31, 0.04))...) + return PIDController(QT, alg; beta=map(Base.Fix1(convert, QT), (0.64, -0.31, 0.04))) end function default_controller_v7(QT, alg::RDPK3SpFSAL35) - return PIDController(map(Base.Fix1(convert, QT), (0.70, -0.23, 0.00))...) + return PIDController(QT, alg; beta=map(Base.Fix1(convert, QT), (0.70, -0.23, 0.00))) end function default_controller_v7(QT, alg::RDPK3Sp49) - return PIDController(map(Base.Fix1(convert, QT), (0.25, -0.12, 0.00))...) + return PIDController(QT, alg; beta=map(Base.Fix1(convert, QT), (0.25, -0.12, 0.00))) end function default_controller_v7(QT, alg::RDPK3SpFSAL49) - return PIDController(map(Base.Fix1(convert, QT), (0.38, -0.18, 0.01))...) + return PIDController(QT, alg; beta=map(Base.Fix1(convert, QT), (0.38, -0.18, 0.01))) end function default_controller_v7(QT, alg::RDPK3Sp510) - return PIDController(map(Base.Fix1(convert, QT), (0.47, -0.20, 0.06))...) + return PIDController(QT, alg; beta=map(Base.Fix1(convert, QT), (0.47, -0.20, 0.06))) end function default_controller_v7(QT, alg::RDPK3SpFSAL510) - return PIDController(map(Base.Fix1(convert, QT), (0.45, -0.13, 0.00))...) + return PIDController(QT, alg; beta=map(Base.Fix1(convert, QT), (0.45, -0.13, 0.00))) end From afd57ac38c00cb18812e0ce31d92da7d901c09cd Mon Sep 17 00:00:00 2001 From: termi-official Date: Mon, 27 Oct 2025 16:37:02 +0100 Subject: [PATCH 12/21] Recover behavior of Extrapolation solver in new interface --- .../src/OrdinaryDiffEqExtrapolation.jl | 4 +-- .../src/alg_utils.jl | 26 +++++++++++++++++++ 2 files changed, 28 insertions(+), 2 deletions(-) diff --git a/lib/OrdinaryDiffEqExtrapolation/src/OrdinaryDiffEqExtrapolation.jl b/lib/OrdinaryDiffEqExtrapolation/src/OrdinaryDiffEqExtrapolation.jl index 0d3c18a3a9..ec663f67b2 100644 --- a/lib/OrdinaryDiffEqExtrapolation/src/OrdinaryDiffEqExtrapolation.jl +++ b/lib/OrdinaryDiffEqExtrapolation/src/OrdinaryDiffEqExtrapolation.jl @@ -2,10 +2,10 @@ module OrdinaryDiffEqExtrapolation import OrdinaryDiffEqCore: alg_order, alg_maximum_order, get_current_adaptive_order, get_current_alg_order, calculate_residuals!, - accept_step_controller, + accept_step_controller, default_controller_v7, legacy_default_controller, beta2_default, beta1_default, gamma_default, initialize!, perform_step!, @unpack, @cache, unwrap_alg, - isthreaded, + isthreaded, isadaptive, LegacyPIController, step_accept_controller!, calculate_residuals, OrdinaryDiffEqMutableCache, OrdinaryDiffEqConstantCache, reset_alg_dependent_opts!, AbstractController, AbstractLegacyController, diff --git a/lib/OrdinaryDiffEqExtrapolation/src/alg_utils.jl b/lib/OrdinaryDiffEqExtrapolation/src/alg_utils.jl index c607dc6dee..7be2b336d7 100644 --- a/lib/OrdinaryDiffEqExtrapolation/src/alg_utils.jl +++ b/lib/OrdinaryDiffEqExtrapolation/src/alg_utils.jl @@ -46,6 +46,32 @@ function legacy_default_controller( return LegacyExtrapolationController(beta1) end +function default_controller_v7(QT, + alg::Union{ + ExtrapolationMidpointDeuflhard, + ImplicitDeuflhardExtrapolation, + ExtrapolationMidpointHairerWanner, + ImplicitHairerWannerExtrapolation, + ImplicitEulerExtrapolation, + ImplicitEulerBarycentricExtrapolation}) + beta1 = QT(beta1_default(alg, beta2_default(alg))) + return LegacyExtrapolationController(beta1) +end + +# FIXME AitkenNeville is missing integration with the extrapolation controller and picks up the PI controller instead. +function legacy_default_controller(alg::AitkenNeville, + cache, + qoldinit, _beta1 = nothing, _beta2 = nothing) + QT = typeof(qoldinit) + beta1, beta2 = _digest_beta1_beta2(alg, cache, Val(QT), _beta1, _beta2) + return LegacyPIController(beta1, beta2) +end +function default_controller_v7(QT, alg::AitkenNeville) + beta2 = QT(beta2_default(alg)) + beta1 = QT(beta1_default(alg, beta2)) + return LegacyPIController(beta1, beta2) +end + beta2_default(alg::ExtrapolationMidpointDeuflhard) = 0 // 1 beta2_default(alg::ImplicitDeuflhardExtrapolation) = 0 // 1 beta2_default(alg::ExtrapolationMidpointHairerWanner) = 0 // 1 From ca8c0788faa126e8c6f9abac32d8549ccba26c8b Mon Sep 17 00:00:00 2001 From: termi-official Date: Mon, 27 Oct 2025 16:46:18 +0100 Subject: [PATCH 13/21] Change name of legacy controllers back and add New suffix to the controllers in the new interface --- lib/OrdinaryDiffEqCore/src/alg_utils.jl | 12 +- .../src/integrators/controllers.jl | 120 +++++++++--------- .../src/OrdinaryDiffEqExtrapolation.jl | 2 +- .../src/alg_utils.jl | 8 +- .../src/controllers.jl | 4 +- .../src/OrdinaryDiffEqFIRK.jl | 4 +- lib/OrdinaryDiffEqFIRK/src/controllers.jl | 4 +- .../src/OrdinaryDiffEqLowStorageRK.jl | 2 +- .../src/alg_utils.jl | 24 ++-- 9 files changed, 91 insertions(+), 89 deletions(-) diff --git a/lib/OrdinaryDiffEqCore/src/alg_utils.jl b/lib/OrdinaryDiffEqCore/src/alg_utils.jl index d773fc1345..78a488aa1b 100644 --- a/lib/OrdinaryDiffEqCore/src/alg_utils.jl +++ b/lib/OrdinaryDiffEqCore/src/alg_utils.jl @@ -311,23 +311,23 @@ alg_adaptive_order(alg::Union{OrdinaryDiffEqAlgorithm, DAEAlgorithm}) = alg_orde function default_controller_v7(QT, alg) if ispredictive(alg) - return PredictiveController(QT, alg) + return NewPredictiveController(QT, alg) elseif isstandard(alg) - return IController(QT, alg) + return NewIController(QT, alg) else - return PIController(QT, alg) + return NewPIController(QT, alg) end end function legacy_default_controller(alg, cache, qoldinit, _beta1 = nothing, _beta2 = nothing) if ispredictive(alg) - return LegacyPredictiveController() + return PredictiveController() elseif isstandard(alg) - return LegacyIController() + return IController() else # Default is PI-controller QT = typeof(qoldinit) beta1, beta2 = _digest_beta1_beta2(alg, cache, Val(QT), _beta1, _beta2) - return LegacyPIController(beta1, beta2) + return PIController(beta1, beta2) end end diff --git a/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl b/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl index b5d9ace551..f72e68fb1b 100644 --- a/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl +++ b/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl @@ -7,7 +7,7 @@ abstract type AbstractControllerCache end setup_controller_cache(alg_cache, controller::AbstractLegacyController) = controller # checks whether the controller should accept a step based on the error estimate -@inline function accept_step_controller(integrator, controller::Union{<:AbstractLegacyController, <:AbstractControllerCache}) +@inline function accept_step_controller(integrator, controller_or_cache::Union{<:AbstractLegacyController, <:AbstractControllerCache}) return integrator.EEst <= 1 end @@ -50,7 +50,7 @@ end # Standard integral (I) step size controller """ - LegacyIController() + IController() The standard (integral) controller is the most basic step size controller. This controller is usually the first one introduced in numerical analysis classes @@ -79,10 +79,10 @@ the predicted step size. Solving Ordinary Differential Equations I Nonstiff Problems [DOI: 10.1007/978-3-540-78862-1](https://doi.org/10.1007/978-3-540-78862-1) """ -struct LegacyIController <: AbstractLegacyController +struct IController <: AbstractLegacyController end -@inline function stepsize_controller!(integrator, controller::LegacyIController, alg) +@inline function stepsize_controller!(integrator, controller::IController, alg) @unpack qmin, qmax, gamma = integrator.opts EEst = DiffEqBase.value(integrator.EEst) @@ -98,7 +98,7 @@ end q end -function step_accept_controller!(integrator, controller::LegacyIController, alg, q) +function step_accept_controller!(integrator, controller::IController, alg, q) @unpack qsteady_min, qsteady_max = integrator.opts if qsteady_min <= q <= qsteady_max @@ -107,12 +107,12 @@ function step_accept_controller!(integrator, controller::LegacyIController, alg, integrator.dt / q # new dt end -function step_reject_controller!(integrator, controller::LegacyIController, alg) +function step_reject_controller!(integrator, controller::IController, alg) @unpack qold = integrator integrator.dt = qold end -struct IController{T} <: AbstractController +struct NewIController{T} <: AbstractController qmin::T qmax::T gamma::T @@ -120,12 +120,12 @@ struct IController{T} <: AbstractController qsteady_max::T end -function IController(alg; kwargs...) - IController(Float64, alg; kwargs...) +function NewIController(alg; kwargs...) + NewIController(Float64, alg; kwargs...) end -function IController(QT, alg; qmin = nothing, qmax = nothing, gamma = nothing, qsteady_min = nothing, qsteady_max = nothing) - IController{QT}( +function NewIController(QT, alg; qmin = nothing, qmax = nothing, gamma = nothing, qsteady_min = nothing, qsteady_max = nothing) + NewIController{QT}( qmin === nothing ? qmin_default(alg) : qmin, qmax === nothing ? qmax_default(alg) : qmax, gamma === nothing ? gamma_default(alg) : gamma, @@ -142,7 +142,7 @@ mutable struct IControllerCache{C, T} <: AbstractControllerCache # EEst::T end -function setup_controller_cache(alg_cache, controller::IController{T}) where T +function setup_controller_cache(alg_cache, controller::NewIController{T}) where T IControllerCache( controller, T(1), @@ -185,7 +185,7 @@ end # PI step size controller """ - LegacyPIController(beta1, beta2) + PIController(beta1, beta2) The proportional-integral (PI) controller is a widespread step size controller with improved stability properties compared to the [`IController`](@ref). @@ -208,7 +208,7 @@ the predicted step size. !!! note The coefficients `beta1, beta2` are not scaled by the order of the method, - in contrast to the [`LegacyPIDController`](@ref). For the `LegacyPIController`, this + in contrast to the [`PIDController`](@ref). For the `PIController`, this scaling by the order must be done when the controller is constructed. ## References @@ -220,12 +220,12 @@ the predicted step size. Solving Ordinary Differential Equations I Nonstiff Problems [DOI: 10.1007/978-3-540-78862-1](https://doi.org/10.1007/978-3-540-78862-1) """ -mutable struct LegacyPIController{QT} <: AbstractLegacyController +mutable struct PIController{QT} <: AbstractLegacyController beta1::QT beta2::QT end -@inline function stepsize_controller!(integrator, controller::LegacyPIController, alg) +@inline function stepsize_controller!(integrator, controller::PIController, alg) @unpack qold = integrator @unpack qmin, qmax, gamma = integrator.opts @unpack beta1, beta2 = controller @@ -242,7 +242,7 @@ end q end -function step_accept_controller!(integrator, controller::LegacyPIController, alg, q) +function step_accept_controller!(integrator, controller::PIController, alg, q) @unpack qsteady_min, qsteady_max, qoldinit = integrator.opts EEst = DiffEqBase.value(integrator.EEst) @@ -253,13 +253,13 @@ function step_accept_controller!(integrator, controller::LegacyPIController, alg return integrator.dt / q # new dt end -function step_reject_controller!(integrator, controller::LegacyPIController, alg) +function step_reject_controller!(integrator, controller::PIController, alg) @unpack q11 = integrator @unpack qmin, gamma = integrator.opts integrator.dt /= min(inv(qmin), q11 / gamma) end -function reset_alg_dependent_opts!(controller::LegacyPIController, alg1, alg2) +function reset_alg_dependent_opts!(controller::PIController, alg1, alg2) if controller.beta2 == beta2_default(alg1) controller.beta2 = beta2_default(alg2) end @@ -268,7 +268,7 @@ function reset_alg_dependent_opts!(controller::LegacyPIController, alg1, alg2) end end -struct PIController{T} <: AbstractController +struct NewPIController{T} <: AbstractController beta1::T beta2::T qmin::T @@ -279,15 +279,15 @@ struct PIController{T} <: AbstractController qoldinit::T end -function PIController(alg; kwargs...) - PIController(Float64, alg; kwargs...) +function NewPIController(alg; kwargs...) + NewPIController(Float64, alg; kwargs...) end -function PIController(QT, alg; beta1 = nothing, beta2 = nothing, qmin = nothing, qmax = nothing, gamma = nothing, qsteady_min = nothing, qsteady_max = nothing, qoldinit = nothing) +function NewPIController(QT, alg; beta1 = nothing, beta2 = nothing, qmin = nothing, qmax = nothing, gamma = nothing, qsteady_min = nothing, qsteady_max = nothing, qoldinit = nothing) beta2 = beta2 === nothing ? beta2_default(alg) : beta2 beta1 = beta1 === nothing ? beta1_default(alg, beta2) : beta1 qoldinit = qoldinit === nothing ? 1 // 10^4 : qoldinit - PIController{QT}( + NewPIController{QT}( beta1, beta2, qmin === nothing ? qmin_default(alg) : qmin, @@ -300,7 +300,7 @@ function PIController(QT, alg; beta1 = nothing, beta2 = nothing, qmin = nothing, end mutable struct PIControllerCache{T} <: AbstractControllerCache - controller::PIController{T} + controller::NewPIController{T} # Propsoed scaling factor for the time step length q::T # Cached εₙ₊₁^β₁ @@ -310,7 +310,7 @@ mutable struct PIControllerCache{T} <: AbstractControllerCache end -function setup_controller_cache(alg_cache, controller::PIController{T}) where T +function setup_controller_cache(alg_cache, controller::NewPIController{T}) where T PIControllerCache( controller, T(1), @@ -371,12 +371,12 @@ end # PID step size controller """ - LegacyPIDController(beta1, beta2, beta3=zero(beta1); + PIDController(beta1, beta2, beta3=zero(beta1); limiter=default_dt_factor_limiter, accept_safety=0.81) The proportional-integral-derivative (PID) controller is a generalization of the -[`LegacyPIController`](@ref) and can have improved stability and efficiency properties. +[`PIController`](@ref) and can have improved stability and efficiency properties. Construct a PID step size controller adapting the time step based on the formula @@ -409,17 +409,17 @@ Some standard controller parameters suggested in the literature are !!! note - In contrast to the [`LegacyPIController`](@ref), the coefficients `beta1, beta2, beta3` + In contrast to the [`PIController`](@ref), the coefficients `beta1, beta2, beta3` are scaled by the order of the method. Thus, standard controllers such as PI42 can use the same coefficients `beta1, beta2, beta3` for different algorithms. !!! note - In contrast to other controllers, the `LegacyPIDController` does not use the keyword + In contrast to other controllers, the `PIDController` does not use the keyword arguments `qmin, qmax` to limit the step size change or the safety factor `gamma`. These common keyword arguments are replaced by the `limiter` and `accept_safety` to guarantee a smooth behavior (Söderlind and Wang, 2006). - Because of this, a `LegacyPIDController` behaves different from a [`LegacyPIController`](@ref), + Because of this, a `PIDController` behaves different from a [`PIController`](@ref), even if `beta1, beta2` are adapted accordingly and `iszero(beta3)`. ## References @@ -435,7 +435,7 @@ Some standard controller parameters suggested in the literature are Compressible Computational Fluid Dynamics # is bigger than this parameter [arXiv:2104.06836](https://arxiv.org/abs/2104.06836) # limiter of the dt factor (before clipping) """ -struct LegacyPIDController{QT, Limiter} <: AbstractLegacyController +struct PIDController{QT, Limiter} <: AbstractLegacyController beta::MVector{3, QT} # controller coefficients err::MVector{3, QT} # history of the error estimates accept_safety::QT # accept a step if the predicted change of the step size @@ -443,17 +443,17 @@ struct LegacyPIDController{QT, Limiter} <: AbstractLegacyController limiter::Limiter # limiter of the dt factor (before clipping) end -function LegacyPIDController(beta1, beta2, beta3 = zero(beta1); +function PIDController(beta1, beta2, beta3 = zero(beta1); limiter = default_dt_factor_limiter, accept_safety = 0.81) beta = MVector(map(float, promote(beta1, beta2, beta3))...) QT = eltype(beta) err = MVector{3, QT}(true, true, true) - return LegacyPIDController(beta, err, convert(QT, accept_safety), limiter) + return PIDController(beta, err, convert(QT, accept_safety), limiter) end -function Base.show(io::IO, controller::LegacyPIDController) - print(io, "LegacyPIDController(beta=", controller.beta, +function Base.show(io::IO, controller::PIDController) + print(io, "PIDController(beta=", controller.beta, ", accept_safety=", controller.accept_safety, ", limiter=", controller.limiter, ")") @@ -461,7 +461,7 @@ end @inline default_dt_factor_limiter(x) = one(x) + atan(x - one(x)) -@inline function stepsize_controller!(integrator, controller::LegacyPIDController, alg) +@inline function stepsize_controller!(integrator, controller::PIDController, alg) @unpack qmax = integrator.opts beta1, beta2, beta3 = controller.beta @@ -505,11 +505,11 @@ end return dt_factor end -@inline function accept_step_controller(integrator, controller::LegacyPIDController) +@inline function accept_step_controller(integrator, controller::PIDController) return integrator.qold >= controller.accept_safety end -function step_accept_controller!(integrator, controller::LegacyPIDController, alg, dt_factor) +function step_accept_controller!(integrator, controller::PIDController, alg, dt_factor) @unpack qsteady_min, qsteady_max = integrator.opts if qsteady_min <= inv(dt_factor) <= qsteady_max @@ -522,13 +522,13 @@ function step_accept_controller!(integrator, controller::LegacyPIDController, al return integrator.dt * dt_factor # new dt end -function step_reject_controller!(integrator, controller::LegacyPIDController, alg) +function step_reject_controller!(integrator, controller::PIDController, alg) integrator.dt *= integrator.qold end -struct PIDController{T, Limiter} <: AbstractController +struct NewPIDController{T, Limiter} <: AbstractController beta::SVector{3, T} # controller coefficients accept_safety::T # accept a step if the predicted change of the step size # is bigger than this parameter @@ -539,11 +539,11 @@ struct PIDController{T, Limiter} <: AbstractController qsteady_max::T end -function PIDController(alg; kwargs...) - PIDController(Float64, alg; kwargs...) +function NewPIDController(alg; kwargs...) + NewPIDController(Float64, alg; kwargs...) end -function PIDController(QT, alg; beta = nothing, accept_safety = 0.81, limiter = default_dt_factor_limiter, qmin = nothing, qmax = nothing, qsteady_min = nothing, qsteady_max = nothing) +function NewPIDController(QT, alg; beta = nothing, accept_safety = 0.81, limiter = default_dt_factor_limiter, qmin = nothing, qmax = nothing, qsteady_min = nothing, qsteady_max = nothing) if beta === nothing beta2 = QT(beta2_default(alg)) beta1 = QT(beta1_default(alg, beta2)) @@ -552,7 +552,7 @@ function PIDController(QT, alg; beta = nothing, accept_safety = 0.81, limiter = beta1, beta2, beta3 = beta end beta = SVector(map(float, promote(beta1, beta2, beta3))...) - PIDController{QT, typeof(limiter)}( + NewPIDController{QT, typeof(limiter)}( beta, QT(accept_safety), limiter, @@ -563,20 +563,20 @@ function PIDController(QT, alg; beta = nothing, accept_safety = 0.81, limiter = ) end -function Base.show(io::IO, controller::PIDController) - print(io, "PIDController(beta=", controller.beta, +function Base.show(io::IO, controller::NewPIDController) + print(io, "NewPIDController(beta=", controller.beta, ", accept_safety=", controller.accept_safety, ", limiter=", controller.limiter, ")") end mutable struct PIDControllerCache{T, Limiter} <: AbstractControllerCache - controller::PIDController{T, Limiter} + controller::NewPIDController{T, Limiter} err::MVector{3, T} # history of the error estimates dt_factor::T end -function setup_controller_cache(alg_cache, controller::PIDController{QT}) where QT +function setup_controller_cache(alg_cache, controller::NewPIDController{QT}) where QT err = MVector{3, QT}(true, true, true) PIDControllerCache( controller, @@ -709,10 +709,10 @@ else end ``` """ -struct LegacyPredictiveController <: AbstractController +struct PredictiveController <: AbstractController end -@inline function stepsize_controller!(integrator, controller::LegacyPredictiveController, alg) +@inline function stepsize_controller!(integrator, controller::PredictiveController, alg) @unpack qmin, qmax, gamma = integrator.opts EEst = DiffEqBase.value(integrator.EEst) if iszero(EEst) @@ -737,7 +737,7 @@ end q end -function step_accept_controller!(integrator, controller::LegacyPredictiveController, alg, q) +function step_accept_controller!(integrator, controller::PredictiveController, alg, q) @unpack qmin, qmax, gamma, qsteady_min, qsteady_max = integrator.opts EEst = DiffEqBase.value(integrator.EEst) @@ -760,13 +760,13 @@ function step_accept_controller!(integrator, controller::LegacyPredictiveControl return integrator.dt / qacc end -function step_reject_controller!(integrator, controller::LegacyPredictiveController, alg) +function step_reject_controller!(integrator, controller::PredictiveController, alg) @unpack dt, success_iter, qold = integrator integrator.dt = success_iter == 0 ? 0.1 * dt : dt / qold end -struct PredictiveController{T} <: AbstractController +struct NewPredictiveController{T} <: AbstractController qmin::T qmax::T gamma::T @@ -774,12 +774,12 @@ struct PredictiveController{T} <: AbstractController qsteady_max::T end -function PredictiveController(alg; kwargs...) +function NewPredictiveController(alg; kwargs...) PIController(Float64, alg; kwargs...) end -function PredictiveController(QT, alg; qmin = nothing, qmax = nothing, gamma = nothing, qsteady_min = nothing, qsteady_max = nothing) - PredictiveController{QT}( +function NewPredictiveController(QT, alg; qmin = nothing, qmax = nothing, gamma = nothing, qsteady_min = nothing, qsteady_max = nothing) + NewPredictiveController{QT}( qmin === nothing ? qmin_default(alg) : qmin, qmax === nothing ? qmax_default(alg) : qmax, gamma === nothing ? gamma_default(alg) : gamma, @@ -789,16 +789,18 @@ function PredictiveController(QT, alg; qmin = nothing, qmax = nothing, gamma = n end mutable struct PredictiveControllerCache{T} <: AbstractControllerCache - controller::PredictiveController{T} + controller::NewPredictiveController{T} dtacc::T erracc::T + qold::T end -function setup_controller_cache(alg_cache, controller::PredictiveController{T}) where T +function setup_controller_cache(alg_cache, controller::NewPredictiveController{T}) where T PredictiveControllerCache{T}( controller, T(1), T(1), + T(1), ) end diff --git a/lib/OrdinaryDiffEqExtrapolation/src/OrdinaryDiffEqExtrapolation.jl b/lib/OrdinaryDiffEqExtrapolation/src/OrdinaryDiffEqExtrapolation.jl index ec663f67b2..1851a2d05f 100644 --- a/lib/OrdinaryDiffEqExtrapolation/src/OrdinaryDiffEqExtrapolation.jl +++ b/lib/OrdinaryDiffEqExtrapolation/src/OrdinaryDiffEqExtrapolation.jl @@ -5,7 +5,7 @@ import OrdinaryDiffEqCore: alg_order, alg_maximum_order, get_current_adaptive_or accept_step_controller, default_controller_v7, legacy_default_controller, beta2_default, beta1_default, gamma_default, initialize!, perform_step!, @unpack, @cache, unwrap_alg, - isthreaded, isadaptive, LegacyPIController, + isthreaded, isadaptive, PIController, step_accept_controller!, calculate_residuals, OrdinaryDiffEqMutableCache, OrdinaryDiffEqConstantCache, reset_alg_dependent_opts!, AbstractController, AbstractLegacyController, diff --git a/lib/OrdinaryDiffEqExtrapolation/src/alg_utils.jl b/lib/OrdinaryDiffEqExtrapolation/src/alg_utils.jl index 7be2b336d7..287b08458f 100644 --- a/lib/OrdinaryDiffEqExtrapolation/src/alg_utils.jl +++ b/lib/OrdinaryDiffEqExtrapolation/src/alg_utils.jl @@ -43,7 +43,7 @@ function legacy_default_controller( qoldinit, _beta1 = nothing, _beta2 = nothing) QT = typeof(qoldinit) beta1, beta2 = _digest_beta1_beta2(alg, cache, Val(QT), _beta1, _beta2) - return LegacyExtrapolationController(beta1) + return ExtrapolationController(beta1) end function default_controller_v7(QT, @@ -55,7 +55,7 @@ function default_controller_v7(QT, ImplicitEulerExtrapolation, ImplicitEulerBarycentricExtrapolation}) beta1 = QT(beta1_default(alg, beta2_default(alg))) - return LegacyExtrapolationController(beta1) + return ExtrapolationController(beta1) end # FIXME AitkenNeville is missing integration with the extrapolation controller and picks up the PI controller instead. @@ -64,12 +64,12 @@ function legacy_default_controller(alg::AitkenNeville, qoldinit, _beta1 = nothing, _beta2 = nothing) QT = typeof(qoldinit) beta1, beta2 = _digest_beta1_beta2(alg, cache, Val(QT), _beta1, _beta2) - return LegacyPIController(beta1, beta2) + return PIController(beta1, beta2) end function default_controller_v7(QT, alg::AitkenNeville) beta2 = QT(beta2_default(alg)) beta1 = QT(beta1_default(alg, beta2)) - return LegacyPIController(beta1, beta2) + return PIController(beta1, beta2) end beta2_default(alg::ExtrapolationMidpointDeuflhard) = 0 // 1 diff --git a/lib/OrdinaryDiffEqExtrapolation/src/controllers.jl b/lib/OrdinaryDiffEqExtrapolation/src/controllers.jl index e47139b474..a62d61cb77 100644 --- a/lib/OrdinaryDiffEqExtrapolation/src/controllers.jl +++ b/lib/OrdinaryDiffEqExtrapolation/src/controllers.jl @@ -1,9 +1,9 @@ # Extrapolation methods -mutable struct LegacyExtrapolationController{QT} <: AbstractLegacyController +mutable struct ExtrapolationController{QT} <: AbstractController beta1::QT end -function reset_alg_dependent_opts!(controller::LegacyExtrapolationController, alg1, alg2) +function reset_alg_dependent_opts!(controller::ExtrapolationController, alg1, alg2) if controller.beta1 == beta1_default(alg1, beta2_default(alg1)) controller.beta1 = beta1_default(alg2, beta2_default(alg2)) end diff --git a/lib/OrdinaryDiffEqFIRK/src/OrdinaryDiffEqFIRK.jl b/lib/OrdinaryDiffEqFIRK/src/OrdinaryDiffEqFIRK.jl index 5c7860651e..e4693dc3b8 100644 --- a/lib/OrdinaryDiffEqFIRK/src/OrdinaryDiffEqFIRK.jl +++ b/lib/OrdinaryDiffEqFIRK/src/OrdinaryDiffEqFIRK.jl @@ -14,8 +14,8 @@ import OrdinaryDiffEqCore: alg_order, calculate_residuals!, qmax_default, alg_adaptive_order, DEFAULT_PRECS, stepsize_controller!, step_accept_controller!, step_reject_controller!, - LegacyPredictiveController, PredictiveControllerCache, - PredictiveController, alg_can_repeat_jac, NewtonAlgorithm, + PredictiveController, PredictiveControllerCache, + NewPredictiveController, alg_can_repeat_jac, NewtonAlgorithm, fac_default_gamma, get_current_adaptive_order, get_fsalfirstlast, isfirk, generic_solver_docstring, _bool_to_ADType, diff --git a/lib/OrdinaryDiffEqFIRK/src/controllers.jl b/lib/OrdinaryDiffEqFIRK/src/controllers.jl index 480503010f..87fad2e5d6 100644 --- a/lib/OrdinaryDiffEqFIRK/src/controllers.jl +++ b/lib/OrdinaryDiffEqFIRK/src/controllers.jl @@ -1,5 +1,5 @@ function step_accept_controller!( - integrator, controller::LegacyPredictiveController, alg::AdaptiveRadau, q) + integrator, controller::PredictiveController, alg::AdaptiveRadau, q) @unpack qmin, qmax, gamma, qsteady_min, qsteady_max = integrator.opts @unpack cache = integrator @unpack num_stages, step, iter, hist_iter, index = cache @@ -43,7 +43,7 @@ function step_accept_controller!( end function step_reject_controller!( - integrator, controller::LegacyPredictiveController, alg::AdaptiveRadau) + integrator, controller::PredictiveController, alg::AdaptiveRadau) @unpack dt, success_iter, qold = integrator @unpack cache = integrator @unpack num_stages, step, iter, hist_iter = cache diff --git a/lib/OrdinaryDiffEqLowStorageRK/src/OrdinaryDiffEqLowStorageRK.jl b/lib/OrdinaryDiffEqLowStorageRK/src/OrdinaryDiffEqLowStorageRK.jl index 416f52efd3..0028ec16be 100644 --- a/lib/OrdinaryDiffEqLowStorageRK/src/OrdinaryDiffEqLowStorageRK.jl +++ b/lib/OrdinaryDiffEqLowStorageRK/src/OrdinaryDiffEqLowStorageRK.jl @@ -8,7 +8,7 @@ import OrdinaryDiffEqCore: alg_order, alg_adaptive_order, calculate_residuals!, OrdinaryDiffEqMutableCache, OrdinaryDiffEqConstantCache, OrdinaryDiffEqAdaptiveAlgorithm, uses_uprev, default_controller_v7, - legacy_default_controller, PIDController, LegacyPIDController, + legacy_default_controller, NewPIDController, PIDController, alg_cache, _vec, _reshape, @cache, isfsal, full_cache, constvalue, _unwrap_val, trivial_limiter!, perform_step!, initialize!, diff --git a/lib/OrdinaryDiffEqLowStorageRK/src/alg_utils.jl b/lib/OrdinaryDiffEqLowStorageRK/src/alg_utils.jl index e849eb6469..00d39c86be 100644 --- a/lib/OrdinaryDiffEqLowStorageRK/src/alg_utils.jl +++ b/lib/OrdinaryDiffEqLowStorageRK/src/alg_utils.jl @@ -86,54 +86,54 @@ uses_uprev(alg::CKLLSRK75_4M_5R, adaptive::Bool) = adaptive function legacy_default_controller(alg::RDPK3Sp35, cache, qoldinit, args...) QT = typeof(qoldinit) - return LegacyPIDController(map(Base.Fix1(convert, QT), (0.64, -0.31, 0.04))...) + return PIDController(map(Base.Fix1(convert, QT), (0.64, -0.31, 0.04))...) end function legacy_default_controller(alg::RDPK3SpFSAL35, cache, qoldinit, args...) QT = typeof(qoldinit) - return LegacyPIDController(map(Base.Fix1(convert, QT), (0.70, -0.23, 0.00))...) + return PIDController(map(Base.Fix1(convert, QT), (0.70, -0.23, 0.00))...) end function legacy_default_controller(alg::RDPK3Sp49, cache, qoldinit, args...) QT = typeof(qoldinit) - return LegacyPIDController(map(Base.Fix1(convert, QT), (0.25, -0.12, 0.00))...) + return PIDController(map(Base.Fix1(convert, QT), (0.25, -0.12, 0.00))...) end function legacy_default_controller(alg::RDPK3SpFSAL49, cache, qoldinit, args...) QT = typeof(qoldinit) - return LegacyPIDController(map(Base.Fix1(convert, QT), (0.38, -0.18, 0.01))...) + return PIDController(map(Base.Fix1(convert, QT), (0.38, -0.18, 0.01))...) end function legacy_default_controller(alg::RDPK3Sp510, cache, qoldinit, args...) QT = typeof(qoldinit) - return LegacyPIDController(map(Base.Fix1(convert, QT), (0.47, -0.20, 0.06))...) + return PIDController(map(Base.Fix1(convert, QT), (0.47, -0.20, 0.06))...) end function legacy_default_controller(alg::RDPK3SpFSAL510, cache, qoldinit, args...) QT = typeof(qoldinit) - return LegacyPIDController(map(Base.Fix1(convert, QT), (0.45, -0.13, 0.00))...) + return PIDController(map(Base.Fix1(convert, QT), (0.45, -0.13, 0.00))...) end function default_controller_v7(QT, alg::RDPK3Sp35) - return PIDController(QT, alg; beta=map(Base.Fix1(convert, QT), (0.64, -0.31, 0.04))) + return NewPIDController(QT, alg; beta=map(Base.Fix1(convert, QT), (0.64, -0.31, 0.04))) end function default_controller_v7(QT, alg::RDPK3SpFSAL35) - return PIDController(QT, alg; beta=map(Base.Fix1(convert, QT), (0.70, -0.23, 0.00))) + return NewPIDController(QT, alg; beta=map(Base.Fix1(convert, QT), (0.70, -0.23, 0.00))) end function default_controller_v7(QT, alg::RDPK3Sp49) - return PIDController(QT, alg; beta=map(Base.Fix1(convert, QT), (0.25, -0.12, 0.00))) + return NewPIDController(QT, alg; beta=map(Base.Fix1(convert, QT), (0.25, -0.12, 0.00))) end function default_controller_v7(QT, alg::RDPK3SpFSAL49) - return PIDController(QT, alg; beta=map(Base.Fix1(convert, QT), (0.38, -0.18, 0.01))) + return NewPIDController(QT, alg; beta=map(Base.Fix1(convert, QT), (0.38, -0.18, 0.01))) end function default_controller_v7(QT, alg::RDPK3Sp510) - return PIDController(QT, alg; beta=map(Base.Fix1(convert, QT), (0.47, -0.20, 0.06))) + return NewPIDController(QT, alg; beta=map(Base.Fix1(convert, QT), (0.47, -0.20, 0.06))) end function default_controller_v7(QT, alg::RDPK3SpFSAL510) - return PIDController(QT, alg; beta=map(Base.Fix1(convert, QT), (0.45, -0.13, 0.00))) + return NewPIDController(QT, alg; beta=map(Base.Fix1(convert, QT), (0.45, -0.13, 0.00))) end From 49a021592cb607b3228b172b5834b5bebc1c0ce0 Mon Sep 17 00:00:00 2001 From: termi-official Date: Mon, 27 Oct 2025 17:32:54 +0100 Subject: [PATCH 14/21] Update Extrapolation alg controller --- .../src/OrdinaryDiffEqExtrapolation.jl | 6 +- .../src/alg_utils.jl | 3 +- .../src/controllers.jl | 156 +++++++++++++++++- 3 files changed, 159 insertions(+), 6 deletions(-) diff --git a/lib/OrdinaryDiffEqExtrapolation/src/OrdinaryDiffEqExtrapolation.jl b/lib/OrdinaryDiffEqExtrapolation/src/OrdinaryDiffEqExtrapolation.jl index 1851a2d05f..ba684b58a1 100644 --- a/lib/OrdinaryDiffEqExtrapolation/src/OrdinaryDiffEqExtrapolation.jl +++ b/lib/OrdinaryDiffEqExtrapolation/src/OrdinaryDiffEqExtrapolation.jl @@ -5,7 +5,7 @@ import OrdinaryDiffEqCore: alg_order, alg_maximum_order, get_current_adaptive_or accept_step_controller, default_controller_v7, legacy_default_controller, beta2_default, beta1_default, gamma_default, initialize!, perform_step!, @unpack, @cache, unwrap_alg, - isthreaded, isadaptive, PIController, + isthreaded, isadaptive, PIController, AbstractControllerCache, step_accept_controller!, calculate_residuals, OrdinaryDiffEqMutableCache, OrdinaryDiffEqConstantCache, reset_alg_dependent_opts!, AbstractController, AbstractLegacyController, @@ -13,12 +13,12 @@ import OrdinaryDiffEqCore: alg_order, alg_maximum_order, get_current_adaptive_or OrdinaryDiffEqAdaptiveAlgorithm, OrdinaryDiffEqAdaptiveImplicitAlgorithm, alg_cache, CompiledFloats, @threaded, stepsize_controller!, - DEFAULT_PRECS, full_cache, qmin_default, + DEFAULT_PRECS, full_cache, qmin_default, qmax_default, constvalue, PolyesterThreads, Sequential, BaseThreads, _digest_beta1_beta2, timedepentdtmin, _unwrap_val, _reshape, _vec, get_fsalfirstlast, generic_solver_docstring, differentiation_rk_docstring, _bool_to_ADType, - _process_AD_choice, LinearAliasSpecifier + _process_AD_choice, LinearAliasSpecifier, setup_controller_cache using FastBroadcast, Polyester, MuladdMacro, RecursiveArrayTools, LinearSolve import OrdinaryDiffEqCore import FastPower diff --git a/lib/OrdinaryDiffEqExtrapolation/src/alg_utils.jl b/lib/OrdinaryDiffEqExtrapolation/src/alg_utils.jl index 287b08458f..425063990c 100644 --- a/lib/OrdinaryDiffEqExtrapolation/src/alg_utils.jl +++ b/lib/OrdinaryDiffEqExtrapolation/src/alg_utils.jl @@ -54,8 +54,7 @@ function default_controller_v7(QT, ImplicitHairerWannerExtrapolation, ImplicitEulerExtrapolation, ImplicitEulerBarycentricExtrapolation}) - beta1 = QT(beta1_default(alg, beta2_default(alg))) - return ExtrapolationController(beta1) + return NewExtrapolationController(QT, alg) end # FIXME AitkenNeville is missing integration with the extrapolation controller and picks up the PI controller instead. diff --git a/lib/OrdinaryDiffEqExtrapolation/src/controllers.jl b/lib/OrdinaryDiffEqExtrapolation/src/controllers.jl index a62d61cb77..a2d292c7bd 100644 --- a/lib/OrdinaryDiffEqExtrapolation/src/controllers.jl +++ b/lib/OrdinaryDiffEqExtrapolation/src/controllers.jl @@ -1,5 +1,5 @@ # Extrapolation methods -mutable struct ExtrapolationController{QT} <: AbstractController +mutable struct ExtrapolationController{QT} <: AbstractLegacyController beta1::QT end @@ -9,7 +9,37 @@ function reset_alg_dependent_opts!(controller::ExtrapolationController, alg1, al end end +struct NewExtrapolationController{QT} <: AbstractController + qmin::QT + qmax::QT +end + +function NewExtrapolationController(QT, alg; qmin = nothing, qmax = nothing, gamma = nothing) + NewExtrapolationController( + QT(qmin === nothing ? qmin_default(alg) : qmin), + QT(qmax === nothing ? qmax_default(alg) : qmax), + ) +end + +mutable struct ExtrapolationControllerCache{QT} <: AbstractControllerCache + controller::NewExtrapolationController{QT} + beta1::QT + gamma::QT +end + +function setup_controller_cache(alg_cache, controller::NewExtrapolationController{T}) where T + ExtrapolationControllerCache{T}( + controller, + T(1), + T(1) + ) +end + +stepsize_controller_internal!(integrator, alg) = stepsize_controller_internal!(integrator, integrator.controller_cache, alg) +stepsize_predictor!(integrator, alg, n_new) = stepsize_predictor!(integrator, integrator.controller_cache, alg, n_new) + @inline function stepsize_controller!(integrator, + controller::ExtrapolationController, alg::Union{ExtrapolationMidpointDeuflhard, ImplicitDeuflhardExtrapolation}) # Dummy function @@ -19,7 +49,19 @@ end zero(typeof(integrator.opts.qmax)) end +@inline function stepsize_controller!(integrator, + cache::ExtrapolationControllerCache, + alg::Union{ExtrapolationMidpointDeuflhard, + ImplicitDeuflhardExtrapolation}) + # Dummy function + # ExtrapolationMidpointDeuflhard's stepsize scaling is stored in the cache; + # it is computed by stepsize_controller_internal! (in perform_step!) resp. stepsize_predictor! + # (in step_accept_controller! and step_reject_controller!) + zero(typeof(cache.controller.qmax)) +end + function stepsize_controller_internal!(integrator, + controller::ExtrapolationController, alg::Union{ExtrapolationMidpointDeuflhard, ImplicitDeuflhardExtrapolation}) # Standard step size controller @@ -41,7 +83,32 @@ function stepsize_controller_internal!(integrator, integrator.cache.Q[integrator.cache.n_curr - alg.min_order + 1] = q end +function stepsize_controller_internal!(integrator, + cache::ExtrapolationControllerCache, + alg::Union{ExtrapolationMidpointDeuflhard, + ImplicitDeuflhardExtrapolation}) + # Standard step size controller + # Compute and save the stepsize scaling based on the latest error estimate of the current order + @unpack controller = cache + + if iszero(integrator.EEst) + q = inv(controller.qmax) + else + # Update gamma and beta1 + cache.beta1 = typeof(cache.beta1)(1 // (2integrator.cache.n_curr + 1)) + cache.gamma = FastPower.fastpower(typeof(cache.gamma)(1 // 4), + cache.beta1) + # Compute new stepsize scaling + qtmp = FastPower.fastpower(integrator.EEst, cache.beta1) / + cache.gamma + @fastmath q = max(inv(controller.qmax), min(inv(controller.qmin), qtmp)) + end + integrator.cache.Q[integrator.cache.n_curr - alg.min_order + 1] = q +end + + function stepsize_predictor!(integrator, + controller::ExtrapolationController, alg::Union{ExtrapolationMidpointDeuflhard, ImplicitDeuflhardExtrapolation}, n_new::Int) # Compute and save the stepsize scaling for order n_new based on the latest error estimate of the current order. @@ -69,6 +136,35 @@ function stepsize_predictor!(integrator, integrator.cache.Q[n_new - alg.min_order + 1] = q end +function stepsize_predictor!(integrator, + cache::ExtrapolationControllerCache, + alg::Union{ExtrapolationMidpointDeuflhard, + ImplicitDeuflhardExtrapolation}, n_new::Int) + # Compute and save the stepsize scaling for order n_new based on the latest error estimate of the current order. + @unpack controller = cache + + if iszero(integrator.EEst) + q = inv(controller.qmax) + else + # Initialize + @unpack t, EEst = integrator + @unpack stage_number = integrator.cache + tol = integrator.opts.internalnorm(integrator.opts.reltol, t) # Deuflhard's approach relies on EEstD ≈ ||relTol|| + s_curr = stage_number[integrator.cache.n_curr - alg.min_order + 1] + s_new = stage_number[n_new - alg.min_order + 1] + # Update gamma and beta1 + cache.beta1 = typeof(cache.beta1)(1 // (2integrator.cache.n_curr + 1)) + cache.gamma = FastPower.fastpower(typeof(cache.gamma)(1 // 4), + cache.beta1) + # Compute new stepsize scaling + qtmp = EEst * + FastPower.fastpower(FastPower.fastpower(tol, (1.0 - s_curr / s_new)), + cache.beta1) / cache.gamma + @fastmath q = max(inv(controller.qmax), min(inv(controller.qmin), qtmp)) + end + integrator.cache.Q[n_new - alg.min_order + 1] = q +end + function step_accept_controller!(integrator, alg::Union{ExtrapolationMidpointDeuflhard, ImplicitDeuflhardExtrapolation}, q) @@ -136,6 +232,7 @@ end end function stepsize_controller_internal!(integrator, + controller::ExtrapolationController, alg::Union{ExtrapolationMidpointHairerWanner, ImplicitHairerWannerExtrapolation, ImplicitEulerExtrapolation, @@ -191,6 +288,63 @@ function stepsize_controller_internal!(integrator, end end +function stepsize_controller_internal!(integrator, + cache::ExtrapolationControllerCache, + alg::Union{ExtrapolationMidpointHairerWanner, + ImplicitHairerWannerExtrapolation, + ImplicitEulerExtrapolation, + ImplicitEulerBarycentricExtrapolation}) + # Standard step size controller + # Compute and save the stepsize scaling based on the latest error estimate of the current order + @unpack controller = cache + + if alg isa + Union{ImplicitEulerExtrapolation, ImplicitEulerBarycentricExtrapolation, + ImplicitHairerWannerExtrapolation} + if iszero(integrator.EEst) + q = inv(controller.qmax) + else + # Update gamma and beta1 + if alg isa ImplicitHairerWannerExtrapolation + cache.beta1 = typeof(cache.beta1)(1 // + (2integrator.cache.n_curr + 1)) + elseif alg isa ImplicitEulerExtrapolation + cache.beta1 = typeof(cache.beta1)(1 // (integrator.cache.n_curr)) + else + cache.beta1 = typeof(cache.beta1)(1 // + (integrator.cache.n_curr - 1)) + end + cache.gamma = FastPower.fastpower( + typeof(cache.gamma)(65 // + 100), + cache.beta1) + # Compute new stepsize scaling + qtmp = FastPower.fastpower(integrator.EEst, cache.beta1) / + (cache.gamma) + @fastmath q = max(inv(controller.qmax), + min(inv(controller.qmin), qtmp)) + end + integrator.cache.Q[integrator.cache.n_curr + 1] = q + else + if iszero(integrator.EEst) + q = inv(controller.qmax) + else + # Update gamma and beta1 + cache.beta1 = typeof(cache.beta1)(1 // (2integrator.cache.n_curr + 1)) + cache.gamma = FastPower.fastpower( + typeof(cache.gamma)(65 // + 100), + cache.beta1) + # Compute new stepsize scaling + qtmp = FastPower.fastpower(integrator.EEst, cache.beta1) / + cache.gamma + @fastmath q = max(inv(controller.qmax), + min(inv(controller.qmin), qtmp)) + end + integrator.cache.Q[integrator.cache.n_curr + 1] = q + end +end + function step_accept_controller!(integrator, alg::Union{ExtrapolationMidpointHairerWanner, ImplicitHairerWannerExtrapolation, From 0698ed2f2158372ceae08632707aaad7bec210da Mon Sep 17 00:00:00 2001 From: termi-official Date: Mon, 27 Oct 2025 17:56:19 +0100 Subject: [PATCH 15/21] Add composite algorithm controller --- lib/OrdinaryDiffEqCore/src/alg_utils.jl | 6 ++++ .../src/integrators/controllers.jl | 35 +++++++++++++++++++ 2 files changed, 41 insertions(+) diff --git a/lib/OrdinaryDiffEqCore/src/alg_utils.jl b/lib/OrdinaryDiffEqCore/src/alg_utils.jl index 78a488aa1b..32f221642a 100644 --- a/lib/OrdinaryDiffEqCore/src/alg_utils.jl +++ b/lib/OrdinaryDiffEqCore/src/alg_utils.jl @@ -319,6 +319,12 @@ function default_controller_v7(QT, alg) end end +function default_controller_v7(QT, alg::CompositeAlgorithm) + return CompositeController( + map(alg->default_controller_v7(QT, alg), alg.algs) + ) +end + function legacy_default_controller(alg, cache, qoldinit, _beta1 = nothing, _beta2 = nothing) if ispredictive(alg) return PredictiveController() diff --git a/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl b/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl index f72e68fb1b..9e2a114696 100644 --- a/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl +++ b/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl @@ -860,3 +860,38 @@ function step_reject_controller!(integrator, cache::PredictiveControllerCache, a @unpack qold = cache integrator.dt = success_iter == 0 ? 0.1 * dt : dt / qold end + +# For a composite algorithm the default strategy is to switch forth and back between the controllers of the individual algorithms +struct CompositeController{T} <: AbstractController + controllers::T +end + +struct CompositeControllerCache{T} <: AbstractControllerCache + caches::T +end + +function setup_controller_cache(alg_cache::CompositeCache, cc::CompositeController) + CompositeControllerCache( + map((cache,controller)->default_controller_v7(cache, controller), alg_cache.caches, cc.controllers) + ) +end + +@inline function accept_step_controller(integrator, cache::CompositeControllerCache) + current_idx = integrator.cache.current + accept_step_controller(integrator, cache.caches[current_idx]) +end + +@inline function stepsize_controller!(integrator, cache::CompositeControllerCache) + current_idx = integrator.cache.current + stepsize_controller!(integrator, cache.caches[current_idx]) +end + +@inline function step_accept_controller!(integrator, cache::CompositeControllerCache, q) + current_idx = integrator.cache.current + step_accept_controller!(integrator, cache.caches[current_idx], q) +end + +@inline function step_reject_controller!(integrator, cache::CompositeControllerCache) + current_idx = integrator.cache.current + step_reject_controller!(integrator, cache.caches[current_idx]) +end From 08d4325f6cb0f3e58d3dc874b002a34c8d75b5eb Mon Sep 17 00:00:00 2001 From: termi-official Date: Mon, 27 Oct 2025 18:09:23 +0100 Subject: [PATCH 16/21] Docstrings --- .../src/integrators/controllers.jl | 47 +++++++++++++++++++ 1 file changed, 47 insertions(+) diff --git a/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl b/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl index 9e2a114696..b7033e469d 100644 --- a/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl +++ b/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl @@ -3,6 +3,53 @@ abstract type AbstractLegacyController <: AbstractController end abstract type AbstractControllerCache end +""" + setup_controller_cache(alg_cache, controller::AbstractController)::AbstractControllerCache + +This function takes a controller together with the cache of the time stepping +algorithm to construct and initialize the respective cache for the controller. +""" +setup_controller_cache + +""" + accept_step_controller(integrator, controller_cache::AbstractControllerCache)::Bool + +This function decides whether the current time step should be accepted or rejected. +A return value of `false` corresponds to a rejection. +""" +accept_step_controller + +""" + stepsize_controller!(integrator, controller_cache::AbstractControllerCache) + +Update the cache to compute the new order of the time marching algorithm and prepare +for an update of the time step. The update can be either due to a rejection or acceptance +of the current time step. +""" +stepsize_controller! + +""" + step_accept_controller!(integrator, controller_cache::AbstractControllerCache, q) + +This function gets called in case of an accepted time step right after [`stepsize_controller!`](@ref). +It returns the proposed new time step length. Please note that the time step length might not be +applied as is and subject to further modification to e.g. match the next time stop. + +!!! warning + The parameter `q` will be removed the next release and will be passed through the + controller cache if needed. +""" +step_accept_controller! + +""" + step_reject_controller!(integrator, controller_cache::AbstractControllerCache) + +This function gets called in case of a rejected time step right after [`stepsize_controller!`](@ref). +It directly sets the time step length (i.e. `integrator.dt`). +""" +step_reject_controller! + + # The legacy controllers do not have this concept. setup_controller_cache(alg_cache, controller::AbstractLegacyController) = controller From aab24b9e64c21634451db19fffbfd8271e255c8d Mon Sep 17 00:00:00 2001 From: termi-official Date: Mon, 27 Oct 2025 19:12:03 +0100 Subject: [PATCH 17/21] Fix default algorithm --- .../src/integrators/controllers.jl | 62 +++++++++---------- lib/OrdinaryDiffEqCore/src/solve.jl | 2 +- .../src/controllers.jl | 2 +- 3 files changed, 31 insertions(+), 35 deletions(-) diff --git a/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl b/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl index b7033e469d..8622af3d68 100644 --- a/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl +++ b/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl @@ -4,15 +4,16 @@ abstract type AbstractLegacyController <: AbstractController end abstract type AbstractControllerCache end """ - setup_controller_cache(alg_cache, controller::AbstractController)::AbstractControllerCache + setup_controller_cache(alg, controller::AbstractController)::AbstractControllerCache -This function takes a controller together with the cache of the time stepping -algorithm to construct and initialize the respective cache for the controller. +This function takes a controller together with the time stepping algorithm to +construct and initialize the respective cache for the controller. """ setup_controller_cache """ - accept_step_controller(integrator, controller_cache::AbstractControllerCache)::Bool + accept_step_controller(integrator, alg)::Bool + accept_step_controller(integrator, controller_cache::AbstractControllerCache, alg)::Bool This function decides whether the current time step should be accepted or rejected. A return value of `false` corresponds to a rejection. @@ -20,7 +21,8 @@ A return value of `false` corresponds to a rejection. accept_step_controller """ - stepsize_controller!(integrator, controller_cache::AbstractControllerCache) + stepsize_controller!(integrator, alg) + stepsize_controller!(integrator, controller_cache::AbstractControllerCache, alg) Update the cache to compute the new order of the time marching algorithm and prepare for an update of the time step. The update can be either due to a rejection or acceptance @@ -29,7 +31,8 @@ of the current time step. stepsize_controller! """ - step_accept_controller!(integrator, controller_cache::AbstractControllerCache, q) + step_accept_controller!(integrator, alg, q) + step_accept_controller!(integrator, controller_cache::AbstractControllerCache, alg, q) This function gets called in case of an accepted time step right after [`stepsize_controller!`](@ref). It returns the proposed new time step length. Please note that the time step length might not be @@ -42,7 +45,8 @@ applied as is and subject to further modification to e.g. match the next time st step_accept_controller! """ - step_reject_controller!(integrator, controller_cache::AbstractControllerCache) + step_reject_controller!(integrator, alg) + step_reject_controller!(integrator, controller_cache::AbstractControllerCache, alg) This function gets called in case of a rejected time step right after [`stepsize_controller!`](@ref). It directly sets the time step length (i.e. `integrator.dt`). @@ -51,7 +55,7 @@ step_reject_controller! # The legacy controllers do not have this concept. -setup_controller_cache(alg_cache, controller::AbstractLegacyController) = controller +setup_controller_cache(alg, controller::AbstractLegacyController) = controller # checks whether the controller should accept a step based on the error estimate @inline function accept_step_controller(integrator, controller_or_cache::Union{<:AbstractLegacyController, <:AbstractControllerCache}) @@ -89,7 +93,7 @@ end struct DummyController <: AbstractController end -setup_controller_cache(alg_cache, controller::DummyController) = controller +setup_controller_cache(alg, controller::DummyController) = controller @inline function accept_step_controller(integrator, controller::DummyController) return integrator.EEst <= 1 @@ -189,7 +193,7 @@ mutable struct IControllerCache{C, T} <: AbstractControllerCache # EEst::T end -function setup_controller_cache(alg_cache, controller::NewIController{T}) where T +function setup_controller_cache(alg, controller::NewIController{T}) where T IControllerCache( controller, T(1), @@ -357,7 +361,7 @@ mutable struct PIControllerCache{T} <: AbstractControllerCache end -function setup_controller_cache(alg_cache, controller::NewPIController{T}) where T +function setup_controller_cache(alg, controller::NewPIController{T}) where T PIControllerCache( controller, T(1), @@ -406,16 +410,6 @@ function step_reject_controller!(integrator, cache::PIControllerCache, alg) integrator.dt /= min(inv(qmin), q11 / gamma) end -# FIXME multi-controller? -# function reset_alg_dependent_opts!(controller::PIControllerCache, alg1, alg2) -# if controller.beta2 == beta2_default(alg1) -# controller.beta2 = beta2_default(alg2) -# end -# if controller.beta1 == beta1_default(alg1, controller.beta2) -# controller.beta1 = beta1_default(alg2, controller.beta2) -# end -# end - # PID step size controller """ PIDController(beta1, beta2, beta3=zero(beta1); @@ -623,7 +617,7 @@ mutable struct PIDControllerCache{T, Limiter} <: AbstractControllerCache dt_factor::T end -function setup_controller_cache(alg_cache, controller::NewPIDController{QT}) where QT +function setup_controller_cache(alg, controller::NewPIDController{QT}) where QT err = MVector{3, QT}(true, true, true) PIDControllerCache( controller, @@ -840,14 +834,16 @@ mutable struct PredictiveControllerCache{T} <: AbstractControllerCache dtacc::T erracc::T qold::T + q::T end -function setup_controller_cache(alg_cache, controller::NewPredictiveController{T}) where T +function setup_controller_cache(alg, controller::NewPredictiveController{T}) where T PredictiveControllerCache{T}( controller, T(1), T(1), T(1), + T(1) ) end @@ -917,28 +913,28 @@ struct CompositeControllerCache{T} <: AbstractControllerCache caches::T end -function setup_controller_cache(alg_cache::CompositeCache, cc::CompositeController) +function setup_controller_cache(alg::CompositeAlgorithm, cc::CompositeController) CompositeControllerCache( - map((cache,controller)->default_controller_v7(cache, controller), alg_cache.caches, cc.controllers) + map((alg, controller)->setup_controller_cache(alg, controller), alg.algs, cc.controllers) ) end -@inline function accept_step_controller(integrator, cache::CompositeControllerCache) +@inline function accept_step_controller(integrator, cache::CompositeControllerCache, alg::CompositeAlgorithm) current_idx = integrator.cache.current - accept_step_controller(integrator, cache.caches[current_idx]) + accept_step_controller(integrator, cache.caches[current_idx], alg.algs[current_idx]) end -@inline function stepsize_controller!(integrator, cache::CompositeControllerCache) +@inline function stepsize_controller!(integrator, cache::CompositeControllerCache, alg::CompositeAlgorithm) current_idx = integrator.cache.current - stepsize_controller!(integrator, cache.caches[current_idx]) + stepsize_controller!(integrator, cache.caches[current_idx], alg.algs[current_idx]) end -@inline function step_accept_controller!(integrator, cache::CompositeControllerCache, q) +@inline function step_accept_controller!(integrator, cache::CompositeControllerCache, alg::CompositeAlgorithm, q) current_idx = integrator.cache.current - step_accept_controller!(integrator, cache.caches[current_idx], q) + step_accept_controller!(integrator, cache.caches[current_idx], alg.algs[current_idx], q) end -@inline function step_reject_controller!(integrator, cache::CompositeControllerCache) +@inline function step_reject_controller!(integrator, cache::CompositeControllerCache, alg::CompositeAlgorithm) current_idx = integrator.cache.current - step_reject_controller!(integrator, cache.caches[current_idx]) + step_reject_controller!(integrator, cache.caches[current_idx], alg.algs[current_idx]) end diff --git a/lib/OrdinaryDiffEqCore/src/solve.jl b/lib/OrdinaryDiffEqCore/src/solve.jl index dffb8fddd0..57956b82f8 100644 --- a/lib/OrdinaryDiffEqCore/src/solve.jl +++ b/lib/OrdinaryDiffEqCore/src/solve.jl @@ -460,7 +460,7 @@ function SciMLBase.__init( qoldinit = hasfield(typeof(controller), :qoldinit) ? controller.qoldinit : (anyadaptive(alg) ? 1 // 10^4 : 0) end - controller_cache = setup_controller_cache(cache, controller) + controller_cache = setup_controller_cache(_alg, controller) save_end_user = save_end save_end = save_end === nothing ? diff --git a/lib/OrdinaryDiffEqExtrapolation/src/controllers.jl b/lib/OrdinaryDiffEqExtrapolation/src/controllers.jl index a2d292c7bd..b1fdd0ce95 100644 --- a/lib/OrdinaryDiffEqExtrapolation/src/controllers.jl +++ b/lib/OrdinaryDiffEqExtrapolation/src/controllers.jl @@ -27,7 +27,7 @@ mutable struct ExtrapolationControllerCache{QT} <: AbstractControllerCache gamma::QT end -function setup_controller_cache(alg_cache, controller::NewExtrapolationController{T}) where T +function setup_controller_cache(alg, controller::NewExtrapolationController{T}) where T ExtrapolationControllerCache{T}( controller, T(1), From eb1e04bee7f65dc5a53b101d7a024fe32676e20a Mon Sep 17 00:00:00 2001 From: termi-official Date: Mon, 27 Oct 2025 19:33:54 +0100 Subject: [PATCH 18/21] Use existing variables to trick downstream packages --- lib/OrdinaryDiffEqCore/src/alg_utils.jl | 2 +- .../src/integrators/controllers.jl | 12 +++++++++--- lib/OrdinaryDiffEqExtrapolation/src/controllers.jl | 7 +++++-- 3 files changed, 15 insertions(+), 6 deletions(-) diff --git a/lib/OrdinaryDiffEqCore/src/alg_utils.jl b/lib/OrdinaryDiffEqCore/src/alg_utils.jl index 32f221642a..6c2282f609 100644 --- a/lib/OrdinaryDiffEqCore/src/alg_utils.jl +++ b/lib/OrdinaryDiffEqCore/src/alg_utils.jl @@ -338,7 +338,7 @@ function legacy_default_controller(alg, cache, qoldinit, _beta1 = nothing, _beta end # TODO remove this when done -default_controller = legacy_default_controller +default_controller(args...) = legacy_default_controller(args...) function _digest_beta1_beta2(alg, cache, ::Val{QT}, _beta1, _beta2) where {QT} if alg isa OrdinaryDiffEqCompositeAlgorithm diff --git a/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl b/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl index 8622af3d68..983b1942dd 100644 --- a/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl +++ b/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl @@ -63,15 +63,21 @@ setup_controller_cache(alg, controller::AbstractLegacyController) = controller end @inline function stepsize_controller!(integrator, alg) - stepsize_controller!(integrator, integrator.controller_cache, alg) + # TODO replace this when done - right now this holds the controller cache! + stepsize_controller!(integrator, integrator.opts.controller, alg) + # stepsize_controller!(integrator, integrator.controller_cache, alg) end @inline function step_accept_controller!(integrator, alg, q) - step_accept_controller!(integrator, integrator.controller_cache, alg, q) + # TODO replace this when done - right now this holds the controller cache! + step_accept_controller!(integrator, integrator.opts.controller, alg, q) + # step_accept_controller!(integrator, integrator.controller_cache, alg, q) end @inline function step_reject_controller!(integrator, alg) - step_reject_controller!(integrator, integrator.controller_cache, alg) + # TODO replace this when done - right now this holds the controller cache! + step_reject_controller!(integrator, integrator.opts.controller, alg) + # step_reject_controller!(integrator, integrator.controller_cache, alg) cache = integrator.cache if hasfield(typeof(cache), :nlsolve) nlsolve = cache.nlsolve diff --git a/lib/OrdinaryDiffEqExtrapolation/src/controllers.jl b/lib/OrdinaryDiffEqExtrapolation/src/controllers.jl index b1fdd0ce95..2f2b25032b 100644 --- a/lib/OrdinaryDiffEqExtrapolation/src/controllers.jl +++ b/lib/OrdinaryDiffEqExtrapolation/src/controllers.jl @@ -35,8 +35,11 @@ function setup_controller_cache(alg, controller::NewExtrapolationController{T}) ) end -stepsize_controller_internal!(integrator, alg) = stepsize_controller_internal!(integrator, integrator.controller_cache, alg) -stepsize_predictor!(integrator, alg, n_new) = stepsize_predictor!(integrator, integrator.controller_cache, alg, n_new) +# TODO replace this when done - right now this holds the controller cache! +stepsize_controller_internal!(integrator, alg) = stepsize_controller_internal!(integrator, integrator.opts.controller, alg) +stepsize_predictor!(integrator, alg, n_new) = stepsize_predictor!(integrator, integrator.opts.controller, alg, n_new) +# stepsize_controller_internal!(integrator, alg) = stepsize_controller_internal!(integrator, integrator.controller_cache, alg) +# stepsize_predictor!(integrator, alg, n_new) = stepsize_predictor!(integrator, integrator.controller_cache, alg, n_new) @inline function stepsize_controller!(integrator, controller::ExtrapolationController, From c49bb50e4500f6775d3835c1552c96fe4af114d2 Mon Sep 17 00:00:00 2001 From: termi-official Date: Mon, 27 Oct 2025 19:46:25 +0100 Subject: [PATCH 19/21] Fix composite algorithm choice --- .../src/perform_step/composite_perform_step.jl | 9 --------- 1 file changed, 9 deletions(-) diff --git a/lib/OrdinaryDiffEqCore/src/perform_step/composite_perform_step.jl b/lib/OrdinaryDiffEqCore/src/perform_step/composite_perform_step.jl index c3c101a47f..4ed9ff39ef 100644 --- a/lib/OrdinaryDiffEqCore/src/perform_step/composite_perform_step.jl +++ b/lib/OrdinaryDiffEqCore/src/perform_step/composite_perform_step.jl @@ -256,9 +256,6 @@ function choose_algorithm!(integrator, cache::DefaultCache) old_cache = cache.cache6 end - integrator.opts.controller.beta2 = beta2 = beta2_default(algs[new_current]) - integrator.opts.controller.beta1 = beta1_default(algs[new_current], beta2) - reset_alg_dependent_opts!(integrator, algs[old_current], algs[new_current]) transfer_cache!(integrator, old_cache, new_cache) end @@ -275,12 +272,6 @@ function reset_alg_dependent_opts!(integrator, alg1, alg2) if integrator.opts.adaptive == isadaptive(alg1) integrator.opts.adaptive = isadaptive(alg2) end - if integrator.opts.qmin == qmin_default(alg1) - integrator.opts.qmin = qmin_default(alg2) - end - if integrator.opts.qmax == qmax_default(alg1) - integrator.opts.qmax == qmax_default(alg2) - end reset_alg_dependent_opts!(integrator.opts.controller, alg1, alg2) nothing end From 1ce97c198897195131c9a4a80ab897539f6536cc Mon Sep 17 00:00:00 2001 From: termi-official Date: Mon, 27 Oct 2025 20:03:23 +0100 Subject: [PATCH 20/21] Missing dispatch --- .../src/perform_step/composite_perform_step.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/lib/OrdinaryDiffEqCore/src/perform_step/composite_perform_step.jl b/lib/OrdinaryDiffEqCore/src/perform_step/composite_perform_step.jl index 4ed9ff39ef..a0159c843f 100644 --- a/lib/OrdinaryDiffEqCore/src/perform_step/composite_perform_step.jl +++ b/lib/OrdinaryDiffEqCore/src/perform_step/composite_perform_step.jl @@ -280,3 +280,4 @@ end # Example: send the history variables from one multistep method to another transfer_cache!(integrator, alg1, alg2) = nothing +reset_alg_dependent_opts!(::Union{AbstractController, AbstractControllerCache}, alg1, alg2) = nothing From 54db5fb0d2b2dc23bec436b45312dbe2b7581d90 Mon Sep 17 00:00:00 2001 From: termi-official Date: Mon, 27 Oct 2025 20:27:35 +0100 Subject: [PATCH 21/21] Fix Verner JET --- lib/OrdinaryDiffEqVerner/src/OrdinaryDiffEqVerner.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/OrdinaryDiffEqVerner/src/OrdinaryDiffEqVerner.jl b/lib/OrdinaryDiffEqVerner/src/OrdinaryDiffEqVerner.jl index ee48edce7e..f9b7223400 100644 --- a/lib/OrdinaryDiffEqVerner/src/OrdinaryDiffEqVerner.jl +++ b/lib/OrdinaryDiffEqVerner/src/OrdinaryDiffEqVerner.jl @@ -4,7 +4,7 @@ import OrdinaryDiffEqCore: alg_order, calculate_residuals!, initialize!, perform_step!, @unpack, unwrap_alg, calculate_residuals, alg_stability_size, OrdinaryDiffEqAlgorithm, - CompositeAlgorithm, + CompositeAlgorithm, accept_step_controller, OrdinaryDiffEqMutableCache, OrdinaryDiffEqConstantCache, OrdinaryDiffEqAdaptiveAlgorithm, CompiledFloats, uses_uprev, alg_cache, _vec, _reshape, @cache, isfsal, full_cache,