Broadcast over `VectorOfArray` doesn't respect dimensions of non-vector parent arrays

[jlchan]

Is your feature request related to a problem? Please describe.

VectorOfArray was extended to multi-dimensional parent arrays in https://github.com/SciML/RecursiveArrayTools.jl/pull/357. However, these arrays don't work with all solvers in OrdinaryDiffEq.jl (for example, SSPRK43).

The error happens during the call to integrator = init(ode, SSPRK43()), specifically this line: https://github.com/SciML/OrdinaryDiffEq.jl/blob/1f2e058c14d2a5c997bf7d5f055cebfb406cbc1b/src/initdt.jl#L51. An MWE is

using OrdinaryDiffEq, RecursiveArrayTools

function rhs!(duu::VectorOfArray, uu::VectorOfArray, p, t)
    du = parent(duu)
    u = parent(uu)
    du .= u
end

u = fill(SVector{2}(ones(2)), 2, 3)
ode = ODEProblem(rhs!, VectorOfArray(u), (0.0, 1.0))
integrator = init(ode, SSPRK43())

The manually truncated stacktrace is

ERROR: DimensionMismatch: array could not be broadcast to match destination
Stacktrace:
  [1] check_broadcast_shape
    @ ./broadcast.jl:579 [inlined]
  [2] check_broadcast_axes
    @ ./broadcast.jl:582 [inlined]
  [3] instantiate
    @ ./broadcast.jl:309 [inlined]
  [4] materialize!
    @ ./broadcast.jl:914 [inlined]
  [5] materialize!
    @ ./broadcast.jl:911 [inlined]
  [6] rhs!(duu::VectorOfArray{…}, uu::VectorOfArray{…}, p::SciMLBase.NullParameters, t::Float64)
    @ Main ./Untitled-1:7
  [7] (::ODEFunction{…})(::VectorOfArray{…}, ::Vararg{…})
    @ SciMLBase ~/.julia/packages/SciMLBase/hSv8d/src/scimlfunctions.jl:2296 [inlined]
  [8] ode_determine_initdt(u0::VectorOfArray{…}, t::Float64, tdir::Float64, dtmax::Float64, abstol::Float64, reltol::Float64, internalnorm::typeof(DiffEqBase.ODE_DEFAULT_NORM), prob::ODEProblem{…}, integrator::OrdinaryDiffEq.ODEIntegrator{…})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/0tf1M/src/initdt.jl:53
  [9] auto_dt_reset!(integrator::OrdinaryDiffEq.ODEIntegrator)
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/0tf1M/src/integrators/integrator_interface.jl:453 [inlined]
...
...

Describe the solution you’d like

The error above is due to the fact that for u::VectorOfArray, zero.(u) returns a VectorOfArray whose parent array is a Vector. For VectorOfArray with multi-dimensional parents, the broadcast should ideally respect the dimension of the parent array.

Additional context

I'm looking at this now, but I'm not the most experienced with broadcast. If anyone has any suggestions on where to start and/or how to implement, I'd welcome them.

[ChrisRackauckas]

For VectorOfArray with multi-dimensional parents, the broadcast should ideally respect the dimension of the parent array.

I think that's the proper way for it to be handled.

[jlchan]

Closed by #374