Slowing of certain orbits with rebound version 4.0.3

[KatKiker]

Since the rebound update to 4.0.3, we noticed some of our orbits were taking a long time to run, or in some cases hanging indefinitely. Pinning the version to 4.0.2 resolves the issue. The following code recreates the behavior:

import rebound
import assist

this_breaks = True

print(f"rebound version: {rebound.__version__}")

ephem = assist.Ephem("../data/linux_p1550p2650.440", "../data/sb441-n16.bsp")

sim = rebound.Simulation()

sim.t = 14658.412205805536

if this_breaks:
    sim.add(
        x=-0.9716337587843651,
        y=0.5557887162674847,
        z=0.15913299510349357,
        vx=-0.004071385625678586,
        vy=-0.016941109480633685,
        vz=-0.004809720968068568,
    )
else:
    sim.add(
        x=-0.971634,
        y=0.555789,
        z=0.159133,
        vx=-0.004071,
        vy=-0.016941,
        vz=-0.00481,
    )

ax = assist.Extras(sim, ephem)

sim.integrate(14718.412205805536)

I don't know if it's helpful or not, but the truncated values do not produce the same error, which is why I was having trouble recreating the issue.

[hannorein]

Try setting sim.ri_ias15.adaptive_mode = 1 to reproduce the old time stepping behaviour.

[hannorein]

Any success with the above change?

[KatKiker]

Yes! It looks like adding that line fixes the slow-down

[hannorein]

Good to hear.

Some background: We've changed the way the time stepping works in IAS15. It should be better, not worse in all cases. So I'd be curious to know more about when this happens. Is there any chance that this is happening during a very close encounter? If so, it might be interesting to look at just how close the encounter was - specifically is the closest distance closer than the physical size of the objects involved in the encounter.

[KatKiker]

We're looking at impacting orbits, so I suspect that is the case. And in fact in the above example, it looks like the issue occurs right around its expected impact date of 66233. Propagating to 14688 (taking into account the jd_ref) finishes quickly, while 14689 hangs.

It does seem like in most other cases the new time stepping method is faster.

[KatKiker]

Setting the min_dt seems to prevent the hanging while still taking advantage of the speed improvements in the new adaptive mode. Do you know what the 'time unit' is in this case? Would setting this to a small default value cause any issues?

[matthewholman]

Glad to hear that works! min_dt would be in units of days for ASSIST.

On Tue, Feb 6, 2024 at 12:33 PM Kathleen Kiker @.***> wrote:

Setting the min_dt seems to prevent the hanging while still taking advantage of the speed improvements in the new adaptive mode. Do you know what the 'time unit' is in this case? Would setting this to a small default value cause any issues?

— Reply to this email directly, view it on GitHub https://github.com/matthewholman/assist/issues/100#issuecomment-1930437097, or unsubscribe https://github.com/notifications/unsubscribe-auth/AA5CVWMTMATJ6DSINNHOFJTYSJSN5AVCNFSM6AAAAABCCVKIG6VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTSMZQGQZTOMBZG4 . You are receiving this because you are subscribed to this thread.Message ID: @.***>

-- Matthew J. Holman, PhD Senior Astrophysicist Center for Astrophysics | Harvard & Smithsonian 60 Garden Street, MS #51 Cambridge, MA 02138 (617) 496-7775

[hannorein]

Thanks for the details, @KatKiker! I'll try to look into it and see if this is expected behaviour (which I suspect) or something is going wrong...

[hannorein]

I can confirm that this is due to a close encounter. The initial conditions which work, get only within 150km of the Earth's centre (note even that is within the Earth). The initial conditions that don't work get to within 1cm of the Earth's centre!? That is super close! I assume you set it up that way on purpose. But the bottom line is, the integrator is working as well as it can. In the old REBOUND version, it might have not slowed down, but the result was probably not very accurate.

The solution is to include physical collisions and stop whenever two particles collide. Right now you need to check for collisions manually. @matthewholman we could consider including this into assist.

[matthewholman]

Thanks for looking into this, @hannorein . I think it makes sense, for this kind of application, to include internal checks for collisions and close approaches. That said, I'm not sure what the best approach is for IAS15. Would we monitor the distances within a time step or only at the beginning or end?

[hannorein]

I would do it at the end of each timestep. One can use the instantaneous positions and velocities and the assumption that particles travel mostly along straight lines to check if there was a collision during the last timestep. This is implemented in REBOUND as REB_COLLISION_LINE. But I'm not sure what the best approach is to use this for ASSIST because the planet particles are not actually part of the simulation.

[bhanson10]

Wanted to follow up on this to see if there been any updates to ASSIST to include automatic collision detection. I'm running a Monte Carlo simulation using ASSIST, and I would like to check how many of my random initial conditions collide with a certain celestial body. The propagation period is quite long, so I use large time steps initially, but as the samples come to the close approach I have to manually change the size of the integration step in order to not "jump over" the epoch where the collision occurs. Is there any way to do this automatically? Could I get the "exact" expected epoch where the collision occurs without changing the time steps manually?

[hannorein]

@bhanson10 IAS15 is an adaptive method. You should not adjust the timestep manually.

There is currently no way to get the "exact" epoch where a collisions occurs. This would be a nice feature to have but it hasn't been implemented. The way to do this would be to make the LINE algorithm in REBOUND compatible with ASSIST.

[bhanson10]

I was interested in plotting the propagated trajectory, so using a fixed timestep, I integrate a step with ASSIST, then output the particle's position, and repeat. The issue is that as I get closer to collision, this fixed timestep may jump over the collision. Thus, I ended up using a quadratic adaptive timestep that outputs the trajectory at smaller steps closer to a collision. Below is a screenshot of the equation used to calculate an adaptive timestep for outputting the trajectory:

0 < Eta < 1 is a safety parameter that can refine or coarsen the timestep (I used 0.2 arbitrarily). As I was only interested in collisions with the Earth and Moon, I only calculated the Delta t for these bodies. Is there a more efficient/streamlined way of carrying out this procedure with ASSIST?

[hannorein]

I'm not sure I completely follow! IAS15 is adaptive to begin with. We tested the timestepping algorithm extensively, see https://arxiv.org/abs/2401.02849. So you should not need to do it yourself.

Also, I don't understand the units in your equations. r^2/v does not have units of time.

[bhanson10]

I'm probably being unclear. Take, for instance, this snippet of code from "assist/examples/asteroid/problem.c":

// Integrate until we reach tend or an error occurs
while (r->t < tend && r->status <= 0){
   // Generate output every 20 days
   reb_simulation_integrate(r, r->t + 20.0);

   // Output test particle positions
   printf("t = %.1f\n", r->t + ephem->jd_ref);
   for(int i=0; i<r->N; i++){
       struct reb_particle p = r->particles[i];
       printf("particles[%d]:  \tx = %.12f \ty = %.12f \tz = %.12f\n", i, p.x, p.y, p.z);
   }
}

This snippet is outputting the position every 20 days, so the Delta t I am referring to would be 20 days for this case (under the hood I know the IAS15 is doing an adaptive time step to get to that epoch, but this is what I mean when I refer to Delta t). I would like to output the trajectory (more precisely, save the trajectory to a .txt file) of the asteroid as it collides with either the Earth or the Moon. Using a large fixed Delta t like 20 would result in missing the epoch where the collision occurs, but using a very small fixed Delta t would result in unnecessary refinement far away from the collision. Therefore, I was using the procedure from my previous comment to get an "adaptive" output rate that is useful for visualizing the trajectory near collisions.

As for the equation itself, this was just a heuristic I was using. I tried a linear relationship r/v first (which would have units of time) but this did not refine the step enough close to collision, and making eta smaller resulted in inefficiencies far away from collision. Using the quadratic relationship r^2/v worked much better, and since this is just heuristic, eta could have units of 1/length to make everything work out.

I was interested in knowing if there is a better, built-in way of doing this using ASSIST. As the IAS15 is refining the step size near a collision, I was wondering if there was anyway I could call the true timestep size and use that in some way.

[hannorein]

I see. You could use reb_simulation_step(r) which does one timestep of whatever stepsize IAS15 thinks is appropriate.