Skip to content

Commit e152e20

Browse files
move BDF controller
1 parent 67485ac commit e152e20

File tree

3 files changed

+43
-42
lines changed

3 files changed

+43
-42
lines changed

lib/OrdinaryDiffEqBDF/src/controllers.jl

Lines changed: 41 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -100,6 +100,47 @@ function step_reject_controller!(integrator, ::FBDF)
100100
bdf_step_reject_controller!(integrator, integrator.cache.terkm1)
101101
end
102102

103+
function bdf_step_reject_controller!(integrator, EEst1)
104+
k = integrator.cache.order
105+
h = integrator.dt
106+
integrator.cache.consfailcnt += 1
107+
integrator.cache.nconsteps = 0
108+
if integrator.cache.consfailcnt > 1
109+
h = h / 2
110+
end
111+
zₛ = 1.2 # equivalent to integrator.opts.gamma
112+
expo = 1 / (k + 1)
113+
z = zₛ * ((integrator.EEst)^expo)
114+
F = inv(z)
115+
if z <= 10
116+
hₖ = F * h
117+
else # z > 10
118+
hₖ = 0.1 * h
119+
end
120+
hₙ = hₖ
121+
kₙ = k
122+
if k > 1
123+
expo = 1 / k
124+
zₖ₋₁ = 1.3 * (EEst1^expo)
125+
Fₖ₋₁ = inv(zₖ₋₁)
126+
if zₖ₋₁ <= 10
127+
hₖ₋₁ = Fₖ₋₁ * h
128+
elseif zₖ₋₁ > 10
129+
hₖ₋₁ = 0.1 * h
130+
end
131+
if integrator.cache.consfailcnt > 2 || hₖ₋₁ > hₖ
132+
hₙ = min(h, hₖ₋₁)
133+
kₙ = k - 1
134+
end
135+
end
136+
# Restart BDf (clear history) when we failed repeatedly
137+
if kₙ == 1 && integrator.cache.consfailcnt > 3
138+
u_modified!(integrator, true)
139+
end
140+
integrator.dt = hₙ
141+
integrator.cache.order = kₙ
142+
end
143+
103144
function post_newton_controller!(integrator, alg::FBDF)
104145
@unpack cache = integrator
105146
if cache.order > 1 && cache.nlsolver.nfails >= 3

src/OrdinaryDiffEq.jl

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -285,7 +285,8 @@ using ..OrdinaryDiffEqDefault
285285
export DefaultODEAlgorithm
286286

287287
using ..OrdinaryDiffEqBDF: reinitFBDF!, error_constant, estimate_terk!, calc_Lagrange_interp!,
288-
calc_finite_difference_weights, estimate_terk, calc_Lagrange_interp
288+
calc_finite_difference_weights, estimate_terk, calc_Lagrange_interp,
289+
bdf_step_reject_controller!
289290
include("nlsolve/newton.jl")
290291
include("perform_step/dae_perform_step.jl")
291292

src/integrators/controllers.jl

Lines changed: 0 additions & 41 deletions
Original file line numberDiff line numberDiff line change
@@ -479,47 +479,6 @@ function step_reject_controller!(integrator, ::DFBDF)
479479
bdf_step_reject_controller!(integrator, integrator.cache.terkm1)
480480
end
481481

482-
function bdf_step_reject_controller!(integrator, EEst1)
483-
k = integrator.cache.order
484-
h = integrator.dt
485-
integrator.cache.consfailcnt += 1
486-
integrator.cache.nconsteps = 0
487-
if integrator.cache.consfailcnt > 1
488-
h = h / 2
489-
end
490-
zₛ = 1.2 # equivalent to integrator.opts.gamma
491-
expo = 1 / (k + 1)
492-
z = zₛ * ((integrator.EEst)^expo)
493-
F = inv(z)
494-
if z <= 10
495-
hₖ = F * h
496-
else # z > 10
497-
hₖ = 0.1 * h
498-
end
499-
hₙ = hₖ
500-
kₙ = k
501-
if k > 1
502-
expo = 1 / k
503-
zₖ₋₁ = 1.3 * (EEst1^expo)
504-
Fₖ₋₁ = inv(zₖ₋₁)
505-
if zₖ₋₁ <= 10
506-
hₖ₋₁ = Fₖ₋₁ * h
507-
elseif zₖ₋₁ > 10
508-
hₖ₋₁ = 0.1 * h
509-
end
510-
if integrator.cache.consfailcnt > 2 || hₖ₋₁ > hₖ
511-
hₙ = min(h, hₖ₋₁)
512-
kₙ = k - 1
513-
end
514-
end
515-
# Restart BDf (clear history) when we failed repeatedly
516-
if kₙ == 1 && integrator.cache.consfailcnt > 3
517-
u_modified!(integrator, true)
518-
end
519-
integrator.dt = hₙ
520-
integrator.cache.order = kₙ
521-
end
522-
523482
function post_newton_controller!(integrator, alg)
524483
integrator.dt = integrator.dt / integrator.opts.failfactor
525484
nothing

0 commit comments

Comments
 (0)