Open
Description
Environment
Which version of REBOUND are you using and on what operating system?
- REBOUND Version: 4.4.5
- API interface: Python
- Operating System (including version): Ubuntu 22.04.5 LTS
Describe the bug
I'm simulating a self-gravitating disc of particles using tree gravity and leapfrog. The simulation runs but afterwards, if I try to load snapshots towards the end, I get a seg fault. The behaviour also appears unphysical throughout the simulation.
To Reproduce
Minimum working example, which runs the sim then afterwards tries to load the final snapshot and seg faults:
'''Minimum working example'''
############################### Libraries ###############################
import rebound
import numpy as np
import time as systime
############################## User inputs ##############################
# Star parameters
mStar_mSun = 1.92
# Disc parameters
mDisc_mSun = 3.000e-2 * 3e-6
aIn_au = 130.
aOut_au = 140.
rSDIndex = 1.5
nDebris = 30
# Sim parameters
endTime_yr = 1e5
simBoxSize_au = 3000
numberOfSnapshots = 10
integratorName = 'leapfrog'
binFilePath = 'simTest.bin'
################################ Program ################################
print('\nSetting up sim...')
# Set up new simulation
sim = rebound.Simulation()
sim.integrator = integratorName
sim.units = ('yr', 'AU', 'Msun')
# Add the star
sim.add(m = mStar_mSun, hash='Star')
# Calculate the number of debris bodies
mDebris_mSun = mDisc_mSun / nDebris
# Semimajor axis distribution index
aIndex = rSDIndex - 1
# Add the debris
for debrisIndex in range(nDebris):
# Add the particle
sim.add(m = mDebris_mSun,
a = (np.random.rand()*(aOut_au**(1-aIndex) - aIn_au**(1-aIndex)) + aIn_au**(1-aIndex))**(1.0/(1-aIndex)),
e = np.random.uniform(0, 0.2),
inc = np.random.uniform(0, 0.2),
Omega = np.random.uniform(0, 2*np.pi),
omega = np.random.uniform(0, 2*np.pi),
M = np.random.uniform(0, 2*np.pi),
hash = 'Debris%s' % (debrisIndex+1))
# Set number of active (big) particles
sim.N_active = 1 + nDebris
# Set the box size and boundary conditions
sim.configure_box(simBoxSize_au)
sim.boundary = 'none'
sim.gravity = 'tree'
# Set the collision prescription
sim.collision = 'linetree'
sim.collision_resolve = 'hardsphere'
# Set timesteps from the smallest orbital period
minOrbitalPeriod_yr = 1e99
for particleIndex, particle in enumerate(sim.particles):
# Omit the primary, since has no orbit
if particleIndex != 0:
minOrbitalPeriod_yr = min(particle.P, minOrbitalPeriod_yr)
sim.dt = 1e-3 * minOrbitalPeriod_yr
# Set the start time
startTime_yr = 0
# Set the output times (log scale)
times_yr = np.geomspace(1e-3*minOrbitalPeriod_yr, endTime_yr, numberOfSnapshots-1)
# Run the simulaton
outputNumber = 0
pcStartTime_s = systime.time()
print('\nStarting sim...')
for i, time_yr in enumerate(times_yr):
outputNumber += 1
sim.integrate(time_yr)
# Save the simulation
sim.save_to_file(binFilePath)
# Print run time and number of bodies remaining
runTimeSoFar_s = systime.time() - pcStartTime_s
print(' %.2E yr (step %s/%s) complete. Runtime: %.2E sec. Bodies remaining: %s/%s' % (sim.t, outputNumber, len(times_yr), runTimeSoFar_s, sim.N, 1+nDebris))
print('\nSimulation complete. Took %.2E seconds\n' % runTimeSoFar_s)
# Try to load the final snapshot
print('Loading sim archive...')
simArchive = rebound.Simulationarchive(binFilePath)
print('Accessing final snapshot...')
finalSnapshot = simArchive[-1]
print('Complete')
########################################################################
Additional context
I am unsure what box size is recommended for tree, but this does not affect the issue.
Name/Affiliation
Tim Pearce, University of Warwick