mercurius testparticles strange behavior

[r-zolotarev]

Hello!

I integrate problem with massless testpartcles and find strange behavior of particles. When IAS15 or WHFast is used there is no significant difference in results with different testparticle_type, when MERCURIUS is used the difference occurs. The masses for testparticles is set to m=0. What this difference may be related with?

Here is an example. The system consist of star with m=1, Jupiter-mass planet and testparticle with m=0.

import matplotlib.pyplot as plt
import numpy as np
import rebound

# system with star, planet (Jupiter) and massless test particle

M = 1
m = 1.e-3
r = 5
R = r*m/M
G = 1
v = np.sqrt(G*(M+m)/r)
V = v*m/M
T = 2.*np.pi*r/v

# test particle
r1 = 4
v1 = np.sqrt(G*(M+m)/r1)

timestep = 2*np.pi/365.25 * 3  # 3 days

sim_ias15 = rebound.Simulation()
sim_mercurius = rebound.Simulation()
sim_whfast = rebound.Simulation()

sim_ias15.integrator = 'ias15'
sim_mercurius.integrator = 'mercurius'
sim_whfast.integrator = 'whfast'

for sim in sim_ias15, sim_mercurius, sim_whfast:
    sim.add(m=M, x=-R, y=0, z=0, vx=0, vy=-V, vz=0)
    sim.add(m=m, x=r, y=0, z=0, vx=0, vy=v, vz=0)
    sim.add(m=0, x=0, y=r1, z=0, vx=-v1, vy=0, vz=0)
    sim.exact_finish_time = 1
    # sim.testparticle_type = 1
    sim.N_active = 2
    sim.move_to_com()

sim_ias15.testparticle_type = 1
sim_mercurius.testparticle_type = 0  # !!!
sim_whfast.testparticle_type = 0

sim_mercurius.dt = timestep
sim_whfast.dt = timestep

sim_mercurius.ri_mercurius.hillfac = 3.

dt_out = T/100
t_tot = 5*T
t_int = np.arange(0, t_tot+dt_out/2, dt_out)

t_ias15, t_mercurius, t_whfast = [], [], []
a_planet = []
a_testparticle = []

for i in range(0, len(t_int)):
    orb = []
    for sim in sim_ias15, sim_mercurius, sim_whfast:
        sim.integrate(t_int[i])
        orb.append(sim.calculate_orbits(primary=sim.particles[0]))

    a_planet.append([orb[0][0].a, orb[1][0].a, orb[2][0].a])
    a_testparticle.append([orb[0][1].a, orb[1][1].a, orb[2][1].a])

    t_ias15.append(sim_ias15.t)
    t_mercurius.append(sim_mercurius.t)
    t_whfast.append(sim_whfast.t)

t_ias15 = np.asarray(t_ias15)
t_mercurius = np.asarray(t_mercurius)
t_whfast = np.asarray(t_whfast)
a_planet = np.asarray(a_planet)
a_testparticle = np.asarray(a_testparticle)

fig, ax = plt.subplots(1, 2)
fig.suptitle(f'sim_mercurius.testparticle_type = {sim_mercurius.testparticle_type}\n'
             f'sim_whfast.testparticle_type = {sim_whfast.testparticle_type}\n'
             f'sim_ias15.testparticle_type = {sim_ias15.testparticle_type}')

ax[0].set_title('planet')
ax[0].plot(t_ias15/T, a_planet[:, 0], c='orange', lw=3, label='ias15')
ax[0].plot(t_mercurius/T, a_planet[:, 1], c='red', lw=2, label='mercurius')
ax[0].plot(t_whfast/T, a_planet[:, 2], c='blue', lw=1, label='whfast')
ax[0].set_xlabel('t/T')
ax[0].set_ylabel('a (AU)')
ax[0].legend()

ax[1].set_title('test particle')
ax[1].plot(t_ias15/T, a_testparticle[:, 0], c='orange', lw=3, label='ias15')
ax[1].plot(t_mercurius/T, a_testparticle[:, 1], c='red', lw=2, label='mercurius')
ax[1].plot(t_whfast/T, a_testparticle[:, 2], c='blue', lw=1, label='whfast')
ax[1].set_xlabel('t/T')
ax[1].set_ylabel('a (AU)')
ax[1].legend()

plt.tight_layout()
plt.show()

REBOUND version: 3.18.1 (installed via pip (Ubuntu 21.04))

Sorry if I miss something important and my question is incorrect

Thanks, Roman

[hannorein]

Hi Roman,

Thanks for the great example to reproduce the issue. I agree, something is weird. I'll investigate.

Hanno

[hannorein]

I think this was indeed a bug. The jump step was not executed for test particles. It should be fixed in version 2.19.2 (648c2a73987712345ce68de4a5bf8ce0a6db3490). Thanks for spotting this.

[r-zolotarev]

Thanks a lot!