Biased error in the Mercurius integrator?

[rafanw72]

Hi Hanno, I would like to report a possible bias in the simulations using Mercurius if of course, I'm not missing anything in my implementation. I did a simulation with a thousand test particles and the Sun and it seems to me that the particles away from the sun are not keeping their eccentricities comparing with WHFAST or IAS15. See my implementation and a graph attached, please. Thanks. program.tar.gz graph.pdf

[hannorein]

Hello,

It's quite possible that there is something to what you're describing. Your code is a quite long and a bit hard to look through. Maybe you can simplify it a bit further by only keeping the minimum to reproduce the issue? And could you describe a little more what the setup is, what you're measuring, and what outcome you're expecting?

Cheers, Hanno

[rafanw72]

Hi Hanno, Ok, I did a simplified version of the codes that reproduce the issue. I think that I just found out what was going on here. Even using r->testparticle_type=0, the particles do influence each other in the simulations. Please look at the first code (problema.c), there are two particles interacting only with the Sun initialized with no-null masses but using r->testparticle_type=0. Note in graphs 1 and 2 that their semimajor axes change with the time, it is not what I expected if the particles dot not influence each other gravitationally. In the second code (problemb.c) I integrated 1000 thousand particles interacting with the Sun using again r->testparticle_type=0 but with not-null masses. Please note in graph 3 that represents a snapshot of 1 My of the integration. The eccentricity and semimajor axis of the particles changed from their initial value (green curve). Again, this result is not expecting if the particles dot not influence each other. Once, you initialized the particles with null-masses then the particles do not influence anymore and the issue goes away. What do you think?

If you need anything more let me know please,

Thank you,

Cheers. problemb.tar.gz problema.tar.gz

[hannorein]

Hello stranger,

Just to let you know, I'll look into it. I just didn't get around to it yet.

Hanno

[rafanw72]

Hi Hanno

Thanks for the feedback. I'm sorry I did not sign my last two messages and I just realize it now. My name is Rafael Ribeiro.

Cheers,

Rafael Ribeiro

[hannorein]

Hi Rafael,

Thanks for the clarification but I have to say it's still not really "simple". 😉

Correct me if I'm wrong, but I think you're saying that when integrating test particles in a Keplerian orbit around a single massive particle, then the eccentricity does not remain constant. If so, then that's indeed a major bug.

One thing that you should definitely change is the way you initialize the particles. Using

struct reb_particle sun;

is dangerous and can lead to non-deterministic behaviour. Do this instead:

struct reb_particle sun = {0};

It would be great if you could cut down your code by another factor of ~4. Get rid off all the units, hashes, restart files, output files, etc. Just focus on the bare minimum to run the code and output the eccentricity.

Hanno

[rafanw72]

Hi Hanno,

Thanks for the clarification as well. Yes, indeed when I'm integrating test particles in a Keplerian orbit around a single massive particle, the eccentricity does not remain constant (the semimajor axis and the mean longitude also do not remain constants).

Ok, thank you! I used to initialize the particles like that "struct reb particle sun = {0};" to be sure that other elements of the structure are null. I'm sorry I did not do this time.

I think now the code is simple and shows the problem. I cut down my code by a factor of 4 or more and I did all the changes on it which you suggested. Please check it out. I'm integrating two test particles and a massive star. The output in your screen will show the eccentricity variation of the test particle b. The output files (part-a.dat and part-b.dat) show the keplerian elements of the two test particles. The columns of the files are (time, semimajor axis, eccentricity and eccentricity variation).

The graph shows the eccentricity variation for the two particles in a function of the time.

If you need anything let me know, simple_code.tar.gz

Cheers,

Rafael Ribeiro

[acpetit]

Hi Hanno, Rafael,

I got intrigued by this issue and I think this comes from the indirect kick in Mercurius.

Looking at the function reb_integrator_mercurius_jump_step you sum the momenta over all particles instead of only the active ones. The eccentricities evolution shown by Rafael are too small for direct interactions given the particles initial conditions.

Cheers, Antoine Petit

[hannorein]

Rafael, thanks for the simpler code! Much better. Looking into it now...

Hi Antoine! Thanks for jumping in! I need to think about it a bit more. But can you elaborate on the argument for not summing over all particles in the jump step? With these semi-active particles, we're just ignoring the interactions between them, but not between them and the star. The Hamiltonian for the jump step is just a sum of the planets' momenta squared. Why would ignoring some of the interaction terms affect it?

Hanno

[hannorein]

For reference, here is a python version of your code, Rafael. I think it should do exactly the same as yours.

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

sim = rebound.Simulation()
sim.add(m=1)
sim.add(m=5e-4,a=1,e=0.1)
sim.add(m=1e-4,a=1,e=0.1)

sim.integrator = "mercurius"
sim.dt = 1e-4*2.*np.pi
sim.N_active = 1
sim.testparticle_type = 0

Nsamples = 1000
times = np.linspace(0.,3e5,Nsamples)
es = np.zeros((Nsamples,2))
for i,t in enumerate(times):
    sim.integrate(t, exact_finish_time = 0)
    es[i] = sim.particles[1].e,sim.particles[2].e

fig,ax = plt.subplots(1,1)
ax.plot(times,es[:,0]-0.1)
ax.plot(times,es[:,1]-0.1)

Plotting the change in eccentricity, I do not see the large oscillation that you see:

I'm at a bit of a loss to what there could be going on. Which version of REBOUND are you using (ideally which commit)?

[acpetit]

Hey,

Maybe I misunderstood what you call particles of type 0. If they are "test bodies" in the sense that they do not affect the massive one i.e. the star in this case, they should not be taken into account in the reflex motion. I would say that it is a bit weird that they even have a mass...

However if they have a mass and interact with the massive bodies but not between each other (what I thought type 1 was) then your code is correct. But in this case, the reported behaviour is expected because the two bodies will necessarily interact through their effect on the star.

Your code do not reproduce Rafael's example because you chose the same semi-major axes for the particles. If you initialize your simulation with 1.1 for the smaller one, you get an oscillation because the two reflex motions are not synchronized.

[rafanw72]

Hi Hanno and Antoine,

Thanks for the feedback. I'm using the most recent version (v3.12.1) of REBOUND (attached). Like Antoine, I also do not think that it comes from the direct gravitational interaction between the two test particles, yes, the eccentricities should be much bigger than I've shown. Somehow, for the other integrators also is happening the same behavior of the eccentricities. Maybe Antoine is right, the test particles type 0 should not interact with the massive particles, right? If so, they do not feel the reflex motion of the star. Or I did not understand right and the particles type 0 are still semi-active? In this case, is everything right but the particles type 0 and 1 are the same types of semi-active particles. Let me know if I can help somehow.

Thank you guys,

Rafael Ribeiro
rebound_last.tar.gz

[hannorein]

Hi Antoine,

I see! I got confused - you got it right!

I think for testparticle_type=1 I would expect this type of oscillation. But it should indeed not be there for testparticle_type=0. In the later case, the mass of the particle should never be used.

I'll see if fixing the jump step fixes the problem...

Thanks so much! Hanno

[hannorein]

I pushed these changes to the fixjump branch (5acf4a9532701). Note that the same problem applies to WHFAST, not just MERCURIUS. The patch fixes the problem encountered by Rafael, but all kind of other tests now fail. I'll continue to look into it...

[hannorein]

Never mind, that was just an error in my patch. Hurray for units test! This seems commit seems to fix everything: f7a34d49

I'll think about it a bit more, at a unit test or two, and will then merge it into the master branch.

Thanks to both of you for pointing this out, providing all the examples, and, best of all, find the solution!

[Shirui-peng]

Hi all,

Just found this thread -- this is great! Konstantin B. and I have been looking into the consequences of this exact effect (i.e., indirect interactions among massive, non-interacting particles) on simulations of planetary accretion. Perhaps unsurprisingly, this effect also appears in MERCURY. The good news is that the artificial excitation of the velocity dispersion in planetesimal disks is very minor in real systems.

Best, Shirui

[hannorein]

This is awesome! Thanks for letting us know. I'm somewhat surprised the same issue exists in MERCURY, but good to know. I guess in Rafael's case it mattered more, because the momentum of even small test particles can be quite large during close encounters.

[rafanw72]

Yeah for my case is really important. This is really great that we have found a solution to this problem! I'm happy to help with it! It surprises me that it appears also in Mercury. Thank you guys that inspire me a lot doing this great work! Cheers,

Rafael Ribeiro

[rafanw72]

Hi Hanno,

I'm sorry for the return to this issue. I do not know if the problem is solved yet. Using the same simple code I plotted again the eccentricity variation of the two test particles after your last commit. Still, there are indirect interactions among massive and non-interacting particles (the same for IAS15). I got different results but still with a large variation of the eccentricity. Interesting that now the large variation of eccentricity is about part-a. Please, look at the code once more and the graphs attached. Graph1 shows the result (eccentricity variation in function of the time) with the Mercurius integrator, in graph2 shows the same but using the IAS15 integrator, in graph3 shows the result that I expected when the particles really have null-masses. Am I missing something? If not I think we would think a lit bit more about that.

simple_code_2.tar.gz

Cheers,

Rafael Ribeiro

[hannorein]

This is because the routines initializing the orbit from orbital elements and the routines calculating the orbital elements still take the particle mass into account. You can verify this by manually setting the particle mass to 0 before calculating the orbital elements, then afterwards setting it back to the finite value. There will not no difference in the actual integration, but you don't see the oscillations anymore.

[hannorein]

(I think this is not a bug, but the expected behaviour. It's just a somewhat unphysical case, having testparticles which do not interact gravitationally in the integrations, but nevertheless have a mass)

[rafanw72]

Hi Hanno,

Thank you! So if you're saying the routines which are initializing the orbit and the routines calculating the orbital elements take the particle mass into account. So, if I integrate, using those same routines, one single particle (part-a) and a Star but now massive and interacting particle (active) I should expect the same variation of eccentricity. It is right? But look at the graph below, The graph shows two different integrations: in green (particle massive and interacting particle or active type and a star) and purple (part-a massive and non-interacting particle or test particle and a star). I would expect that they have had the same eccentricity variation because the routines are taking into account the same masses. Do you agree?

Code:

Green: struct reb_simulation r = reb_create_simulation(); r = reb_create_simulation(); r->integrator = REB_INTEGRATOR_IAS15; r->testparticle_type = 0; // test particles will not feel the gravity from other test particles. r->dt = 0.00012.*M_PI; struct reb_particle sun={0}; sun.m=1.; sun.r=1e-2; reb_add(r, sun); r->N_active = 2; struct reb_particle part1 ={0}; part1= reb_tools_orbit_to_particle(r->G, sun, 5e-4, 1.0, 0.1, 0.001, M_PI, M_PI, M_PI); part1.r=1e-6; reb_add(r, part1); reb_move_to_com(r); r->heartbeat = heartbeat; reb_integrate(r, 3e5);

Purple:

struct reb_simulation r = reb_create_simulation(); r = reb_create_simulation(); r->integrator = REB_INTEGRATOR_IAS15; r->testparticle_type = 0; // test particles will not feel the gravity from other test particles. r->dt = 0.00012.*M_PI; struct reb_particle sun={0}; sun.m=1.; sun.r=1e-2; reb_add(r, sun); r->N_active = 1; struct reb_particle part1 ={0}; part1= reb_tools_orbit_to_particle(r->G, sun, 5e-4, 1.0, 0.1, 0.001, M_PI, M_PI, M_PI); part1.r=1e-6; reb_add(r, part1); reb_move_to_com(r); r->heartbeat = heartbeat; reb_integrate(r, 3e5);

Thank you,

Cheers,

Rafael Ribeiro

[hannorein]

If you have the exact same initial conditions in both cases, then the eccentricities should be different because the star feels the testparticle in one case, but not in the other. Don't get too hung up on this! I'm not sure if there is really something here worth thinking about. If you want to use a test particle, I suggest you just set the mass to 0 and you should always get what you expect.

[rafanw72]

Ok, thank you very much for clarifying it to me! Now, I understand, I will try to adapt all my problems, where the mass of the test particles matters, setting the particle mass to 0 before calculating the orbital elements. I was doing some calculations in real systems of long-period comets and I notice that particles far away (100, 000 au) from the star are presenting large variations in eccentricities (even setting the particles with zero mass). However, I do not want you to be worried about it. I'll do more tests here to be sure that I'm not missing anything. Thanks a lot!

[hannorein]

Ok! Feel free to open another issue if you think there is still a problem.

[hannorein]

@Shirui-peng noticed that there might still be an issue even after the above fixes. Here is a short code to illustrate the issue:

import numpy as np
import rebound
import rebound.data
%matplotlib inline
import matplotlib.pyplot as plt
print(rebound.__build__)

sim = rebound.Simulation()
rebound.data.add_outer_solar_system(sim)
sim.integrator = "whfast"
sim.ri_whfast.coordinates = "democraticheliocentric"
sim.N_active = sim.N - 1
sim.testparticle_type = 1

sim.dt = 0.01

sim.move_to_com()
E0 = sim.calculate_energy()
G0 = sim.calculate_angular_momentum()

times = np.linspace(0.,20000.,100)
errors = np.zeros(len(times))
amerrs = np.zeros(len(times))

for i,t in enumerate(times):
  sim.integrate(t,exact_finish_time=0)
  errors[i] = abs((sim.calculate_energy() - E0)/E0)
  G = sim.calculate_angular_momentum()
  amerrs[i] = abs((G[2] - G0[2])/G0[2])

plt.yscale("log")
plt.plot(times/(2e3*np.pi),errors,'r',lw=2,label='$|\Delta E/E_0|$')
plt.plot(times/(2e3*np.pi),amerrs,'b',lw=2,label='$|\Delta G/G_0|$')
plt.xlabel("Time (kyr)")
plt.ylabel("Errors")
plt.legend()
plt.show()

When testparticle_type=1 is set, then the angular momentum error is significantly larger compared to when all particles are active, or when testparticle_type=0. I think this is due to a bug in the coordinate transformations. I think the commit 75ece6da should fix this.

[rafanw72]

Hi Hanno,

Sorry for being late to answer! Thank you very much for let me know about it and to fix it! Cheers,

Rafael Ribeiro