We can reproduce REBOUNDx simulations with extra effects bit by bit with the SimulationArchive (see Rein & Tamayo 2017 for details), and make sure to read SimulationArchive.ipynb and SimulationArchiveRestart.ipynb first.
However, this will only work under some conditions. In particular, for bit-by-bit reproducibility, currently one requires that:
REBOUNDx effects in the simulation are machine independent (can't have, e.g., trig functions, pow or exp in implementation)To use the SimulationArchive with REBOUNDx, we need to save a REBOUNDx binary SavingAndLoadingSimulations.ipynb. Since we can't change effects or particle parameters, it doesn't matter at what point we save this binary (we will just need it to load the SimulationArchive). Here we do a WHFAST integration with the symplectic gr_potential implementation for general relativity corrections. We set up the simulation like we usually would:
In [1]:
import rebound
import reboundx
from reboundx import constants
sim = rebound.Simulation()
sim.add(m=1.)
sim.add(m=1e-3, a=1., e=0.2)
sim.add(m=1e-3, a=1.9)
sim.move_to_com()
sim.dt = sim.particles[1].P*0.05 # timestep is 5% of orbital period
sim.integrator = "whfast"
rebx = reboundx.Extras(sim)
gr = rebx.load_force("gr_potential")
rebx.add_force(gr)
gr.params["c"] = constants.C
rebx.save("rebxarchive.bin")
We now set up the SimulationArchive and integrate like we normally would (SimulationArchive.ipynb):
In [2]:
sim.automateSimulationArchive("archive.bin", interval=1e3, deletefile=True)
sim.integrate(1.e6)
Once we're ready to inspect our simulation, we use the reboundx.SimulationArchive wrapper that additionally takes a REBOUNDx binary:
In [3]:
sa = reboundx.SimulationArchive("archive.bin", rebxfilename = "rebxarchive.bin")
We now create a different simulation from a snapshot in the SimulationArchive halfway through:
In [4]:
sim2, rebx = sa[500]
sim2.t
Out[4]:
We now integrate our loaded simulation to the same time as above (1.e6):
In [5]:
sim2.integrate(1.e6)
and see that we obtain exactly the same particle positions in the original and reloaded simulations:
In [6]:
sim.status()
In [7]:
sim2.status()