From 483d62a88af0fb761072ba57856957e5ea1c3f73 Mon Sep 17 00:00:00 2001 From: Qingyu Qu <2283984853@qq.com> Date: Sun, 11 May 2025 15:35:26 +0800 Subject: [PATCH 1/2] Fix incorrect interpolant evaluation --- lib/BoundaryValueDiffEqMIRK/src/adaptivity.jl | 29 ++++++++++--------- .../test/mirk_basic_tests.jl | 20 +++++++++++++ 2 files changed, 36 insertions(+), 13 deletions(-) diff --git a/lib/BoundaryValueDiffEqMIRK/src/adaptivity.jl b/lib/BoundaryValueDiffEqMIRK/src/adaptivity.jl index c81b8332f..62d94455f 100644 --- a/lib/BoundaryValueDiffEqMIRK/src/adaptivity.jl +++ b/lib/BoundaryValueDiffEqMIRK/src/adaptivity.jl @@ -7,7 +7,7 @@ After we construct an interpolant, we use interp_eval to evaluate it. i = interval(mesh, t) dt = mesh_dt[i] τ = (t - mesh[i]) / dt - w, w′ = interp_weights(τ, cache.alg) + w, _ = interp_weights(τ, cache.alg) sum_stages!(y, cache, w, i, dt) return y end @@ -626,32 +626,35 @@ function sum_stages!(cache::MIRKCache{iip, T, use_both, NoDiffCacheNeeded}, w, sum_stages!(cache.fᵢ_cache, cache.fᵢ₂_cache, cache, w, w′, i, dt) end +# Here we should not directly in-place change z in several steps +# because in final step we actually need to use the original z(which is cache.y₀.u[i]) +# we use fᵢ₂_cache to avoid additional allocations. function sum_stages!(z::AbstractArray, cache::MIRKCache{iip, T, use_both, DiffCacheNeeded}, w, i::Int, dt = cache.mesh_dt[i]) where {iip, T, use_both} - (; stage, k_discrete, k_interp) = cache + (; stage, k_discrete, k_interp, fᵢ₂_cache) = cache (; s_star) = cache.ITU - z .= zero(z) - __maybe_matmul!(z, k_discrete[i].du[:, 1:stage], w[1:stage]) + fᵢ₂_cache .= zero(z) + __maybe_matmul!(fᵢ₂_cache, k_discrete[i].du[:, 1:stage], w[1:stage]) __maybe_matmul!( - z, k_interp.u[i][:, 1:(s_star - stage)], w[(stage + 1):s_star], true, true) - z .= z .* dt .+ cache.y₀.u[i] + fᵢ₂_cache, k_interp.u[i][:, 1:(s_star - stage)], w[(stage + 1):s_star], true, true) + z .= fᵢ₂_cache .* dt .+ cache.y₀.u[i] - return z + return nothing end function sum_stages!( z::AbstractArray, cache::MIRKCache{iip, T, use_both, NoDiffCacheNeeded}, w, i::Int, dt = cache.mesh_dt[i]) where {iip, T, use_both} - (; stage, k_discrete, k_interp) = cache + (; stage, k_discrete, k_interp, fᵢ₂_cache) = cache (; s_star) = cache.ITU - z .= zero(z) - __maybe_matmul!(z, k_discrete[i][:, 1:stage], w[1:stage]) + fᵢ₂_cache .= zero(z) + __maybe_matmul!(fᵢ₂_cache, k_discrete[i][:, 1:stage], w[1:stage]) __maybe_matmul!( - z, k_interp.u[i][:, 1:(s_star - stage)], w[(stage + 1):s_star], true, true) - z .= z .* dt .+ cache.y₀.u[i] + fᵢ₂_cache, k_interp.u[i][:, 1:(s_star - stage)], w[(stage + 1):s_star], true, true) + z .= fᵢ₂_cache .* dt .+ cache.y₀.u[i] - return z + return nothing end @views function sum_stages!(z, z′, cache::MIRKCache{iip, T, use_both, DiffCacheNeeded}, w, diff --git a/lib/BoundaryValueDiffEqMIRK/test/mirk_basic_tests.jl b/lib/BoundaryValueDiffEqMIRK/test/mirk_basic_tests.jl index 70b0c2e45..ab72a287f 100644 --- a/lib/BoundaryValueDiffEqMIRK/test/mirk_basic_tests.jl +++ b/lib/BoundaryValueDiffEqMIRK/test/mirk_basic_tests.jl @@ -381,3 +381,23 @@ end prob, MIRK4(), dt = 0.05, abstol = 1e-5, controller = error_control) end end + +# https://github.com/SciML/BoundaryValueDiffEq.jl/issues/319 +@testitem "Test interpolant evaluation with big defect" setup=[MIRKConvergenceTests] begin + function lotka!(du, u, p, t) + x = u[1] + y = u[2] + du[1] = p[1] * x - p[2] * x * y + du[2] = -p[3] * y + p[4] * x * y + end + + function bc!(resid, u, p, t) + resid[1] = u(0.0)[1] - 1.0 + resid[2] = u(0.0)[2] - 2.0 + end + bvp = BVProblem(lotka!, bc!, [1.0, 2.0], (0, 10), [7.5, 4.0, 8, 5]) + for order in (4, 5, 6) + sol = solve(bvp, mirk_solver(Val(order)), dt = 0.01) + @test SciMLBase.successful_retcode(sol) + end +end From c689ac59016056123b19f7461d5d69f86c544ee7 Mon Sep 17 00:00:00 2001 From: Qingyu Qu <2283984853@qq.com> Date: Sun, 11 May 2025 23:45:36 +0800 Subject: [PATCH 2/2] Bump version --- lib/BoundaryValueDiffEqMIRK/Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/BoundaryValueDiffEqMIRK/Project.toml b/lib/BoundaryValueDiffEqMIRK/Project.toml index 4a53bc965..74d250551 100644 --- a/lib/BoundaryValueDiffEqMIRK/Project.toml +++ b/lib/BoundaryValueDiffEqMIRK/Project.toml @@ -1,6 +1,6 @@ name = "BoundaryValueDiffEqMIRK" uuid = "1a22d4ce-7765-49ea-b6f2-13c8438986a6" -version = "1.6.1" +version = "1.6.2" [deps] ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"