MIRK 1.15+ adaptive mesh refinement: UndefRefError on nonlinear BVPs
## Summary
`BoundaryValueDiffEqMIRK` 1.15.0 and 1.16.0 throw `UndefRefError` from inside `__perform_mirk_iteration` (mesh refinement code) when solving sufficiently nonlinear BVPs with default kwargs. Setting `adaptive = false` is a clean workaround. MIRK 1.14.0 works correctly on the same problem (with old SciMLBase 2 stack).
## Versions
- `BoundaryValueDiffEqMIRK` v1.15.0, v1.16.0 — affected
- `BoundaryValueDiffEqMIRK` v1.14.0 — not affected
- `BoundaryValueDiffEqCore` 2.5.0
- `SciMLBase` 3.7.1
- `RecursiveArrayTools` 4.3.0
- Julia 1.12.6
## MWE
Standalone — no Manifolds.jl needed. The ODE is the torus geodesic system on a `(R=3, r=2)` torus (smooth and nonlinear enough that the solver decides to refine the mesh).
```julia
using BoundaryValueDiffEqMIRK
const R = 3.0
const r = 2.0
function f!(du, u, p, t)
sinθ, cosθ = sincos(u[1])
R_θ = R + r * cosθ
du[1] = u[3]
du[2] = u[4]
du[3] = -u[4]^2 * R_θ * sinθ / r
du[4] = 2 * r * sinθ / R_θ * u[3] * u[4]
end
a1 = [0.5, -1.2]
a2 = [-0.5, 0.3]
function bc!(residual, u, p, t)
ua = u(0.0)
ub = u(1.0)
residual[1:2] = ua[1:2] - a1
residual[3:4] = ub[1:2] - a2
end
prob = BVProblem(f!, bc!, vcat(a1, zero(a1)), (0.0, 1.0))
solve(prob, MIRK4(); dt = 0.05) # ERROR: UndefRefError
solve(prob, MIRK4(); dt = 0.05, adaptive = false) # Success
```
## Failing trace
```
ERROR: UndefRefError: access to undefined reference
Stacktrace:
[1] getindex
@ ./essentials.jl:920 [inlined]
[2] __perform_mirk_iteration(cache::BoundaryValueDiffEqMIRK.MIRKCache{...})
@ BoundaryValueDiffEqMIRK ~/.julia/packages/BoundaryValueDiffEqMIRK/dR6jB/src/mirk.jl:342
[3] solve!(cache::BoundaryValueDiffEqMIRK.MIRKCache{...})
@ BoundaryValueDiffEqMIRK ~/.julia/packages/BoundaryValueDiffEqMIRK/dR6jB/src/mirk.jl:275
...
```
## Likely culprit
`mirk.jl:342` in MIRK 1.16 (https://github.com/SciML/BoundaryValueDiffEq.jl/blob/master/lib/BoundaryValueDiffEqMIRK/src/mirk.jl#L334-L346):
```julia
if info == ReturnCode.Success # Nonlinear Solve Successful and defect norm is acceptable
if error_norm > abstol
# We construct a new mesh to equidistribute the defect
mesh, mesh_dt, _, info = mesh_selector!(cache, controller)
if info == ReturnCode.Success
(length(mesh) < length(cache.mesh)) &&
__resize!(cache.y₀, length(cache.mesh), cache.M)
for (i, m) in enumerate(cache.mesh)
interp_eval!(cache.y₀.u[i], cache, m, mesh, mesh_dt)
end
__expand_cache!(cache)
end
end
end
```
`__resize!` only fires when the new `mesh` is strictly shorter than `cache.mesh`. The subsequent loop iterates `cache.mesh` and writes into `cache.y₀.u[i]`, so it relies on `length(cache.y₀.u) >= length(cache.mesh)`. If `cache.y₀.u` was allocated with placeholder slots that were never filled (the `safe_similar(u0, ifelse(adaptive, N, 0))` path in `__init`, where the inner vectors are length-0), then `cache.y₀.u[i]` may be an undef-ref `Vector{Float64}` and the `getindex` inside `interp_eval!` blows up.
## Impact
This blocks every Manifolds.jl user who hits `solve_chart_log_bvp` / `estimate_distance_from_bvp` on the SciMLBase 3 stack (https://github.com/JuliaManifolds/Manifolds.jl/pull/887). The Manifolds extension currently works around it with `adaptive = false`, but ideally users should be able to opt back into adaptive refinement once this is fixed.
## Workaround
Pass `adaptive = false` to `solve` when using MIRK on a nonlinear BVP that would otherwise trigger refinement.
0 条评论