Merge pull request #6 from Shreyas-Ekanathan/upstream

Upstream

Author: Shreyas-Ekanathan

lib/OrdinaryDiffEqFIRK/src/OrdinaryDiffEqFIRK.jl

@@ -30,6 +30,6 @@ include("firk_tableaus.jl")
 include("firk_perform_step.jl")
 include("integrator_interface.jl")
 
-export RadauIIA3, RadauIIA5, RadauIIA7
+export RadauIIA3, RadauIIA5, RadauIIA9
 
 end

lib/OrdinaryDiffEqFIRK/src/alg_utils.jl

@@ -1,13 +1,13 @@
-qmax_default(alg::Union{RadauIIA3, RadauIIA5, RadauIIA7}) = 8
+qmax_default(alg::Union{RadauIIA3, RadauIIA5, RadauIIA9}) = 8
 
 alg_order(alg::RadauIIA3) = 3
 alg_order(alg::RadauIIA5) = 5
-alg_order(alg::RadauIIA7) = 7
+alg_order(alg::RadauIIA9) = 9
 
 isfirk(alg::RadauIIA3) = true
 isfirk(alg::RadauIIA5) = true
-isfirk(alg::RadauIIA7) = true
+isfirk(alg::RadauIIA9) = true
 
 alg_adaptive_order(alg::RadauIIA3) = 1
 alg_adaptive_order(alg::RadauIIA5) = 3
-alg_adaptive_order(alg::RadauIIA7) = 5
+alg_adaptive_order(alg::RadauIIA9) = 7

lib/OrdinaryDiffEqFIRK/src/algorithms.jl

@@ -114,10 +114,10 @@ year={1999},
 publisher={Elsevier}
 }
 
-RadauIIA7: Fully-Implicit Runge-Kutta Method
+RadauIIA9: 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, StepLimiter} <:
+struct RadauIIA9{CS, AD, F, P, FDT, ST, CJ, Tol, C1, C2, StepLimiter} <:
        OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ}
     linsolve::F
     precs::P
@@ -131,15 +131,15 @@ struct RadauIIA7{CS, AD, F, P, FDT, ST, CJ, Tol, C1, C2, StepLimiter} <:
     step_limiter!::StepLimiter
 end
 
-function RadauIIA7(; chunk_size = Val{0}(), autodiff = Val{true}(),
+function RadauIIA9(; 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,
         step_limiter! = trivial_limiter!)
-    RadauIIA7{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve),
+    RadauIIA9{_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,

lib/OrdinaryDiffEqFIRK/src/firk_caches.jl

@@ -273,7 +273,7 @@ function alg_cache(alg::RadauIIA5, u, rate_prototype, ::Type{uEltypeNoUnits},
         Convergence, alg.step_limiter!)
 end
 
-mutable struct RadauIIA7ConstantCache{F, Tab, Tol, Dt, U, JType} <:
+mutable struct RadauIIA9ConstantCache{F, Tab, Tol, Dt, U, JType} <:
                OrdinaryDiffEqConstantCache
     uf::F
     tab::Tab
@@ -290,22 +290,22 @@ mutable struct RadauIIA7ConstantCache{F, Tab, Tol, Dt, U, JType} <:
     J::JType
 end
 
-function alg_cache(alg::RadauIIA7, u, rate_prototype, ::Type{uEltypeNoUnits},
+function alg_cache(alg::RadauIIA9, u, rate_prototype, ::Type{uEltypeNoUnits},
         ::Type{uBottomEltypeNoUnits},
         ::Type{tTypeNoUnits}, uprev, uprev2, f, t, dt, reltol, p, calck,
         ::Val{false}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
     uf = UDerivativeWrapper(f, t, p)
     uToltype = constvalue(uBottomEltypeNoUnits)
-    tab = RadauIIA7Tableau(uToltype, constvalue(tTypeNoUnits))
+    tab = RadauIIA9Tableau(uToltype, constvalue(tTypeNoUnits))
 
     κ = 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,
+    RadauIIA9ConstantCache(uf, tab, κ, one(uToltype), 10000, u, u, u, u, dt, dt,
         Convergence, J)
 end
 
-mutable struct RadauIIA7Cache{uType, cuType, uNoUnitsType, rateType, JType, W1Type, W2Type,
+mutable struct RadauIIA9Cache{uType, cuType, uNoUnitsType, rateType, JType, W1Type, W2Type,
     UF, JC, F1, F2, Tab, Tol, Dt, rTol, aTol, StepLimiter} <:
                OrdinaryDiffEqMutableCache
     u::uType
@@ -370,13 +370,13 @@ mutable struct RadauIIA7Cache{uType, cuType, uNoUnitsType, rateType, JType, W1Ty
     step_limiter!::StepLimiter
 end
 
-function alg_cache(alg::RadauIIA7, u, rate_prototype, ::Type{uEltypeNoUnits},
+function alg_cache(alg::RadauIIA9, u, rate_prototype, ::Type{uEltypeNoUnits},
         ::Type{uBottomEltypeNoUnits},
         ::Type{tTypeNoUnits}, uprev, uprev2, f, t, dt, reltol, p, calck,
         ::Val{true}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
     uf = UJacobianWrapper(f, t, p)
     uToltype = constvalue(uBottomEltypeNoUnits)
-    tab = RadauIIA7Tableau(uToltype, constvalue(tTypeNoUnits))
+    tab = RadauIIA9Tableau(uToltype, constvalue(tTypeNoUnits))
 
     κ = alg.κ !== nothing ? convert(uToltype, alg.κ) : convert(uToltype, 1 // 100)
 
@@ -457,7 +457,7 @@ function alg_cache(alg::RadauIIA7, u, rate_prototype, ::Type{uEltypeNoUnits},
     rtol = reltol isa Number ? reltol : zero(reltol)
     atol = reltol isa Number ? reltol : zero(reltol)
 
-    RadauIIA7Cache(u, uprev,
+    RadauIIA9Cache(u, uprev,
         z1, z2, z3, z4, z5, w1, w2, w3, w4, w5,
         dw1, ubuff, dw23, dw45, cubuff1, cubuff2, cont1, cont2, cont3, cont4,
         du1, fsalfirst, k, k2, k3, k4, k5, fw1, fw2, fw3, fw4, fw5,

lib/OrdinaryDiffEqFIRK/src/firk_perform_step.jl

@@ -51,7 +51,7 @@ function initialize!(integrator, cache::RadauIIA5ConstantCache)
     nothing
 end
 
-function initialize!(integrator, cache::RadauIIA7ConstantCache)
+function initialize!(integrator, cache::RadauIIA9ConstantCache)
     integrator.kshortsize = 2
     integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
     integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t) # Pre-start fsal
@@ -98,7 +98,7 @@ function initialize!(integrator, cache::RadauIIA5Cache)
     nothing
 end
 
-function initialize!(integrator, cache::RadauIIA7Cache)
+function initialize!(integrator, cache::RadauIIA9Cache)
     integrator.kshortsize = 2
     integrator.fsalfirst = cache.fsalfirst
     integrator.fsallast = cache.k
@@ -784,7 +784,7 @@ end
     return
 end
 
-@muladd function perform_step!(integrator, cache::RadauIIA7ConstantCache,
+@muladd function perform_step!(integrator, cache::RadauIIA9ConstantCache,
         repeat_step = false)
     @unpack t, dt, uprev, u, f, p = integrator
     @unpack T11, T12, T13, T14, T15, T21, T22, T23, T24, T25, T31, T32, T33, T34, T35, T41, T42, T43, T44, T45, T51 = cache.tab #= T52 = 1, T53 = 0, T54 = 1, T55 = 0=#
@@ -1014,7 +1014,7 @@ end
     return
 end
 
-@muladd function perform_step!(integrator, cache::RadauIIA7Cache, repeat_step = false)
+@muladd function perform_step!(integrator, cache::RadauIIA9Cache, repeat_step = false)
     @unpack t, dt, uprev, u, f, p, fsallast, fsalfirst = integrator
     @unpack T11, T12, T13, T14, T15, T21, T22, T23, T24, T25, T31, T32, T33, T34, T35, T41, T42, T43, T44, T45, T51 = cache.tab #= T52 = 1, T53 = 0, T54 = 1, T55 = 0=#
     @unpack TI11, TI12, TI13, TI14, TI15, TI21, TI22, TI23, TI24, TI25, TI31, TI32, TI33, TI34, TI35, TI41, TI42, TI43, TI44, TI45, TI51, TI52, TI53, TI54, TI55 = cache.tab

lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl

@@ -111,7 +111,7 @@ function RadauIIA5Tableau(T, T2)
         e1, e2, e3)
 end
 
-struct RadauIIA7Tableau{T, T2}
+struct RadauIIA9Tableau{T, T2}
     T11::T
     T12::T
     T13::T
@@ -178,7 +178,7 @@ struct RadauIIA7Tableau{T, T2}
     e5::T
 end
 
-function RadauIIA7Tableau(T, T2)
+function RadauIIA9Tableau(T, T2)
     T11 = convert(T, -1.251758622050104589014e-2)
     T12 = convert(T, -1.024204781790882707009e-2)
     T13 = convert(T, 4.767387729029572386318e-2)
@@ -248,7 +248,7 @@ function RadauIIA7Tableau(T, T2)
     e4 = convert(T, 5.920031671845428725662e-1)
     e5 = convert(T, -2.000000000000000000000e-1)
 
-    RadauIIA7Tableau{T, T2}(T11, T12, T13, T14, T15,
+    RadauIIA9Tableau{T, T2}(T11, T12, T13, T14, T15,
         T21, T22, T23, T24, T25, T31, T32, T33, T34, T35,
         T41, T42, T43, T44, T45, T51, #=T52, T53, T54, T55=#
         TI11, TI12, TI13, TI14, TI15, TI21, TI22, TI23, TI24, TI25,

lib/OrdinaryDiffEqFIRK/src/integrator_interface.jl

@@ -1,4 +1,4 @@
-@inline function DiffEqBase.get_tmp_cache(integrator, alg::Union{RadauIIA3, RadauIIA5, RadauIIA7},
+@inline function DiffEqBase.get_tmp_cache(integrator, alg::Union{RadauIIA3, RadauIIA5, RadauIIA9},
         cache::OrdinaryDiffEqMutableCache)
     (cache.tmp, cache.atmp)
 end

lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl

@@ -8,10 +8,10 @@ for prob in [prob_ode_linear, prob_ode_2Dlinear]
     @test sim21.𝒪est[:final]≈5 atol=testTol
 end
 
-sim21 = test_convergence(1 ./ 2 .^ (2.777:-1:0.777), prob_ode_linear, RadauIIA7())
-@test sim21.𝒪est[:final]≈7 atol=testTol
+sim21 = test_convergence(1 ./ 2 .^ (2.5:-1:0.5), prob_ode_linear, RadauIIA9())
+@test sim21.𝒪est[:final]≈8 atol=testTol
 
-sim21 = test_convergence(1 ./ 2 .^ (2.777:-1:0.777), prob_ode_2Dlinear, RadauIIA7())
+sim21 = test_convergence(1 ./ 2 .^ (2.5:-1:0.5), prob_ode_2Dlinear, RadauIIA9())
 @test sim21.𝒪est[:final]≈8 atol=testTol
 
 # test adaptivity
@@ -57,4 +57,4 @@ for iip in (true, false)
         @test sol.stats.njacs < sol.stats.nw # W reuse
     end
     @test length(sol) < 5000 # the error estimate is not very good
-end
+end