diff --git a/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl b/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl index 9a0bb53267..a38ce1566a 100644 --- a/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl +++ b/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl @@ -400,3 +400,56 @@ function post_newton_controller!(integrator, alg) integrator.dt = integrator.dt / integrator.opts.failfactor nothing end + +@inline function stepsize_controller!(integrator, controller::PredictiveController, 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.fastpower(EEst, expo) / fac + @fastmath q = DiffEqBase.value(max(inv(qmax), min(inv(qmin), qtmp))) + integrator.qold = q + end + q +end + +function step_accept_controller!(integrator, controller::PredictiveController, alg, q) + @unpack qmin, qmax, gamma, qsteady_min, qsteady_max = integrator.opts + + EEst = DiffEqBase.value(integrator.EEst) + + if integrator.success_iter > 0 + expo = 1 / (get_current_adaptive_order(alg, integrator.cache) + 1) + qgus = (integrator.dtacc / integrator.dt) * + FastPower.fastpower((EEst^2) / integrator.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 + integrator.dtacc = integrator.dt + integrator.erracc = max(1e-2, EEst) + + return integrator.dt / qacc +end + +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 \ No newline at end of file diff --git a/lib/OrdinaryDiffEqFIRK/src/controllers.jl b/lib/OrdinaryDiffEqFIRK/src/controllers.jl index c8cb86f195..6de40914d0 100644 --- a/lib/OrdinaryDiffEqFIRK/src/controllers.jl +++ b/lib/OrdinaryDiffEqFIRK/src/controllers.jl @@ -1,52 +1,3 @@ -@inline function stepsize_controller!(integrator, controller::PredictiveController, 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.fastpower(EEst, expo) / fac - @fastmath q = DiffEqBase.value(max(inv(qmax), min(inv(qmin), qtmp))) - integrator.qold = q - end - q -end - -function step_accept_controller!(integrator, controller::PredictiveController, alg, q) - @unpack qmin, qmax, gamma, qsteady_min, qsteady_max = integrator.opts - - EEst = DiffEqBase.value(integrator.EEst) - - if integrator.success_iter > 0 - expo = 1 / (get_current_adaptive_order(alg, integrator.cache) + 1) - qgus = (integrator.dtacc / integrator.dt) * - FastPower.fastpower((EEst^2) / integrator.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 - integrator.dtacc = integrator.dt - integrator.erracc = max(1e-2, EEst) - - return integrator.dt / qacc -end - - function step_accept_controller!(integrator, controller::PredictiveController, alg::AdaptiveRadau, q) @unpack qmin, qmax, gamma, qsteady_min, qsteady_max = integrator.opts @unpack cache = integrator @@ -85,11 +36,6 @@ function step_accept_controller!(integrator, controller::PredictiveController, a return integrator.dt / qacc end -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 - function step_reject_controller!(integrator, controller::PredictiveController, alg::AdaptiveRadau) @unpack dt, success_iter, qold = integrator @unpack cache = integrator