Collision resolve can get messed up if multiple mergers per timestep

[joshuawallace]

When a collision is detected, collision.c currently stores the indices of the particles involved in the collision, specifically the indices correspond to their positions in the particles array. In the case of

-multiple collisions in a single timestep, and -particle mergers

I think I have found an issue. I give a specific example to highlight this. Let's say we have four particles, referred to by their original indices in the particles array: 0,1, 2, and 3. Let's say after a time step there is a collision between 0 and 1 as well as a collision between 2 and 3. The collision between 0 and 1 is resolved first and results in a merger. The information for particle 0 is updated accordingly and particle 1 is removed. If keepSorted=1 in reb_remove(), particles 2 and 3 will move down to fill the gap left by particle 1 (final order: 0, 2, 3. Now it's time to resolve the collision between particles 2 and 3. The collision resolve function will now be looking in the memory at where the particles used to be located in the particles array, not where they are now. Where particle 3 used to be in the array is no longer technically included in the simulation particle, but the particles array still has memory allocated for particle 3's old spot and particle 3's information is still stored there. The collision resolve function will now be trying to resolve a collision between particle 3 (at particle 2's old spot) and particle 3 (at particle 3's old spot, now off the portion of the array that the simulation cares about).

So, it seems that making keepSorted=0 in reb_remove() for the initial removal of particle 1 would fix this. However, it will not. If the collision between particle 2 and 3 leads to a merger, then particle 2's information will be updated based on the merger information, and then reb_remove() will be called on particle 3's old position, which is outside of the array currently recognized by the simulation. reb_remove() will print out the following error: "\nIndex %d passed to particles_remove was out of range (N=%d). Did not remove particle.\n" and particle 3 (which is at particle 1's old spot) won't actually be removed from the simulation.

My proposal to fix this:

-Assign each particle a unique ID and rewrite the reb_collision_search() function to give particle ID's to the collision struct. These ID's will follow the particles wherever they go in the array.

The simple fix I propose above would make it so particle ID's would have to be assigned. REBOUND currently does not require particle ID's to be assigned (as far as I know) so the fix proposed above would probably not be ideal for a general fix. Any other ideas?

[hannorein]

I think you're right.

I would suggest a different solution, however. When you resolve the collisions found in a given timestep and one collision ends up to be a merger, then I would go through the list of all other collisions that occurred during that timestep and update the particle references in that list accordingly.

A few reasons for why I think this might be better:

• It let's the user use whatever they want for the ID. For example, some people might want to have all planets have the index 1 and all planetesimals the index 2.
• One only needs to change the code in one place.
• Multiple mergers per timestep should be rare, so it will not negatively impact performance.

Which routine do you use for the merging part?

[joshuawallace]

The routine I use for merging is basically identical to the collision_resolve_merger() function you have defined in the problem.c file in the example/mergers directory. I've replaced these two lines:

r->N--;
particles[c.p2] = particles[r->N];

with

reb_remove(r,c.p2,1);

For what you propose, it would be good for the resolve() function to know both the values of collisions_N and what collision number it is, so resolve() can know whether it needs to update the index values in other collisions or if it is the last/only one and doesn't need to bother.

An alternative method would be to let reb_collision_search() itself update the index values when looping over the collisions. In this case, it would be useful for reb_collision_search() to know whether keepSorted is set to 1 or 0.

I don't see a clear way to implement either of the methods above without editing the input parameters for resolve(), or setting a value for keepSorted with a sufficiently large scope that reb_collision_search would be able to see it (perhaps part of the reb_simulation struct). Any ideas/comments?

[hannorein]

Yeah. You're right. This would indeed involve changing the arguments for resolve(). My hope was to find a solution that doesn't require us to change the API.

Two simple work-arounds I could imagine:

1. Limit the number of mergers to one per timestep. If, for your case, mergers are rare and you are only interested in the evolution in a statistical sense, that would be fine.
2. Implement your own collision routine. If you're using the N^2 algorithm, it's really simple and you can just copy and paste the code from collision.c. This would also allow you to optimize the search a bit more (regarding issue #84).

Just FYI, I never thought much about merging particles when writing the collision routine. I mainly used it for Saturn's rings myself, where particles bounce off each other.

[hannorein]

Have a look at the above suggestion for a workaround.

[7ayushgupta]

@hannorein I am simulating the decay of the planetismals around a planet, and hence there are multiple collisions of the particles with the central planet. I want to remove the planetismals after such a collision. I cannot limit the number of mergers as the workaround shows. What should I actually do then?

[hannorein]

Welcome to the 4 year old thread ;-).

I'm not quite sure I understand the issue. You should be able to do something like in the example here: https://rebound.readthedocs.io/en/latest/c_example_mergers.html

If you can post a short snipped of code that shows the issue, that would help!

[7ayushgupta]

No problem @hannorein. I solved it by actually maintaining which particles are colliding within a particular timestep and then modifying the reb_remove code, so that it removes the particles sequentially.

[hannorein]

I'll close this for now.

I think the default (only one collision per particle per timestep for merging collisions) does make sense. If that's not the desired behaviour, one can easily implement a new merger routine. To do that, just copy the reb_collision_resolve_merge function in collision.c and change the first line of code.