Here we show how to add force steps before and after a REBOUND timestep. This can be important if the force is velocity-dependent and you are using WHFast, since it assumes that additional forces are only position-dependent (see Sec. 5.1 of Tamayo et al. 2019, where an ultra-short period planet is spuriously ejected after 3 million years by velocity-dependent GR forces). Velocity-dependent forces should be integrated across the timestep, especially if they're non-dissipative, like GR.
In [1]:
import rebound
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
def system():
sim = rebound.Simulation()
sim.G = 4*np.pi**2
sim.add(m=0.93)
sim.add(m=4.5*3.e-7, P=0.571/365.25, e=0.01)
sim.add(m=41.*3.e-7, P=13.34/365.25, e=0.01)
sim.move_to_com()
sim.dt = 0.07*sim.particles[1].P
return sim
If we simply add the GR force and integrate as normal with WHFast, the energy error will grow secularly:
In [2]:
sim = system()
sim.integrator = "whfast"
In [3]:
import reboundx
rebx = reboundx.Extras(sim)
gr = rebx.load_force("gr")
rebx.add_force(gr)
gr.params["c"] = 63197.8# AU/yr
In [4]:
Nout = 1000
times = np.linspace(0, 1e4*sim.particles[1].P, Nout)
E0 = rebx.gr_hamiltonian(gr)
Eerr = np.zeros(Nout)
for i, time in enumerate(times):
sim.integrate(time, exact_finish_time=0)
E = rebx.gr_hamiltonian(gr)
Eerr[i] = np.abs((E-E0)/E0)
In [5]:
fig, ax = plt.subplots(figsize=(9,6))
ax.loglog(times/sim.particles[1].P, Eerr, '.')
ax.set_xlabel('Time (Inner Planet Orbits)', fontsize=24)
ax.set_ylabel('Relative Energy Error', fontsize=24)
Out[5]:
If instead we add an operator step before and after the timestep that integrates the GR force across it, we fix the problem. We again load the gr force, but instead of adding it to rebx, we instead add an 'integrate_force' operator, and add the gr force as a parameter to that operator step:
In [6]:
sim = system()
sim.integrator = "whfast"
rebx = reboundx.Extras(sim)
gr = rebx.load_force("gr")
gr.params["c"] = 63197.8# AU/yr
intforce = rebx.load_operator('integrate_force')
rebx.add_operator(intforce)
intforce.params['force'] = gr
intforce.params['integrator'] = reboundx.integrators['implicit_midpoint']
Above we have also chosen which integrator to use across the timestep, here implicit midpoint. See the REBOUNDx paper for more discussion
In [7]:
Nout = 1000
times = np.logspace(0, np.log10(1e4*sim.particles[1].P), Nout)
E0 = rebx.gr_hamiltonian(gr)
Eerr = np.zeros(Nout)
for i, time in enumerate(times):
sim.integrate(time, exact_finish_time=0)
E = rebx.gr_hamiltonian(gr)
Eerr[i] = np.abs((E-E0)/E0)
In [8]:
fig, ax = plt.subplots(figsize=(9,6))
ax.loglog(times/sim.particles[1].P, Eerr, '.')
ax.set_xlabel('Time (Inner Planet Orbits)', fontsize=24)
ax.set_ylabel('Relative Energy Error', fontsize=24)
Out[8]:
In [ ]: