Skip to content

Using ArrayPartition of VectorOfArrays gives wrong results #2728

Open
@JoshuaLampert

Description

@JoshuaLampert

Describe the bug 🐞

I want to use an ArrayPartition of two VectorOfArrays as underlying data structure (I specifically want a vector of vector structure instead of a matrix structure because the inner vectors can have different lengths). However, the result overwrites the initial condition and all intermediate saved time steps by the result of the last time step (which is correct in this MRE, but also see the Additional context below).

Expected behavior

The result should not have the same solution for each time step, but rather the correct ones.

Minimal Reproducible Example 👇

using OrdinaryDiffEqLowOrderRK
using RecursiveArrayTools: ArrayPartition, VectorOfArray

function rhs!(dω, ω, p, t)
    du, dv =.x
    du .= 1.0
    dv .= 1.0
end

u_ini = VectorOfArray([[1.0, 2.0], [3.0, 4.0]])
v_ini = VectorOfArray([[1.0, 2.0], [3.0, 4.0]])
# u_ini = [1.0, 2.0, 3.0, 4.0] # This works as expected
# v_ini = [3.0, 4.0, 3.0, 4.0]
ω_ini = ArrayPartition(u_ini, v_ini)

tspan = (0.0, 1.0)
ode = ODEProblem(rhs!, ω_ini, tspan)
sol = solve(ode, Euler(), dt = 0.2)

Error & Stacktrace ⚠️

The result is

retcode: Success
Interpolation: 3rd order Hermite
t: 6-element Vector{Float64}:
 0.0
 0.2
 0.4
 0.6000000000000001
 0.8
 1.0
u: 6-element Vector{ArrayPartition{Float64, Tuple{VectorOfArray{Float64, 2, Vector{Vector{Float64}}}, VectorOfArray{Float64, 2, Vector{Vector{Float64}}}}}}:
 (VectorOfArray{Float64, 2, Vector{Vector{Float64}}}([[1.9999999999999998, 3.000000000000001], [4.000000000000001, 5.000000000000001]]), VectorOfArray{Float64, 2, Vector{Vector{Float64}}}([[1.9999999999999998, 3.000000000000001], [4.000000000000001, 5.000000000000001]]))
 (VectorOfArray{Float64, 2, Vector{Vector{Float64}}}([[1.9999999999999998, 3.000000000000001], [4.000000000000001, 5.000000000000001]]), VectorOfArray{Float64, 2, Vector{Vector{Float64}}}([[1.9999999999999998, 3.000000000000001], [4.000000000000001, 5.000000000000001]]))
 (VectorOfArray{Float64, 2, Vector{Vector{Float64}}}([[1.9999999999999998, 3.000000000000001], [4.000000000000001, 5.000000000000001]]), VectorOfArray{Float64, 2, Vector{Vector{Float64}}}([[1.9999999999999998, 3.000000000000001], [4.000000000000001, 5.000000000000001]]))
 (VectorOfArray{Float64, 2, Vector{Vector{Float64}}}([[1.9999999999999998, 3.000000000000001], [4.000000000000001, 5.000000000000001]]), VectorOfArray{Float64, 2, Vector{Vector{Float64}}}([[1.9999999999999998, 3.000000000000001], [4.000000000000001, 5.000000000000001]]))
 (VectorOfArray{Float64, 2, Vector{Vector{Float64}}}([[1.9999999999999998, 3.000000000000001], [4.000000000000001, 5.000000000000001]]), VectorOfArray{Float64, 2, Vector{Vector{Float64}}}([[1.9999999999999998, 3.000000000000001], [4.000000000000001, 5.000000000000001]]))
 (VectorOfArray{Float64, 2, Vector{Vector{Float64}}}([[1.9999999999999998, 3.000000000000001], [4.000000000000001, 5.000000000000001]]), VectorOfArray{Float64, 2, Vector{Vector{Float64}}}([[1.9999999999999998, 3.000000000000001], [4.000000000000001, 5.000000000000001]]))

which is clearly wrong. The result should not be the same at each time step.

Environment (please complete the following information):

  • Output of using Pkg; Pkg.status()
Status `/tmp/jl_AoO3TL/Project.toml`
  [1344f307] OrdinaryDiffEqLowOrderRK v1.2.0
  [731186ca] RecursiveArrayTools v3.33.0
  • Output of using Pkg; Pkg.status(; mode = PKGMODE_MANIFEST)
  [47edcb42] ADTypes v1.14.0
  [7d9f7c33] Accessors v0.1.42
  [79e6a3ab] Adapt v4.3.0
  [4fba245c] ArrayInterface v7.19.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
  [a8cc5b0e] Crayons v4.1.1
  [9a962f9c] DataAPI v1.16.0
  [864edb3b] DataStructures v0.18.22
  [e2d170a0] DataValueInterfaces v1.0.0
  [2b5f629d] DiffEqBase v6.175.0
  [ffbed154] DocStringExtensions v0.9.4
  [4e289a0a] EnumX v1.0.5
  [f151be2c] EnzymeCore v0.8.9
  [e2ba6199] ExprTools v0.1.10
  [55351af7] ExproniconLite v0.10.14
  [7034ab61] FastBroadcast v0.3.5
  [9aa1b823] FastClosures v0.3.2
  [a4df4552] FastPower v1.1.2
  [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
  [b964fa9f] LaTeXStrings v1.4.0
  [10f19ff3] LayoutPointers v0.1.17
  [1914dd2f] MacroTools v0.5.16
  [d125e4d3] ManualMemory v0.1.8
  [2e0e35c7] Moshi v0.3.5
  [46d2c3a1] MuladdMacro v0.2.4
  [bac558e1] OrderedCollections v1.8.1
  [bbf590c4] OrdinaryDiffEqCore v1.26.0
  [1344f307] OrdinaryDiffEqLowOrderRK v1.2.0
  [d96e819e] Parameters v0.12.3
  [f517fe37] Polyester v0.7.17
  [1d0040c9] PolyesterWeave v0.2.2
⌅ [aea7be01] PrecompileTools v1.2.1
  [21216c6a] Preferences v1.4.3
  [08abe8d2] PrettyTables v2.4.0
  [3cdcf5f2] RecipesBase v1.3.4
  [731186ca] RecursiveArrayTools v3.33.0
  [189a3867] Reexport v1.2.2
  [ae029012] Requires v1.3.1
  [7e49a35a] RuntimeGeneratedFunctions v0.5.14
  [94e857df] SIMDTypes v0.1.0
  [0bca4576] SciMLBase v2.95.0
  [c0aeaf25] SciMLOperators v1.3.0
  [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
  [1e83bf80] StaticArraysCore v1.4.3
  [10745b16] Statistics v1.11.1
  [7792a7ef] StrideArraysCore v0.5.7
  [892a3eda] StringManipulation v0.4.1
  [2efcf032] SymbolicIndexingInterface v0.3.40
  [3783bdb8] TableTraits v1.0.1
  [bd369af6] Tables v1.12.0
  [8290d209] ThreadingUtilities v0.5.4
  [781d530d] TruncatedStacktraces v1.4.0
  [3a884ed6] UnPack v1.0.2
  [56f22d72] Artifacts v1.11.0
  [2a0f44e3] Base64 v1.11.0
  [ade2ca70] Dates v1.11.0
  [8ba89e20] Distributed v1.11.0
  [9fa8497b] Future v1.11.0
  [b77e0a4c] InteractiveUtils v1.11.0
  [8f399da3] Libdl v1.11.0
  [37e2e46d] LinearAlgebra v1.11.0
  [56ddb016] Logging v1.11.0
  [d6f4376e] Markdown v1.11.0
  [de0858da] Printf v1.11.0
  [9a3f8284] Random v1.11.0
  [ea8e919c] SHA v0.7.0
  [9e88b42a] Serialization v1.11.0
  [6462fe0b] Sockets v1.11.0
  [fa267f1f] TOML v1.0.3
  [cf7118a7] UUIDs v1.11.0
  [4ec0a83e] Unicode v1.11.0
  [e66e0078] CompilerSupportLibraries_jll v1.1.1+0
  [4536629a] OpenBLAS_jll v0.3.27+1
  [8e850b90] libblastrampoline_jll v5.11.0+0
  • Output of versioninfo()
Julia Version 1.11.5
Commit 760b2e5b739 (2025-04-14 06:53 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 12 × 13th Gen Intel(R) Core(TM) i5-1345U
  WORD_SIZE: 64
  LLVM: libLLVM-16.0.6 (ORCJIT, goldmont)
Threads: 1 default, 0 interactive, 1 GC (on 12 virtual cores)

Additional context

I also tried other integrators, which have similar, yet slightly different issues. E.g., with

using OrdinaryDiffEqLowStorageRK
using RecursiveArrayTools: ArrayPartition, VectorOfArray

function rhs!(dω, ω, p, t)
    du, dv =.x
    du .= 1.0
    dv .= 1.0
end

u_ini = VectorOfArray([[1.0, 2.0], [3.0, 4.0]])
v_ini = VectorOfArray([[1.0, 2.0], [3.0, 4.0]])
# u_ini = [1.0, 2.0, 3.0, 4.0] # This works
# v_ini = [3.0, 4.0, 3.0, 4.0]
ω_ini = ArrayPartition(u_ini, v_ini)

tspan = (0.0, 1.0)
ode = ODEProblem(rhs!, ω_ini, tspan)
sol = solve(ode, RDPK3SpFSAL49(), dt = 0.2, adaptive = false)

I get

retcode: Success
Interpolation: 3rd order Hermite
t: 6-element Vector{Float64}:
 0.0
 0.2
 0.4
 0.6000000000000001
 0.8
 1.0
u: 6-element Vector{ArrayPartition{Float64, Tuple{VectorOfArray{Float64, 2, Vector{Vector{Float64}}}, VectorOfArray{Float64, 2, Vector{Vector{Float64}}}}}}:
 (VectorOfArray{Float64, 2, Vector{Vector{Float64}}}([[1.8980136154766194, 2.8980136154766165], [3.898013615476608, 4.8980136154766125]]), VectorOfArray{Float64, 2, Vector{Vector{Float64}}}([[1.8980136154766194, 2.8980136154766165], [3.898013615476608, 4.8980136154766125]]))
 (VectorOfArray{Float64, 2, Vector{Vector{Float64}}}([[1.8980136154766194, 2.8980136154766165], [3.898013615476608, 4.8980136154766125]]), VectorOfArray{Float64, 2, Vector{Vector{Float64}}}([[1.8980136154766194, 2.8980136154766165], [3.898013615476608, 4.8980136154766125]]))
 (VectorOfArray{Float64, 2, Vector{Vector{Float64}}}([[1.8980136154766194, 2.8980136154766165], [3.898013615476608, 4.8980136154766125]]), VectorOfArray{Float64, 2, Vector{Vector{Float64}}}([[1.8980136154766194, 2.8980136154766165], [3.898013615476608, 4.8980136154766125]]))
 (VectorOfArray{Float64, 2, Vector{Vector{Float64}}}([[1.8980136154766194, 2.8980136154766165], [3.898013615476608, 4.8980136154766125]]), VectorOfArray{Float64, 2, Vector{Vector{Float64}}}([[1.8980136154766194, 2.8980136154766165], [3.898013615476608, 4.8980136154766125]]))
 (VectorOfArray{Float64, 2, Vector{Vector{Float64}}}([[1.8980136154766194, 2.8980136154766165], [3.898013615476608, 4.8980136154766125]]), VectorOfArray{Float64, 2, Vector{Vector{Float64}}}([[1.8980136154766194, 2.8980136154766165], [3.898013615476608, 4.8980136154766125]]))
 (VectorOfArray{Float64, 2, Vector{Vector{Float64}}}([[1.8980136154766194, 2.8980136154766165], [3.898013615476608, 4.8980136154766125]]), VectorOfArray{Float64, 2, Vector{Vector{Float64}}}([[1.8980136154766194, 2.8980136154766165], [3.898013615476608, 4.8980136154766125]]))

Again, every time step is the same, but this time it's also not giving the correct final solution.
Using adaptive time stepping instead gives similar results.

Metadata

Metadata

Assignees

Labels

Type

No type

Projects

No projects

Milestone

No milestone

Relationships

None yet

Development

No branches or pull requests

Issue actions