From 2d6b5f6156a46c2368da9620bf7a781a0f2b417b Mon Sep 17 00:00:00 2001 From: Shreyas Ekanathan Date: Sun, 25 Aug 2024 19:16:50 -0400 Subject: [PATCH] fixes --- lib/OrdinaryDiffEqFIRK/src/algorithms.jl | 2 +- lib/OrdinaryDiffEqFIRK/src/firk_caches.jl | 28 ++++++++----------- .../src/firk_perform_step.jl | 1 + lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl | 12 ++++++-- 4 files changed, 22 insertions(+), 21 deletions(-) diff --git a/lib/OrdinaryDiffEqFIRK/src/algorithms.jl b/lib/OrdinaryDiffEqFIRK/src/algorithms.jl index ba4c51c0d0..4429fb78b6 100644 --- 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 index 39e70e5831..5d7ccba8c4 100644 --- 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 index 87d4de366e..81740fbf0b 100644 --- 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 index bbea52cbf2..bead50a768 100644 --- a/lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl +++ b/lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl @@ -9,9 +9,9 @@ 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 @@ -19,6 +19,12 @@ sim21 = test_convergence(1 ./ 2 .^ (2.5:-1:0.5), prob_ode_linear, RadauIIA9()) 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