Add Order Adaptivity to Radau

Project.toml

@@ -121,8 +121,8 @@ OrdinaryDiffEqQPRK = "1"
 OrdinaryDiffEqRKN = "1"
 OrdinaryDiffEqRosenbrock = "1"
 OrdinaryDiffEqSDIRK = "1"
-OrdinaryDiffEqStabilizedIRK = "1"
 OrdinaryDiffEqSSPRK = "1"
+OrdinaryDiffEqStabilizedIRK = "1"
 OrdinaryDiffEqStabilizedRK = "1"
 OrdinaryDiffEqSymplecticRK = "1"
 OrdinaryDiffEqTsit5 = "1"

lib/OrdinaryDiffEqFIRK/src/algorithms.jl

@@ -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 = 7, 
      linsolve = nothing, precs = DEFAULT_PRECS,
      extrapolant = :dense, fast_convergence_cutoff = 1 // 5,
      new_W_γdt_cutoff = 1 // 5,
@@ -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
 

lib/OrdinaryDiffEqFIRK/src/controllers.jl

@@ -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
@@ -42,6 +43,44 @@ 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 = 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
+    @show cache.num_stages
+    if (step > 10)
+        Ψ = θ * θprev
+        if (Ψ <= 0.001 && num_stages < alg.max_num_stages)
+            cache.num_stages += 2
+        elseif ((Ψ >= 0.1 || 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
 

lib/OrdinaryDiffEqFIRK/src/firk_caches.jl

@@ -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
@@ -486,6 +486,10 @@ mutable struct AdaptiveRadauConstantCache{F, Tab, Tol, Dt, U, JType} <:
     W_γdt::Dt
     status::NLStatus
     J::JType
+    num_stages::Int
+    step::Int
+    θ::BigFloat
+    θprev::BigFloat
 end
 
 function alg_cache(alg::AdaptiveRadau, u, rate_prototype, ::Type{uEltypeNoUnits},
@@ -494,30 +498,24 @@ 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
-
-    if (num_stages == 3)
-        tab = BigRadauIIA5Tableau(uToltype, constvalue(tTypeNoUnits))
-    elseif (num_stages == 5)
-        tab = BigRadauIIA9Tableau(uToltype, constvalue(tTypeNoUnits))
-    elseif (num_stages == 7)
-        tab = BigRadauIIA13Tableau(uToltype, constvalue(tTypeNoUnits))
-    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)
+    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
-
-    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)
+    AdaptiveRadauConstantCache(uf, tabs, κ, one(uToltype), 10000, cont, dt, dt,
+        Convergence, J, num_stages, 1, big"1.0", big"1.0")
 end
 
 mutable struct AdaptiveRadauCache{uType, cuType, tType, uNoUnitsType, rateType, JType, W1Type, W2Type,
@@ -544,7 +542,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
@@ -559,6 +557,10 @@ mutable struct AdaptiveRadauCache{uType, cuType, tType, uNoUnitsType, rateType,
     W_γdt::Dt
     status::NLStatus
     step_limiter!::StepLimiter
+    num_stages::Int
+    step::Int
+    θ::BigFloat
+    θprev::BigFloat
 end
 
 function alg_cache(alg::AdaptiveRadau, u, rate_prototype, ::Type{uEltypeNoUnits},
@@ -567,52 +569,51 @@ 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
-
-    if (num_stages == 3)
-        tab = BigRadauIIA5Tableau(uToltype, constvalue(tTypeNoUnits))
-    elseif (num_stages == 5)
-        tab = BigRadauIIA9Tableau(uToltype, constvalue(tTypeNoUnits))
-    elseif (num_stages == 7)
-        tab = BigRadauIIA13Tableau(uToltype, constvalue(tTypeNoUnits))
-    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)
+
+    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
 
     κ = 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]
@@ -622,7 +623,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)
@@ -640,7 +641,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)
@@ -649,9 +650,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")
 end