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.