MERCURIUS vs. IAS15 speed

[nchoksi4456]

Hi all,

I've been running some simulations of closely packed systems (10 ~jupiter mass planets from ~0.03 to 10 au). Since mergers are common in this setup, I've tried integrating with both MERCURIUS and IAS15 only. I was surprised to find MERCURIUS requiring runtimes a factor of a ~5 longer than IAS15 alone. Below I give a sketch of my code for the MERCURIUS integration, and highlight any lines that differ between the MERCURIUS and IAS15 runs.

sim = rebound.Simulation()
*sim.integrator = "mercurius"*;  # Obviously for IAS15-only integration this would change
sim.G = cgs.G; sim.collision = 'direct'; sim.collision_resolve = 'merge'
sim.move_to_com()

# Omitting a few lines here to randomly generate planet masses, radii, orbital elements...
for mv, rv, av, pv, fv, Ov in zip(mlist, rlist, alist, pomegalist, flist):
      sim.add(m = mv, r = rv, a = av,pomega = pv, f = fv)
   
# Extra lines just for MERCURIUS
**sim.dt = 1e-3*sim.particles[1].P
sim.ri_ias15.min_dt = 1e-4*sim.dt
sim.ri_mercurius.hillfac = 3.0**
for time in times:
    sim.integrate(time, exact_finish_time = 0)
    particles = sim.particles
    Plist = []
    for particle in particles[1::]:
	    Plist.append(particle.P)
    # For MERCURIUS only, update the timestep every 1 yr (spacing in times array)
    **sim.dt = 1e-3*np.amin(Plist)
    sim.ri_ias15.min_dt = 1e-4*sim.dt**

Is this behavior expected, or I am using MERCURIUS incorrectly?

Thanks! PS: out of curiosity, are there significant speed-ups associated with using "natural units" (where semi-major axes, periods, etc. are ~1) instead of, e.g., CGS?

[hannorein]

Have a look at the post at the very end of this thread #401. I think that should answer all your questions!

[nchoksi4456]

Thanks Hanno. I think my question is a little different though: I'm interested in the runtime of REBOUND's hybrid integrator MERCURIUS (WHFast + IAS15) vs. using just IAS15 alone (both algorithms implemented in REBOUND), whereas the post you referred to discusses MERCURIUS vs. mercury6. Wouldn't we expect MERCURIUS to run faster than IAS15 alone for the setup I described? Sorry if I am missing something basic.

[hannorein]

I see. Well, MERCURIUS is not automatically faster than IAS15. It depends on the system you're integrating, what you're interested in, what kind of accuracy you're aiming for, etc. There is no simple solution. You need to play around with some of the parameters and try to find a sweet spot that is fast but still accurate enough for what you want to achieve.

A few hints:

• Your timestep seems rather small.
• Check that the planets are not constantly having close encounters with each other. This depends on the hillfac parameter, but also on what kind of orbits your planets are on.
• Look into the safe_mode option.

[nchoksi4456]

Okay, thanks Hanno!

[hannorein]

Feel free to reopen the issue again later. Especially if you can post a short example that I can run myself, then we can brainstorm on how to make it faster...

[nchoksi4456]

Thanks Hanno. Before sharing any working example scripts with you, I am trying to turn off safe_mode as you suggested. But I occasionally have ejections from the system so I am removing those ejected particles. I think this means I need to recalculate Jacobi coordinates upon removing a particle. Could you confirm that I am doing so correctly in the code below?

for time in times:
	sim.integrate(time, exact_finish_time = 0)
	particles = sim.particles

	Plist = [] 
	for j in reversed(range(1, sim.N)):
		particle = particles[j]
		dist = np.sqrt(particle.x**2 + particle.y**2)
		if(dist > 1000*cgsAU): # Remove particles if they get too far from the star. For hybrid integrator, will need to recalculate Jacobi coordinates
			sim.remove(j)
			sim.ri_whfast.recalculate_coordinates_this_timestep = 1
		else:
			particle.r = calc_Rplanet(particle.m, particle.a) # For my problem I am letting planet sizes evolve w/ time
			Plist.append(particle.P) # Record periods of non-ejected particles so that timestep can be updated accordingly
			sim.ri_whfast.recalculate_coordinates_this_timestep = 0

	sim.dt = 1e-3*np.amin(Plist) # Can fiddle with precise choice later
	sim.ri_ias15.min_dt = 1e-3*sim.dt

[hannorein]

• Mercurius doesn't use Jacobi coordinates (they don't make much sense if particles have close encounters) but democratic heliocentric coordinates.
• You do not need to recalculate coordinates manually when you remove a particle. Only when you manually reshuffle them, or add new ones (in a weird way).
• You use a timestep of 1e-3 times the inner planet's period. That's seems small. For example, in a typical solar system integration, I'd expect ~5% of Mercury's period. In particular, when there are close encounters, you're probably not interested in accurately predicting the motion of each individual planet, you just want to get it right in a statistical sense. So I'd go with something like 1e-2. But this obviously depends on your application.

[nchoksi4456]

Thanks. So my setup looks like this:

	sim = rebound.Simulation()
	sim.integrator = "mercurius"; sim.G = cgsG; sim.collision = 'direct'; sim.collision_resolve = 'merge'
	sim.collision_resolve_keep_sorted = 1; sim.track_energy_offset = 1

        # Omitted lines here adding particles to simulation

	sim.ri_mercurius.safe_mode = 0	# Safe mode OFF
        sim.dt = 1e-2*np.amin([particle.P for particle in sim.particles[1::]])
	sim.ri_ias15.min_dt = 1e-1*sim.dt  # Just a guess
	sim.ri_mercurius.hillfac = 3.0
	
	sim.automateSimulationArchive("test.bin", interval=10*cgsyr,deletefile=True)
        sim.move_to_com()
	for time in times:
		sim.integrate(time, exact_finish_time = 0)
		particles = sim.particles

		Plist = []
		for j in reversed(range(1, sim.N)):
			particle = particles[j]
			dist = np.sqrt(particle.x**2 + particle.y**2)
			if(dist > 1000*cgsAU): 
				sim.remove(j)
			else:
				particle.r = calc_Rplanet(particle.m, particle.a) 
				Plist.append(particle.P)
		sim.dt = 1e-2*np.amin(Plist)
		sim.ri_ias15.min_dt = 1e-1*sim.dt # Just a guess

and is performing much faster than before! A couple small follow-up questions: -- Do you have a rule of thumb for how short the IAS15 timestep ought to be? With track_energy_offset = 1 I suppose I can just make a guess for ias15.min_dt and check afterward to ensure that energy is conserved well enough. -- I am not recalculating the switchover distance after each collision, for simplicity. Am I right in thinking that this is OK so long as masses aren't changing by orders of magnitude through collisions?

[hannorein]

Glad to hear it's faster.

Yes, regarding ias15.min_dt. You could estimate it by dividing the planet radius by its escape speed. But it depends what you want to do. You might be fine with something much larger if you're just interested in some statistical properties.

In these "hybrid" simulations you're always going to get something wrong, but in many cases the physical results remain relatively robust. So a convergence study is never a bad idea. I would run the simulations, measure whichever physical quantity you're interested in, then change one of the numerical fudge factors (preferably in the direction where the simulation runs faster) and see if your physical quantity remains the same. If so, I'd be pretty confident in the result. Energy conservation alone might not be the most telling quantity to monitor!

If two ~equal mass planet merge, I think I would recalculate the switchover radii. Or account for the increase by inflating the radii initially already.

[nchoksi4456]

Hi all,

Sorry to revive this old thread. I had a couple other questions about MERCURIUS.

In a merger between two bodies m1 and m2, how does REBOUND/MERCURIUS assign dcrit of the merger product m3? I ran a few merger test cases myself and found that dcrit(m_new) equaled (to machine precision) either dcrit(m1) or dcrit(m2). In some runs dcrit(m_new) = dcrit(m1) and in some dcrit(m_new) = dcrit(m2). I couldn't find the pattern.

I am running simulations where semi-major axis and mass of bodies may change substantially over time (due to scattering / mergers). So I thought the simplest thing to do would be to just have REBOUND recalculate dcrit (as some multiple of current Rhill) at every timestep. I went into reb_integrator_mercurius_part1 in integrator_mercurius.c and changed the following line (line 413 in my version):

if (rim->recalculate_dcrit_this_timestep){ rim->recalculate_dcrit_this_timestep = 0;

to

if (rim->recalculate_dcrit_this_timestep){ rim->recalculate_dcrit_this_timestep = 1;

This seemed to have the desired effect, but I was finding runtimes that were on average ~25% longer, which I thought was a pretty surprising slowdown for a small extra step. Did I make this change correctly? Is this level of slowdown expected?

Thanks! Nick

[hannorein]

Hi Nick!

The merge function doesn't do anything to dcrit. It's just the dcrit of one of the original particles (which one of the two is random, on purpose). That's the default but you can easily change that by using your own merge function. For example, you might want to always merge the smaller particle into the larger particle.

I agree that the default is not great but I'm not sure what the expected behaviour should be. Depending on what merges with what, one might either want to recalculate all dcrits (you could set recalculate_dcrit_this_timestep=1 in the collision routine to do that), or just calculate a single new dcrit for the one post-merger particle (you could do that by calling reb_integrator_mercurius_calculate_dcrit_for_particle()).

Hanno

[nchoksi4456]

Hi Hanno, thanks for the quick response! Just to follow-up, is recalculating dcrit supposed to be computationally expensive? I had experimented with just recalculating dcrit EVERY timestep, regardless of whether or not there is a collision (that way I don't have to guess as to when $a$ and $m$ have changed enough to warrant re-calculation). I was surprised to find that when I set rim->recalculate_dcrit_this_timestep = 1 (instead of 0) in reb_integrator_mercurius_part1 -- thus forcing the function to recalculate dcrit at every timestep -- that the simulations took ~25% longer (see also the latter half of my last post). I expected a minimal slowdown, since its just recomputing (m/Mstar)^(1/3)*a. ** After staring at the code some more, I'm guessing the slowdown comes not from actually recalculating dcrit itself, but from the synchronization step which happens before calculating dcrit:

    if (rim->recalculate_dcrit_this_timestep){
        // Force recalculation of dcrit at every timestep
        //rim->recalculate_dcrit_this_timestep = 0; 
        rim->recalculate_dcrit_this_timestep        = 1;

        if (rim->is_synchronized==0){
            reb_integrator_mercurius_synchronize(r);
            reb_integrator_mercurius_inertial_to_dh(r);
            rim->recalculate_coordinates_this_timestep = 0;
            reb_warning(r,"MERCURIUS: Recalculating dcrit but pos/vel were not synchronized before.");
        }
        rim->dcrit[0] = 2.*r->particles[0].r; // central object only uses physical radius
        for (int i=1;idcrit[i] = reb_integrator_mercurius_calculate_dcrit_for_particle(r, i);
        }
    }

I think by recalculating dcrit I'm effectively wiping out the speed-up I was getting from having set safe_mode =0, since I have to keep synchronizing. Is that right? If so, is synchronizing before calculating dcrit actually necessary? The $a$ that would be calculated would be slightly incorrect (half a Kepler step off?), but I would guess this is a negligible error for the simple purpose of calculating Rhill (I see there are other criteria for dcrit like "average velocity" or "current velocity" but I am just interested in "hillfac").

Thanks again! Sorry to keep pestering you :) Nick

[hannorein]

I think you're correct with everything (slowdown is primarily due to the synchronization, ~25% sounds reasonable). But it's tricky. I'm not sure you really want to change dcrit every timestep. It's similar to the other problem where you changed the timestep frequently. MERCURIUS is smoothly switching between integrators. It's not a hard switch at precisely dcrit. By changing dcrit for a particle currently near dcrit, you introduce an error. If these small changes are frequent and add up then it could lead to very bad long term behaviour (like in the other case).

I think the better option is to set dcrit manually in your merge routine. It's not very difficult (copy and paste the original merge routine, then add the logic to calculate dcrit for the remaining particle). But let me know if you want me to elaborate.

[nchoksi4456]

Gotcha, thanks. Mergers is one concern, which your suggestion to edit the collision routine addresses. My other concern was just scattering: if a planet's semi-major axis changes substantially, then its dcrit should be updated to reflect the new $a$, right? That was my other motivation for updating dcrit (\propto $a$) at every timestep, since I don't really know a priori when $a$ has changed by enough to warrant recalculation of dcrit. Do you have an idea how to address this?

[hannorein]

You're asking exactly the right questions. But I don't have good answers. In fact I don't think there exists a satisfying solution in the framework of a hybrid integrator.

If the orbits change significantly (because of scattering, etc), then one has to recalculate the dcrit. But it's impossible to know when this needs to be done ahead of time. You could do it on regular intervals, making sure it's not too often (which would introduce errors), but also not too rarely (if you miss a collision that also introduces errors).

A hybrid integrator only performs well when there is some sort of "order" (think of a planet migrating semi-smoothly through a planetesimal disk). If everything is just colliding and scattering with each other all the time (think of star cluster), then the best one can do is some brute force integration (e.g. IAS15). It seems that your case falls somewhere in-between.

(Happy to brain storm more...)

[dmhernan]

@nchoksi4456 If your orbits vary substantially, it's not really what a code like Mercurius is built for. However, people working on stellar dynamics have thought about critical radii for more chaotic type problems, where there isn't even necessarily any dominant mass. The code PeTAR, for instance, is also a hybrid algorithm, like Mercurius, with a transition function. Your problem may be more suited to such a code.

[nchoksi4456]

Thanks Hanno and David, this has been super informative! Right now I'm thinking my best bet might be some combination of IAS15 at early times, when the system is most dynamically active, and MERCURIUS for the longer-term integration, when things have settled down a bit.

[hannorein]

I think that's probably what I would do.

[dmhernan]

No problem, combining Mercurius and IAS15 sounds reasonable! I forgot to mention, if the orbits are unpredictable, another solution for keeping the integration symplectic and fast-ish is our NbodyGradient code in Julia: https://github.com/ericagol/NbodyGradient.jl

[nchoksi4456]

Hi Hanno, I made the following update to the reb_collision_resolve_merge function in REBOUND to only update dcrit when there is a merger. Would you mind confirming that I'm not accidentally breaking anything? I did some tests and it seemed to work as desired, but I thought I would check with you before waiting days for the simulations to run :)

    // Merge by conserving mass, volume and momentum ...
    pi->lastcollision = r->t; // No change here just including for context
    // Update dcrit for merger product
    if (r->integrator == REB_INTEGRATOR_MERCURIUS) {
        struct reb_simulation_integrator_mercurius* const rim = &(r->ri_mercurius);
        rim->dcrit[i] = reb_integrator_mercurius_calculate_dcrit_for_particle(r, i); 
    }

I didn't synchronize before recalculating dcrit, since you seemed to agree earlier in this thread that it should make no appreciable difference, especially since I'm now only updating dcrit rarely (when there is a merger).

Thanks!

[hannorein]

Yes. That looks good. As you say, it's not formally correct, and there might be special cases where this gives the wrong result, but I'd say this is a good workable solution. You could add an if statement to make sure dcrit only increases, but never decreases. I think that would make sense from a physical point of view.