Skip to content

Commit 19312aa

Browse files
Refactor predictive controller
It shouldn't be in Radau, Radau has its own.
1 parent f65bcca commit 19312aa

File tree

2 files changed

+53
-54
lines changed

2 files changed

+53
-54
lines changed

lib/OrdinaryDiffEqCore/src/integrators/controllers.jl

Lines changed: 53 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -400,3 +400,56 @@ function post_newton_controller!(integrator, alg)
400400
integrator.dt = integrator.dt / integrator.opts.failfactor
401401
nothing
402402
end
403+
404+
@inline function stepsize_controller!(integrator, controller::PredictiveController, alg)
405+
@unpack qmin, qmax, gamma = integrator.opts
406+
EEst = DiffEqBase.value(integrator.EEst)
407+
if iszero(EEst)
408+
q = inv(qmax)
409+
else
410+
if fac_default_gamma(alg)
411+
fac = gamma
412+
else
413+
if isfirk(alg)
414+
@unpack iter = integrator.cache
415+
@unpack maxiters = alg
416+
else
417+
@unpack iter, maxiters = integrator.cache.nlsolver
418+
end
419+
fac = min(gamma, (1 + 2 * maxiters) * gamma / (iter + 2 * maxiters))
420+
end
421+
expo = 1 / (get_current_adaptive_order(alg, integrator.cache) + 1)
422+
qtmp = FastPower.fastpower(EEst, expo) / fac
423+
@fastmath q = DiffEqBase.value(max(inv(qmax), min(inv(qmin), qtmp)))
424+
integrator.qold = q
425+
end
426+
q
427+
end
428+
429+
function step_accept_controller!(integrator, controller::PredictiveController, alg, q)
430+
@unpack qmin, qmax, gamma, qsteady_min, qsteady_max = integrator.opts
431+
432+
EEst = DiffEqBase.value(integrator.EEst)
433+
434+
if integrator.success_iter > 0
435+
expo = 1 / (get_current_adaptive_order(alg, integrator.cache) + 1)
436+
qgus = (integrator.dtacc / integrator.dt) *
437+
FastPower.fastpower((EEst^2) / integrator.erracc, expo)
438+
qgus = max(inv(qmax), min(inv(qmin), qgus / gamma))
439+
qacc = max(q, qgus)
440+
else
441+
qacc = q
442+
end
443+
if qsteady_min <= qacc <= qsteady_max
444+
qacc = one(qacc)
445+
end
446+
integrator.dtacc = integrator.dt
447+
integrator.erracc = max(1e-2, EEst)
448+
449+
return integrator.dt / qacc
450+
end
451+
452+
function step_reject_controller!(integrator, controller::PredictiveController, alg)
453+
@unpack dt, success_iter, qold = integrator
454+
integrator.dt = success_iter == 0 ? 0.1 * dt : dt / qold
455+
end

lib/OrdinaryDiffEqFIRK/src/controllers.jl

Lines changed: 0 additions & 54 deletions
Original file line numberDiff line numberDiff line change
@@ -1,52 +1,3 @@
1-
@inline function stepsize_controller!(integrator, controller::PredictiveController, alg)
2-
@unpack qmin, qmax, gamma = integrator.opts
3-
EEst = DiffEqBase.value(integrator.EEst)
4-
if iszero(EEst)
5-
q = inv(qmax)
6-
else
7-
if fac_default_gamma(alg)
8-
fac = gamma
9-
else
10-
if isfirk(alg)
11-
@unpack iter = integrator.cache
12-
@unpack maxiters = alg
13-
else
14-
@unpack iter, maxiters = integrator.cache.nlsolver
15-
end
16-
fac = min(gamma, (1 + 2 * maxiters) * gamma / (iter + 2 * maxiters))
17-
end
18-
expo = 1 / (get_current_adaptive_order(alg, integrator.cache) + 1)
19-
qtmp = FastPower.fastpower(EEst, expo) / fac
20-
@fastmath q = DiffEqBase.value(max(inv(qmax), min(inv(qmin), qtmp)))
21-
integrator.qold = q
22-
end
23-
q
24-
end
25-
26-
function step_accept_controller!(integrator, controller::PredictiveController, alg, q)
27-
@unpack qmin, qmax, gamma, qsteady_min, qsteady_max = integrator.opts
28-
29-
EEst = DiffEqBase.value(integrator.EEst)
30-
31-
if integrator.success_iter > 0
32-
expo = 1 / (get_current_adaptive_order(alg, integrator.cache) + 1)
33-
qgus = (integrator.dtacc / integrator.dt) *
34-
FastPower.fastpower((EEst^2) / integrator.erracc, expo)
35-
qgus = max(inv(qmax), min(inv(qmin), qgus / gamma))
36-
qacc = max(q, qgus)
37-
else
38-
qacc = q
39-
end
40-
if qsteady_min <= qacc <= qsteady_max
41-
qacc = one(qacc)
42-
end
43-
integrator.dtacc = integrator.dt
44-
integrator.erracc = max(1e-2, EEst)
45-
46-
return integrator.dt / qacc
47-
end
48-
49-
501
function step_accept_controller!(integrator, controller::PredictiveController, alg::AdaptiveRadau, q)
512
@unpack qmin, qmax, gamma, qsteady_min, qsteady_max = integrator.opts
523
@unpack cache = integrator
@@ -85,11 +36,6 @@ function step_accept_controller!(integrator, controller::PredictiveController, a
8536
return integrator.dt / qacc
8637
end
8738

88-
function step_reject_controller!(integrator, controller::PredictiveController, alg)
89-
@unpack dt, success_iter, qold = integrator
90-
integrator.dt = success_iter == 0 ? 0.1 * dt : dt / qold
91-
end
92-
9339
function step_reject_controller!(integrator, controller::PredictiveController, alg::AdaptiveRadau)
9440
@unpack dt, success_iter, qold = integrator
9541
@unpack cache = integrator

0 commit comments

Comments
 (0)