Number of collision exceptions thrown varies with integration step

[blackhaz]

Dear All,

Initially address this to Hanno Rein via e-mail, he suggested I post this here for collective wisdom. I am simulating collisions in the inner Solar System using IAS15. I catch the rebound.Collision exception to detect collisions. I use Hill radius to get the exception thrown then calculate collisional radius per each incoming projectile to see if there was a collision. If a particle is inside the collisional radius then I count this as a collision.

The problem I have is that the number of exceptions thrown by Rebound is different depending on the step between simulation outputs. I have a program structure similar to this example:

https://rebound.readthedocs.io/en/latest/ipython/Churyumov-Gerasimenko.html

If I set Noutputs to, say, 1000 then I get about 20 collision exceptions thrown. If I set it to, say, 300000 then I have 166000 exceptions thrown. My expectation was that each sim.integrate() call is interrupted in the middle by an exception, so I should be getting collision exceptions in the middle of the integration, so should be getting the same amount of exceptions without respect of the Noutputs.

I have tried setting exact_finish_time to 0 but what it did result in is that the amount of exceptions thrown stayed about the same, the collisions reported has decreased significantly. I can still get the amount of collisions up by bumping up Noutputs, which shouldn't happen. Ideally, I should be getting the same amount of exceptions even if I do the integration in one go, is that correct?

I am trying to catch collisions with Earth at particles coming in at 30 km/s, so would have to set Noutputs really high to capture data every few hours in the simulation, which would slow down the whole integration of ~ 50 kyr significantly.

Any ideas what could be wrong? Would appreciate any pointers. Thank you!

[hannorein]

What do you do when you get an exception? If you do nothing and continue to run, you will get another exception the following timestep.

[blackhaz]

Each exception is handled as following: I first calculate relative velocity between the planet and all the incoming particles. Based on that I calculate the collisional radius for the planet - per each particle. If a particle is inside the radius, it is removed from the simulation, and I increase the collision counter.

[hannorein]

What is a collisional radius?

[blackhaz]

I adopt the one from Rickman et al. [1], please see equation 5 in their paper. I am still not sure if this is the right thing to do, as gravitational focusing, in principle, should be handled by the integrator.

REFERENCES [1] https://www.aanda.org/articles/aa/abs/2014/09/aa23966-14/aa23966-14.html

[blackhaz]

Addendum: I have just removed the collisional radius whatsoever from my script, and I get wildly different number of collisions now depending on Noutput. I am removing particles from simulation immediately when the exception is thrown without checking for the distance. The number of collision exceptions thrown (and particles removed) still decreases with Noutput. This is with exact_finish_time=0.

UPDATE: If I use sim.collision_resolve = "merge" then I end up with the same amount of particles in the system left after simulation, without respect of Noutput setting, so this issue has probably something to do with the way collisions are handled at the exception. In the collision handler I do a check of which particles have collided depending on the last collision time, as in the documentation example:

for p in sim.particles: if p.lastcollision == sim.t: ...

Perhaps, this has to be done somehow else?

UPDATE #2: I have tried a little different approach. Instead of catching exceptions, I set my own collision resolver:

sim.collision_resolve = handle_encounter_v2

In my resolver function, I follow a strategy similar to reb_collision_resolve_merge(), whereby I remove the particle with the higher index. Now I get the same amount of collisions without respect of the Noutput setting, so everything works fine this way. This still does not explain why I get different number of collisions while catching exceptions, which is the way described in the documentation.

[hannorein]

I'm still not sure what you consider a "collisional radius". If you want to resolve collisions, then using sim.collision_resolve is clearly the way to go.

[blackhaz]

Hanno, when calculating impact probabilities using collisional radius instead of target's radius effectively increases target's cross-sectional area due to the gravitational focusing, e.g. Earth attracting incoming projectile. Thinking this all over, the integrator should probably account for that change of trajectory of the projectile, shouldn't it?

[hannorein]

Sorry for being slow here, but since you haven't posted any code, I still don't know what you're doing with this collisional radius. Maybe this helps:

REBOUND allows you to set one radius per particle. You can set it to whatever you like, but in general it makes sense to think of this as the physical radius of a particle. When two particles overlap, then it is considered a physical collision. Whether you catch every physical collision depends on things like the integrator used, the collision detection method used, the timestep used, etc. It's ultimately a tradeoff between speed and precision. Every application has a different goal.

Complementary to detecting physical collisions, the sim.exit_min_distance variable provides a quick way to stop simulations when particles get close to each other in situations where one does not care about the details of encounters, for example when checking the stability of planetary systems. After every timestep, it checks if any two particles are closer than a sim.exit_min_distance from each other. It's easy to come up with scenarios where one misses an encounter.

[blackhaz]

I apologize, I didn't realize it would be useful. Attached is part of my old exception-based collision handler which we now know does not work, as we should handle collisions via sim.collision_resolve and not by catching exceptions.

Basically, I tried to catch an exception when the particle entered a planet's Hill sphere. Then I calculated the collisional radius (Rcoll) based on the particle's speed relative to the planet (U). If a particle was found inside the collisional radius then I assume this was a physical collision and the particle got removed. Still unsure if that was the right strategy, as the integrator should, in principle, account for this gravitational focusing.

handle_encounter.txt

[hannorein]

Ok. Got it! If you're capturing "collision" with radii of about the Hill radius, make sure your timestep is small enough. Although IAS15 is an adaptive scheme and makes the timestep smaller during close encounters, a Hill sphere is pretty large and the timestep might not get very small. Thus there is a chance that you jump over an "collision", especially if you only scratch the surface of the Hill sphere. This might be what you saw. Given that you don't actually want to resolve very deep close encounters, you might not need IAS15. You could get away with a fixed timestep integrator such as WHFast. Then you can do a normal convergence study, i.e. making the timestep smaller and smaller and checking that your results do not depend on the timestep.

[blackhaz]

Thank you, Hanno. I appreciate your time. I am not sure I can ask you this here, but does using this collisional radius make any sense at all? The integrator should itself simulate the gravitational attraction of an incoming particle by the planet, right? I've tried searching but couldn't find any mention of using increasing radii due to gravitational focusing in their simulations.

[hannorein]

Yes, the integrator does simulate the interactions between the incoming particle and the planet (just try it out 😉) . But there are cases where you might want to ignore these interactions to speed up the simulation. Only in that case would it make sense to artificially increase the radii to emulate gravitational focusing.

[blackhaz]

Thank you very much, Hanno. Appreciate your input.

[hannorein]

You're welcome. I'll close this for now, but feel free to reopen an issue if you have more questions.