Fix __resize! on AbstractVectorOfArray under RecursiveArrayTools v4

Closes #484.

`__resize!(::AbstractVectorOfArray, n, M)` was missed by the
SciMLBase-v3 / RecursiveArrayTools-v4 compatibility patch
(bea82b8). Under RAT v4, `AbstractVectorOfArray` subtypes
`AbstractArray` with column-major linear-indexing semantics, so:

- `length(x)`  is now the total scalar count (`prod(size(x))`),
  not the number of timesteps.
- `last(x)`    is now the last scalar, not the last timestep vector.
- `append!(x, ::AbstractVectorOfArray)` falls back to the generic
  AbstractArray `append!`, which iterates scalars rather than
  timestep vectors.

The old code in `__resize!`:

    N = n - length(x)
    N > 0 ? append!(x, VectorOfArray([safe_similar(last(x)) for _ in 1:N])) :
        resize!(x, n)

silently miscomputes `N` (target timesteps minus scalar count)
and never grows `x.u`, so subsequent `cache.y₀.u[i]` reads on
the post-mesh-refinement loop in MIRK / FIRK hit `#undef` and
throw `UndefRefError` from `mirk.jl:342` (`firk.jl` analog).

Switch to `x.u` for both `length` and `last`, and `append!`
directly into `x.u` instead of constructing a wrapper
`VectorOfArray`. `resize!(x, n)` already delegates to
`resize!(x.u, n)` in RAT v4 and is left as-is.

The same pattern fix applies to the FIRK-specific overloads
in `BoundaryValueDiffEqFIRK/src/utils.jl` for both
`FIRKTableau{false}` and `FIRKTableau{true}`.

Adds a regression test based on the issue MWE (torus geodesic
BVP that triggers adaptive mesh refinement).

Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>

Author: ChrisRackauckas

lib/BoundaryValueDiffEqCore/Project.toml

@@ -1,6 +1,6 @@
 name = "BoundaryValueDiffEqCore"
 uuid = "56b672f2-a5fe-4263-ab2d-da677488eb3a"
-version = "2.5.0"
+version = "2.5.1"
 authors = ["Qingyu Qu <erikqqy123@gmail.com>"]
 
 [deps]

lib/BoundaryValueDiffEqCore/src/utils.jl

@@ -262,9 +262,9 @@ function __resize!(x::AbstractVector{<:DiffCache}, n, M)
 end
 
 function __resize!(x::AbstractVectorOfArray, n, M)
-    N = n - length(x)
+    N = n - length(x.u)
     N == 0 && return x
-    N > 0 ? append!(x, VectorOfArray([safe_similar(last(x)) for _ in 1:N])) : resize!(x, n)
+    N > 0 ? append!(x.u, [safe_similar(last(x.u)) for _ in 1:N]) : resize!(x, n)
     return x
 end
 

lib/BoundaryValueDiffEqFIRK/src/utils.jl

@@ -24,18 +24,18 @@ end
 
 function BoundaryValueDiffEqCore.__resize!(x::AbstractVectorOfArray, n, M, TU::FIRKTableau{false})
     (; s) = TU
-    N = (n - 1) * (s + 1) + 1 - length(x)
+    N = (n - 1) * (s + 1) + 1 - length(x.u)
     N == 0 && return x
-    N > 0 ? append!(x, VectorOfArray([safe_similar(last(x)) for _ in 1:N])) :
+    N > 0 ? append!(x.u, [safe_similar(last(x.u)) for _ in 1:N]) :
         resize!(x, (n - 1) * (s + 1) + 1)
     return x
 end
 
 function BoundaryValueDiffEqCore.__resize!(x::AbstractVectorOfArray, n, M, TU::FIRKTableau{true})
     (; s) = TU
-    N = n - length(x)
+    N = n - length(x.u)
     N == 0 && return x
-    N > 0 ? append!(x, VectorOfArray([safe_similar(last(x)) for _ in 1:N])) : resize!(x, n)
+    N > 0 ? append!(x.u, [safe_similar(last(x.u)) for _ in 1:N]) : resize!(x, n)
     return x
 end
 @inline __K0_on_u0(

lib/BoundaryValueDiffEqMIRK/test/mirk_basic_tests.jl

@@ -687,3 +687,39 @@ end
     sol2 = solve(bvp2, MIRK4(), dt = 0.1, adaptive = false, nlsolve_kwargs = (; maxiters = 0))
     @test sol2.u == u_guess
 end
+
+# https://github.com/SciML/BoundaryValueDiffEq.jl/issues/484
+# Adaptive mesh refinement on a sufficiently nonlinear BVP must not throw
+# UndefRefError. The torus geodesic system is smooth and nonlinear enough
+# that the solver decides to refine the mesh.
+@testitem "Adaptive mesh refinement (issue 484)" begin
+    using BoundaryValueDiffEqMIRK
+
+    R = 3.0
+    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))
+    sol = solve(prob, MIRK4(); dt = 0.05)
+    @test SciMLBase.successful_retcode(sol)
+    @test sol.u[1][1:2] ≈ a1 atol = 1.0e-8
+    @test sol.u[end][1:2] ≈ a2 atol = 1.0e-8
+end