Merge pull request #4 from Shreyas-Ekanathan/master

update

src/OrdinaryDiffEq.jl

@@ -428,7 +428,7 @@ export FunctionMap, Euler, Heun, Ralston, Midpoint, RK4, ExplicitRK, OwrenZen3,
        FRK65, PFRK87,
        RKM, MSRK5, MSRK6, Stepanov5, SIR54, QPRK98, PSRK4p7q6, PSRK3p6q5, PSRK3p5q4
 
-export RadauIIA3, RadauIIA5, RadauIIA7
+export RadauIIA3, RadauIIA5, RadauIIA9
 
 export MagnusMidpoint, LinearExponential, MagnusLeapfrog, LieEuler, CayleyEuler,
        MagnusGauss4, MagnusNC6, MagnusGL6, MagnusGL8, MagnusNC8, MagnusGL4,

src/alg_utils.jl

@@ -175,7 +175,7 @@ qmin_default(alg::DP8) = 1 // 3
 qmax_default(alg::Union{OrdinaryDiffEqAlgorithm, DAEAlgorithm}) = 10
 qmax_default(alg::CompositeAlgorithm) = minimum(qmax_default.(alg.algs))
 qmax_default(alg::DP8) = 6
-qmax_default(alg::Union{RadauIIA3, RadauIIA5, RadauIIA7}) = 8
+qmax_default(alg::Union{RadauIIA3, RadauIIA5, RadauIIA9}) = 8
 
 function has_chunksize(alg::OrdinaryDiffEqAlgorithm)
     return alg isa Union{OrdinaryDiffEqExponentialAlgorithm,
@@ -430,7 +430,8 @@ alg_order(alg::TanYam7) = 7
 alg_order(alg::TsitPap8) = 8
 alg_order(alg::RadauIIA3) = 3
 alg_order(alg::RadauIIA5) = 5
-alg_order(alg::RadauIIA7) = 7
+alg_order(alg::RadauIIA9) = 9
+alg_order(alg::ImplicitEuler) = 1
 alg_order(alg::RKMK2) = 2
 alg_order(alg::RKMK4) = 4
 alg_order(alg::LieRK4) = 4
@@ -530,6 +531,7 @@ alg_adaptive_order(alg::Rosenbrock32) = 2
 
 alg_adaptive_order(alg::RadauIIA3) = 1
 alg_adaptive_order(alg::RadauIIA5) = 3
+alg_adaptive_order(alg::RadauIIA9) = 7
 
 # this is actually incorrect and is purposefully decreased as this tends
 # to track the real error much better

src/algorithms.jl

@@ -673,10 +673,10 @@ year={1999},
 publisher={Elsevier}
 }
 
-RadauIIA7: Fully-Implicit Runge-Kutta Method
+RadauII97: 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
@@ -690,15 +690,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,
@@ -712,7 +712,7 @@ function RadauIIA7(; chunk_size = Val{0}(), autodiff = Val{true}(),
         controller,
         step_limiter!)
 end
-TruncatedStacktraces.@truncate_stacktrace RadauIIA7
+TruncatedStacktraces.@truncate_stacktrace RadauIIA9
 
 ################################################################################
 

src/caches/firk_caches.jl

@@ -274,7 +274,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
@@ -291,22 +291,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,15 +370,15 @@ mutable struct RadauIIA7Cache{uType, cuType, uNoUnitsType, rateType, JType, W1Ty
     status::NLStatus
     step_limiter!::StepLimiter
 end
-TruncatedStacktraces.@truncate_stacktrace RadauIIA7Cache 1
+TruncatedStacktraces.@truncate_stacktrace RadauIIA9Cache 1
 
-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)
 
@@ -459,7 +459,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,

src/integrators/controllers.jl

@@ -406,7 +406,7 @@ end
         if fac_default_gamma(alg)
             fac = gamma
         else
-            if alg isa Union{RadauIIA3, RadauIIA5, RadauIIA7}
+            if alg isa Union{RadauIIA3, RadauIIA5, RadauIIA9}
                 @unpack iter = integrator.cache
                 @unpack maxiters = alg
             else

src/integrators/integrator_interface.jl

@@ -113,7 +113,7 @@ end
 end
 
 # avoid method ambiguity
-# for typ in (Union{RadauIIA3, RadauIIA5, RadauIIA7})
+# for typ in (Union{RadauIIA3, RadauIIA5, RadauIIA9})
 #     @eval @inline function DiffEqBase.get_tmp_cache(integrator, alg::$typ,
 #             cache::OrdinaryDiffEqConstantCache)
 #         nothing
@@ -126,7 +126,7 @@ end
     (cache.tmp,)
 end
 @inline function DiffEqBase.get_tmp_cache(
-        integrator, alg::Union{RadauIIA3, RadauIIA5, RadauIIA7},
+        integrator, alg::Union{RadauIIA3, RadauIIA5, RadauIIA9},
         cache::OrdinaryDiffEqMutableCache)
     (cache.tmp, cache.atmp)
 end

src/perform_step/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

src/tableaus/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,
@@ -258,3 +258,5 @@ function RadauIIA7Tableau(T, T2)
         γ, α1, β1, α2, β2,
         e1, e2, e3, e4, e5)
 end
+
+

test/algconvergence/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.75:-0.5:0.25), 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.75:-0.5:0.25), prob_ode_2Dlinear, RadauIIA0())
 @test sim21.𝒪est[:final]≈8 atol=testTol
 
 # test adaptivity