Fix similar and fill! for mixed nested VectorOfArray

[JoshuaLampert]

Checklist

• Appropriate tests were added
• Any code changes were done in a way that does not break public API
• All documentation related to code changes were updated
• The new code follows the
  contributor guidelines, in particular the SciML Style Guide and
  COLPRAC.
• Any new documentation only uses public API

Additional context

This is similar to #542. The same MWE from there, but now with SSPRK43():

using RecursiveArrayTools
using OrdinaryDiffEqSSPRK

function rhs!(du, u, p, t)
    du .= u
    return nothing
end

n = 2
left = ones(n)
right = left

u_flat = VectorOfArray([copy(left), copy(right)])
u_nested = VectorOfArray([copy(left),
                          VectorOfArray([right[e] for e in 1:n])])

tspan = (0.0, 1.0)
prob_flat = ODEProblem(rhs!, u_flat, tspan)
prob_nested = ODEProblem(rhs!, u_nested, tspan)

solver = SSPRK43()
kwargs = (; dt = 0.1, adaptive = false, save_everystep = false)

sol_flat = solve(prob_flat, solver; kwargs...)
sol_nested = solve(prob_nested, solver; kwargs...)

previously it errored:

ERROR: LoadError: MethodError: no method matching similar(::Float64, ::Type{Float64})
The function `similar` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  similar(::JuliaInterpreter.NonRecursiveInterpreter, ::Any)
   @ JuliaInterpreter ~/.julia/packages/JuliaInterpreter/KPJTC/src/types.jl:40
  similar(::Type{A}, ::Type{T}, ::StaticArraysCore.Size{S}) where {A<:Array, T, S}
   @ StaticArrays ~/.julia/packages/StaticArrays/mI3IC/src/abstractarray.jl:136
  similar(::Type{SA}, ::Type{T}, ::StaticArraysCore.Size{S}) where {SA<:StaticArraysCore.SizedArray, T, S}
   @ StaticArrays ~/.julia/packages/StaticArrays/mI3IC/src/abstractarray.jl:135
  ...

Stacktrace:
  [1] _broadcast_getindex_evalf
    @ ./broadcast.jl:699 [inlined]
  [2] _broadcast_getindex
    @ ./broadcast.jl:687 [inlined]
  [3] _getindex
    @ ./broadcast.jl:620 [inlined]
  [4] getindex
    @ ./broadcast.jl:616 [inlined]
  [5] copy
    @ ./broadcast.jl:933 [inlined]
  [6] materialize
    @ ./broadcast.jl:894 [inlined]
  [7] similar
    @ ~/.julia/packages/RecursiveArrayTools/AvG3i/src/vector_of_array.jl:1407 [inlined]
  [8] _broadcast_getindex_evalf
    @ ./broadcast.jl:699 [inlined]
  [9] _broadcast_getindex
    @ ./broadcast.jl:687 [inlined]
 [10] _getindex
    @ ./broadcast.jl:620 [inlined]
 [11] getindex
    @ ./broadcast.jl:616 [inlined]
 [12] copyto_nonleaf!(dest::Vector{…}, bc::Base.Broadcast.Broadcasted{…}, iter::Base.OneTo{…}, state::Int64, count::Int64)
    @ Base.Broadcast ./broadcast.jl:1104
 [13] copy
    @ ./broadcast.jl:941 [inlined]
 [14] materialize(bc::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1}, Nothing, typeof(similar), Tuple{Vector{Union{VectorOfArray{…}, Vector{…}}}, Base.RefValue{Type{Float64}}}})
    @ Base.Broadcast ./broadcast.jl:894
 [15] similar(VA::VectorOfArray{Float64, 2, Vector{Union{VectorOfArray{Float64, 1, Vector{Float64}}, Vector{Float64}}}}, ::Type{Float64})
    @ RecursiveArrayTools ~/.julia/packages/RecursiveArrayTools/AvG3i/src/vector_of_array.jl:1407
 [16] alg_cache(alg::SSPRK43{…}, u::VectorOfArray{…}, rate_prototype::VectorOfArray{…}, ::Type{…}, ::Type{…}, ::Type{…}, uprev::VectorOfArray{…}, uprev2::VectorOfArray{…}, f::ODEFunction{…}, t::Float64, dt::Float64, reltol::Float64, p::SciMLBase.NullParameters, calck::Bool, ::Val{…}, verbose::OrdinaryDiffEqCore.DEVerbosity{…})
    @ OrdinaryDiffEqSSPRK ~/.julia/packages/OrdinaryDiffEqSSPRK/4XvhT/src/ssprk_caches.jl:816
 [17] __init(prob::ODEProblem{…}, alg::SSPRK43{…}, timeseries_init::Tuple{}, ts_init::Tuple{}, ks_init::Tuple{}; saveat::Tuple{}, tstops::Tuple{}, d_discontinuities::Tuple{}, save_idxs::Nothing, save_everystep::Bool, save_on::Bool, save_discretes::Bool, save_start::Bool, save_end::Nothing, callback::Nothing, dense::Bool, calck::Bool, dt::Float64, dtmin::Float64, dtmax::Float64, force_dtmin::Bool, adaptive::Bool, abstol::Nothing, reltol::Nothing, gamma::Nothing, qmin::Nothing, qmax::Nothing, qsteady_min::Nothing, qsteady_max::Nothing, beta1::Nothing, beta2::Nothing, qoldinit::Nothing, fullnormalize::Bool, failfactor::Int64, maxiters::Int64, internalnorm::typeof(DiffEqBase.ODE_DEFAULT_NORM), internalopnorm::typeof(LinearAlgebra.opnorm), isoutofdomain::typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), unstable_check::typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), verbose::SciMLLogging.Standard, controller::Nothing, timeseries_errors::Bool, dense_errors::Bool, advance_to_tstop::Bool, stop_at_next_tstop::Bool, initialize_save::Bool, progress::Bool, progress_steps::Int64, progress_name::String, progress_message::typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), progress_id::Symbol, userdata::Nothing, allow_extrapolation::Bool, initialize_integrator::Bool, alias::ODEAliasSpecifier, initializealg::DiffEqBase.DefaultInit, rng::Nothing, save_noise::Bool, delta::Nothing, W::Nothing, P::Nothing, sqdt::Nothing, noise::Nothing, c::Nothing, rate_constants::Nothing, _cache::Nothing, seed::UInt64, kwargs::@Kwargs{})
    @ OrdinaryDiffEqCore ~/.julia/packages/OrdinaryDiffEqCore/tusSd/src/solve.jl:495
 [18] __init (repeats 2 times)
    @ ~/.julia/packages/OrdinaryDiffEqCore/tusSd/src/solve.jl:19 [inlined]
 [19] __solve(::ODEProblem{…}, ::SSPRK43{…}; kwargs::@Kwargs{…})
    @ OrdinaryDiffEqCore ~/.julia/packages/OrdinaryDiffEqCore/tusSd/src/solve.jl:9
 [20] __solve
    @ ~/.julia/packages/OrdinaryDiffEqCore/tusSd/src/solve.jl:1 [inlined]
 [21] #solve_call#22
    @ ~/.julia/packages/DiffEqBase/4RFWW/src/solve.jl:172 [inlined]
 [22] solve_call
    @ ~/.julia/packages/DiffEqBase/4RFWW/src/solve.jl:137 [inlined]
 [23] #solve_up#29
    @ ~/.julia/packages/DiffEqBase/4RFWW/src/solve.jl:630 [inlined]
 [24] solve_up
    @ ~/.julia/packages/DiffEqBase/4RFWW/src/solve.jl:603 [inlined]
 [25] #solve#28
    @ ~/.julia/packages/DiffEqBase/4RFWW/src/solve.jl:587 [inlined]
 [26] top-level scope
    @ ~/.julia/dev/RecursiveArrayTools/run/mwe.jl:28
 [27] include(mapexpr::Function, mod::Module, _path::String)
    @ Base ./Base.jl:307
 [28] top-level scope
    @ REPL[22]:1
in expression starting at /home/lampert/.julia/dev/RecursiveArrayTools/run/mwe.jl:28
Some type information was truncated. Use `show(err)` to see complete types.

with this fix it runs through with:

julia> sol_nested.u[end]
VectorOfArray{Float64,2}:
2-element Vector{Union{VectorOfArray{Float64, 1, Vector{Float64}}, Vector{Float64}}}:
 [6.400934199221854, 6.400934199221854]
 VectorOfArray{Float64, 1, Vector{Float64}}([2.718228502873718, 2.718228502873718])

and on top of #542 it returns the correct:

julia> sol_nested.u[end]
VectorOfArray{Float64,2}:
2-element Vector{Union{VectorOfArray{Float64, 1, Vector{Float64}}, Vector{Float64}}}:
 [2.718228502873718, 2.718228502873718]
 VectorOfArray{Float64, 1, Vector{Float64}}([2.718228502873718, 2.718228502873718])

For transparency, I used AI tools to help me come up with this solution.