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 <#100 (comment)>, 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?