Particles in system stop moving using Mercurius integrator
alelatt opened this issue · 8 comments
Environment
Which version of REBOUND are you using and on what operating system?
- REBOUND Version: 4.3.1-4.3.2
- API interface: Python
- Operating System (including version): Linux on WSL
Describe the bug
During a simulation using Mercurius integrator two of the planets in the system get "stuck" in place.
To Reproduce
The system was created as part of a test to learn how to use parallelization and Simulation Archives using
def Simulation(niter):
print(niter)
t_init = time.time()
sim = rebound.Simulation()
sim.units = ['mearth', 'year', 'AU']
sim.integrator = "mercurius"
sim.ri_mercurius.r_crit_hill = 3.
sim.ri_mercurius.L = "infinity"
sim.add(m = 333000)
np.random.seed(niter)
sim.add(m = 1, a = abs(np.random.normal(1, 0.2)), e = abs(np.random.normal(0.1, 0.05)), Omega = np.random.uniform(0, 2*np.pi))
sim.add(m = 1, a = abs(np.random.normal(2, 0.2)), e = abs(np.random.normal(0.1, 0.05)), Omega = np.random.uniform(0, 2*np.pi))
sim.add(m = 1, a = abs(np.random.normal(3, 0.2)), e = abs(np.random.normal(0.1, 0.05)), Omega = np.random.uniform(0, 2*np.pi))
sim.add(m = 1, a = abs(np.random.normal(4, 0.2)), e = abs(np.random.normal(0.1, 0.05)), Omega = np.random.uniform(0, 2*np.pi))
sim.add(m = 1, a = abs(np.random.normal(5, 0.2)), e = abs(np.random.normal(0.1, 0.05)), Omega = np.random.uniform(0, 2*np.pi))
sim.move_to_com()
times = np.linspace(0, -1e6, 100)
for timestep in times:
sim.integrate(timestep, exact_finish_time=0)
sim.save_to_file("./simarchive/archive{0}.bin".format(niter))
if np.any(Compute_Energy(sim) > 0):
return [1, time.time() - t_init, sim.t]
return [0, time.time() - t_init, sim.t]
if __name__ == '__main__':
Nsym = 40
pool = Pool()
results = np.array(pool.map(Simulation, range(Nsym)))
where the simulation stops after some time or if one of the objects becomes unbound at the end of a step (checked using the function "ComputeEnergy").
The systems that have objects that become unbound all present the same problem: after some time two planets that are next to each other get "stuck" in place whilst the rest of the planets show no abnormal behaviour.
I attach the bin file for the archive of one of the systems that got stuck (before t = -40404) and to which the images are referred to
problem.zip
The planets that get stuck stop moving even though their velocities seem to grow (thus the energy grows over zero and the simulation stops) as shown in the plots for one of such systems
Reducing or enlarging the timestep using sim.dt doesn't change the outcome.
Additional context
I first downloaded and started using Rebound a couple of weeks ago to use it for my Master's thesis and chose to use Mercurius as the integrator.
When doing the first tests with mercurius I got the error
File "/home/alelatt/Thesis_2024/Test_Examples/test_expulsion.py", line 40, in
sim.ri_mercurius.L = "infinity"
File "/home/alelatt/.local/lib/python3.10/site-packages/rebound/integrators/mercurius.py", line 59, in L
self._L = cast(clibrebound.reb_integrator_mercurius_L_minfinity,MERCURIUSLF)
NameError: name 'cast' is not defined
which I solved by changing
def L(self, func):
if func == "mercury":
self._L = cast(clibrebound.reb_integrator_mercurius_L_mercury,MERCURIUSLF)
elif func == "C4":
self._L = cast(clibrebound.reb_integrator_mercurius_L_C4,MERCURIUSLF)
elif func == "C5":
self._L = cast(clibrebound.reb_integrator_mercurius_L_C5,MERCURIUSLF)
elif func == "infinity":
self._L = cast(clibrebound.reb_integrator_mercurius_L_minfinity,MERCURIUSLF)
to
def L(self, func):
if func == "mercury":
self._L = ctypes.cast(clibrebound.reb_integrator_mercurius_L_mercury,MERCURIUSLF)
elif func == "C4":
self._L = ctypes.cast(clibrebound.reb_integrator_mercurius_L_C4,MERCURIUSLF)
elif func == "C5":
self._L = ctypes.cast(clibrebound.reb_integrator_mercurius_L_C5,MERCURIUSLF)
elif func == "infinity":
self._L = ctypes.cast(clibrebound.reb_integrator_mercurius_L_infinity,MERCURIUSLF)
in "mercurius.py" and got the first basic snippets of code to work .
I don't know if this is a separate bug (in which case tell me if I should open another bug report) or if this is the source of my problems.
The set of simulations was first computed on a pc using rebound 4.3.2. The same result is obtained on a pc using rebound 4.3.1.
A test was conducted on the above simulation by rescaling the semimajor axes of all planets by a factor 1/10 and leaving the same other starting conditions, giving the exact same end result.
Another test was conducted by removing from the above system one by one the planets that were not effected by the issue, but this didn't solve the issue.
It's a bit hard to disentangle what's going on here. Can you provide a set of initial conditions that:
- does not depend on random numbers
- leads to this issue within a short amount of time (a few seconds CPU time)
That way I could reproduce it and look into it.
The other issue with ctypes.cast is a bug not related to this, but thanks for letting me know.
I've found one example which solves in a few seconds
sim = rebound.Simulation()
sim.units = ['mearth', 'year', 'AU']
sim.add(m=333000)
sim.add(m = 1, a = 1.0882454973770086, e = 0.08345649240529567, Omega = 1.2988547595412578)
sim.add(m = 1, a = 1.9953392507100292, e = 0.14208888335059075, Omega = 3.843700051147186)
sim.add(m = 1, a = 3.0219219683156386, e = 0.17912405585307878, Omega = 1.8648525506672495)
sim.add(m = 1, a = 3.919381337633177, e = 0.18188142448613187, Omega = 0.9946902347936728)
sim.add(m = 1, a = 4.928234210599751, e = 0.13017358013047473, Omega = 2.602715385609317)
sim.integrator = "mercurius"
sim.ri_mercurius.r_crit_hill = 3.
sim.ri_mercurius.L = "infinity"
sim.move_to_com()
sim.integrate(-700)
in which the problem appears between -600 and -700 years.
The orbital elements are taken from those of one of the failed simulations created using random values, this is why they are so precise.
Perfect! Thank you. I'll look into it...
I can reproduce the issue. I suspect this might be bug that occurs when you're integrating backwards in time. Still need to fix it...
I think this should be fixed in 7a6f596. Please let me know if this fixes all the issues for you. And thanks again for reporting it and especially for providing such a nice example that allowed me to quickly reproduce the issue!
Thank you very much for the quick fix! I will check and get back to you as soon as possible.
I checked all the systems that showed the problem and it seems to be solved. Thanks again!
Glad to hear. And thanks again for reporting this issue!