Encounter Prediction Switching Condition for TRACE

[tigerchenlu98]

@hannorein @dmhernan

This is in response to a few comments/suggestions we've gotten - my first try at implementing a switching condition that accounts for close encounters mid-timestep. In brief, we crudely calculate positions and velocities at dt/2 behind and ahead of the check and use the interpolation scheme from MERCURIUS encounter_predict. From preliminary testing it catches many more close encounters in large N problems, where relative velocities may be quite high compared to rcrit.

At the moment this is quite a bit slower than the default condition without a huge accuracy boost, so I don't think it's worth making the default, but having it as an option may be useful. There's also certainly room for optimization, let me know if you have immediate thoughts.

Also not quite sure why the checks are not passing... everything is fine on my machine (python and C) for what it's worth.

[hannorein]

If the timestep is this large, then aren't you basically ignoring planet-planet interactions during conjunctions unless there is an encounter.

The runs on windows sometimes fail because of a timing unit test which is probably too strict (you should be able to click on the failed tests and see the error).

[tigerchenlu98]

Yes, you are. Increasing the switching radius/decreasing the timestep could fix this in principle, but for these large N problems this seems more robust.

Hmm, I'm seeing a different error in my tests:

ERROR: ERROR: Failed to build installable wheels for some pyproject.toml based projects (rebound)

[hannorein]

Ah sorry, I just assumed it's the usual timing thing, but it's a syntax error. My bad. If you look at the output, the compiler error says that the following and similar expressions are not correct:

struct reb_particle pi_pre = {};

You probably want to have

struct reb_particle pi_pre = {0};

A few more things:

• I don't understand what this algorithm does. Sorry! Can you explain it in a bit more detail?
• Is this something like the LINE collision detection algorithm? If not, wouldn't something like that be better?
• Why do you need to estimate coordinates half a timestep ahead AND behind of the current timestep? What exactly is stored in these new arrays?
• Do you need to safe and load these new arrays to bit-wise reproduce the results after restarting a simulation?
• Does this break the time symmetry of the algorithm in new ways?
• You might want to add some basic units tests.

[dmhernan]

Hello, to add some comments (we can also discuss via email):
--- the idea in extrapolating forward and backwards in time was to maintain a time-symmetric criteria that doesn't depend on the sign of velocity. Otherwise, if we only extrapolate forwards (or backwards), the condition gives a different result/close encounter prediction if we switch the sign of velocity. This is solved if we simply extrapolate both forwards and backwards.
---Tiger and I also discussed doing a simple linear extrapolation rather than the cubic Mercurius interpolation, which should avoid all overhead (Not sure if you tried it, Tiger). If r_n and v_n are the current planet positions and velocities during step n, we simply predict r_{n \pm 1/2} = r_n \pm t v_n where t \in [0,h/2]. Then we take the separations as the difference of all the r_{n \pm 1/2}. We take the absolute value of separations and find the minimum, which is used in the switching criteria, rather than using r_n as is done currently. This probably would catch something.

[dmhernan]

And yes, after checking, I guess this linear extrapolation could reuse something from LINE collision detection?

[tigerchenlu98]

I'll keep things on the PR in case others want to follow along:

• The algorithm is exactly the same as the MERCURIUS encounter prediction algorithm: you have your beginning and endpoints and approximate the motion between the two as a cubic, then take the derivative to find the minimum distance within the timestep. A couple differences, all regarding the choice of pre-state and post-state. MERCURIUS uses the present simulation as the pre-state and uses a Kepler step to find the post-state, which is convenient with the way MERCURIUS is set up. We use times at t-h/2 and t+h/2 as the pre- and post- states for reversibility reasons as @dmhernan mentioned. It's not convenient for us to use a Kepler step since we don't get it for "free" computationally as MERCURIUS does, so I just use a simple Euler scheme. This can probably be improved, both in accuracy and efficiency.

• Looking briefly at LINE, yes this looks useful, I'll dig a bit more into it

• I don't think we need to store the pre- and -post arrays. They should be able to be constructed from scratch at the beginning of each timestep.

• Yes, unit tests would be great, I'll add a few.

[dmhernan]

OK, so if I understood everything, here's the difference between what LINE is doing and Tiger's change. He's doing a prediction with Euler as
r_{n \pm 1/2} = r_n \pm t v_n (1)
v_{n \pm 1/2} = v_n \pm t a_n, (2)
with a_n the acceleration. Using the four r and v (r_n, r_{n \pm 1/2}, v_n, v_{n \pm 1/2}) we construct a cubic describing the position change. However, according to Euler above, r is changing only linearly in time so I'm not sure what the meaning of the cubic minimum is in this case.

LINE would do no velocity update (because velocities are constant, yes?), and eq (2) above would be gone. And definitely no cubic minimum would be needed. I'm guessing all the slowdown is coming from these cubic solvers right?

[hannorein]

Thanks for the explanations. That helps!

• I'm still confused as to how this is working. I see one function which uses the pre/post arrays and one which sets them. So there are two switching functions? How can I set two switching functions?
• I'm also not sure what cubic refers to here. Do you really solve a cubic equation? My gut feeling is that doing the simplest thing (assuming the velocities are constant = "LINE") might work best.

[tigerchenlu98]

Ah, my mistake! I meant to implement

r_{n \pm 1/2} = r_n \pm t v_n \pm 1/2 t^2 a_n

but forgot the acceleration term. It's now in (so now the cubic makes sense) and seems to be meaningfully better than the default case. On the accretion example I tested we catch 50% more close encounters and go from 1.3e-2 to 5.3e-3 relative error.

@hannorein The two switching functions and arrays are a hack, this can certainly be optimized still. My logic was this: we need the pre- and -post positions/velocities for all particles. But the way the switching function syntax is set up right now, we calculate these positions/velocities N^2 times (since the switching function is called for every pair of particles), rather than the N which is all we need. So I just computed them in the pericenter switching condition loop (since we get the approximate accelerations "for free" for the PRS condition anyway). The final implementation will have to move them out of here, and I'm not sure what the best way to do this is -- we shouldn't demand every switching function does this, since the default doesn't need the encounter prediction.

[dmhernan]

Hello, by cubic I just mean the Chambers 1999 cubic interpolation polynomial, eq. (11) and pasted here. Given r_n, r_{n \pm 1/2}, v_n, v_{n \pm 1/2}, we calculate the minimum r. But @tigerchenlu98 , note that in your new modification, I believe r is described by a quadratic polynomial, so the cubic machinery is still not needed. We just need the minimum of a quadratic equation now.

And I'd also encourage what @hannorein mentions as well that we try velocities constant (so then it's only the minimum of a linear function!)

[dmhernan]

Also, note I think we need this,
r_{n \pm 1/2} = r_n \pm t v_n + 1/2 t^2 a_n (acceleration doesn't change sign)

[tigerchenlu98]

You're both totally right, the cubic is completely unnecessary -- don't know what I was thinking

[dmhernan]

Might be good, once the new switching criteria is implemented, to verify it's still reversible and there is still no energy drift in say, chaotic exchange problem.

[hannorein]

@tigerchenlu98 The unit test fails because of the {0} issue. Otherwise your c/python interface looks right. But let me know if you think there is another issue...

[hannorein]

The PR looks good to me. Just one thing regarding the naming: right now we have reb_integrator_trace_switch_default but there is not longer a default function for the pericenter. Would it make sense to rename reb_integrator_trace_switch_default to a more descriptive name as well?

[tigerchenlu98]

Hi Hanno,

There is a default pericenter condition -- reb_integrator_switch_peri_default. Unless you're proposing changing both? We could do switch_hill and switch_peri_prs24. I lean slightly towards just keeping the defaults since I think prs24 may be confusing, but no strong preference. Let me know what you think and I can make the appropriate changes.

[hannorein]

It looks like you're removing reb_integrator_trace_switch_peri_default?

[tigerchenlu98]

Ah I didn't mean to do that -- it's added back now.

[hannorein]

Ok. Now it makes sense. Thanks!