fixes

Shreyas-Ekanathan · Aug 25, 2024 · 2d6b5f6 · 2d6b5f6
1 parent c69f1c1
commit 2d6b5f6
Show file tree

Hide file tree

Showing 4 changed files with 22 additions and 21 deletions.
diff --git a/lib/OrdinaryDiffEqFIRK/src/algorithms.jl b/lib/OrdinaryDiffEqFIRK/src/algorithms.jl
@@ -171,7 +171,7 @@ end
 
 function AdaptiveRadau(; chunk_size = Val{0}(), autodiff = Val{true}(),
      standardtag = Val{true}(), concrete_jac = nothing,
-     diff_type = Val{:forward}, num_stages = 5, 
+     diff_type = Val{:forward}, num_stages = 3, 
      linsolve = nothing, precs = DEFAULT_PRECS,
      extrapolant = :dense, fast_convergence_cutoff = 1 // 5,
      new_W_γdt_cutoff = 1 // 5,

diff --git a/lib/OrdinaryDiffEqFIRK/src/firk_caches.jl b/lib/OrdinaryDiffEqFIRK/src/firk_caches.jl
@@ -511,7 +511,7 @@ end
 
 mutable struct AdaptiveRadauCache{uType, cuType, uNoUnitsType, rateType, JType, W1Type, W2Type,
     UF, JC, F1, F2, Tab, Tol, Dt, rTol, aTol, StepLimiter} <:
-               OrdinaryDiffEqMutableCache
+               FIRKMutableCache
     u::uType
     uprev::uType
     z::AbstractVector{uType}
@@ -533,7 +533,6 @@ mutable struct AdaptiveRadauCache{uType, cuType, uNoUnitsType, rateType, JType,
     κ::Tol
     ηold::Tol
     iter::Int
-    tmp::AbstractVector{uType}
     atmp::uNoUnitsType
     jac_config::JC
     linsolve1::F1 #real
@@ -565,12 +564,12 @@ function alg_cache(alg::AdaptiveRadau, u, rate_prototype, ::Type{uEltypeNoUnits}
 
     dw1 = zero(u)
     ubuff = zero(u)
-    dw2 = Vector{typeof(u)}(undef, floor(Int, num_stages / 2))
+    dw2 = Vector{Any}(undef, floor(Int, num_stages / 2))
     for i in 1 : floor(Int, num_stages / 2)
         dw2[i] = similar(u, Complex{eltype(u)})
-        recursivefill!(dw[i], false)
+        recursivefill!(dw2[i], false)
     end
-    cubuff = Vector{typeof(u)}(undef, floor(Int, num_stages / 2))
+    cubuff = Vector{Any}(undef, floor(Int, num_stages / 2))
     for i in 1 : floor(Int, num_stages / 2)
         cubuff[i] = similar(u, Complex{eltype(u)})
         recursivefill!(cubuff[i], false)
@@ -592,32 +591,27 @@ function alg_cache(alg::AdaptiveRadau, u, rate_prototype, ::Type{uEltypeNoUnits}
     if J isa AbstractSciMLOperator
         error("Non-concrete Jacobian not yet supported by RadauIIA5.")
     end
-    W2 = vector{typeof(Complex{W1})}(undef, floor(Int, num_stages/2))
+    W2 = Vector{Any}(undef, floor(Int, num_stages/2))
     for i in 1 : floor(Int, num_stages / 2)
         W2[i] = similar(J, Complex{eltype(W1)})
-        recursivefill!(w2[i], false)
+        recursivefill!(W2[i], false)
     end
 
     du1 = zero(rate_prototype)
 
-    tmp = Vector{typeof(u)}(undef, binomial(num_stages , 2))
-    for i in 1 : binomial(num_stages , 2)
-        tmp[i] = zero(u)
-    end
-
     atmp = similar(u, uEltypeNoUnits)
     recursivefill!(atmp, false)
 
-    jac_config = build_jac_config(alg, f, uf, du1, uprev, u, tmp, dw1)
+    jac_config = build_jac_config(alg, f, uf, du1, uprev, u, zero(u), dw1)
 
     linprob = LinearProblem(W1, _vec(ubuff); u0 = _vec(dw1))
     linsolve1 = init(linprob, alg.linsolve, alias_A = true, alias_b = true,
         assumptions = LinearSolve.OperatorAssumptions(true))
 
-    linsolve2 = Vector{typeof(linsolve1)}(undef, floor(Int, num_stages / 2))
-    for i in 1 : floor(int, num_stages / 2)
+    linsolve2 = Vector{Any}(undef, floor(Int, num_stages / 2))
+    for i in 1 : floor(Int, num_stages / 2)
         linprob = LinearProblem(W2[i], _vec(cubuff[i]); u0 = _vec(dw2[i]))
-        linsolve2 = init(linprob, alg.linsolve, alias_A = true, alias_b = true,
+        linsolve2[i] = init(linprob, alg.linsolve, alias_A = true, alias_b = true,
         assumptions = LinearSolve.OperatorAssumptions(true))
     end
 
@@ -629,7 +623,7 @@ function alg_cache(alg::AdaptiveRadau, u, rate_prototype, ::Type{uEltypeNoUnits}
         du1, fsalfirst, k, fw,
         J, W1, W2,
         uf, tab, κ, one(uToltype), 10000,
-        tmp, atmp, jac_config,
+        atmp, jac_config,
         linsolve1, linsolve2, rtol, atol, dt, dt,
         Convergence, alg.step_limiter!)
 end

diff --git a/lib/OrdinaryDiffEqFIRK/src/firk_perform_step.jl b/lib/OrdinaryDiffEqFIRK/src/firk_perform_step.jl
@@ -90,6 +90,7 @@ function initialize!(integrator, cache::RadauIIA9Cache)
 end
 
 function initialize!(integrator, cache::AdaptiveRadauCache)
+    println("here")
     integrator.kshortsize = 2
     resize!(integrator.k, integrator.kshortsize)
     integrator.k[1] = integrator.fsalfirst

diff --git a/lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl b/lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl
@@ -9,16 +9,22 @@ for prob in [prob_ode_linear, prob_ode_2Dlinear]
 end
 
 
-sol = solve(prob_ode_linear, AdaptiveRadau(), adaptive = false, dt = 1e-2)
-sol = solve(prob_ode_linear, RadauIIA9(), adaptive = false, dt = 1e-2)
-sol = solve(prob_ode_linear, RadauIIA5(), adaptive = false, dt = 1e-2)
+sol = solve(prob_ode_2Dlinear, AdaptiveRadau(), adaptive = false, dt = 1e-2)
+sol = solve(prob_ode_2Dlinear, RadauIIA9(), adaptive = false, dt = 1e-2)
+sol = solve(prob_ode_2Dlinear, RadauIIA5(), adaptive = false, dt = 1e-2)
 
 sim21 = test_convergence(1 ./ 2 .^ (2.5:-1:0.5), prob_ode_linear, RadauIIA9())
 @test sim21.𝒪est[:final]≈9 atol=testTol
 
 sim21 = test_convergence(1 ./ 2 .^ (2.5:-1:0.5), prob_ode_2Dlinear, RadauIIA9())
 @test sim21.𝒪est[:final]≈9 atol=testTol
 
+sim21 = test_convergence(1 ./ 2 .^ (2.25:-1:0.25), prob_ode_linear, AdaptiveRadau())
+@test sim21.𝒪est[:final]≈9 atol=testTol
+
+sim21 = test_convergence(1 ./ 2 .^(2.25:-1:0.25), prod_ode_2Dlinear, AdaptiveRadau())
+@test sim21.𝒪est[:final]≈9 atol=testTol
+
 # test adaptivity
 for iip in (true, false)
     if iip