small tweaks

Author: Shreyas-Ekanathan

src/algorithms.jl

@@ -1724,8 +1724,8 @@ publisher={Elsevier}
 RadauIIA7: Fully-Implicit Runge-Kutta Method
 An A-B-L stable fully implicit Runge-Kutta method with internal tableau complex basis transform for efficiency.
 """
-struct RadauIIA7{CS, AD, F, P, FDT, ST, CJ, Tol, C1, C2} <:
-       OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ}
+struct RadauIIA7{CS, AD, F, P, FDT, ST, CJ, Tol, C1, C2, StepLimiter} <:
+    OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ}
     linsolve::F
     precs::P
     smooth_est::Bool
@@ -1735,26 +1735,30 @@ struct RadauIIA7{CS, AD, F, P, FDT, ST, CJ, Tol, C1, C2} <:
     fast_convergence_cutoff::C1
     new_W_γdt_cutoff::C2
     controller::Symbol
+    step_limiter!::StepLimiter
 end
 
 function RadauIIA7(; chunk_size = Val{0}(), autodiff = Val{true}(),
-        standardtag = Val{true}(), concrete_jac = nothing,
-        diff_type = Val{:forward},
-        linsolve = nothing, precs = DEFAULT_PRECS,
-        extrapolant = :dense, fast_convergence_cutoff = 1 // 5,
-        new_W_γdt_cutoff = 1 // 5,
-        controller = :Predictive, κ = nothing, maxiters = 10, smooth_est = true)
-    RadauIIA7{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve),
-        typeof(precs), diff_type, _unwrap_val(standardtag), _unwrap_val(concrete_jac),
-        typeof(κ), typeof(fast_convergence_cutoff), typeof(new_W_γdt_cutoff)}(linsolve,
-        precs,
-        smooth_est,
-        extrapolant,
-        κ,
-        maxiters,
-        fast_convergence_cutoff,
-        new_W_γdt_cutoff,
-        controller)
+    standardtag = Val{true}(), concrete_jac = nothing,
+    diff_type = Val{:forward},
+    linsolve = nothing, precs = DEFAULT_PRECS,
+    extrapolant = :dense, fast_convergence_cutoff = 1 // 5,
+    new_W_γdt_cutoff = 1 // 5,
+    controller = :Predictive, κ = nothing, maxiters = 10, smooth_est = true,
+    step_limiter! = trivial_limiter!)
+RadauIIA7{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve),
+    typeof(precs), diff_type, _unwrap_val(standardtag), _unwrap_val(concrete_jac),
+    typeof(κ), typeof(fast_convergence_cutoff),
+    typeof(new_W_γdt_cutoff), typeof(step_limiter!)}(linsolve,
+    precs,
+    smooth_est,
+    extrapolant,
+    κ,
+    maxiters,
+    fast_convergence_cutoff,
+    new_W_γdt_cutoff,
+    controller,
+    step_limiter!)
 end
 TruncatedStacktraces.@truncate_stacktrace RadauIIA7
 

src/caches/firk_caches.jl

@@ -299,7 +299,7 @@ uf = UDerivativeWrapper(f, t, p)
 uToltype = constvalue(uBottomEltypeNoUnits)
 tab = RadauIIA7Tableau(uToltype, constvalue(tTypeNoUnits))
 
-κ = convert(uToltype, 1 // 100)
+κ = alg.κ !== nothing ? convert(uToltype, alg.κ) : convert(uToltype, 1 // 100)
 J = false .* _vec(rate_prototype) .* _vec(rate_prototype)'
 
 RadauIIA7ConstantCache(uf, tab, κ, one(uToltype), 10000, u, u, u, u, dt, dt,

src/integrators/controllers.jl

@@ -397,7 +397,7 @@ end
         if alg isa Union{RKC, IRKC, SERK2}
             fac = gamma
         else
-            if alg isa Union{RadauIIA3, RadauIIA5}
+            if alg isa Union{RadauIIA3, RadauIIA5, RadauIIA7}
                 @unpack iter = integrator.cache
                 @unpack maxiters = alg
             else

src/perform_step/firk_perform_step.jl

@@ -788,8 +788,8 @@ end
     @unpack maxiters = alg
     mass_matrix = integrator.f.mass_matrix
 
-    # precalculations
-    rtol = @.. broadcast=false reltol^(2 / 3)/10
+    # precalculations rtol pow is (num stages + 1)/(2*num stages)
+    rtol = @.. broadcast=false reltol#^(5/8)/10
     atol = @.. broadcast=false rtol*(abstol / reltol)
     c1m1 = c1 - 1
     c2m1 = c2 - 1
@@ -927,7 +927,7 @@ end
         z2 = @.. broadcast=false T21*w1 + T22*w2 + T23*w3 + T24*w4 + T25*w5
         z3 = @.. broadcast=false T31*w1 + T32*w2 + T33*w3 + T34*w4 + T35*w5
         z4 = @.. broadcast=false T41*w1 + T42*w2 + T43*w3 + T44*w4 + T45*w5
-        z5 = @.. broadcast=false T51*w1 + w2 + w4 #= T52=1, T53 = 0, T54 = 1, T55 = 0=#
+        z5 = @.. broadcast=false T51*w1 + w2 + w4 #= T52=1, T53=0, T54=1, T55=0 =#
 
         # check stopping criterion
         iter > 1 && (η = θ / (1 - θ))
@@ -974,16 +974,16 @@ end
     if integrator.EEst <= oneunit(integrator.EEst)
         cache.dtprev = dt
         if alg.extrapolant != :constant
-            cache.cont1 = @.. broadcast=false (z4 - z5) / c4m1 # first derivative on [c4, 1]
-            tmp1 = @.. broadcast=false (z3 - z4) / c3mc4 # first derivative on [c3, c4]
-            cache.cont2 = @.. broadcast=false (tmp1 - cache.cont1) / c3m1 # second derivative on [c3, 1]
-            tmp2 = @.. broadcast=false (z2 - z3) / c2mc3 # first derivative on [c2, c3]
-            tmp3 = @.. broadcast=false (tmp2 - tmp1) / c2mc4 # second derivative on [c2, c4]
-            cache.cont3 = @.. broadcast=false (tmp3 - cache.cont2) / c2m1 # third derivative on [c2, 1]
-            tmp4 = @.. broadcast=false (z1 - z2) / c1mc2 # first derivative on [c1, c2]
-            tmp5 = @.. broadcast=false (tmp4 - tmp2) / c1mc3 # second derivative on [c1, c3]
-            tmp6 = @.. broadcast=false (tmp5 - tmp3) / c1mc4 # third derivative on [c1, c4]
-            cache.cont4 = @.. broadcast=false (tmp6 - cache.cont3) / c1m1 #fourth derivative on [c1, 1]
+            cache.cont1 = @..  (z4 - z5) / c4m1 # first derivative on [c4, 1]
+            tmp1 = @..  (z3 - z4) / c3mc4 # first derivative on [c3, c4]
+            cache.cont2 = @..  (tmp1 - cache.cont1) / c3m1 # second derivative on [c3, 1]
+            tmp2 = @..  (z2 - z3) / c2mc3 # first derivative on [c2, c3]
+            tmp3 = @..  (tmp2 - tmp1) / c2mc4 # second derivative on [c2, c4]
+            cache.cont3 = @..  (tmp3 - cache.cont2) / c2m1 # third derivative on [c2, 1]
+            tmp4 = @..  (z1 - z2) / c1mc2 # first derivative on [c1, c2]
+            tmp5 = @..  (tmp4 - tmp2) / c1mc3 # second derivative on [c1, c3]
+            tmp6 = @..  (tmp5 - tmp3) / c1mc4 # third derivative on [c1, c4]
+            cache.cont4 = @..  (tmp6 - cache.cont3) / c1m1 #fourth derivative on [c1, 1]
         end
     end