Segfault in interpolation

[antoine-levitt]

I get a segfault when getting results out using VelocityVerlet, in a getindex. I'm guessing overagressive inbounds somewhere? I'm too lazy to do an MWE but the below should reproduce the bug. Changing the integrator to something else fixes it. I'm filing against NBS since that's what I'm using, but I can file to DifferentialEquations.jl if that is more appropriate.

using SPICE
using Downloads: download
using NBodySimulator, StaticArrays

# units: kilometers, seconds, kilograms

const LSK = "https://naif.jpl.nasa.gov/pub/naif/generic_kernels/lsk/naif0012.tls"
const SPK = "https://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/planets/de440.bsp"
# Download kernels
# time kernel
download(LSK, "naif0012.tls")
download(SPK, "de440.bsp")
download("https://naif.jpl.nasa.gov/pub/naif/generic_kernels/pck/Gravity.tpc", "Gravity.tpc")

furnsh("naif0012.tls")
# ephemeris kernel
furnsh("de440.bsp")
# masses kernel
furnsh("Gravity.tpc")

# Convert the calendar date to ephemeris seconds past J2000
start_date = str2et("Jan 1, 2020")

years_to_observe = 1
dates = collect(range(start_date, start_date + 60*60*24*365*years_to_observe, length=365*years_to_observe))

ref_pos = spkpos.("earth", dates, "J2000", "none", "sun")
ref_pos = [p[1] for p in ref_pos]
ref_pos = reduce(hcat, ref_pos)'

const grav_constant = 6.6743e-11 / 1e9 # m^3 kg^-1 s^-2 -> km^3

planets = ["earth_barycenter", "sun"]
# planets = ["earth_barycenter", "Mercury_barycenter", "Venus_barycenter", "Mars_barycenter", "Jupiter_barycenter", "Saturn_barycenter", "Uranus_barycenter", "Neptune_barycenter", "Sun"]
# planets = ("earth_barycenter", "Mercury_barycenter", "Venus_barycenter", "Mars_barycenter", "Jupiter_barycenter", "Saturn_barycenter", "Uranus_barycenter", "Sun")
planets_to_observe = [1]

masses = bodvrd.(planets, "GM") ./ grav_constant
posvel = [spkezr(p, start_date, "J2000", "none", "sun") for p in planets]

bodies = [MassBody(SVector(posvel[i][1][1:3]...), SVector(posvel[i][1][4:6]...), masses[i][1]) for i = 1:length(planets)]
system = GravitationalSystem(bodies, grav_constant)
tspan = (dates[1], dates[end])
simulation = NBodySimulation(system, tspan)
result = run_simulation(simulation, VelocityVerlet(), dt=60*60*24).solution
pos = [result(t)[:, planets_to_observe] .- result(t)[:, end] for t in dates] # observe wrt sun
pos = cat(pos...; dims=3) # α x p x t

function get_ref_pos(t)
    hcat([spkpos(p, t, "J2000", "none", "sun")[1] for p in planets[planets_to_observe]]...)
end
ref_pos = cat(get_ref_pos.(dates)...; dims=3)

println(maximum(abs, pos-ref_pos))

# using PyPlot
# plot3D(pos[1, 1, :], pos[2, 1, :], pos[3, 1, :])

[ChrisRackauckas]

If you run this with inbounds checking forced does it still segfault? I can't run this without the data files.

[antoine-levitt]

The data files should be downloaded automatically by the script

[oscardssmith]

For me this runs properly until the pos = [result(t)[:, planets_to_observe] .- result(t)[:, end] for t in dates] # observe wrt sun line

[antoine-levitt]

Good debugging tips, check-bounds=yes gets me

antoine@epsilon:~$ julia --check-bounds=yes scratch.jl 
ERROR: LoadError: BoundsError: attempt to access Tuple{Matrix{Float64}, Matrix{Float64}} at index [3]
Stacktrace:
  [1] getindex
    @ ~/.julia/packages/RecursiveArrayTools/gH9Kb/src/array_partition.jl:216 [inlined]
  [2] macro expansion
    @ ./multidimensional.jl:903 [inlined]
  [3] macro expansion
    @ ./cartesian.jl:64 [inlined]
  [4] _unsafe_getindex!
    @ ./multidimensional.jl:898 [inlined]
  [5] _unsafe_getindex(::IndexCartesian, ::ArrayPartition{Float64, Tuple{Matrix{Float64}, Matrix{Float64}}}, ::Base.Slice{Base.OneTo{Int64}}, ::Vector{Int64})
    @ RecursiveArrayTools ~/.julia/packages/RecursiveArrayTools/gH9Kb/src/array_partition.jl:245
  [6] _getindex
    @ ./multidimensional.jl:875 [inlined]
  [7] getindex(::ArrayPartition{Float64, Tuple{Matrix{Float64}, Matrix{Float64}}}, ::Function, ::Vector{Int64})
    @ Base ./abstractarray.jl:1241
  [8] (::var"#7#8")(t::Float64)
    @ Main ./none:0
  [9] iterate
    @ ./generator.jl:47 [inlined]
 [10] collect(itr::Base.Generator{Vector{Float64}, var"#7#8"})
    @ Base ./array.jl:787
 [11] top-level scope
    @ ~/scratch.jl:46
in expression starting at /home/antoine/scratch.jl:46

[ChrisRackauckas]

Latest versions? Show ]st on Julia v1.8?

[antoine-levitt]

(@v1.8) pkg> st
Status `~/.julia/environments/v1.8/Project.toml`
  [ac637c84] AbbreviatedStackTraces v0.1.5
  [6e4b80f9] BenchmarkTools v1.3.2
  [acf6eb54] DFTK v0.5.15 `~/.julia/dev/DFTK`
  [56d4f2e9] Gridap v0.17.15 `~/.julia/dev/Gridap`
  [3025c34a] GridapGmsh v0.6.1 `~/.julia/dev/GridapGmsh`
  [7073ff75] IJulia v1.23.3
  [5903a43b] Infiltrator v1.6.3
⌃ [5078a376] LazyArrays v0.22.6
  [0e6f8da7] NBodySimulator v1.9.0
  [d330b81b] PyPlot v2.11.0
  [295af30f] Revise v3.4.0
  [5bab7191] SPICE v0.2.3
  [90137ffa] StaticArrays v1.5.11
  [64499a7a] WriteVTK v1.16.0
  [f43a241f] Downloads v1.6.0
Info Packages marked with ⌃ have new versions available and may be upgradable.

[antoine-levitt]

(julia 1.8.3)

[ChrisRackauckas]

Line 46 is after the integration, so it's not to do with the integration. Something is odd with the interpolation after the solution. Are the dates within the bounds of the solution? I can dig in in a bit, just running out of computers right now.

[antoine-levitt]

Yes, it's the interpolation.

tspan = (dates[1], dates[end])
simulation = NBodySimulation(system, tspan)
result = run_simulation(simulation, VelocityVerlet(), dt=60*60*24).solution
pos = [result(t)[:, planets_to_observe] .- result(t)[:, end] for t in dates] # observe wrt sun

is the relevant part, dates should definitely be in range.

[ChrisRackauckas]

The interpolation is fine, but ArrayPartitions are AbstractVector so I'm not sure what the matrix indexing is doing here. Is this what you meant?

julia> result(0.0)
([246892.39098961183 -0.7507504327994997; 44191.712533044185 -0.1343762957324868; 19156.371340334987 -0.058249886089935206], [-1.560440141564781e14 4.744985087272809e8; -2.6788610085964547e13 8.14584384415191e7; -1.1612398946939709e13 3.531082356110646e7])

julia> result(dates[1])
([-29.84582650955506 0.0; -4.726462025631025 0.0; -2.0488332065817096 0.0], [-2.4882295399688877e7 0.0; 1.3301623108934854e8 0.0; 5.766241108251575e7 0.0])

julia> result(dates[1])[:, 1]
ERROR: BoundsError: attempt to access Tuple{Matrix{Float64}, Matrix{Float64}} at index [3]
Stacktrace:
 [1] getindex
   @ C:\Users\accou\.julia\packages\RecursiveArrayTools\gH9Kb\src\array_partition.jl:216 [inlined]
 [2] macro expansion
   @ .\multidimensional.jl:903 [inlined]
 [3] macro expansion
   @ .\cartesian.jl:64 [inlined]
 [4] _unsafe_getindex!
   @ .\multidimensional.jl:898 [inlined]
 [5] _unsafe_getindex(::IndexCartesian, ::ArrayPartition{Float64, Tuple{Matrix{Float64}, Matrix{Float64}}}, ::Base.Slice{Base.OneTo{Int64}}, ::Int64)
   @ RecursiveArrayTools C:\Users\accou\.julia\packages\RecursiveArrayTools\gH9Kb\src\array_partition.jl:245
 [6] _getindex
   @ .\multidimensional.jl:875 [inlined]
 [7] getindex(::ArrayPartition{Float64, Tuple{Matrix{Float64}, Matrix{Float64}}}, ::Function, ::Int64)
   @ Base .\abstractarray.jl:1241
 [8] top-level scope
   @ REPL[47]:1

julia> result(dates[1])[1]
-29.84582650955506

julia> pos = [result(t)[planets_to_observe] .- result(t)[end] for t in dates] # observe wrt sun      
365-element Vector{Vector{Float64}}:
 [-29.84582650955506]
 [-29.7787085265913]
 [-29.756971049702006]
 [-29.780454567120422]
 [-29.848984986504647]
 [-29.96237369489716]
 [-30.120417624435966]
 [-30.32289932378238]
 ⋮
 [165.28898095767556]
 [165.64466104997476]
 [165.9542939696105]
 [166.21794639967436]
 [166.43569984785287]
 [166.607650623168]
 [166.73390980684837]

[antoine-levitt]

This is weird, because the same code with this works: result = run_simulation(simulation, alg=Vern9(), abstol=tol, reltol=tol).solution. Anyway, it should not segfault even if invalid.

[ChrisRackauckas]

The segfault is fixed in SciML/RecursiveArrayTools.jl#239 .