Skip to content

Commit

Permalink
prelim adaptivity
Browse files Browse the repository at this point in the history
  • Loading branch information
@Shreyas-Ekanathan
Shreyas-Ekanathan committed Sep 24, 2024
1 parent 733bfa5 commit 4ac990c
Show file tree
Hide file tree
Showing 6 changed files with 214 additions and 76 deletions.
7 changes: 4 additions & 3 deletions lib/OrdinaryDiffEqFIRK/src/algorithms.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -163,12 +163,13 @@ struct AdaptiveRadau{CS, AD, F, P, FDT, ST, CJ, Tol, C1, C2, StepLimiter} <:
new_W_γdt_cutoff::C2
controller::Symbol
step_limiter!::StepLimiter
num_stages::Int
min_num_stages::Int
max_num_stages::Int
end

function AdaptiveRadau(; chunk_size = Val{0}(), autodiff = Val{true}(),
standardtag = Val{true}(), concrete_jac = nothing,
diff_type = Val{:forward}, num_stages = 3,
diff_type = Val{:forward}, min_num_stages = 3, max_num_stages = 3,
linsolve = nothing, precs = DEFAULT_PRECS,
extrapolant = :dense, fast_convergence_cutoff = 1 // 5,
new_W_γdt_cutoff = 1 // 5,
Expand All @@ -186,6 +187,6 @@ function AdaptiveRadau(; chunk_size = Val{0}(), autodiff = Val{true}(),
fast_convergence_cutoff,
new_W_γdt_cutoff,
controller,
step_limiter!, num_stages)
step_limiter!, min_num_stages, max_num_stages)
end

40 changes: 39 additions & 1 deletion lib/OrdinaryDiffEqFIRK/src/controllers.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,8 +24,9 @@
q
end

function step_accept_controller!(integrator, controller::PredictiveController, alg, q)
function step_accept_controller!(integrator, controller::PredictiveController, alg::Union{RadauIIA3, RadauIIA5, RadauIIA9}, q)
@unpack qmin, qmax, gamma, qsteady_min, qsteady_max = integrator.opts

EEst = DiffEqBase.value(integrator.EEst)

if integrator.success_iter > 0
Expand All @@ -42,6 +43,43 @@ function step_accept_controller!(integrator, controller::PredictiveController, a
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
@unpack num_stages, step, θ, θprev, orders = cache

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) *
DiffEqBase.fastpow((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)

cache.step = step + 1
if (step > 10)
Ψ = θ * θprev
if (Ψ <= 0.002 && num_stages < alg.max_num_stages)
cache.num_stages += 2
elseif ((Ψ >= 0.8 || cache.status == VerySlowConvergence || cache.status == Divergence) && num_stages > alg.min_num_stages)
cache.num_stages -= 2
cache.step = 1
end
end
return integrator.dt / qacc
end

Expand Down
90 changes: 60 additions & 30 deletions lib/OrdinaryDiffEqFIRK/src/firk_caches.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -477,7 +477,7 @@ end
mutable struct AdaptiveRadauConstantCache{F, Tab, Tol, Dt, U, JType} <:
OrdinaryDiffEqConstantCache
uf::F
tab::Tab
tabs::Vector{Tab}
κ::Tol
ηold::Tol
iter::Int
Expand All @@ -486,6 +486,11 @@ mutable struct AdaptiveRadauConstantCache{F, Tab, Tol, Dt, U, JType} <:
W_γdt::Dt
status::NLStatus
J::JType
num_stages::Int
step::Int
θ::BigFloat
θprev::BigFloat
orders::Vector{Int}
end

function alg_cache(alg::AdaptiveRadau, u, rate_prototype, ::Type{uEltypeNoUnits},
Expand All @@ -494,8 +499,16 @@ function alg_cache(alg::AdaptiveRadau, u, rate_prototype, ::Type{uEltypeNoUnits}
::Val{false}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
uf = UDerivativeWrapper(f, t, p)
uToltype = constvalue(uBottomEltypeNoUnits)
num_stages = alg.num_stages

num_stages = alg.min_num_stages
max = alg.max_num_stages
tabs = [BigRadauIIA5Tableau(uToltype, constvalue(tTypeNoUnits)), BigRadauIIA9Tableau(uToltype, constvalue(tTypeNoUnits)), BigRadauIIA13Tableau(uToltype, constvalue(tTypeNoUnits))]

i = 9
while i <= alg.max_num_stages
push!(tabs, adaptiveRadauTableau(uToltype, constvalue(tTypeNoUnits), i))
i += 2
end
#=
if (num_stages == 3)
tab = BigRadauIIA5Tableau(uToltype, constvalue(tTypeNoUnits))
elseif (num_stages == 5)
Expand All @@ -505,19 +518,19 @@ function alg_cache(alg::AdaptiveRadau, u, rate_prototype, ::Type{uEltypeNoUnits}
elseif iseven(num_stages) || num_stages <3
error("num_stages must be odd and 3 or greater")
else
tab = adaptiveRadauTableau(uToltype, constvalue(tTypeNoUnits), num_stages)
tab = adaptiveRadauTableau(uToltype, constvalue(tTypeNoUnits), min_num_stages)
end

cont = Vector{typeof(u)}(undef, num_stages)
for i in 1: num_stages
=#
cont = Vector{typeof(u)}(undef, max)
for i in 1: max
cont[i] = zero(u)
end

κ = alg.κ !== nothing ? convert(uToltype, alg.κ) : convert(uToltype, 1 // 100)
J = false .* _vec(rate_prototype) .* _vec(rate_prototype)'

AdaptiveRadauConstantCache(uf, tab, κ, one(uToltype), 10000, cont, dt, dt,
Convergence, J)
orders = [0,0,0,0]
AdaptiveRadauConstantCache(uf, tabs, κ, one(uToltype), 10000, cont, dt, dt,
Convergence, J, num_stages, 1, big"1.0", big"1.0", orders)
end

mutable struct AdaptiveRadauCache{uType, cuType, tType, uNoUnitsType, rateType, JType, W1Type, W2Type,
Expand All @@ -544,7 +557,7 @@ mutable struct AdaptiveRadauCache{uType, cuType, tType, uNoUnitsType, rateType,
W1::W1Type #real
W2::Vector{W2Type} #complex
uf::UF
tab::Tab
tabs::Vector{Tab}
κ::Tol
ηold::Tol
iter::Int
Expand All @@ -559,6 +572,11 @@ mutable struct AdaptiveRadauCache{uType, cuType, tType, uNoUnitsType, rateType,
W_γdt::Dt
status::NLStatus
step_limiter!::StepLimiter
num_stages::Int
step::Int
θ::BigFloat
θprev::BigFloat
orders::Vector{Int}
end

function alg_cache(alg::AdaptiveRadau, u, rate_prototype, ::Type{uEltypeNoUnits},
Expand All @@ -567,8 +585,19 @@ function alg_cache(alg::AdaptiveRadau, u, rate_prototype, ::Type{uEltypeNoUnits}
::Val{true}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
uf = UJacobianWrapper(f, t, p)
uToltype = constvalue(uBottomEltypeNoUnits)
num_stages = alg.num_stages

min = alg.min_num_stages
max = alg.max_num_stages

num_stages = min

tabs = [BigRadauIIA5Tableau(uToltype, constvalue(tTypeNoUnits)), BigRadauIIA9Tableau(uToltype, constvalue(tTypeNoUnits)), BigRadauIIA13Tableau(uToltype, constvalue(tTypeNoUnits))]
i = 9
while i <= max
push!(tabs, adaptiveRadauTableau(uToltype, constvalue(tTypeNoUnits), i))
i += 2
end
#=
if (num_stages == 3)
tab = BigRadauIIA5Tableau(uToltype, constvalue(tTypeNoUnits))
elseif (num_stages == 5)
Expand All @@ -580,39 +609,40 @@ function alg_cache(alg::AdaptiveRadau, u, rate_prototype, ::Type{uEltypeNoUnits}
else
tab = adaptiveRadauTableau(uToltype, constvalue(tTypeNoUnits), num_stages)
end
=#

κ = alg.κ !== nothing ? convert(uToltype, alg.κ) : convert(uToltype, 1 // 100)

z = Vector{typeof(u)}(undef, num_stages)
w = Vector{typeof(u)}(undef, num_stages)
for i in 1 : num_stages
z = Vector{typeof(u)}(undef, max)
w = Vector{typeof(u)}(undef, max)
for i in 1 : max
z[i] = w[i] = zero(u)
end

c_prime = Vector{typeof(t)}(undef, num_stages) #time stepping
c_prime = Vector{typeof(t)}(undef, max) #time stepping

dw1 = zero(u)
ubuff = zero(u)
dw2 = [similar(u, Complex{eltype(u)}) for _ in 1 : (num_stages - 1) ÷ 2]
dw2 = [similar(u, Complex{eltype(u)}) for _ in 1 : (max - 1) ÷ 2]
recursivefill!.(dw2, false)
cubuff = [similar(u, Complex{eltype(u)}) for _ in 1 : (num_stages - 1) ÷ 2]
cubuff = [similar(u, Complex{eltype(u)}) for _ in 1 : (max - 1) ÷ 2]
recursivefill!.(cubuff, false)
dw = Vector{typeof(u)}(undef, num_stages - 1)
dw = Vector{typeof(u)}(undef, max - 1)

cont = Vector{typeof(u)}(undef, num_stages)
for i in 1 : num_stages
cont = Vector{typeof(u)}(undef, max)
for i in 1 : max
cont[i] = zero(u)
end

derivatives = Matrix{typeof(u)}(undef, num_stages, num_stages)
for i in 1 : num_stages, j in 1 : num_stages
derivatives = Matrix{typeof(u)}(undef, max, max)
for i in 1 : max, j in 1 : max
derivatives[i, j] = zero(u)
end

fsalfirst = zero(rate_prototype)
fw = Vector{typeof(rate_prototype)}(undef, num_stages)
ks = Vector{typeof(rate_prototype)}(undef, num_stages)
for i in 1: num_stages
fw = Vector{typeof(rate_prototype)}(undef, max)
ks = Vector{typeof(rate_prototype)}(undef, max)
for i in 1: max
ks[i] = fw[i] = zero(rate_prototype)
end
k = ks[1]
Expand All @@ -622,7 +652,7 @@ function alg_cache(alg::AdaptiveRadau, u, rate_prototype, ::Type{uEltypeNoUnits}
error("Non-concrete Jacobian not yet supported by AdaptiveRadau.")
end

W2 = [similar(J, Complex{eltype(W1)}) for _ in 1 : (num_stages - 1) ÷ 2]
W2 = [similar(J, Complex{eltype(W1)}) for _ in 1 : (max - 1) ÷ 2]
recursivefill!.(W2, false)

du1 = zero(rate_prototype)
Expand All @@ -640,7 +670,7 @@ function alg_cache(alg::AdaptiveRadau, u, rate_prototype, ::Type{uEltypeNoUnits}

linsolve2 = [
init(LinearProblem(W2[i], _vec(cubuff[i]); u0 = _vec(dw2[i])), alg.linsolve, alias_A = true, alias_b = true,
assumptions = LinearSolve.OperatorAssumptions(true)) for i in 1 : (num_stages - 1) ÷ 2]
assumptions = LinearSolve.OperatorAssumptions(true)) for i in 1 : (max - 1) ÷ 2]

rtol = reltol isa Number ? reltol : zero(reltol)
atol = reltol isa Number ? reltol : zero(reltol)
Expand All @@ -649,9 +679,9 @@ function alg_cache(alg::AdaptiveRadau, u, rate_prototype, ::Type{uEltypeNoUnits}
z, w, c_prime, dw1, ubuff, dw2, cubuff, dw, cont, derivatives,
du1, fsalfirst, ks, k, fw,
J, W1, W2,
uf, tab, κ, one(uToltype), 10000, tmp,
uf, tabs, κ, one(uToltype), 10000, tmp,
atmp, jac_config,
linsolve1, linsolve2, rtol, atol, dt, dt,
Convergence, alg.step_limiter!)
Convergence, alg.step_limiter!, num_stages, 1, big"1.0", big"1.0", [0,0,0,0])
end

Loading

0 comments on commit 4ac990c

Please sign in to comment.