Mercurius for violent systems

[hannorein]

Here is a planetary system which is unstable and violently disrupts itself over short timescale (thanks to @nchoksi4456 for providing this example). MERCURIUS fails at this problem and has energy errors of order unity. I can't see a way to improve this without choosing a very restrictive timestep or switching radius. I suspect that this might not be a bug and that MERCURIUS just can't handle this level of non-Keplerian motion. But I'm wondering if I'm missing something. @dmhernan, I've already tried out if a smoother switching function helps, but it doesn't seem that alone will do the trick. Nevertheless, I wonder if this might be a good test problem since you were looking for one. Note that one can generate different initial conditions by changing the random seed.

Some diagnostics:

Full Code:

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

random.seed(1)

sim = rebound.Simulation()
sim.add(m=1, r=0.005)
for i in range(10):
    loga = random.uniform(np.log10(0.3), np.log10(10.))
    sim.add(m=0.005, r=0.0005, a=10**loga, f="uniform", Omega="uniform", inc=i*np.pi/180.)
sim.move_to_com()    

sim.integrator = "mercurius"
sim.ri_mercurius.hillfac = 3
sim.dt = min([o.P for o in sim.calculate_orbits()]) * 0.05 # 5% timestep
#sim.ri_mercurius.L = "C4"

times = np.arange(0, 5e2, step=1)*2.*np.pi
Es = np.zeros(len(times))
As = np.zeros((len(times),sim.N-1))
Eccs = np.zeros((len(times),sim.N-1))
for j, time in enumerate(times):
    sim.integrate(time)
    Es[j] = sim.calculate_energy()
    os = sim.calculate_orbits()
    for i,o in enumerate(os):
        As[j,i] = o.a
        Eccs[j,i] = o.e   


fig, axs = plt.subplots(1,3,figsize=(15,4))
axs[0].set_yscale("log")
axs[0].set_ylabel("Energy error")
axs[1].set_ylim([0.1,20])
axs[1].set_yscale("log")
axs[1].set_ylabel("semi-major axis")
axs[2].set_ylim([0.,1])
axs[2].set_ylabel("eccentricity")

axs[0].plot(times/(np.pi*2.),np.abs((Es-Es[0])/Es[0]))
for i in range(As.shape[1]):
    axs[1].plot(times/(np.pi*2.),As[:,i])
    axs[2].plot(times/(np.pi*2.),Eccs[:,i])

[acpetit]

Hey,

This is the kind of system I tested my adaptive integrator on so I can provide you a bit of feedback here. The main sources of large energy error were due to unresolved pericenter passages and this is why I implemented Mikkola (1997) regularization (that has been tested on Symba too IIRC).
I just looked at your test case and there is a trend between the minimum periapsis and the energy jump between snapshots (seed 3 here).

I looked rapidly at Mercuruius code and I don't think you switch to IAS15 for planet star close encounters, so it might be a direction to explore.

fig, axs = plt.subplots(1,3,figsize=(15,4))
axs[0].set_yscale("log")
axs[0].set_ylabel("Energy variation $|E_{k+1}-E_k|/E_0$")
axs[1].set_ylim([0.01,20])
axs[1].set_yscale("log")
axs[2].set_xscale("log")
axs[2].set_yscale("log")
axs[1].set_ylabel("Periapsis")
axs[2].set_ylim([0.001,1])
axs[2].set_ylabel("Minimum periapsis")
axs[0].set_xlabel("Energy variation $|E_{k+1}-E_k|/E_0$")

axs[0].plot(times[:-1]/(np.pi*2.),np.abs(np.diff(Es)/Es[0]))

for i in range(As.shape[1]):
    axs[1].plot(times/(np.pi*2.),As[:,i]*(1-Eccs[:,i]))
axs[2].plot(np.abs(np.diff(Es)/Es[0]),np.min(As[:-1]*(1-Eccs[:-1]),axis=1),'.',ms=2)

[hannorein]

Thanks for jumping in @acpetit! I agree, it has most likely something to do with the highly eccentric orbits. I'm not quite sure I understand why regularization helps in this case. Assuming there are no collisions during periapsis, shouldn't the WH part of the code handle this? I guess it's never going to be exact (because of the jump step)...

[acpetit]

I was not thinking about collision regularization but a regularization at the pericentre. In my paper that is the section 6, inspired by this work https://doi.org/10.1023/A:1008217427749 . The problem is that the WH scheme is terrible if the periapsis is unresolved (basically there are not enough perturbations steps when ideally you would like them spread out regularly along the orbit).

When I was testing my code I realized that the big energy errors were happening when a close encounter created very eccentric orbits. However, the error was not induced by the close encounter but by the very next periapsis passage.
This was also the conclusion by Rauch and Holman (https://ui.adsabs.harvard.edu/abs/1999AJ....117.1087R/abstract)
However, I am not sure there is an optimal way to adapt the timestep to always resolve the pericentre in the context of many planet problem.

Anyway, I am happy to exchange more about it if you want some clarifications :)

[acpetit]

Maybe an idea would be to switch to IAS15 whenever a planet enters the region where the pericentre passage cannot be resolved? I just fear that someone may as well directly run in IAS15 in this case because it will lead to tons of switching and I don't think it is great for the energy error too.

[hannorein]

Thanks! FWIW, I don't think it makes sense to switch to IAS15 for every pericenter approach to the sun. Then one might just use IAS15 all the time. But even so, I don't think this would resolve the issue. The WH Kepler solver can solve arbitrarily large eccentricities. It's the perturbations to the Keplerian orbit and the jump step that cause the energy error. And those effects are still there when using IAS15, because of the Hamiltonian splitting of the hybrid integration scheme. Getting rid of that would mean getting rid of the hybrid integration scheme. Then one can just use an arbitrary adaptive integrator (with or without regularization).

[dmhernan]

Hi, thanks for this cool problem and code description! This is definitely a challenging problem. Here are some thoughts, also in response to @acpetit's nice ideas.

1. @hannorein have you experimented by chance with keeping the switching radii fixed, to some large, reasonable values? I wondered how much the energy jumps can be attributed to the switching radius shifting, which then breaks symplecticity.
2. @acpetit if the issue is the resolution of periapsis, that's really interesting! Of course, the regularization approach has the limitation, as I understand it, that it depends on a global adaptive timestep. Ideally we'd want more efficient solutions.
3. If resolving the peripasis is an issue, modifying Mercurius as suggested by Duncan and Levison 2000 (https://iopscience.iop.org/article/10.1086/301553/fulltext/ ) I think could work. The idea would be to keep Mercurius exactly as it is, except to modify according to their eq. (5) and (6), which introduces a new, extra switching function which moves the jump Hamiltonian around. One can make this new switching activation as rare or frequent as one wants.

[hannorein]

Just a quick response

1. The switching radii are fixed by default. Increasing the radii will in general improve things a little, but I don't think there is a reasonable value that would fix the problem completely and still be a hybrid integrator.
2. I agree, Duncan and Levison 2000 might offer a good solution. I can't quite foresee how difficult the implementation would be (it's always harder than I think).

[acpetit]

@dmhernan : Thanks for the paper, I believe I never looked at the details they discuss. So if I understand correctly the pericentre problem in the WH splitting mostly comes from the indirect part?
Now that I think about it, it completely makes sense although I never looked at it that way.

Regarding implementation, one thing to keep in mind is that H_ind and H_int no longer commute in democratic heliocentric coordinates so one cannot combine the two half steps together for both terms.

[dmhernan]

Hi @acpetit, I guess there are different issues here. A specific problem of the Mercurius splitting is being able to apply the jump Hamiltonian frequently enough. This is not a problem for other splittings like WH with Jacobi coordinates. A solution to this problem would be the LD2000 approach, which is to change the Hamiltonian, and is quite effective, as shown in their Fig. 1. I think another solution would be to reduce the timestep, which would achieve the same goal.

A separate issue is that of the stepsize resonances with the orbital period. With or without the LD2000 solution, I guess you'd be subject to this, and the solution is only smaller timesteps, as far as I know.

[hannorein]

FWIW, I'm going to attempt to implement the LD2000 approach. It seems simple enough (famous last words). If you want to follow along (or contribute!) watch out for the LD2000 branch...

[hannorein]

One issue might be that IAS15 currently solves only second order differential equations. But the LD2000 EoM are a set of coupled first order differential equations. I'm not sure how to handle those without changing a lot in IAS15.

[dmhernan]

@hannorein Sounds good. Yeah, the new EoM will probably look messy and very different from the Mercurius ODEs when they converted to a second order system. Another option would be to use a different numerical method than IAS15. Since switching is not meant to be used for long time scales where Brouwer's Law matters, I would think using some other conventional method like Bulirsch--Stoer would work for this new method or even Mercurius. Maybe you've looked into this.

[hannorein]

I agree. Any suggestion for a nice BS implementation that I could use as a starting point? All the ones I have seen are very messy: not scale invariant, lots of magic numbers and goto statements, etc.

[dmhernan]

Hmm, I have used versions from Numerical Recipes in C, or the one in MERCURY in fortran, or one in Matlab, https://www.mathworks.com/matlabcentral/fileexchange/55528-bulirsch-stoer, and maybe others. I don't know how good they are, exactly.

[danielk333]

@hannorein You probably already know these but maybe the links are useful for others reading this thread in the future:

• (free to read + examples) Numerical Recipes C++-edition
• (free to read + examples) Numerical Recipes C-edition
• (free to read + examples) Numerical Recipes FORTRAN-edition
• SWIFT also has a FORTRAN implantation of BS, don't know how clean it is
• I know that Orekit has a BS integrator but I think that is just sub-classed on Apache commons BS integrator, so if one can suffer trough the bleeding eyes from all the Java code maybe that one is actually conceptually cleaned up 😅
• AMUSE has a BRUTUS implementation in C++ and it might have other BS implementations as well

[hannorein]

Thanks! The cleanest implementation I've seen is based on a code originally from Hairer: https://github.com/Hipparchus-Math/hipparchus/blob/master/hipparchus-ode/src/main/java/org/hipparchus/ode/nonstiff/GraggBulirschStoerIntegrator.java.

That might be worth porting into REBOUND (not just for this purpose)...

[danielk333]

Aha, now that I actually looked at the Orekit code it is using the Hipparchus one!

I have found it very useful to have a BS laying around for various purposes so 👍 from me to see it included in rebound.

[hannorein]

I'm slowly making progress. The BS implementation seems to work nicely already (but still needs some polishing before I merge it into the main branch)!

I've been staring a bit more at the equations in the Levison and Duncan 2000 paper. Are the equations of motions coming from Hamiltonian in Eq. 6 well defined for particle $i$ in the limit of $m_i -> 0$? It looks like there is a term $\ddot r_i = 1/m_i ...$ I'm not sure how to interpret that. Does the splitting just not work for test particles? Or am I missing something here...

[acpetit]

I don't think there is a problem here. This is the same idea that when you need to write the EOM of a test particle from the planetary Hamiltonian. Formally, let's consider $1\leq i \leq N$ planets of mass $m_i$ and DH coordinates $(p_i,q_i)$ and a soon to be test particle of coordinates $(p,q)$ and mass $m$ (that will go to zero afterward).

You make the change of coordinates $(p',q') = (p/m,mq)$ and write the Hamiltonian and the equation for these coordinates. By doing so,you can safely take m->0 in the planets EOM and they are not affected anymore by the test particle (except in F, but I'll come back to this).

The EOM for the test particles are now
$\dot q' = m[p' + (\sum (p_i)+mp')/m_0F(q_i,q))]$ and $\dot p'= [standard acceleration term] -(\sum p_i+mp')^2/2m_0 dF/dq'$

The second equation is now without singularity if you take m->0 and in the first one, you can simplify by m on both sides and get the standard equation too.

The only question is what to do with F(r)? The thing is you don't need to include a factor for the test particle there since it will do nothing to the reflex motion. So you might as well not have it there. The case of semi-active TPs can be treated as planets I would say (but semi-active test particles already do weird things so I don't think it matters).

[hannorein]

Hm. I think it's different than for normal Hamiltonian. The problem comes from the F(r) factor. If I implement it naively that will introduce a 1/m_i term. Test particles can be treated differently, but it still seems unphysical to have this divergence for small masses.

[hannorein]

One could add a factor of m_i/m_0 in F(r). That way, massless particles having close encounters with the star would not trigger a changeover.

[dmhernan]

Hmm, I'll check the details later.

[hannorein]

@dmhernan, any more thoughts on this? If this splitting really behaves badly for small masses, it might not be worth implementing.

[dmhernan]

@hannorein Hi, sorry that I'm a bit tied up. Taking a glance, I think I see the problem with \dot{v}~1/m, which blows up for test particles. But I can see in Fig. 7 from LD2000 that their method worked well for test particles. (By the way, is the end goal here a paper or just a new Rebound code without paper?)

[hannorein]

No worries. I agree that it works for test particles. But there might be an issue with particles that have a very small but finite mass.

I'm not sure where this is leading. I'd like to make Mercurius more robust for "violent systems". But if it's just Mercurius + LD2000, then it's not worth writing up.

[dmhernan]

Hmm, I think I see. So I think the issue is avoiding rewriting the Hamiltonian just for the special case of test particles (and also avoiding the small masses issue)? I agree that's unsatisfying and isn't necessary for regular MERCURY... interesting.

[dmhernan]

Something that I'd also be curious about is if changing the Hamiltonian splitting could help in this problem; e.g., what Dehnen and I called WHDS, https://arxiv.org/abs/1612.05329

[hannorein]

I'm not sure. It could be that WHDS solves the problem of planet-sun encounters. But the flip side might be that it leads to larger errors for planet-planet encounters (because the relative position between planets changes during the jump step).

[dmhernan]

Hi, to update this thread, I think I have a simple solution (for sure simpler than LD2000) for this, although I haven't tested on this particular problem yet. Would it be useful to send you some comparison results I have so far, @hannorein ?

[nchoksi4456]

Hi all - just pinging on this to see if there's been any progress. Thanks!

[dmhernan]

It seems this thread gets a new comment each year :D. If this regards my comment, I'm glad to discuss with whomever would like. Others may be working on this problem too, I don't know. If you have ideas, @nchoksi4456, I'm certainly happy to hear about it.

[tigerchenlu98]

Hi @nchoksi4456, to update this thread we have a solution! The new TRACE integrator has been added to REBOUND, which should be a strict improvement to MERCURIUS and does correctly resolve the periapsis (as suspected, this was the issue). We specifically tested TRACE for violent scattering systems, and it performs very well.

The documentation and examples are now on the REBOUND website, and the paper describing it's implementation will be on the arXiv tonight. Please don't hesitate to reach out if you have any questions about it!