From b1c691586410a2aed32fe7d12e7633a174c0ad1c Mon Sep 17 00:00:00 2001
From: CompatHelper Julia <compathelper_noreply@julialang.org>
Date: Fri, 20 Sep 2024 12:33:05 +0000
Subject: [PATCH 01/10] CompatHelper: add new compat entry for OrdinaryDiffEq
 at version 6 for package downstream, (keep existing compat)

---
 test/downstream/Project.toml | 1 +
 1 file changed, 1 insertion(+)

diff --git a/test/downstream/Project.toml b/test/downstream/Project.toml
index 345d90f382..4ba47ff201 100644
--- a/test/downstream/Project.toml
+++ b/test/downstream/Project.toml
@@ -12,6 +12,7 @@ Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 DDEProblemLibrary = "0.1"
 DelayDiffEq = "5.42"
 Measurements = "2.9"
+OrdinaryDiffEq = "6"
 SciMLSensitivity = "7.30"
 StaticArrays = "1"
 StochasticDiffEq = "6.60.1"

From 4ac990ca40e826660d75a615332bbe7800173e71 Mon Sep 17 00:00:00 2001
From: Shreyas Ekanathan <shreyase39@gmail.com>
Date: Tue, 24 Sep 2024 19:00:44 -0400
Subject: [PATCH 02/10] prelim adaptivity

---
 lib/OrdinaryDiffEqFIRK/src/algorithms.jl      |  7 +-
 lib/OrdinaryDiffEqFIRK/src/controllers.jl     | 40 ++++++++-
 lib/OrdinaryDiffEqFIRK/src/firk_caches.jl     | 90 ++++++++++++-------
 .../src/firk_perform_step.jl                  | 60 ++++++++-----
 lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl   | 39 ++++----
 lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl | 54 ++++++++++-
 6 files changed, 214 insertions(+), 76 deletions(-)

diff --git a/lib/OrdinaryDiffEqFIRK/src/algorithms.jl b/lib/OrdinaryDiffEqFIRK/src/algorithms.jl
index 25165422c8..a209eb8799 100644
--- a/lib/OrdinaryDiffEqFIRK/src/algorithms.jl
+++ b/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 = 3, 
      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
 
diff --git a/lib/OrdinaryDiffEqFIRK/src/controllers.jl b/lib/OrdinaryDiffEqFIRK/src/controllers.jl
index df6cf153ab..d4f721e60c 100644
--- a/lib/OrdinaryDiffEqFIRK/src/controllers.jl
+++ b/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,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
 
diff --git a/lib/OrdinaryDiffEqFIRK/src/firk_caches.jl b/lib/OrdinaryDiffEqFIRK/src/firk_caches.jl
index af6da2a282..380695be0d 100644
--- a/lib/OrdinaryDiffEqFIRK/src/firk_caches.jl
+++ b/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,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},
@@ -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)
@@ -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,
@@ -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
@@ -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},
@@ -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)
@@ -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]
@@ -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)
@@ -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)
@@ -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
 
diff --git a/lib/OrdinaryDiffEqFIRK/src/firk_perform_step.jl b/lib/OrdinaryDiffEqFIRK/src/firk_perform_step.jl
index 044257afad..9b936fd5df 100644
--- a/lib/OrdinaryDiffEqFIRK/src/firk_perform_step.jl
+++ b/lib/OrdinaryDiffEqFIRK/src/firk_perform_step.jl
@@ -90,7 +90,7 @@ function initialize!(integrator, cache::RadauIIA9Cache)
 end
 
 function initialize!(integrator, cache::AdaptiveRadauCache)
-    @unpack num_stages = cache.tab
+    @unpack num_stages = cache
     integrator.kshortsize = 2
     resize!(integrator.k, integrator.kshortsize)
     integrator.k[1] = integrator.fsalfirst
@@ -1354,8 +1354,10 @@ end
 @muladd function perform_step!(integrator, cache::AdaptiveRadauConstantCache,
     repeat_step = false)
     @unpack t, dt, uprev, u, f, p = integrator
-    @unpack T, TI, γ, α, β, c, e, num_stages = cache.tab
-    @unpack κ, cont = cache
+    @unpack tabs, num_stages = cache
+    tab = tabs[(num_stages - 1) ÷ 2]
+    @unpack T, TI, γ, α, β, c, e = tab
+    @unpack κ, cont, θ, θprev = cache
     @unpack internalnorm, abstol, reltol, adaptive = integrator.opts
     alg = unwrap_alg(integrator, true)
     @unpack maxiters = alg
@@ -1398,7 +1400,7 @@ end
             cache.cont[i] = @.. map(zero, u)
         end
     else
-        c_prime = Vector{typeof(u)}(undef, num_stages) #time stepping
+        c_prime = Vector{typeof(dt)}(undef, num_stages) #time stepping
         c_prime[num_stages] = @.. dt / cache.dtprev
         for i in 1 : num_stages - 1
             c_prime[i] = @.. c[i] * c_prime[num_stages]
@@ -1435,7 +1437,7 @@ end
         for i in 1 : num_stages
             ff[i] = f(uprev + z[i], p, t + c[i] * dt)
         end
-        integrator.stats.nf += num_stages
+        OrdinaryDiffEqCore.increment_nf!(integrator.stats, 5)
 
         #fw = TI * ff
         fw = Vector{typeof(u)}(undef, num_stages)
@@ -1483,7 +1485,9 @@ end
         # check divergence (not in initial step)
 
         if iter > 1
+            cache.θprev = θ
             θ = ndw / ndwprev
+            cache.θ = θ
             (diverge = θ > 1) && (cache.status = Divergence)
             (veryslowconvergence = ndw * θ^(maxiters - iter) > κ * (1 - θ)) &&
                 (cache.status = VerySlowConvergence)
@@ -1535,9 +1539,13 @@ end
     
     if adaptive
         edt = e ./ dt
-        tmp = dot(edt, z)
+        tmp = 0
+        for i in 1 : num_stages
+            tmp = @.. tmp + edt[i] * z[i]
+        end
         mass_matrix != I && (tmp = mass_matrix * tmp)
-        utilde = @.. broadcast=false 1 / γ * dt * integrator.fsalfirst+tmp
+        #utilde = @.. broadcast=false 1 / γ * dt * integrator.fsalfirst + tmp
+        utilde = @.. broadcast=false integrator.fsalfirst + tmp
         if alg.smooth_est
             utilde = _reshape(LU1 \ _vec(utilde), axes(u))
             integrator.stats.nsolve += 1
@@ -1548,8 +1556,9 @@ end
         if !(integrator.EEst < oneunit(integrator.EEst)) && integrator.iter == 1 ||
         integrator.u_modified
             f0 = f(uprev .+ utilde, p, t)
-            integrator.stats.nf += 1
-            utilde = @.. broadcast=false 1 / γ * dt * f0 + tmp
+            OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1)
+            #utilde = @.. broadcast=false 1 / γ * dt * f0 + tmp
+            utilde = @.. broadcast=false f0+tmp
             if alg.smooth_est
                 utilde = _reshape(LU1 \ _vec(utilde), axes(u))
                 integrator.stats.nsolve += 1
@@ -1589,11 +1598,13 @@ end
 
 @muladd function perform_step!(integrator, cache::AdaptiveRadauCache, repeat_step = false)
     @unpack t, dt, uprev, u, f, p, fsallast, fsalfirst = integrator
-    @unpack T, TI, γ, α, β, c, e, num_stages = cache.tab
+    @unpack num_stages, tabs = cache
+    tab = tabs[(num_stages - 1) ÷ 2]
+    @unpack T, TI, γ, α, β, c, e = tab
     @unpack κ, cont, derivatives, z, w, c_prime = cache
     @unpack dw1, ubuff, dw2, cubuff, dw = cache
     @unpack ks, k, fw, J, W1, W2 = cache
-    @unpack tmp, atmp, jac_config, linsolve1, linsolve2, rtol, atol, step_limiter! = cache
+    @unpack tmp, atmp, jac_config, linsolve1, linsolve2, rtol, atol, step_limiter!, θ, θprev = cache
     @unpack internalnorm, abstol, reltol, adaptive = integrator.opts
     alg = unwrap_alg(integrator, true)
     @unpack maxiters = alg
@@ -1659,7 +1670,7 @@ end
             @.. tmp = uprev + z[i]
             f(ks[i], tmp, p, t + c[i] * dt)
         end
-        integrator.stats.nf += num_stages
+        OrdinaryDiffEqCore.increment_nf!(integrator.stats, 5)
 
         #mul!(fw, TI, ks)
         for i in 1:num_stages
@@ -1729,7 +1740,9 @@ end
         # check divergence (not in initial step)
 
         if iter > 1
+            cache.θprev = θ
             θ = ndw / ndwprev
+            cache.θ = θ
             (diverge = θ > 1) && (cache.status = Divergence)
             (veryslowconvergence = ndw * θ^(maxiters - iter) > κ * (1 - θ)) &&
                 (cache.status = VerySlowConvergence)
@@ -1782,16 +1795,20 @@ end
     step_limiter!(u, integrator, p, t + dt)
     
     if adaptive
-        utilde = w2
+        utilde = w[2]
         edt = e./dt
-        @.. tmp = dot(edt, z) + 1 / γ * dt * fsalfirst
-        mass_matrix != I && (mul!(w1, mass_matrix, tmp); copyto!(tmp, w1))
-        @.. ubuff=integrator.fsalfirst + tmp
+        @.. tmp = 0
+        for i in 1 : num_stages
+            @.. tmp += edt[i] * z[i]
+        end
+        mass_matrix != I && (mul!(w[1], mass_matrix, tmp); copyto!(tmp, w[1]))
+        #@.. ubuff=1 / γ * dt * integrator.fsalfirst + tmp
+        @.. broadcast=false ubuff=integrator.fsalfirst + tmp
 
         if alg.smooth_est
-            linres1 = dolinsolve(integrator, linres1.cache; b = _vec(ubuff),
+            linres = dolinsolve(integrator, linres.cache; b = _vec(ubuff),
                 linu = _vec(utilde))
-            cache.linsolve1 = linres1.cache
+            cache.linsolve1 = linres.cache
             integrator.stats.nsolve += 1
         end
 
@@ -1804,13 +1821,14 @@ end
            integrator.u_modified
             @.. broadcast=false utilde=uprev + utilde
             f(fsallast, utilde, p, t)
-            integrator.stats.nf += 1
+            OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1)
+            #@.. broadcast=false ubuff = 1 / γ * dt * fsallast + tmp
             @.. broadcast=false ubuff=fsallast + tmp
 
             if alg.smooth_est
-                linres1 = dolinsolve(integrator, linres1.cache; b = _vec(ubuff),
+                linres = dolinsolve(integrator, linres.cache; b = _vec(ubuff),
                     linu = _vec(utilde))
-                cache.linsolve1 = linres1.cache
+                cache.linsolve1 = linres.cache
                 integrator.stats.nsolve += 1
             end
 
diff --git a/lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl b/lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl
index 750a9a0b34..28d94854d2 100644
--- a/lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl
+++ b/lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl
@@ -119,7 +119,6 @@ struct RadauIIATableau{T1, T2}
     α::Vector{T1}
     β::Vector{T1}
     e::Vector{T1}
-    num_stages::Int
 end
 
 function BigRadauIIA5Tableau(T1, T2)
@@ -135,9 +134,9 @@ function BigRadauIIA5Tableau(T1, T2)
     c[3] = big"1"
 
     e = Vector{T1}(undef, 3)
-    e[1] = big"-0.428298294115368098113417591057340987284723986878723769598436582629558514867595819404541110575367601847354683220540647741034052880983125451920841851066713815"
-    e[2] = big"0.245039074384916438547779982042524963814308131639809054920681599475247366387530452311400517264433163448821319862243292031101333835037542080594323438762615605"
-    e[3] = big"-0.0916296098652259003910911359018870983677439465358993542342383692664256286871842771687905889462502859518430848371143253189423645209875225522045944109790661043"
+    e[1] = big"-10.048809399827415944628228317014873027801513671875"
+    e[2] = big"1.3821427331607492039466933420044369995594024658203125"
+    e[3] = big"-0.333333333333333314829616256247390992939472198486328125"
 
     TI = Matrix{T1}(undef, 3, 3)
     TI[1, 1] = big"4.32557989006315535102435095295614882731995158490590784287320458848019483341979047442263696495019938973156007686663488090615420049217658854859024016717169837"
@@ -161,7 +160,7 @@ function BigRadauIIA5Tableau(T1, T2)
     T[3, 2] = big"1.0"
     T[3, 3] = big"0.0"
     RadauIIATableau{T1, T2}(T, TI,  
-    c, γ, α, β, e, 3)
+    c, γ, α, β, e)
 end
 
 function BigRadauIIA9Tableau(T1, T2)
@@ -181,11 +180,11 @@ function BigRadauIIA9Tableau(T1, T2)
     c[5] = big"1.0"
 
     e = Vector{T1}(undef, 5)
-    e[1] = big"-0.239909571163200476817707991076962793618683493830916562279975225042872448414819259070978815977101189851237591634144816820425592764432710089981892023441549743"
-    e[2] = big"0.125293484229223300606887443525929336197638450194973243323835481816625682669684634271160851996902445722310139152840852195000603355301430153561918341655137084"
-    e[3] = big"-0.0688288849083782089370741375422279772873399871312158026536514369967825732836040693366396751988019623495452730460577336089848105495733304937016477629621990433"
-    e[4] = big"0.0372433600301293198267284480468961585078575499935477902539140092558239369583143189611737274891328175087350382605569395937945872776839987254373869550943049195"
-    e[5] = big"-0.012863950751139890646895902730137465239579845437088427263654957605488133042638147254426913683751171160691603649073170415735165315443538347036196755548109703"
+    e[1] = big"27.78093394406463730479"
+    e[2] = big"3.641478498049213152712"
+    e[3] = big"-1.252547721169118720491"
+    e[4] = big"0.5920031671845428725662"
+    e[5] = big"-0.2000000000000000000000"
 
     TI = Matrix{T1}(undef, 5, 5)
     TI[1, 1] = big"30.0415677215444016277146611632467970747634862837368422955138470463852339244593400023985957753164599415374157317627305099177616927640413043608408838747985125"
@@ -242,7 +241,7 @@ function BigRadauIIA9Tableau(T1, T2)
     T[5, 5] = big"0.0"
 
     RadauIIATableau{T1, T2}(T, TI,  
-    c, γ, α, β, e, 5)
+    c, γ, α, β, e)
 end
 
 
@@ -415,13 +414,13 @@ function BigRadauIIA13Tableau(T1, T2)
     c[7] = big"1.0"
 
     e = Vector{T1}(undef, 7)
-    e[1] = big"-0.171003707892600662399316289094713451418682228802554671094543075316220821099884263713032871578607576486535539488632407200766379971279557791263375813632287147"
-    e[2] = big"0.0934967172358652400317534533028674569308657324394316331629203486361371292312231403973668315582870547753526899857449840409175610009437530537219068836721686211"
-    e[3] = big"-0.0538908303114758775848180855518003793385120454908028879947132475960997563222416509683480350817114056356343433378213334826358824351525243402758389616051172681"
-    e[4] = big"0.03036786965048439581219923157368590250090409822952169396779792168990510618756404452728392288892998298088147691907128776816545685599760715439221674418662785"
-    e[5] = big"-0.0169792974425458224044481617230998766694942329759644144506911809530904808476835739189122151558601434810772520378036293579816345384682687602578758514350075723"
-    e[6] = big"0.00942688256820236884916415666439281573527695349338346787075550606528435808025071461023926432308932314282041885090975780812273515256740094297311541275151861292"
-    e[7] = big"-0.00331409873565629283448601269346047459594635696619041493081994712789731442072563377354737487903843138987115421498455722938021358621090485566426506726181005806"
+    e[1] = big"-54.37443689412861451458"
+    e[2] = big"7.000024004259186512041"
+    e[3] = big"-2.355661091987557192256"
+    e[4] = big"1.132289066106134386384"
+    e[5] = big"-0.6468913267673587118673"
+    e[6] = big"0.3875333853753523774248"
+    e[7] = big"-0.1428571428571428571429"
 
     TI = Matrix{T1}(undef, 7, 7)
     TI[1, 1] = big"258.131926319982229276108947425184471333411128774462923076434633414645220927977539758484670571338176678808837829326061674950321562391576244286310404028770676"
@@ -526,7 +525,7 @@ function BigRadauIIA13Tableau(T1, T2)
     T[7, 7] = big"0.0"
 
     RadauIIATableau{T1, T2}(T, TI,  
-    c, γ, α, β, e, 7)
+    c, γ, α, β, e)
 end
 
 using Polynomials, LinearAlgebra, GenericSchur, RootedTrees, Symbolics
@@ -629,5 +628,5 @@ function adaptiveRadauTableau(T1, T2, num_stages::Int)
         b_hat = Symbolics.expand.((AA \ -bb))
         e = [Symbolics.symbolic_to_float(b_hat[i] - b[i]) for i in 1 : num_stages]    
     end
-    RadauIIATableau{T1, T2}(T, TI, c, γ, α, β, e, num_stages)
+    RadauIIATableau{T1, T2}(T, TI, c, γ, α, β, e)
 end
diff --git a/lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl b/lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl
index 2eb857827d..8ad32f5e80 100644
--- a/lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl
+++ b/lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl
@@ -19,10 +19,62 @@ prob_ode_2Dlinear_big = remake(prob_ode_2Dlinear, u0 = big.(prob_ode_2Dlinear.u0
 
 for i in [3, 5, 7, 9], prob in [prob_ode_linear_big, prob_ode_2Dlinear_big]
     dts = 1 ./ 2 .^ (4.25:-1:0.25)
-    sim21 = test_convergence(dts, prob, AdaptiveRadau(num_stages = i))
+    sim21 = test_convergence(dts, prob, AdaptiveRadau(min_num_stages = i, max_num_stages = i))
     @test sim21.𝒪est[:final]≈ (2 * i - 1) atol=testTol
 end
 
+solve(prob_ode_linear, AdaptiveRadau(min_num_stages = 3, max_num_stages = 3))
+solve(prob_ode_linear, RadauIIA5())
+
+using OrdinaryDiffEq, StaticArrays, LinearSolve, ParameterizedFunctions
+
+hires = @ode_def Hires begin
+    dy1 = -1.71 * y1 + 0.43 * y2 + 8.32 * y3 + 0.0007
+    dy2 = 1.71 * y1 - 8.75 * y2
+    dy3 = -10.03 * y3 + 0.43 * y4 + 0.035 * y5
+    dy4 = 8.32 * y2 + 1.71 * y3 - 1.12 * y4
+    dy5 = -1.745 * y5 + 0.43 * y6 + 0.43 * y7
+    dy6 = -280.0 * y6 * y8 + 0.69 * y4 + 1.71 * y5 - 0.43 * y6 + 0.69 * y7
+    dy7 = 280.0 * y6 * y8 - 1.81 * y7
+    dy8 = -280.0 * y6 * y8 + 1.81 * y7
+end
+
+u0 = SA[1, 0, 0, 0, 0, 0, 0, 0.0057]
+probiip = ODEProblem{true}(hires, Vector(u0), (0.0, 10.0))
+proboop = ODEProblem{false}(hires, Vector(u0), (0.0, 10.0))
+proboop = ODEProblem{false}(hires, u0, (0.0, 10.0))
+
+#=@btime =# sol = solve(proboop, AdaptiveRadau(min_num_stages = 3, max_num_stages = 7), reltol = 1e-1)
+#=@btime =# sol = solve(proboop, RadauIIA5())
+
+function rober!(du, u, p, t)
+    y₁, y₂, y₃ = u
+    k₁, k₂, k₃ = p
+    du[1] = -k₁ * y₁ + k₃ * y₂ * y₃
+    du[2] = k₁ * y₁ - k₂ * y₂^2 - k₃ * y₂ * y₃
+    du[3] = k₂ * y₂^2
+    nothing
+end
+prob = ODEProblem(rober!, [1.0, 0.0, 0.0], (0.0, 1e5), [0.04, 3e7, 1e4])
+sol = solve(prob, AdaptiveRadau(min_num_stages = 3, max_num_stages =7))
+sol2 = solve(prob, RadauIIA5())
+
+using BenchmarkTools
+#oop
+@btime solve(prob_ode_linear, RadauIIA5(), adaptive = false, dt = 1e-2)
+@btime solve(prob_ode_linear, AdaptiveRadau(num_stages = 3), adaptive = false, dt = 1e-2)
+
+#ip
+@btime solve(prob_ode_2Dlinear, RadauIIA5(), adaptive = false, dt = 1e-2)
+@btime solve(prob_ode_2Dlinear, AdaptiveRadau(num_stages = 3), adaptive = false, dt = 1e-2)
+
+VSCodeServer.@profview solve(prob_ode_linear, AdaptiveRadau(num_stages = 3), adaptive = false, dt = 1e-5)
+VSCodeServer.@profview solve(prob_ode_linear, RadauIIA5(), adaptive = false, dt = 1e-5)
+
+
+solve(prob_ode_linear, AdaptiveRadau(num_stages = 3), adaptive = false, dt = 1e-2)
+
+
 # test adaptivity
 for iip in (true, false)
     if iip

From 7a58a9e4d5cd2a313bdfd30a84047b559c1cab2e Mon Sep 17 00:00:00 2001
From: Shreyas Ekanathan <shreyase39@gmail.com>
Date: Thu, 26 Sep 2024 21:17:04 -0400
Subject: [PATCH 03/10] tweaks

---
 Project.toml                                  |  2 +-
 lib/OrdinaryDiffEqFIRK/src/algorithms.jl      |  2 +-
 lib/OrdinaryDiffEqFIRK/src/controllers.jl     |  7 +--
 lib/OrdinaryDiffEqFIRK/src/firk_caches.jl     | 33 +-----------
 lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl   |  4 +-
 lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl | 52 +------------------
 6 files changed, 11 insertions(+), 89 deletions(-)

diff --git a/Project.toml b/Project.toml
index db7de365a2..69e498bc76 100644
--- a/Project.toml
+++ b/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"
diff --git a/lib/OrdinaryDiffEqFIRK/src/algorithms.jl b/lib/OrdinaryDiffEqFIRK/src/algorithms.jl
index a209eb8799..552130eca3 100644
--- a/lib/OrdinaryDiffEqFIRK/src/algorithms.jl
+++ b/lib/OrdinaryDiffEqFIRK/src/algorithms.jl
@@ -169,7 +169,7 @@ end
 
 function AdaptiveRadau(; chunk_size = Val{0}(), autodiff = Val{true}(),
      standardtag = Val{true}(), concrete_jac = nothing,
-     diff_type = Val{:forward}, min_num_stages = 3, max_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,
diff --git a/lib/OrdinaryDiffEqFIRK/src/controllers.jl b/lib/OrdinaryDiffEqFIRK/src/controllers.jl
index d4f721e60c..cc4fb097c2 100644
--- a/lib/OrdinaryDiffEqFIRK/src/controllers.jl
+++ b/lib/OrdinaryDiffEqFIRK/src/controllers.jl
@@ -51,7 +51,7 @@ 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
+    @unpack num_stages, step, θ, θprev = cache
  
     EEst = DiffEqBase.value(integrator.EEst)
 
@@ -71,11 +71,12 @@ function step_accept_controller!(integrator, controller::PredictiveController, a
     integrator.erracc = max(1e-2, EEst)
     
     cache.step = step + 1
+    @show cache.num_stages
     if (step > 10)
         Ψ = θ * θprev
-        if (Ψ <= 0.002 && num_stages < alg.max_num_stages)
+        if (Ψ <= 0.001 && 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)
+        elseif ((Ψ >= 0.1 || cache.status == VerySlowConvergence || cache.status == Divergence) && num_stages > alg.min_num_stages)
             cache.num_stages -= 2 
             cache.step = 1
         end
diff --git a/lib/OrdinaryDiffEqFIRK/src/firk_caches.jl b/lib/OrdinaryDiffEqFIRK/src/firk_caches.jl
index 380695be0d..5e15d40b29 100644
--- a/lib/OrdinaryDiffEqFIRK/src/firk_caches.jl
+++ b/lib/OrdinaryDiffEqFIRK/src/firk_caches.jl
@@ -490,7 +490,6 @@ mutable struct AdaptiveRadauConstantCache{F, Tab, Tol, Dt, U, JType} <:
     step::Int
     θ::BigFloat
     θprev::BigFloat
-    orders::Vector{Int}
 end
 
 function alg_cache(alg::AdaptiveRadau, u, rate_prototype, ::Type{uEltypeNoUnits},
@@ -508,19 +507,6 @@ function alg_cache(alg::AdaptiveRadau, u, rate_prototype, ::Type{uEltypeNoUnits}
         push!(tabs, adaptiveRadauTableau(uToltype, constvalue(tTypeNoUnits), i))
         i += 2
     end
-    #=
-    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), min_num_stages)
-    end
-    =#
     cont = Vector{typeof(u)}(undef, max)
     for i in 1: max
         cont[i] = zero(u)
@@ -528,9 +514,8 @@ function alg_cache(alg::AdaptiveRadau, u, rate_prototype, ::Type{uEltypeNoUnits}
 
     κ = alg.κ !== nothing ? convert(uToltype, alg.κ) : convert(uToltype, 1 // 100)
     J = false .* _vec(rate_prototype) .* _vec(rate_prototype)'
-    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)
+        Convergence, J, num_stages, 1, big"1.0", big"1.0")
 end
 
 mutable struct AdaptiveRadauCache{uType, cuType, tType, uNoUnitsType, rateType, JType, W1Type, W2Type,
@@ -576,7 +561,6 @@ mutable struct AdaptiveRadauCache{uType, cuType, tType, uNoUnitsType, rateType,
     step::Int
     θ::BigFloat
     θprev::BigFloat
-    orders::Vector{Int}
 end
 
 function alg_cache(alg::AdaptiveRadau, u, rate_prototype, ::Type{uEltypeNoUnits},
@@ -597,19 +581,6 @@ function alg_cache(alg::AdaptiveRadau, u, rate_prototype, ::Type{uEltypeNoUnits}
         push!(tabs, adaptiveRadauTableau(uToltype, constvalue(tTypeNoUnits), i))
         i += 2
     end
-    #=
-    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)
-    end
-    =#
 
     κ = alg.κ !== nothing ? convert(uToltype, alg.κ) : convert(uToltype, 1 // 100)
 
@@ -682,6 +653,6 @@ function alg_cache(alg::AdaptiveRadau, u, rate_prototype, ::Type{uEltypeNoUnits}
         uf, tabs, κ, one(uToltype), 10000, tmp,
         atmp, jac_config,
         linsolve1, linsolve2, rtol, atol, dt, dt,
-        Convergence, alg.step_limiter!, num_stages, 1, big"1.0", big"1.0", [0,0,0,0])
+        Convergence, alg.step_limiter!, num_stages, 1, big"1.0", big"1.0")
 end
 
diff --git a/lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl b/lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl
index 28d94854d2..49b51d298a 100644
--- a/lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl
+++ b/lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl
@@ -136,7 +136,7 @@ function BigRadauIIA5Tableau(T1, T2)
     e = Vector{T1}(undef, 3)
     e[1] = big"-10.048809399827415944628228317014873027801513671875"
     e[2] = big"1.3821427331607492039466933420044369995594024658203125"
-    e[3] = big"-0.333333333333333314829616256247390992939472198486328125"
+    e[3] = big"-0.333333333333333333333333333333333333333333333333333"
 
     TI = Matrix{T1}(undef, 3, 3)
     TI[1, 1] = big"4.32557989006315535102435095295614882731995158490590784287320458848019483341979047442263696495019938973156007686663488090615420049217658854859024016717169837"
@@ -180,7 +180,7 @@ function BigRadauIIA9Tableau(T1, T2)
     c[5] = big"1.0"
 
     e = Vector{T1}(undef, 5)
-    e[1] = big"27.78093394406463730479"
+    e[1] = big"-27.78093394406463730479"
     e[2] = big"3.641478498049213152712"
     e[3] = big"-1.252547721169118720491"
     e[4] = big"0.5920031671845428725662"
diff --git a/lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl b/lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl
index 8ad32f5e80..ca4edc9ff9 100644
--- a/lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl
+++ b/lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl
@@ -23,57 +23,7 @@ for i in [3, 5, 7, 9], prob in [prob_ode_linear_big, prob_ode_2Dlinear_big]
     @test sim21.𝒪est[:final]≈ (2 * i - 1) atol=testTol
 end
 
-solve(prob_ode_linear, AdaptiveRadau(min_num_stages = 3, max_num_stages = 3))
-solve(prob_ode_linear, RadauIIA5())
-
-using OrdinaryDiffEq, StaticArrays, LinearSolve, ParameterizedFunctions
-
-hires = @ode_def Hires begin
-    dy1 = -1.71 * y1 + 0.43 * y2 + 8.32 * y3 + 0.0007
-    dy2 = 1.71 * y1 - 8.75 * y2
-    dy3 = -10.03 * y3 + 0.43 * y4 + 0.035 * y5
-    dy4 = 8.32 * y2 + 1.71 * y3 - 1.12 * y4
-    dy5 = -1.745 * y5 + 0.43 * y6 + 0.43 * y7
-    dy6 = -280.0 * y6 * y8 + 0.69 * y4 + 1.71 * y5 - 0.43 * y6 + 0.69 * y7
-    dy7 = 280.0 * y6 * y8 - 1.81 * y7
-    dy8 = -280.0 * y6 * y8 + 1.81 * y7
-end
-
-u0 = SA[1, 0, 0, 0, 0, 0, 0, 0.0057]
-probiip = ODEProblem{true}(hires, Vector(u0), (0.0, 10.0))
-proboop = ODEProblem{false}(hires, Vector(u0), (0.0, 10.0))
-proboop = ODEProblem{false}(hires, u0, (0.0, 10.0))
-
-#=@btime =# sol = solve(proboop, AdaptiveRadau(min_num_stages = 3, max_num_stages = 7), reltol = 1e-1)
-#=@btime =# sol = solve(proboop, RadauIIA5())
-
-function rober!(du, u, p, t)
-    y₁, y₂, y₃ = u
-    k₁, k₂, k₃ = p
-    du[1] = -k₁ * y₁ + k₃ * y₂ * y₃
-    du[2] = k₁ * y₁ - k₂ * y₂^2 - k₃ * y₂ * y₃
-    du[3] = k₂ * y₂^2
-    nothing
-end
-prob = ODEProblem(rober!, [1.0, 0.0, 0.0], (0.0, 1e5), [0.04, 3e7, 1e4])
-sol = solve(prob, AdaptiveRadau(min_num_stages = 3, max_num_stages =7))
-sol2 = solve(prob, RadauIIA5())
-
-using BenchmarkTools
-#oop
-@btime solve(prob_ode_linear, RadauIIA5(), adaptive = false, dt = 1e-2)
-@btime solve(prob_ode_linear, AdaptiveRadau(num_stages = 3), adaptive = false, dt = 1e-2)
-
-#ip
-@btime solve(prob_ode_2Dlinear, RadauIIA5(), adaptive = false, dt = 1e-2)
-@btime solve(prob_ode_2Dlinear, AdaptiveRadau(num_stages = 3), adaptive = false, dt = 1e-2)
-
-VSCodeServer.@profview solve(prob_ode_linear, AdaptiveRadau(num_stages = 3), adaptive = false, dt = 1e-5)
-VSCodeServer.@profview solve(prob_ode_linear, RadauIIA5(), adaptive = false, dt = 1e-5)
-
-
-solve(prob_ode_linear, AdaptiveRadau(num_stages = 3), adaptive = false, dt = 1e-2)
-
+solve(prob_ode_2Dlinear, AdaptiveRadau())
 
 # test adaptivity
 for iip in (true, false)

From 775d6ba56abe38c913595cab817dd1c9e7ee0aba Mon Sep 17 00:00:00 2001
From: Shreyas Ekanathan <shreyase39@gmail.com>
Date: Thu, 26 Sep 2024 21:20:47 -0400
Subject: [PATCH 04/10] Update ode_firk_tests.jl

---
 lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl | 2 --
 1 file changed, 2 deletions(-)

diff --git a/lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl b/lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl
index ca4edc9ff9..1dc482ffa3 100644
--- a/lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl
+++ b/lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl
@@ -23,8 +23,6 @@ for i in [3, 5, 7, 9], prob in [prob_ode_linear_big, prob_ode_2Dlinear_big]
     @test sim21.𝒪est[:final]≈ (2 * i - 1) atol=testTol
 end
 
-solve(prob_ode_2Dlinear, AdaptiveRadau())
-
 # test adaptivity
 for iip in (true, false)
     if iip

From 5ff81ed26a673f5f7a8f2c15ff3deb5ff7cea64e Mon Sep 17 00:00:00 2001
From: Shreyas-Ekanathan <142109039+Shreyas-Ekanathan@users.noreply.github.com>
Date: Thu, 26 Sep 2024 21:22:41 -0400
Subject: [PATCH 05/10] Update Project.toml

---
 test/downstream/Project.toml | 1 -
 1 file changed, 1 deletion(-)

diff --git a/test/downstream/Project.toml b/test/downstream/Project.toml
index 4ba47ff201..abd479bb5f 100644
--- a/test/downstream/Project.toml
+++ b/test/downstream/Project.toml
@@ -2,7 +2,6 @@
 DDEProblemLibrary = "f42792ee-6ffc-4e2a-ae83-8ee2f22de800"
 DelayDiffEq = "bcd4f6db-9728-5f36-b5f7-82caef46ccdb"
 Measurements = "eff96d63-e80a-5855-80a2-b1b0885c5ab7"
-OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 SciMLSensitivity = "1ed8b502-d754-442c-8d5d-10ac956f44a1"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0"

From 125e88a04d690dcc0b876ccda5d304daf1a77828 Mon Sep 17 00:00:00 2001
From: Shreyas-Ekanathan <142109039+Shreyas-Ekanathan@users.noreply.github.com>
Date: Thu, 26 Sep 2024 21:23:20 -0400
Subject: [PATCH 06/10] Update Project.toml

---
 test/downstream/Project.toml | 1 +
 1 file changed, 1 insertion(+)

diff --git a/test/downstream/Project.toml b/test/downstream/Project.toml
index abd479bb5f..6f934b345f 100644
--- a/test/downstream/Project.toml
+++ b/test/downstream/Project.toml
@@ -2,6 +2,7 @@
 DDEProblemLibrary = "f42792ee-6ffc-4e2a-ae83-8ee2f22de800"
 DelayDiffEq = "bcd4f6db-9728-5f36-b5f7-82caef46ccdb"
 Measurements = "eff96d63-e80a-5855-80a2-b1b0885c5ab7"
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"	
 SciMLSensitivity = "1ed8b502-d754-442c-8d5d-10ac956f44a1"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0"

From 377e8ed0cd1ca5b62e84c1b4148cd9b93ace1a41 Mon Sep 17 00:00:00 2001
From: Shreyas Ekanathan <shreyase39@gmail.com>
Date: Fri, 27 Sep 2024 07:23:56 -0400
Subject: [PATCH 07/10] fixes

---
 lib/OrdinaryDiffEqFIRK/src/algorithms.jl      |  8 ++++----
 lib/OrdinaryDiffEqFIRK/src/controllers.jl     |  6 +++---
 lib/OrdinaryDiffEqFIRK/src/firk_caches.jl     | 10 +++++-----
 lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl |  2 +-
 4 files changed, 13 insertions(+), 13 deletions(-)

diff --git a/lib/OrdinaryDiffEqFIRK/src/algorithms.jl b/lib/OrdinaryDiffEqFIRK/src/algorithms.jl
index 552130eca3..75ccae073c 100644
--- a/lib/OrdinaryDiffEqFIRK/src/algorithms.jl
+++ b/lib/OrdinaryDiffEqFIRK/src/algorithms.jl
@@ -163,13 +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
- min_num_stages::Int
- max_num_stages::Int
+ min_stages::Int
+ max_stages::Int
 end
 
 function AdaptiveRadau(; chunk_size = Val{0}(), autodiff = Val{true}(),
      standardtag = Val{true}(), concrete_jac = nothing,
-     diff_type = Val{:forward}, min_num_stages = 3, max_num_stages = 7, 
+     diff_type = Val{:forward}, min_stages = 3, max_stages = 7, 
      linsolve = nothing, precs = DEFAULT_PRECS,
      extrapolant = :dense, fast_convergence_cutoff = 1 // 5,
      new_W_γdt_cutoff = 1 // 5,
@@ -187,6 +187,6 @@ function AdaptiveRadau(; chunk_size = Val{0}(), autodiff = Val{true}(),
      fast_convergence_cutoff,
      new_W_γdt_cutoff,
      controller,
-     step_limiter!, min_num_stages, max_num_stages)
+     step_limiter!, min_stages, max_stages)
 end
 
diff --git a/lib/OrdinaryDiffEqFIRK/src/controllers.jl b/lib/OrdinaryDiffEqFIRK/src/controllers.jl
index cc4fb097c2..404f51dad6 100644
--- a/lib/OrdinaryDiffEqFIRK/src/controllers.jl
+++ b/lib/OrdinaryDiffEqFIRK/src/controllers.jl
@@ -24,7 +24,7 @@
     q
 end
 
-function step_accept_controller!(integrator, controller::PredictiveController, alg::Union{RadauIIA3, RadauIIA5, RadauIIA9}, q)
+function step_accept_controller!(integrator, controller::PredictiveController, alg, q)
     @unpack qmin, qmax, gamma, qsteady_min, qsteady_max = integrator.opts
  
     EEst = DiffEqBase.value(integrator.EEst)
@@ -74,9 +74,9 @@ function step_accept_controller!(integrator, controller::PredictiveController, a
     @show cache.num_stages
     if (step > 10)
         Ψ = θ * θprev
-        if (Ψ <= 0.001 && num_stages < alg.max_num_stages)
+        if (Ψ <= 0.001 && num_stages < alg.max_stages)
             cache.num_stages += 2
-        elseif ((Ψ >= 0.1 || cache.status == VerySlowConvergence || cache.status == Divergence) && num_stages > alg.min_num_stages)
+        elseif ((Ψ >= 0.1 || cache.status == VerySlowConvergence || cache.status == Divergence) && num_stages > alg.min_stages)
             cache.num_stages -= 2 
             cache.step = 1
         end
diff --git a/lib/OrdinaryDiffEqFIRK/src/firk_caches.jl b/lib/OrdinaryDiffEqFIRK/src/firk_caches.jl
index 5e15d40b29..d31fc7c7e9 100644
--- a/lib/OrdinaryDiffEqFIRK/src/firk_caches.jl
+++ b/lib/OrdinaryDiffEqFIRK/src/firk_caches.jl
@@ -498,12 +498,12 @@ 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.min_num_stages
-    max = alg.max_num_stages
+    num_stages = alg.min_stages
+    max = alg.max_stages
     tabs = [BigRadauIIA5Tableau(uToltype, constvalue(tTypeNoUnits)), BigRadauIIA9Tableau(uToltype, constvalue(tTypeNoUnits)), BigRadauIIA13Tableau(uToltype, constvalue(tTypeNoUnits))]
     
     i = 9
-    while i <= alg.max_num_stages
+    while i <= alg.max_stages
         push!(tabs, adaptiveRadauTableau(uToltype, constvalue(tTypeNoUnits), i))
         i += 2
     end
@@ -570,8 +570,8 @@ function alg_cache(alg::AdaptiveRadau, u, rate_prototype, ::Type{uEltypeNoUnits}
     uf = UJacobianWrapper(f, t, p)
     uToltype = constvalue(uBottomEltypeNoUnits)
 
-    min = alg.min_num_stages
-    max = alg.max_num_stages
+    min = alg.min_stages
+    max = alg.max_stages
 
     num_stages = min
 
diff --git a/lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl b/lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl
index 1dc482ffa3..eac935ec5f 100644
--- a/lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl
+++ b/lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl
@@ -19,7 +19,7 @@ prob_ode_2Dlinear_big = remake(prob_ode_2Dlinear, u0 = big.(prob_ode_2Dlinear.u0
 
 for i in [3, 5, 7, 9], prob in [prob_ode_linear_big, prob_ode_2Dlinear_big]
     dts = 1 ./ 2 .^ (4.25:-1:0.25)
-    sim21 = test_convergence(dts, prob, AdaptiveRadau(min_num_stages = i, max_num_stages = i))
+    sim21 = test_convergence(dts, prob, AdaptiveRadau(min_stages = i, max_stages = i))
     @test sim21.𝒪est[:final]≈ (2 * i - 1) atol=testTol
 end
 

From 58baa983b28f94f2e650463b2d357210b91c8810 Mon Sep 17 00:00:00 2001
From: Shreyas Ekanathan <shreyase39@gmail.com>
Date: Sun, 20 Oct 2024 12:05:55 -0400
Subject: [PATCH 08/10] edits, adaptivity fixed

---
 lib/OrdinaryDiffEqFIRK/src/alg_utils.jl       |  7 +-
 lib/OrdinaryDiffEqFIRK/src/controllers.jl     | 33 ++++++--
 lib/OrdinaryDiffEqFIRK/src/firk_caches.jl     | 39 ++++-----
 .../src/firk_perform_step.jl                  | 79 +++++++++----------
 lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl   |  4 +-
 5 files changed, 89 insertions(+), 73 deletions(-)

diff --git a/lib/OrdinaryDiffEqFIRK/src/alg_utils.jl b/lib/OrdinaryDiffEqFIRK/src/alg_utils.jl
index 428393e0bf..ef91903e63 100644
--- a/lib/OrdinaryDiffEqFIRK/src/alg_utils.jl
+++ b/lib/OrdinaryDiffEqFIRK/src/alg_utils.jl
@@ -1,9 +1,9 @@
-qmax_default(alg::Union{RadauIIA3, RadauIIA5, RadauIIA9}) = 8
+qmax_default(alg::Union{RadauIIA3, RadauIIA5, RadauIIA9, AdaptiveRadau}) = 8
 
 alg_order(alg::RadauIIA3) = 3
 alg_order(alg::RadauIIA5) = 5
 alg_order(alg::RadauIIA9) = 9
-alg_order(alg::AdaptiveRadau) = 5
+alg_order(alg::AdaptiveRadau) = 5 #dummy value
 
 isfirk(alg::RadauIIA3) = true
 isfirk(alg::RadauIIA5) = true
@@ -13,3 +13,6 @@ isfirk(alg::AdaptiveRadau) = true
 alg_adaptive_order(alg::RadauIIA3) = 1
 alg_adaptive_order(alg::RadauIIA5) = 3
 alg_adaptive_order(alg::RadauIIA9) = 5
+
+get_current_alg_order(alg::AdaptiveRadau, cache) = cache.num_stages * 2 - 1
+get_current_adaptive_order(alg::AdaptiveRadau, cache) = cache.num_stages
diff --git a/lib/OrdinaryDiffEqFIRK/src/controllers.jl b/lib/OrdinaryDiffEqFIRK/src/controllers.jl
index 404f51dad6..b91856d25d 100644
--- a/lib/OrdinaryDiffEqFIRK/src/controllers.jl
+++ b/lib/OrdinaryDiffEqFIRK/src/controllers.jl
@@ -1,7 +1,6 @@
 @inline function stepsize_controller!(integrator, controller::PredictiveController, alg)
     @unpack qmin, qmax, gamma = integrator.opts
     EEst = DiffEqBase.value(integrator.EEst)
-
     if iszero(EEst)
         q = inv(qmax)
     else
@@ -51,7 +50,7 @@ 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
+    @unpack num_stages, step, iter, hist_iter = cache
  
     EEst = DiffEqBase.value(integrator.EEst)
 
@@ -69,16 +68,18 @@ function step_accept_controller!(integrator, controller::PredictiveController, a
     end
     integrator.dtacc = integrator.dt
     integrator.erracc = max(1e-2, EEst)
-    
     cache.step = step + 1
-    @show cache.num_stages
+    hist_iter = hist_iter * 0.8 + iter * 0.2
+    cache.hist_iter = hist_iter
     if (step > 10)
-        Ψ = θ * θprev
-        if (Ψ <= 0.001 && num_stages < alg.max_stages)
+        if (hist_iter < 2.6 && num_stages < alg.max_stages)
             cache.num_stages += 2
-        elseif ((Ψ >= 0.1 || cache.status == VerySlowConvergence || cache.status == Divergence) && num_stages > alg.min_stages)
+            cache.step = 1
+            cache.hist_iter = iter
+        elseif ((hist_iter > 8 || cache.status == VerySlowConvergence || cache.status == Divergence) && num_stages > alg.min_stages)
             cache.num_stages -= 2 
             cache.step = 1
+            cache.hist_iter = iter
         end
     end
     return integrator.dt / qacc
@@ -88,3 +89,21 @@ function step_reject_controller!(integrator, controller::PredictiveController, a
     @unpack dt, success_iter, qold = integrator
     integrator.dt = success_iter == 0 ? 0.1 * dt : dt / qold
 end
+
+function step_reject_controller!(integrator, controller::PredictiveController, alg::AdaptiveRadau)
+    @unpack dt, success_iter, qold = integrator
+    @unpack cache = integrator 
+    @unpack num_stages, step, iter, hist_iter = cache
+    integrator.dt = success_iter == 0 ? 0.1 * dt : dt / qold
+    cache.step = step + 1
+    hist_iter = hist_iter * 0.8 + iter * 0.2
+    cache.hist_iter = hist_iter
+    if (step > 10)
+        if ((hist_iter > 8 || cache.status == VerySlowConvergence || cache.status == Divergence) && num_stages > alg.min_stages)
+            cache.num_stages -= 2 
+            cache.step = 1
+            cache.hist_iter = iter
+        end
+    end
+end
+
diff --git a/lib/OrdinaryDiffEqFIRK/src/firk_caches.jl b/lib/OrdinaryDiffEqFIRK/src/firk_caches.jl
index d31fc7c7e9..4f4045a897 100644
--- a/lib/OrdinaryDiffEqFIRK/src/firk_caches.jl
+++ b/lib/OrdinaryDiffEqFIRK/src/firk_caches.jl
@@ -488,8 +488,7 @@ mutable struct AdaptiveRadauConstantCache{F, Tab, Tol, Dt, U, JType} <:
     J::JType
     num_stages::Int
     step::Int
-    θ::BigFloat
-    θprev::BigFloat
+    hist_iter::Float64
 end
 
 function alg_cache(alg::AdaptiveRadau, u, rate_prototype, ::Type{uEltypeNoUnits},
@@ -515,11 +514,11 @@ function alg_cache(alg::AdaptiveRadau, u, rate_prototype, ::Type{uEltypeNoUnits}
     κ = alg.κ !== nothing ? convert(uToltype, alg.κ) : convert(uToltype, 1 // 100)
     J = false .* _vec(rate_prototype) .* _vec(rate_prototype)'
     AdaptiveRadauConstantCache(uf, tabs, κ, one(uToltype), 10000, cont, dt, dt,
-        Convergence, J, num_stages, 1, big"1.0", big"1.0")
+        Convergence, J, num_stages, 1, 0.0)
 end
 
 mutable struct AdaptiveRadauCache{uType, cuType, tType, uNoUnitsType, rateType, JType, W1Type, W2Type,
-    UF, JC, F1, F2, Tab, Tol, Dt, rTol, aTol, StepLimiter} <:
+    UF, JC, F1, F2, #=F3,=# Tab, Tol, Dt, rTol, aTol, StepLimiter} <:
                FIRKMutableCache
     u::uType
     uprev::uType
@@ -551,6 +550,7 @@ mutable struct AdaptiveRadauCache{uType, cuType, tType, uNoUnitsType, rateType,
     jac_config::JC
     linsolve1::F1 #real
     linsolve2::Vector{F2} #complex
+    #linres2::Vector{F3} 
     rtol::rTol
     atol::aTol
     dtprev::Dt
@@ -559,8 +559,7 @@ mutable struct AdaptiveRadauCache{uType, cuType, tType, uNoUnitsType, rateType,
     step_limiter!::StepLimiter
     num_stages::Int
     step::Int
-    θ::BigFloat
-    θprev::BigFloat
+    hist_iter::Float64
 end
 
 function alg_cache(alg::AdaptiveRadau, u, rate_prototype, ::Type{uEltypeNoUnits},
@@ -598,12 +597,9 @@ function alg_cache(alg::AdaptiveRadau, u, rate_prototype, ::Type{uEltypeNoUnits}
     recursivefill!.(dw2, false)
     cubuff = [similar(u, Complex{eltype(u)}) for _ in 1 : (max - 1) ÷ 2]
     recursivefill!.(cubuff, false)
-    dw = Vector{typeof(u)}(undef, max - 1)
+    dw = [zero(u) for i in 1 : max]
 
-    cont = Vector{typeof(u)}(undef, max)
-    for i in 1 : max
-        cont[i] = zero(u)
-    end
+    cont = [zero(u) for i in 1:max]
 
     derivatives = Matrix{typeof(u)}(undef, max, max)
     for i in 1 : max, j in 1 : max
@@ -611,11 +607,9 @@ function alg_cache(alg::AdaptiveRadau, u, rate_prototype, ::Type{uEltypeNoUnits}
     end
 
     fsalfirst = zero(rate_prototype)
-    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
+    fw = [zero(rate_prototype) for i in 1 : max]
+    ks = [zero(rate_prototype) for i in 1 : max]
+    
     k = ks[1]
 
     J, W1 = build_J_W(alg, u, uprev, p, t, dt, f, uEltypeNoUnits, Val(true))
@@ -642,7 +636,14 @@ 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 : (max - 1) ÷ 2]
-
+    #=
+    linres_tmp = dolinsolve(nothing, linsolve2[1]; A = W2[1], b = _vec(cubuff[1]), linu = _vec(dw2[1]))
+    linres2 = Vector{typeof(linres_tmp)}(undef , (max - 1) ÷ 2)
+    linres2[1] = linres_tmp
+    for i in 2 : (num_stages - 1) ÷ 2
+        linres2[i] = dolinsolve(nothing, linsolve2[1]; A = W2[1], b = _vec(cubuff[i]), linu = _vec(dw2[i]))
+    end
+    =#
     rtol = reltol isa Number ? reltol : zero(reltol)
     atol = reltol isa Number ? reltol : zero(reltol)
 
@@ -652,7 +653,7 @@ function alg_cache(alg::AdaptiveRadau, u, rate_prototype, ::Type{uEltypeNoUnits}
         J, W1, W2,
         uf, tabs, κ, one(uToltype), 10000, tmp,
         atmp, jac_config,
-        linsolve1, linsolve2, rtol, atol, dt, dt,
-        Convergence, alg.step_limiter!, num_stages, 1, big"1.0", big"1.0")
+        linsolve1, linsolve2, #=linres2,=# rtol, atol, dt, dt,
+        Convergence, alg.step_limiter!, num_stages, 1, 0.0)
 end
 
diff --git a/lib/OrdinaryDiffEqFIRK/src/firk_perform_step.jl b/lib/OrdinaryDiffEqFIRK/src/firk_perform_step.jl
index 9b936fd5df..2138a60f0c 100644
--- a/lib/OrdinaryDiffEqFIRK/src/firk_perform_step.jl
+++ b/lib/OrdinaryDiffEqFIRK/src/firk_perform_step.jl
@@ -1357,7 +1357,7 @@ end
     @unpack tabs, num_stages = cache
     tab = tabs[(num_stages - 1) ÷ 2]
     @unpack T, TI, γ, α, β, c, e = tab
-    @unpack κ, cont, θ, θprev = cache
+    @unpack κ, cont = cache
     @unpack internalnorm, abstol, reltol, adaptive = integrator.opts
     alg = unwrap_alg(integrator, true)
     @unpack maxiters = alg
@@ -1401,15 +1401,15 @@ end
         end
     else
         c_prime = Vector{typeof(dt)}(undef, num_stages) #time stepping
-        c_prime[num_stages] = @.. dt / cache.dtprev
+        c_prime[num_stages] = dt / cache.dtprev
         for i in 1 : num_stages - 1
-            c_prime[i] = @.. c[i] * c_prime[num_stages]
+            c_prime[i] = c[i] * c_prime[num_stages]
         end
         for i in 1 : num_stages # collocation polynomial
-            z[i] = @.. cont[num_stages] * (c_prime[i] - c[1] + 1) + cont[num_stages - 1]
+            z[i] = @.. cont[num_stages - 1] + cont[num_stages] * (c_prime[i] - c[1] + 1)
             j = num_stages - 2
             while j > 0
-                z[i] = @.. z[i] * (c_prime[i] - c[num_stages - j] + 1) + cont[j]
+                z[i] = @.. cont[j] + z[i] * (c_prime[i] - c[num_stages - j] + 1)
                 j = j - 1
             end
             z[i] = @.. z[i] * c_prime[i]
@@ -1433,9 +1433,9 @@ end
         integrator.stats.nnonliniter += 1
 
         # evaluate function
-        ff = Vector{typeof(u)}(undef, num_stages)
+        #ff = Vector{typeof(u)}(undef, num_stages)
         for i in 1 : num_stages
-            ff[i] = f(uprev + z[i], p, t + c[i] * dt)
+            z[i] = f(uprev + z[i], p, t + c[i] * dt)
         end
         OrdinaryDiffEqCore.increment_nf!(integrator.stats, 5)
 
@@ -1444,34 +1444,34 @@ end
         for i in 1 : num_stages
             fw[i] = @.. zero(u)
             for j in 1:num_stages
-                fw[i] = @.. fw[i] + TI[i,j] * ff[j]
+                fw[i] = @.. fw[i] + TI[i,j] * z[j]
             end
         end
 
-        Mw = Vector{typeof(u)}(undef, num_stages)
+        #Mw = Vector{typeof(u)}(undef, num_stages)
         if mass_matrix isa UniformScaling # `UniformScaling` doesn't play nicely with broadcast
             for i in 1 : num_stages
-                Mw[i] = @.. mass_matrix.λ * w[i] #scaling by eigenvalue
+                z[i] = @.. mass_matrix.λ * w[i] #scaling by eigenvalue
             end
         else
-            Mw = mass_matrix * w #standard multiplication
+            z = mass_matrix * w #standard multiplication
         end
 
         rhs = Vector{typeof(u)}(undef, num_stages)
-        rhs[1] = @.. fw[1] - γdt * Mw[1]
+        rhs[1] = @.. fw[1] - γdt * z[1]
         i = 2
         while i <= num_stages #block by block multiplication
-            rhs[i] = @.. fw[i] - αdt[i ÷ 2] * Mw[i] + βdt[i ÷ 2] * Mw[i + 1]
-            rhs[i + 1] = @.. fw[i + 1] - βdt[i ÷ 2] * Mw[i] - αdt[i ÷ 2] * Mw[i + 1]
+            rhs[i] = @.. fw[i] - αdt[i ÷ 2] * z[i] + βdt[i ÷ 2] * z[i + 1]
+            rhs[i + 1] = @.. fw[i + 1] - βdt[i ÷ 2] * z[i] - αdt[i ÷ 2] * z[i + 1]
             i += 2
         end
 
-        dw = Vector{typeof(u)}(undef, num_stages)
-        dw[1] = _reshape(LU1 \ _vec(rhs[1]), axes(u))
+        #dw = Vector{typeof(u)}(undef, num_stages)
+        z[1] = _reshape(LU1 \ _vec(rhs[1]), axes(u))
         for i in 2 :(num_stages + 1) ÷ 2
             tmp = _reshape(LU2[i - 1] \ _vec(@.. rhs[2 * i - 2] + rhs[2 * i - 1] * im), axes(u))
-            dw[2 * i - 2] = @.. real(tmp)
-            dw[2 * i - 1] = @.. imag(tmp)
+            z[2 * i - 2] = @.. real(tmp)
+            z[2 * i - 1] = @.. imag(tmp)
         end
         integrator.stats.nsolve +=(num_stages + 1) ÷ 2
 
@@ -1479,15 +1479,13 @@ end
         iter > 1 && (ndwprev = ndw)
         ndw = 0.0
         for i in 1 : num_stages
-            ndw += internalnorm(calculate_residuals(dw[i], uprev, u, atol, rtol, internalnorm, t), t)
+            ndw += internalnorm(calculate_residuals(z[i], uprev, u, atol, rtol, internalnorm, t), t)
         end
 
         # check divergence (not in initial step)
 
         if iter > 1
-            cache.θprev = θ
             θ = ndw / ndwprev
-            cache.θ = θ
             (diverge = θ > 1) && (cache.status = Divergence)
             (veryslowconvergence = ndw * θ^(maxiters - iter) > κ * (1 - θ)) &&
                 (cache.status = VerySlowConvergence)
@@ -1497,7 +1495,7 @@ end
         end
         
         for i in 1 : num_stages
-            w[i] = @.. w[i] - dw[i]
+            w[i] = @.. w[i] - z[i]
         end
 
         # transform `w` to `z`
@@ -1538,10 +1536,9 @@ end
     u = @.. uprev + z[num_stages]
     
     if adaptive
-        edt = e ./ dt
         tmp = 0
         for i in 1 : num_stages
-            tmp = @.. tmp + edt[i] * z[i]
+            tmp = @.. tmp + e[i]/dt * z[i]
         end
         mass_matrix != I && (tmp = mass_matrix * tmp)
         #utilde = @.. broadcast=false 1 / γ * dt * integrator.fsalfirst + tmp
@@ -1604,7 +1601,7 @@ end
     @unpack κ, cont, derivatives, z, w, c_prime = cache
     @unpack dw1, ubuff, dw2, cubuff, dw = cache
     @unpack ks, k, fw, J, W1, W2 = cache
-    @unpack tmp, atmp, jac_config, linsolve1, linsolve2, rtol, atol, step_limiter!, θ, θprev = cache
+    @unpack tmp, atmp, jac_config, linsolve1, linsolve2, #=linres2,=# rtol, atol, step_limiter! = cache
     @unpack internalnorm, abstol, reltol, adaptive = integrator.opts
     alg = unwrap_alg(integrator, true)
     @unpack maxiters = alg
@@ -1638,13 +1635,14 @@ end
             c_prime[i] = c[i] * c_prime[num_stages]
         end
         for i in 1 : num_stages # collocation polynomial
-            @.. z[i] = cont[num_stages] * (c_prime[i] - c[1] + 1) + cont[num_stages - 1]
+            z[i] = cont[num_stages] * (c_prime[i] - c[1] + 1) + cont[num_stages - 1]
             j = num_stages - 2
             while j > 0
-                @.. z[i] *= (c_prime[i] - c[num_stages - j] + 1) + cont[j]
+                @.. z[i] *= (c_prime[i] - c[num_stages - j] + 1) 
+                @.. z[i] += cont[j]
                 j = j - 1
             end
-            @.. z[i] *= c_prime[i]
+            z[i] = z[i] * c_prime[i]
         end
         #mul!(w, TI, z)
         for i in 1:num_stages
@@ -1674,9 +1672,9 @@ end
 
         #mul!(fw, TI, ks)
         for i in 1:num_stages
-            fw[i] = @.. zero(u)
+            fw[i] =  zero(u)
             for j in 1:num_stages
-                fw[i] = @.. fw[i] + TI[i,j] * ks[j]
+                @.. fw[i] += TI[i,j] * ks[j]
             end
         end
 
@@ -1722,10 +1720,8 @@ end
         integrator.stats.nsolve += (num_stages + 1) / 2
 
         for i in 1 : (num_stages - 1) ÷ 2
-            dw[2 * i - 1] = z[2 * i - 1]
-            dw[2 * i] = z[2 * i]
-            dw[2 * i - 1] = real(dw2[i])
-            dw[2 * i] = imag(dw2[i])
+            @.. dw[2 * i - 1] = real(dw2[i])
+            @.. dw[2 * i] = imag(dw2[i])
         end
 
         # compute norm of residuals
@@ -1740,9 +1736,7 @@ end
         # check divergence (not in initial step)
 
         if iter > 1
-            cache.θprev = θ
             θ = ndw / ndwprev
-            cache.θ = θ
             (diverge = θ > 1) && (cache.status = Divergence)
             (veryslowconvergence = ndw * θ^(maxiters - iter) > κ * (1 - θ)) &&
                 (cache.status = VerySlowConvergence)
@@ -1759,15 +1753,15 @@ end
         # transform `w` to `z`
         #mul!(z, T, w)
         for i in 1:num_stages - 1
-            z[i] = @.. zero(u)
+            z[i] = zero(u)
             for j in 1:num_stages
-                z[i] = @.. z[i] + T[i,j] * w[j]
+                @.. z[i] += T[i,j] * w[j]
             end
         end
-        z[num_stages] = @.. T[num_stages, 1] * w[1]
+        z[num_stages] = T[num_stages, 1] * w[1]
         i = 2
         while i < num_stages
-            z[num_stages] = @.. z[num_stages] + w[i]
+            @.. z[num_stages] += w[i]
             i += 2
         end
 
@@ -1796,10 +1790,9 @@ end
     
     if adaptive
         utilde = w[2]
-        edt = e./dt
         @.. tmp = 0
         for i in 1 : num_stages
-            @.. tmp += edt[i] * z[i]
+            @.. tmp += e[i]/dt * z[i]
         end
         mass_matrix != I && (mul!(w[1], mass_matrix, tmp); copyto!(tmp, w[1]))
         #@.. ubuff=1 / γ * dt * integrator.fsalfirst + tmp
@@ -1851,7 +1844,7 @@ end
                 end
             end
             for i in 1 : num_stages
-                cache.cont[i] = derivatives[i, num_stages]
+                @.. cache.cont[i] = derivatives[i, num_stages]
             end
         end
     end
diff --git a/lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl b/lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl
index 49b51d298a..50a06c059f 100644
--- a/lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl
+++ b/lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl
@@ -134,8 +134,8 @@ function BigRadauIIA5Tableau(T1, T2)
     c[3] = big"1"
 
     e = Vector{T1}(undef, 3)
-    e[1] = big"-10.048809399827415944628228317014873027801513671875"
-    e[2] = big"1.3821427331607492039466933420044369995594024658203125"
+    e[1] = big"-10.0488093998274155624603295076470799145872107881988"
+    e[2] = big"1.38214273316074889579366284098041324792054412153223"
     e[3] = big"-0.333333333333333333333333333333333333333333333333333"
 
     TI = Matrix{T1}(undef, 3, 3)

From 7afc20ce49253e70236c22d3a58925a7caf3e28f Mon Sep 17 00:00:00 2001
From: Shreyas Ekanathan <shreyase39@gmail.com>
Date: Sat, 9 Nov 2024 14:33:41 -0500
Subject: [PATCH 09/10] add adaptivity set and misc fixes

---
 lib/OrdinaryDiffEqFIRK/src/firk_caches.jl     |  21 ++--
 .../src/firk_perform_step.jl                  |  15 ++-
 lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl   | 102 +++++++++++-------
 3 files changed, 74 insertions(+), 64 deletions(-)

diff --git a/lib/OrdinaryDiffEqFIRK/src/firk_caches.jl b/lib/OrdinaryDiffEqFIRK/src/firk_caches.jl
index 4f4045a897..8dbcd690c7 100644
--- a/lib/OrdinaryDiffEqFIRK/src/firk_caches.jl
+++ b/lib/OrdinaryDiffEqFIRK/src/firk_caches.jl
@@ -518,7 +518,7 @@ function alg_cache(alg::AdaptiveRadau, u, rate_prototype, ::Type{uEltypeNoUnits}
 end
 
 mutable struct AdaptiveRadauCache{uType, cuType, tType, uNoUnitsType, rateType, JType, W1Type, W2Type,
-    UF, JC, F1, F2, #=F3,=# Tab, Tol, Dt, rTol, aTol, StepLimiter} <:
+    UF, JC, F1, F2, Tab, Tol, Dt, rTol, aTol, StepLimiter} <:
                FIRKMutableCache
     u::uType
     uprev::uType
@@ -550,7 +550,6 @@ mutable struct AdaptiveRadauCache{uType, cuType, tType, uNoUnitsType, rateType,
     jac_config::JC
     linsolve1::F1 #real
     linsolve2::Vector{F2} #complex
-    #linres2::Vector{F3} 
     rtol::rTol
     atol::aTol
     dtprev::Dt
@@ -569,10 +568,8 @@ function alg_cache(alg::AdaptiveRadau, u, rate_prototype, ::Type{uEltypeNoUnits}
     uf = UJacobianWrapper(f, t, p)
     uToltype = constvalue(uBottomEltypeNoUnits)
 
-    min = alg.min_stages
     max = alg.max_stages
-
-    num_stages = min
+    num_stages = alg.min_stages
 
     tabs = [BigRadauIIA5Tableau(uToltype, constvalue(tTypeNoUnits)), BigRadauIIA9Tableau(uToltype, constvalue(tTypeNoUnits)), BigRadauIIA13Tableau(uToltype, constvalue(tTypeNoUnits))]
     i = 9
@@ -590,6 +587,9 @@ function alg_cache(alg::AdaptiveRadau, u, rate_prototype, ::Type{uEltypeNoUnits}
     end
 
     c_prime = Vector{typeof(t)}(undef, max) #time stepping
+    for i in 1 : max
+        c_prime[i] = zero(t)
+    end
 
     dw1 = zero(u)
     ubuff = zero(u)
@@ -636,14 +636,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 : (max - 1) ÷ 2]
-    #=
-    linres_tmp = dolinsolve(nothing, linsolve2[1]; A = W2[1], b = _vec(cubuff[1]), linu = _vec(dw2[1]))
-    linres2 = Vector{typeof(linres_tmp)}(undef , (max - 1) ÷ 2)
-    linres2[1] = linres_tmp
-    for i in 2 : (num_stages - 1) ÷ 2
-        linres2[i] = dolinsolve(nothing, linsolve2[1]; A = W2[1], b = _vec(cubuff[i]), linu = _vec(dw2[i]))
-    end
-    =#
+
     rtol = reltol isa Number ? reltol : zero(reltol)
     atol = reltol isa Number ? reltol : zero(reltol)
 
@@ -653,7 +646,7 @@ function alg_cache(alg::AdaptiveRadau, u, rate_prototype, ::Type{uEltypeNoUnits}
         J, W1, W2,
         uf, tabs, κ, one(uToltype), 10000, tmp,
         atmp, jac_config,
-        linsolve1, linsolve2, #=linres2,=# rtol, atol, dt, dt,
+        linsolve1, linsolve2, rtol, atol, dt, dt,
         Convergence, alg.step_limiter!, num_stages, 1, 0.0)
 end
 
diff --git a/lib/OrdinaryDiffEqFIRK/src/firk_perform_step.jl b/lib/OrdinaryDiffEqFIRK/src/firk_perform_step.jl
index 2138a60f0c..5e42aa5b94 100644
--- a/lib/OrdinaryDiffEqFIRK/src/firk_perform_step.jl
+++ b/lib/OrdinaryDiffEqFIRK/src/firk_perform_step.jl
@@ -1601,7 +1601,7 @@ end
     @unpack κ, cont, derivatives, z, w, c_prime = cache
     @unpack dw1, ubuff, dw2, cubuff, dw = cache
     @unpack ks, k, fw, J, W1, W2 = cache
-    @unpack tmp, atmp, jac_config, linsolve1, linsolve2, #=linres2,=# rtol, atol, step_limiter! = cache
+    @unpack tmp, atmp, jac_config, linsolve1, linsolve2, rtol, atol, step_limiter! = cache
     @unpack internalnorm, abstol, reltol, adaptive = integrator.opts
     alg = unwrap_alg(integrator, true)
     @unpack maxiters = alg
@@ -1635,14 +1635,14 @@ end
             c_prime[i] = c[i] * c_prime[num_stages]
         end
         for i in 1 : num_stages # collocation polynomial
-            z[i] = cont[num_stages] * (c_prime[i] - c[1] + 1) + cont[num_stages - 1]
+            @.. z[i] = cont[num_stages] * (c_prime[i] - c[1] + 1) + cont[num_stages - 1]
             j = num_stages - 2
             while j > 0
                 @.. z[i] *= (c_prime[i] - c[num_stages - j] + 1) 
                 @.. z[i] += cont[j]
                 j = j - 1
             end
-            z[i] = z[i] * c_prime[i]
+            @.. z[i] *= c_prime[i]
         end
         #mul!(w, TI, z)
         for i in 1:num_stages
@@ -1672,7 +1672,7 @@ end
 
         #mul!(fw, TI, ks)
         for i in 1:num_stages
-            fw[i] =  zero(u)
+            @.. fw[i] = zero(u)
             for j in 1:num_stages
                 @.. fw[i] += TI[i,j] * ks[j]
             end
@@ -1704,17 +1704,14 @@ end
 
         cache.linsolve1 = linres.cache
 
-        linres2 = Vector{Any}(undef,(num_stages - 1) ÷ 2)
-
         for i in 1 :(num_stages - 1) ÷ 2
             @.. cubuff[i]=complex(
             fw[2 * i] - αdt[i] * Mw[2 * i] + βdt[i] * Mw[2 * i + 1], fw[2 * i + 1] - βdt[i] * Mw[2 * i] - αdt[i] * Mw[2 * i + 1])
             if needfactor
-                linres2[i] = dolinsolve(integrator, linsolve2[i]; A = W2[i], b = _vec(cubuff[i]), linu = _vec(dw2[i]))
+                cache.linsolve2[i] = dolinsolve(integrator, linsolve2[i]; A = W2[i], b = _vec(cubuff[i]), linu = _vec(dw2[i])).cache
             else
-                linres2[i] = dolinsolve(integrator, linsolve2[i]; A = nothing, b = _vec(cubuff[i]), linu = _vec(dw2[i]))
+                cache.linsolve2[i] = dolinsolve(integrator, linsolve2[i]; A = nothing, b = _vec(cubuff[i]), linu = _vec(dw2[i])).cache
             end
-            cache.linsolve2[i] = linres2[i].cache
         end
 
         integrator.stats.nsolve += (num_stages + 1) / 2
diff --git a/lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl b/lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl
index 50a06c059f..a0e05f28fd 100644
--- a/lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl
+++ b/lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl
@@ -134,9 +134,9 @@ function BigRadauIIA5Tableau(T1, T2)
     c[3] = big"1"
 
     e = Vector{T1}(undef, 3)
-    e[1] = big"-10.0488093998274155624603295076470799145872107881988"
-    e[2] = big"1.38214273316074889579366284098041324792054412153223"
-    e[3] = big"-0.333333333333333333333333333333333333333333333333333"
+    e[1] = big"-10.0488093998274155624603295076470799145872107881988969663429493235855742140670683952596720105774938812433874028620997746246706860729547671304601625528869782"
+    e[2] = big"1.38214273316074889579366284098041324792054412153223029967628265691890754740040172859300534391082721457672073619543310795800401940628810046379349588622031217"
+    e[3] = big"-0.333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333644"
 
     TI = Matrix{T1}(undef, 3, 3)
     TI[1, 1] = big"4.32557989006315535102435095295614882731995158490590784287320458848019483341979047442263696495019938973156007686663488090615420049217658854859024016717169837"
@@ -180,11 +180,11 @@ function BigRadauIIA9Tableau(T1, T2)
     c[5] = big"1.0"
 
     e = Vector{T1}(undef, 5)
-    e[1] = big"-27.78093394406463730479"
-    e[2] = big"3.641478498049213152712"
-    e[3] = big"-1.252547721169118720491"
-    e[4] = big"0.5920031671845428725662"
-    e[5] = big"-0.2000000000000000000000"
+    e[1] = big"-27.7809339440646373047872078172168798923674228687740760060378492475924178050505976287227228556471699142365371740120443650701118024100678675823465762727483305"
+    e[2] = big"3.64147849804921315271165508774289722904088750334220956841022786858917594981395319605788667956024462601802006251583142928630101075351336314632135787805261686"
+    e[3] = big"-1.25254772116911872049065249430114914889315244289570569309128740586057170336299694248256681515155624683225624015343224399700466177251702555220815764199263189"
+    e[4] = big"0.592003167184542872566205223775131812219687808327572130718908784863813558599641375147402991238481535050773351649645179780815453429071529988233376036688329872"
+    e[5] = big"-0.199999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999997076"
 
     TI = Matrix{T1}(undef, 5, 5)
     TI[1, 1] = big"30.0415677215444016277146611632467970747634862837368422955138470463852339244593400023985957753164599415374157317627305099177616927640413043608408838747985125"
@@ -414,13 +414,13 @@ function BigRadauIIA13Tableau(T1, T2)
     c[7] = big"1.0"
 
     e = Vector{T1}(undef, 7)
-    e[1] = big"-54.37443689412861451458"
-    e[2] = big"7.000024004259186512041"
-    e[3] = big"-2.355661091987557192256"
-    e[4] = big"1.132289066106134386384"
-    e[5] = big"-0.6468913267673587118673"
-    e[6] = big"0.3875333853753523774248"
-    e[7] = big"-0.1428571428571428571429"
+    e[1] = big"-54.374436894128614514583710369683221528326818668136315170227649609831132483812209590903458627819914413600703287942266678601263304348350182019714004102122958"
+    e[2] = big"7.00002400425918651204068363735192307633403862621907697222411411256593188888314837387690262103761082115674234000933589934965063951414231971808906314491204573"
+    e[3] = big"-2.35566109198755719225604586775720723211163199654640573606711168106849118084357027539414093812951288166804790294091903523762277368547775099880612390898224076"
+    e[4] = big"1.13228906610613438638449290827978318662460499026070073842612187085281352278780837966549916347601259689966925986653914736463076068138934273474363230390185871"
+    e[5] = big"-0.646891326767358711867345222439989069591870662562921671446738173180691199552327090727940249497816198076028398716990245669520129053944261569921119452534594627"
+    e[6] = big"0.387533385375352377424782057105854424214534853623007724234120623518712309680007346340280888076477218145510846867158055651267664035097674992751409157682864641"
+    e[7] = big"-0.142857142857142857142857142857142857142857142857142857142857142857142857142857142857142857142857142857142857142857142857142857142857142857142857142857092806"
 
     TI = Matrix{T1}(undef, 7, 7)
     TI[1, 1] = big"258.131926319982229276108947425184471333411128774462923076434633414645220927977539758484670571338176678808837829326061674950321562391576244286310404028770676"
@@ -599,34 +599,54 @@ function adaptiveRadauTableau(T1, T2, num_stages::Int)
 
     if (num_stages == 9)
         e = Vector{BigFloat}(undef, 9)
-        e[1] = big"-0.133101731359431287515066981129913748644705107621439651956220186897253838380345034218235538734967567153163030284540660584311040323114847240173627907922903296"
-        e[2] = big"0.0754476228408557299650196603226967248368445025181771896522057250989188754588885465998346476425502117889420021664297319179240040109156780754680742172762707621"
-        e[3] = big"-0.0458369394236156144604575482137179697005739995740615341890112217655441769701945378217626766299683076189687755618065050383493055018324395934911567207485032988"
-        e[4] = big"0.0271430329153098694457979735602502142083095152399102869109830450899844979409229538982100527256348792152825816553434603418662939944133319974874915933773657075"
-        e[5] = big"-0.0156126300301219212217568535995825232086423550686814635293876744035364259647929167763641353639085929285192729729570945658304937255929114458885296622493040224"
-        e[6] = big"0.00890598154557403928205152521539967562877335780940124672915181111908317890891659158654221736499522823959933517986673010006749138291836676520080172845444352328"
-        e[7] = big"-0.00514824122639241252178399021479378841872099572255461304439292434131750195489022869965968028106854978547414579491205935930595041763060069987112580994637398395"
-        e[8] = big"0.00296533914055503317169967748114188676589522458557982039693426239853498956125735811263087631479968309978854200615027412311940897061471388689986239742919640848"
-        e[9] = big"-0.0010634368308888065260482548541946175520274736959410047497431569257848032902381738362547705844630238841535652230832162703806430112125115777122361837311714267"
+        e[1] = big"-89.8315397040376845865027298766511166861131537901479318008187013574099993398844876573472315778350373191126204142357525815115482293843777624541394691345885716"
+        e[2] = big"11.4742766094687721590222610299234578063148408248968597722844661019124491691448775794163842022854672278004372474682761156236829237591471118886342174262239472"
+        e[3] = big"-3.81419058476042873698615187248837320040477891376179026064712181641592908409919668221598902628694008903410444392769866137859041139561191341971835412426311966"
+        e[4] = big"1.81155300867853110911564243387531599775142729190474576183505286509346678884073482369609308584446518479366940471952219053256362416491879701351428578466580598"
+        e[5] = big"-1.03663781378817415276482837566889343026914084945266083480559060702535168750966084568642219911350874500410428043808038021858812311835772945467924877281164517"
+        e[6] = big"0.660865688193716483757690045578935452512421753840843511309717716369201467579470723336314286637650332622546110594223451602017981477424498704954672224534648119"
+        e[7] = big"-0.444189256280526730087023435911479370800996444567516110958885112499737452734669537494435549195615660656770091500773942469075264796140815048410568498349675229"
+        e[8] = big"0.290973163636905565556251162453264542120491238398561072912173321087011249774042707406397888774630179702057578431394918930648610404108923880955576205699885598"
+        e[9] = big"-0.111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111222795"
+    elseif (num_stages == 11)
+        e = Vector{BigFloat}(undef, 11)
+        e[1] = big"-134.152626015465044063378550835075318643291579891352838474367124350171545245813244797505763447327562609902792066283575334085390478517120485782603677022267543"
+        e[2] = big"17.0660253399060146849212356299749772423073416838121578997449942694355150369717420038613850964748566731121793290881077515821557030349184664685171028112845693"
+        e[3] = big"-5.63464089555106294823267450977601185069165875295372865523759287935369597689662768988715406731927279137711764532851201746616033935275093116699140897901326857"
+        e[4] = big"2.65398285960564394428637524662555134392389271086844331137910389226095922845489762567700560496915255196379049844894623384211693438658842276927416827629120392"
+        e[5] = big"-1.50753272514563441873424939425410006034401178578882643601844794171149654717227697249290904230103304153661631200445957060050895700394738491883951084826421405"
+        e[6] = big"0.960260572218344245935269463733859188992760928707230734981795807797858324380878500135029848170473080912207529262984056182004711806457345405466997261506487216"
+        e[7] = big"-0.658533932484491373507110339620843007350146695468297825313721271556868110859353953892288534787571420691760379406525738632649863532050280264983313133523641674"
+        e[8] = big"0.47189364490739958527881800092758816959227958959727295348380187162217987951960275929676019062173412149363239153353720640122975284789262792027244826613784432"
+        e[9] = big"-0.34181016557091711933253384050957887606039737751222218385118573305954222606860932803075900338195356026497059819558648780544900376040113065955083806288937526"
+        e[10] = big"0.233890408488838371854329668882967402012428680999899584289285425645726546573900943747784263972086087200538161975992991491742449181322441138528940521648041699"
+        e[11] = big"-0.0909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909093788951"
+    elseif (num_stages == 13)
+        e = Vector{BigFloat}(undef, 13)
+        e[1] = big"-187.337806666035250696387113105488477375830948862159770885826492736743460038872636916422100706359786154665214547894636085276885830138994748219148357620227002"
+        e[2] = big"23.775705048946302520021716862887025159493544949407763131913924588605891085865877529749667170060976683489861224477421212170329019074926368036881685518012728"
+        e[3] = big"-7.81823724708755833325842676798052630403951326380926053607036280237871312516353176794790424805918285990907426633641064901501063343970205708057561515795364672"
+        e[4] = big"3.66289388251066047904501665386587373682645522696191680651425553890800106379174431775463608296821504040006089759980653462003322200870566661322334735061646223"
+        e[5] = big"-2.06847094952801462392548700163367193433237251061765813625197254100990426184032443671875204952150187523419743001493620194301209589692419776688692360679336566"
+        e[6] = big"1.31105635982993157063104433803023633257356281733787535204132865785504258558244947718491624714070193102812968996631302993877989767202703509685785407541965509"
+        e[7] = big"-0.897988270828178667954874573865888835427640297795141000639881363403080887358272161865529150995401606679722232843051402663087372891040498351714982629218397165"
+        e[8] = big"0.648958340079591709325028357505725843500310779765000237611355105578356380892509437805732950287939403489669590070670546599339082534053791877148407548785389408"
+        e[9] = big"-0.485906120880156534303797908584178831869407602334908394589833216071089678420073112977712585616439120156658051446412515753614726507868506301824972455936531663"
+        e[10] = big"0.370151313405058266144090771980402238126294149688261261935258556082315591034906662511634673912342573394958760869036835172495369190026354174118335052418701339"
+        e[11] = big"-0.27934271062931554435643589252670994638477019847143394253283050767117135003630906657393675748475838251860910095199485920686192935009874559019443503474805827"
+        e[12] = big"0.195910097140006778096161342733266840441407888950433028972173797170889557600583114422425296743817444283872389581116632280572920821812614435192580036549169031"
+        e[13] = big"-0.0769230769230769230769230769230769230769230769230769230769230769230769230769230769230769230769230769230769230769230769230769230769230769230769230769254590189"
     else
-        p = num_stages
-        eb = variables(:b, 1:num_stages + 1)
-        @variables y
-        zz = zeros(size(a, 1) + 1)
-        zz2 = zeros(size(a, 1))
-        eA = [zz'
-            zz2 a]
-        ec = [0; c]
-        constraints = map(Iterators.flatten(RootedTreeIterator(i) for i in 1:2*p-3)) do t
-            residual_order_condition(t, RungeKuttaMethod(eA, eb, ec))
+        e_sym = variables(:e, 1:num_stages)
+        constraints = map(Iterators.flatten(RootedTreeIterator(i) for i in 1:num_stages)) do t
+            residual_order_condition(t, RungeKuttaMethod(a, e_sym, c))
         end 
-        AA, bb, islinear = Symbolics.linear_expansion(Symbolics.substitute.(constraints, (eb[1]=>1/γ,)), eb[2:end])
-        AA = Float64.(map(unwrap, AA))
-        idxs = qr(AA', ColumnNorm()).p[1:num_stages]
-        @assert rank(AA[idxs, :]) == num_stages
-        @assert islinear
-        b_hat = Symbolics.expand.((AA \ -bb))
-        e = [Symbolics.symbolic_to_float(b_hat[i] - b[i]) for i in 1 : num_stages]    
+        AA, bb, islinear = Symbolics.linear_expansion(constraints, e_sym[1:end])
+        AA = BigFloat.(map(unwrap, AA))
+        bb = BigFloat.(map(unwrap, bb))
+        A = vcat([zeros(num_stages -1); 1]', AA)
+        b_2 = vcat(-1/big(num_stages), -(num_stages)^2, -1, zeros(size(A, 1) - 3))
+        e = A \ b_2
     end
     RadauIIATableau{T1, T2}(T, TI, c, γ, α, β, e)
 end

From 3f07e1d2a5b1d13762bd2975ef0dc36be5c53fc3 Mon Sep 17 00:00:00 2001
From: Shreyas Ekanathan <shreyase39@gmail.com>
Date: Sat, 9 Nov 2024 18:04:21 -0500
Subject: [PATCH 10/10] Update firk_tableaus.jl

---
 lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl b/lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl
index a0e05f28fd..f0684fb259 100644
--- a/lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl
+++ b/lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl
@@ -136,7 +136,7 @@ function BigRadauIIA5Tableau(T1, T2)
     e = Vector{T1}(undef, 3)
     e[1] = big"-10.0488093998274155624603295076470799145872107881988969663429493235855742140670683952596720105774938812433874028620997746246706860729547671304601625528869782"
     e[2] = big"1.38214273316074889579366284098041324792054412153223029967628265691890754740040172859300534391082721457672073619543310795800401940628810046379349588622031217"
-    e[3] = big"-0.333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333644"
+    e[3] = big(-1)/3
 
     TI = Matrix{T1}(undef, 3, 3)
     TI[1, 1] = big"4.32557989006315535102435095295614882731995158490590784287320458848019483341979047442263696495019938973156007686663488090615420049217658854859024016717169837"
@@ -184,7 +184,7 @@ function BigRadauIIA9Tableau(T1, T2)
     e[2] = big"3.64147849804921315271165508774289722904088750334220956841022786858917594981395319605788667956024462601802006251583142928630101075351336314632135787805261686"
     e[3] = big"-1.25254772116911872049065249430114914889315244289570569309128740586057170336299694248256681515155624683225624015343224399700466177251702555220815764199263189"
     e[4] = big"0.592003167184542872566205223775131812219687808327572130718908784863813558599641375147402991238481535050773351649645179780815453429071529988233376036688329872"
-    e[5] = big"-0.199999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999997076"
+    e[5] = big(-1)/5
 
     TI = Matrix{T1}(undef, 5, 5)
     TI[1, 1] = big"30.0415677215444016277146611632467970747634862837368422955138470463852339244593400023985957753164599415374157317627305099177616927640413043608408838747985125"
@@ -420,7 +420,7 @@ function BigRadauIIA13Tableau(T1, T2)
     e[4] = big"1.13228906610613438638449290827978318662460499026070073842612187085281352278780837966549916347601259689966925986653914736463076068138934273474363230390185871"
     e[5] = big"-0.646891326767358711867345222439989069591870662562921671446738173180691199552327090727940249497816198076028398716990245669520129053944261569921119452534594627"
     e[6] = big"0.387533385375352377424782057105854424214534853623007724234120623518712309680007346340280888076477218145510846867158055651267664035097674992751409157682864641"
-    e[7] = big"-0.142857142857142857142857142857142857142857142857142857142857142857142857142857142857142857142857142857142857142857142857142857142857142857142857142857092806"
+    e[7] = big(-1)/7
 
     TI = Matrix{T1}(undef, 7, 7)
     TI[1, 1] = big"258.131926319982229276108947425184471333411128774462923076434633414645220927977539758484670571338176678808837829326061674950321562391576244286310404028770676"