From d8c42744e3ee30590470abad988d425cc8dc75bd Mon Sep 17 00:00:00 2001 From: singhharsh1708 Date: Wed, 28 Oct 2026 05:24:41 +0530 Subject: [PATCH] Reuse in-place FIRK W only at the dt it was factored for --- lib/OrdinaryDiffEqFIRK/Project.toml | 2 +- lib/OrdinaryDiffEqFIRK/src/algorithms.jl | 8 ++++---- lib/OrdinaryDiffEqFIRK/src/firk_perform_step.jl | 4 ++++ lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl | 13 +++++++++++++ 4 files changed, 22 insertions(+), 5 deletions(-) diff --git a/lib/OrdinaryDiffEqFIRK/Project.toml b/lib/OrdinaryDiffEqFIRK/Project.toml index 01e0399e54..dd5beeb6ad 100644 --- a/lib/OrdinaryDiffEqFIRK/Project.toml +++ b/lib/OrdinaryDiffEqFIRK/Project.toml @@ -1,7 +1,7 @@ name = "OrdinaryDiffEqFIRK" uuid = "5960d6e9-dd7a-4743-88e7-cf307b64f125" authors = ["ParamThakkar123 "] -version = "2.8.9" +version = "2.8.10" [deps] FastBroadcast = "7034ab61-46d4-4ed7-9d0f-46aef9175898" diff --git a/lib/OrdinaryDiffEqFIRK/src/algorithms.jl b/lib/OrdinaryDiffEqFIRK/src/algorithms.jl index 5dd5d9afd6..2b604102c3 100644 --- a/lib/OrdinaryDiffEqFIRK/src/algorithms.jl +++ b/lib/OrdinaryDiffEqFIRK/src/algorithms.jl @@ -52,7 +52,7 @@ function RadauIIA3(; concrete_jac = nothing, linsolve = nothing, extrapolant = :dense, fast_convergence_cutoff = 1 // 5, - new_W_γdt_cutoff = 1 // 5, + new_W_γdt_cutoff = 0 // 1, κ = nothing, maxiters = 10, step_limiter! = trivial_limiter! ) @@ -99,7 +99,7 @@ function RadauIIA5(; concrete_jac = nothing, linsolve = nothing, extrapolant = :dense, fast_convergence_cutoff = 1 // 5, - new_W_γdt_cutoff = 1 // 5, + new_W_γdt_cutoff = 0 // 1, κ = nothing, maxiters = 10, smooth_est = true, step_limiter! = trivial_limiter! ) @@ -148,7 +148,7 @@ function RadauIIA9(; concrete_jac = nothing, linsolve = nothing, extrapolant = :dense, fast_convergence_cutoff = 1 // 5, - new_W_γdt_cutoff = 1 // 5, + new_W_γdt_cutoff = 0 // 1, κ = nothing, maxiters = 10, smooth_est = true, step_limiter! = trivial_limiter! ) @@ -200,7 +200,7 @@ function AdaptiveRadau(; min_order = 5, max_order = 13, threading = false, linsolve = nothing, extrapolant = :dense, fast_convergence_cutoff = 1 // 5, - new_W_γdt_cutoff = 1 // 5, + new_W_γdt_cutoff = 0 // 1, κ = nothing, maxiters = 10, smooth_est = true, step_limiter! = trivial_limiter! ) diff --git a/lib/OrdinaryDiffEqFIRK/src/firk_perform_step.jl b/lib/OrdinaryDiffEqFIRK/src/firk_perform_step.jl index 3ccc3865f5..16e42e03be 100644 --- a/lib/OrdinaryDiffEqFIRK/src/firk_perform_step.jl +++ b/lib/OrdinaryDiffEqFIRK/src/firk_perform_step.jl @@ -447,6 +447,7 @@ end (new_jac = do_newJ(integrator, alg, cache, repeat_step)) && (firk_new_J!(J, W1, integrator, cache); cache.W_γdt = dt) new_W = do_newW(integrator, alg, new_jac, cache.W_γdt) + new_W && (cache.W_γdt = dt) if is_lazy_W(W1) set_W_gamma!(W1, αdt + βdt * im) integrator.stats.nw += 1 @@ -825,6 +826,7 @@ end (new_jac = do_newJ(integrator, alg, cache, repeat_step)) && (firk_new_J!(J, W1, integrator, cache, W2); cache.W_γdt = dt) new_W = do_newW(integrator, alg, new_jac, cache.W_γdt) + new_W && (cache.W_γdt = dt) if is_lazy_W(W1) set_W_gamma!(W1, γdt) set_W_gamma!(W2, αdt + βdt * im) @@ -1403,6 +1405,7 @@ end (new_jac = do_newJ(integrator, alg, cache, repeat_step)) && (firk_new_J!(J, W1, integrator, cache, W2, W3); cache.W_γdt = dt) new_W = do_newW(integrator, alg, new_jac, cache.W_γdt) + new_W && (cache.W_γdt = dt) if is_lazy_W(W1) set_W_gamma!(W1, γdt) set_W_gamma!(W2, α1dt + β1dt * im) @@ -2052,6 +2055,7 @@ end (new_jac = do_newJ(integrator, alg, cache, repeat_step)) && (firk_new_J!(J, W1, integrator, cache, W2); cache.W_γdt = dt) new_W = do_newW(integrator, alg, new_jac, cache.W_γdt) + new_W && (cache.W_γdt = dt) if is_lazy_W(W1) set_W_gamma!(W1, γdt) for i in 1:((num_stages - 1) ÷ 2) diff --git a/lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl b/lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl index 0c186fc70b..1a4169e585 100644 --- a/lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl +++ b/lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl @@ -339,3 +339,16 @@ end ) @test @inferred(OrdinaryDiffEqCore.get_current_adaptive_order(integ.alg.algs[2], integ.cache)) isa Int end + +@testset "In-place FIRK stays accurate after a callback" begin + f!(du, u, p, t) = (du[1] = -3u[1] + 1; nothing) + exact(t) = (1 - exp(-3t)) / 3 + callback = DiscreteCallback((u, t, integrator) -> t == 0.5, integrator -> nothing) + for alg in (RadauIIA5(), RadauIIA9(), AdaptiveRadau()) + sol = solve( + ODEProblem(f!, [0.0], (0.0, 1.0)), alg; reltol = 1.0e-8, abstol = 1.0e-11, + callback, tstops = [0.5] + ) + @test maximum(abs(u[1] - exact(t)) for (t, u) in zip(sol.t, sol.u)) < 1.0e-8 + end +end