dtamayo / reboundx

A library for adding additional forces to the REBOUND N-body integration package
GNU General Public License v3.0
80 stars 61 forks source link

safe_mode #36

Closed hannorein closed 5 years ago

hannorein commented 5 years ago

It seems that REBOUNDx overwrites the safe_mode flag in whfast. This leads to a significant slowdown in simulations where I add a GR potential and use symplectic correctors. For example:

import rebound
import reboundx
from rebound.data import add_outer_solar_system

sim = rebound.Simulation()
add_outer_solar_system(sim)
sim.integrator = "whfast"
sim.ri_whfast.corrector = 7
sim.ri_whfast.safe_mode = 0

rebx = reboundx.Extras(sim)
gr = rebx.add("gr_potential")
gr.params['c'] = constants.C

sim.integrate(100)

I should probably remember this, but is this by design or a bug?

hannorein commented 5 years ago

I think I narrowed down the problem.

When calling rebx = reboundx.Extras(sim), not only additional forces, but also function pointers for pre and post timestep modifications get set to the REBOUNDx wrappers. In the above examples, those remain dummies which don't do anything. However, in REBOUND, it will force the integrator to synchronize every timestep. Again, not sure if this is good or bad!

dtamayo commented 5 years ago

Thanks for looking into this. I'm not sure what the best solution is, but this is a clear use case where the current implementation is unusable and you'd just copy-paste your own gr_potential routine.

My hesitation before was that you don't know what effects the user will add ahead of time, but maybe a better default is that whenever a pre or post timestep modification is added in REBOUNDx, it overwrites the REBOUND function pointer to the generic REBOUNDx wrapper. I don't see any obvious downsides to this.

dtamayo commented 5 years ago

I see now that this isn't straightforward in the previous version, because whether a force gets executed as a force or as an additional operator can be changed by the user at run time (there was a force_as_operator field in each effect). The fact that this is straightforward in the new version I think means that the structural changes to make everything more explicit were done well :) So I'll patch this in the lookuptable branch for now.

hannorein commented 5 years ago

Thanks! I agree, I don't see a downside to setting the function pointers later. And in cases where they are getting set, we want to synchronize. It might not be necessary in every case, but the additional complexity to doing something smarter would be too high. I hacked a together a workaround for the simulations I'm currently running, so no rush to do anything!