Solvers using `VectorOfArray` behave differently for different `parent` array types

[jlchan]

Describe the bug 🐞

Solutions represented as nested arrays (e.g., AbstractArray{<:AbstractArray}) are supposed to be wrapped in VectorOfArray. However, different parent array representations result in different solver behavior. For example, wrapping StructArray{SVector{1, Float64}} with VectorOfArray results in a different number of adaptive time-steps than wrapping Vector{SVector{1, Float64}} with VectorOfArray.

Expected behavior

I would expect the behavior of an adaptive time-stepper to be independent of the underlying storage format.

Minimal Reproducible Example 👇

This runs the same solver with two solution representations (which should be equivalent). The number of time-steps is different. Changing the solver to Tsit5() also returns different numbers of time-steps (5 vs 3).

using OrdinaryDiffEqSSPRK, RecursiveArrayTools
using StaticArrays, StructArrays

function rhs!(du_voa, u_voa, params, t)
    du = parent(du_voa)
    u = parent(u_voa)
    du .= u
end

# StructArray storage
u = StructArray{SVector{1, Float64}}(ntuple(_ -> [1.0, 2.0], 1))
ode = ODEProblem(rhs!, VectorOfArray(u), (0, 0.7))
sol_SA = solve(ode, SSPRK43())
@show sol_SA.stats.naccept

# Vector{<:SVector} storage
u = SVector{1, Float64}.([1.0, 2.0])
ode = ODEProblem(rhs!, VectorOfArray(u), (0, 0.7))
sol_SV = solve(ode, SSPRK43())
@show sol_SV.stats.naccept

Error & Stacktrace ⚠️

Not an error; just the output. I expected the two naccept values to coincide.

sol_SA.stats.naccept = 3
sol_SV.stats.naccept = 9

Environment (please complete the following information):

• Output of using Pkg; Pkg.status()

Status `~/Downloads/scratch/Trixi_v0p10/MWE/Project.toml`
  [669c94d9] OrdinaryDiffEqSSPRK v1.2.0
  [731186ca] RecursiveArrayTools v3.31.0
  [90137ffa] StaticArrays v1.9.13
  [09ab397b] StructArrays v0.7.0

• Output of using Pkg; Pkg.status(; mode = PKGMODE_MANIFEST)

Status `~/Downloads/scratch/Trixi_v0p10/MWE/Manifest.toml`
  [47edcb42] ADTypes v1.14.0
  [7d9f7c33] Accessors v0.1.42
  [79e6a3ab] Adapt v4.2.0
  [4fba245c] ArrayInterface v7.18.0
  [62783981] BitTwiddlingConvenienceFunctions v0.1.6
  [2a0fbf3d] CPUSummary v0.2.6
  [fb6a15b2] CloseOpenIntervals v0.1.13
  [38540f10] CommonSolve v0.2.4
  [f70d9fcc] CommonWorldInvalidations v1.0.0
  [34da2185] Compat v4.16.0
  [a33af91c] CompositionsBase v0.1.2
  [2569d6c7] ConcreteStructs v0.2.3
  [187b0558] ConstructionBase v1.5.8
  [adafc99b] CpuId v0.3.1
  [9a962f9c] DataAPI v1.16.0
  [864edb3b] DataStructures v0.18.20
  [e2d170a0] DataValueInterfaces v1.0.0
  [2b5f629d] DiffEqBase v6.164.2
  [ffbed154] DocStringExtensions v0.9.3
  [4e289a0a] EnumX v1.0.4
  [f151be2c] EnzymeCore v0.8.8
  [e2ba6199] ExprTools v0.1.10
  [55351af7] ExproniconLite v0.10.14
  [7034ab61] FastBroadcast v0.3.5
  [9aa1b823] FastClosures v0.3.2
  [a4df4552] FastPower v1.1.1
  [1a297f60] FillArrays v1.13.0
  [069b7b12] FunctionWrappers v1.1.3
  [77dc65aa] FunctionWrappersWrappers v0.1.3
  [46192b85] GPUArraysCore v0.2.0
  [615f187c] IfElse v0.1.1
  [3587e190] InverseFunctions v0.1.17
  [82899510] IteratorInterfaceExtensions v1.0.0
  [ae98c720] Jieko v0.2.1
  [10f19ff3] LayoutPointers v0.1.17
  [1914dd2f] MacroTools v0.5.15
  [d125e4d3] ManualMemory v0.1.8
  [2e0e35c7] Moshi v0.3.5
  [46d2c3a1] MuladdMacro v0.2.4
  [bac558e1] OrderedCollections v1.8.0
  [bbf590c4] OrdinaryDiffEqCore v1.19.0
  [669c94d9] OrdinaryDiffEqSSPRK v1.2.0
  [d96e819e] Parameters v0.12.3
  [f517fe37] Polyester v0.7.16
  [1d0040c9] PolyesterWeave v0.2.2
  [aea7be01] PrecompileTools v1.2.1
  [21216c6a] Preferences v1.4.3
  [3cdcf5f2] RecipesBase v1.3.4
  [731186ca] RecursiveArrayTools v3.31.0
  [189a3867] Reexport v1.2.2
  [ae029012] Requires v1.3.1
  [7e49a35a] RuntimeGeneratedFunctions v0.5.13
  [94e857df] SIMDTypes v0.1.0
  [0bca4576] SciMLBase v2.75.1
  [c0aeaf25] SciMLOperators v0.3.12
  [53ae85a6] SciMLStructures v1.7.0
  [efcf1570] Setfield v1.1.2
  [ce78b400] SimpleUnPack v1.1.0
  [aedffcd0] Static v1.2.0
  [0d7ed370] StaticArrayInterface v1.8.0
  [90137ffa] StaticArrays v1.9.13
  [1e83bf80] StaticArraysCore v1.4.3
  [7792a7ef] StrideArraysCore v0.5.7
  [09ab397b] StructArrays v0.7.0
  [2efcf032] SymbolicIndexingInterface v0.3.38
  [3783bdb8] TableTraits v1.0.1
  [bd369af6] Tables v1.12.0
  [8290d209] ThreadingUtilities v0.5.2
  [781d530d] TruncatedStacktraces v1.4.0
  [3a884ed6] UnPack v1.0.2
  [56f22d72] Artifacts
  [2a0f44e3] Base64
  [ade2ca70] Dates
  [8ba89e20] Distributed
  [9fa8497b] Future
  [b77e0a4c] InteractiveUtils
  [76f85450] LibGit2
  [8f399da3] Libdl
  [37e2e46d] LinearAlgebra
  [56ddb016] Logging
  [d6f4376e] Markdown
  [ca575930] NetworkOptions v1.2.0
  [de0858da] Printf
  [9a3f8284] Random
  [ea8e919c] SHA v0.7.0
  [9e88b42a] Serialization
  [6462fe0b] Sockets
  [2f01184e] SparseArrays v1.10.0
  [10745b16] Statistics v1.10.0
  [fa267f1f] TOML v1.0.3
  [cf7118a7] UUIDs
  [4ec0a83e] Unicode
  [e66e0078] CompilerSupportLibraries_jll v1.1.1+0
  [e37daf67] LibGit2_jll v1.6.4+0
  [29816b5a] LibSSH2_jll v1.11.0+1
  [c8ffd9c3] MbedTLS_jll v2.28.2+1
  [4536629a] OpenBLAS_jll v0.3.23+4
  [bea87d4a] SuiteSparse_jll v7.2.1+1
  [8e850b90] libblastrampoline_jll v5.11.0+0

• Output of versioninfo()

Julia Version 1.10.8
Commit 4c16ff44be8 (2025-01-22 10:06 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: macOS (arm64-apple-darwin24.0.0)
  CPU: 12 × Apple M2 Max
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-15.0.7 (ORCJIT, apple-m1)
Threads: 8 default, 0 interactive, 4 GC (on 8 virtual cores)
Environment:
  JULIA_EDITOR = code
  JULIA_NUM_THREADS = 8

Additional context

If you have a suggestion for what might be causing this discrepancy, I'm happy to try to fix it. I would guess that some norm-like quantity in the adaptive time-stepping algorithm causes this, but I don't know which one.

[jlchan]

@ranocha any ideas what might be causing this discrepancy in SSPRK43?

[ChrisRackauckas]

If I had to guess, I'd check if they are handling the FastBroadcast.jl the same? Are the ODE_DEFAULT_NORM definitions on the arrays giving the same number? Different floating point accumulations could cause at least small differences, so the summation ordering in those two would be the thing to look at.

[ranocha]

I can't think of anything in the ODE solvers themselves that causes this issue. Are the norm computations the same?

[jlchan]

I can't think of anything in the ODE solvers themselves that causes this issue. Are the norm computations the same?

Yes, at least for the 1 and 2 norms.

If I had to guess, I'd check if they are handling the FastBroadcast.jl the same? Are the ODE_DEFAULT_NORM definitions on the arrays giving the same number? Different floating point accumulations could cause at least small differences, so the summation ordering in those two would be the thing to look at.

DiffEqBase.ODE_DEFAULT_NORM(_, 0.0) returns the same thing for both array types (all digits identical).

For FastBroadcast.jl, I'll try to check in more detail tomorrow, but they seem to both behave identically at the moment.

[jlchan]

Looking at the solution history, the two solutions start to differ only at the second time-step. They're exactly the same for the initial condition and first time-step.

Could this be FSAL related?

[ranocha]

Is the second time step size the same?

[jlchan]

Yes - the second time-step size is exactly the same. The time-step size and solution size start to differ in the third time-step.

EDIT: looking more closely, after one step! call, the different parent array types have different values of dtpropose. The values of q11 and EEst are different as well, but they seem to be computed by opts.internalnorm, which seems to be set to ODE_DEFAULT_NORM.

MWE:

using OrdinaryDiffEqSSPRK, RecursiveArrayTools
using StaticArrays, StructArrays

function rhs!(du_voa, u_voa, params, t)
    du = parent(du_voa)
    u = parent(u_voa)
    du .= u
end

# StructArray storage
u_SA = StructArray{SVector{1, Float64}}(ntuple(_ -> [1.0, 2.0], 1))
u_SV = SVector{1, Float64}.([1.0, 2.0])

iSA = init(ODEProblem(rhs!, VectorOfArray(u_SA), (0, 0.7)), SSPRK43())
iSV = init(ODEProblem(rhs!, VectorOfArray(u_SV), (0, 0.7)), SSPRK43())
step!(iSA)
step!(iSV)

@show iSA.dtpropose, iSV.dtpropose

[jlchan]

One more data point - the solver behavior for VectorOfArray{Vector{<:SVector}} is consistent with the solver behavior for a Vector{Float64} solution, so it's probably VectorOfArray{StructArray} that's causing this issue.

[ranocha]

It has to do with the adaptive time stepping since I get the same results for fixed time step sizes:

julia> using OrdinaryDiffEqSSPRK, RecursiveArrayTools, StaticArrays, StructArrays

julia> function rhs!(du_voa, u_voa, params, t)
           du = parent(du_voa)
           u = parent(u_voa)
           du .= u
       end

julia> begin
       # StructArray storage
       u = StructArray{SVector{1, Float64}}(ntuple(_ -> [1.0, 2.0], 1))
       ode = ODEProblem(rhs!, VectorOfArray(u), (0, 0.7))
       sol_SA = solve(ode, SSPRK43(), adaptive = false, dt = 0.1)
       #@show sol_SA.u[end]

       # Vector{<:SVector} storage
       u = SVector{1, Float64}.([1.0, 2.0])
       ode = ODEProblem(rhs!, VectorOfArray(u), (0, 0.7))
       sol_SV = solve(ode, SSPRK43(), adaptive = false, dt = 0.1)
       #@show sol_SV.u[end]

       # Difference
       sol_SA.u[end] - sol_SV.u[end]
       end
VectorOfArray{Float64,2}:
2-element StructArray(::Vector{Float64}) with eltype SVector{1, Float64}:
 [0.0]
 [0.0]

[ranocha]

Do calls to DiffEqBase.calculate_residuals! like in

OrdinaryDiffEq.jl/lib/OrdinaryDiffEqSSPRK/src/ssprk_perform_step.jl

Lines 815 to 817 in 8fa9987

	calculate_residuals!(atmp, utilde, uprev, u, integrator.opts.abstol,
	integrator.opts.reltol, integrator.opts.internalnorm, t,
	thread)

give the same results? The code can be found in https://github.com/SciML/DiffEqBase.jl/blob/master/src/calculate_residuals.jl

[jlchan]

No, and I think it's because of the way StructArrays treats assignment: https://juliaarrays.github.io/StructArrays.jl/stable/counterintuitive/#Modifying-a-field-of-a-struct-element. For example:

u = VectorOfArray(StructArray{MVector{1, Float64}}(ntuple(_ -> [1.0, 2.0], 1)))
u[1,1] = 100
@show u[1,1] # this is still 1.0

I think this can be fixed by special-casing VectorOfArray{StructArray} in RecursiveArrayTools.jl.

[jlchan]

Hm, SciML/RecursiveArrayTools.jl#425 doesn't fix the issue. I'll take a closer look at this over the weekend.

[jlchan]

OK, SciML/RecursiveArrayTools.jl#425 only partially addressed the issue. The problem appears to be in calculate_residuals! where the assignment is broadcasted.

using StaticArrays, StructArrays, RecursiveArrayTools
u = VectorOfArray(StructArray{MVector{1, Float64}}(ntuple(_ -> [1.0, 2.0], 1)))
v = 2 * copy(u)
u[1,1] = 100 # works now
@. u = v     # doesn't work, @.. gives the same behavior too