Fix incorrect interpolant evaluation

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"

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,

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