tweaks

Author: Shreyas-Ekanathan

Project.toml

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

lib/OrdinaryDiffEqFIRK/src/algorithms.jl

@@ -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,

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

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,29 +507,15 @@ 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)
     end
 
     κ = 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
 

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"

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)