test particles influencing active particles during integration

[ryurongliu]

Hi Hanno,

I'm running into a strange issue while simulating test particle-active particle collisions in the TRAPPIST-1 system. If I integrate the planets without any test particles, I get a reproducible final planet configuration after 1000 years, no problem. When I put massless test particles into the system, they somehow influence the orbit of the planets enough to reach a very different (but also reproducible) final planet configuration after 1000 years. I've checked my sim's parameters for N_active = num_planets+1, test_particle_type = 0, and all the masses of the test particles = 0, and I've run multiple tests with varying collision resolution types, and I keep getting very different final planet configurations after integrating for 1000 years. Is there any reason why test particles might influence the behavior of active particles like this?

Thanks,
Raina

[hannorein]

Hi Raina,

I need a bit more information to give you an answer (either a short piece of code which reproduces the effect, or details about the integrator you use, collision model, etc).

There can be various small differences at the floating point level. They should not have an effect on the physical outcomes, but if the system is chaotic then an individual trajectory might be significantly altered.

Of course, there is always the possibility that this could be a bug, or something might be configured incorrectly.

Hanno

[ryurongliu]

Hi Hanno,

Thanks for the quick reply. I am setting up my system like so:

sim = rebound.Simulation()


sim.units = ('yr', 'AU', 'Msun')
sim.G = 39.441677
sim.integrator = "IAS15"
sim.collision = "none"
sim.dt = (1/365.25)/10

I then added my star and 7 planets by orbital elements:

sim.add(m = 0.0898, r = 0.1192*Rsun_to_AU, hash = 'a')
sim.add(m = 1.374*Mearth_to_Msun, a=1.154e-2, e=0, inc=(90-89.728)*deg_to_rad, 
            Omega=1, omega=1, M=320.2328615005874, r = 1.116*Rearth_to_AU, hash='b')
sim.add(m = 1.308*Mearth_to_Msun, a=1.580e-2, e=0, inc=(90-89.778)*deg_to_rad,
            Omega=1, omega=1, M=302.03085481266083, r = 1.097*Rearth_to_AU, hash='c')
sim.add(m = 0.388*Mearth_to_Msun, a=2.227e-2, e=0, inc=(90-89.896)*deg_to_rad,
            Omega=1, omega=1, M=-421.17961121323583, r = 0.788*Rearth_to_AU, hash='d')
sim.add(m = 0.692*Mearth_to_Msun, a=2.925e-2, e=0, inc=(90-89.793)*deg_to_rad,
            Omega=1, omega=1, M=-270.15332650367066, r = 0.920*Rearth_to_AU, hash='e')
sim.add(m = 1.039*Mearth_to_Msun, a=3.849e-2, e=0, inc=(90-89.740)*deg_to_rad,
            Omega=1, omega=1, M=-187.20059131394225, r = 1.045*Rearth_to_AU, hash='f')
sim.add(m = 1.321*Mearth_to_Msun, a=4.683e-2, e=0, inc=(90-89.742)*deg_to_rad,
            Omega=1, omega=1, M=-136.97239841027107, r = 1.129*Rearth_to_AU, hash='g')
sim.add(m = 0.326*Mearth_to_Msun, a=6.189e-2, e=0, inc=(90-89.805)*deg_to_rad,
            Omega=1, omega=1, M=-89.87573477315452, r = 0.755*Rearth_to_AU, hash='h')

Then I finish setting up the system with

sim.N_active = 8
sim.testparticle_type = 0
sim.move_to_com()

I then added in one randomly generated massless test particle with the following, which adds a particle 1km above the surface of TRAPPIST1-d with a radially outwards initial velocity:

random.seed(108938)
planet = sim.particles['d'] 
M = planet.m                
r = planet.r                      
v_esc = math.sqrt(2*G*M/r)       
    
px = planet.x                     
py = planet.y
pz = planet.z
pvx = planet.vx
pvy = planet.vy
pvz = planet.vz
 
dir_v = [random.random()*2-1, random.random()*2-1, random.random()*2-1]
n_v = dir_v/np.linalg.norm(dir_v)
vi = n_v*(v_esc)
        
vx = vi[0]+pvx               
vy = vi[1]+pvy
vz = vi[2]+pvz
        
x = px+n_v[0]*(r + 1*km_to_AU)
y = py+n_v[1]*(r + 1*km_to_AU)   
z = pz+n_v[2]*(r + 1*km_to_AU)
        
sim.add(m=0, x=x, y=y, z=z, vx=vx, vy=vy, vz=vz, hash=1)

I first integrated the system without the particle for 100 years, then restarted the simulation with one particle added in and integrated for 100 years. I would expect the planet configurations to be basically the same across these two tests, but as you can see in the attached graphs, the planets have reached extremely different final configurations after 100 years.

Thanks,
Raina

integrations.pdf
integrations.pdf

[ryurongliu]

Hmm, it seems the graphs did not attach - here they are again:

[hannorein]

Thanks! I think this is because IAS15's timestep changes as soon as you add another particles. Even if it is a test particle, it needs to reduce the timestep if there is a close encounter. That will change the planets' coordinates at the 16th digits because of floating point roundoff errors. Over time, the difference will grow exponentially. There is no way to avoid this with IAS15. This shouldn't make a difference when it comes to interpreting your results. But the non-reproducibility can be a bit annoying to work with.

We've had a similar discussion recently, but with the MERCURIUS integrator (#494). Have a look ,maybe this is helpful for your case too.

[ryurongliu]

Ah, I see now, thanks for the explanation - I've been pulling my hair out over this for a couple days now, hah!
Re: the Mercurius solution, has that been added to the official REBOUND package, and does it work with Python? I noticed in the thread you mentioned some issues with the Python version which is what I'm using.

[hannorein]

Yes, that has been added to the main branch. I think it should be working fine. But I haven't used this feature much myself. So I'd be keen to hear if this works for your problem.