different particles during collision and afterwards

[Findus23]

Okay, I preface this by the fact that this is most likely a stupid mistake by me. But I totally fail to find the issue and maybe this example helps someone in the future as it is more complex than https://rebound.readthedocs.io/en/latest/ipython/User_Defined_Collision_Resolve.html

The issue seems to be that sometimes two objects collide, a new one is formed, but some years later another object collides with one of the two initial ones even though they should no longer exist.

I'd love to provide a 30 line script that reproduces this, but unfortunately I can only reproduce this with my exact data.
I was still able to boil it down to a simpler script that only takes two minutes to run. But it is too large to paste it here, so it can be found here: https://gist.github.com/Findus23/c228493ea62bdd2034b13515a265f4b3

The basic idea would be the same as https://rebound.readthedocs.io/en/latest/ipython/User_Defined_Collision_Resolve.html, but it also tries to name the more massive particle target (it doesn't matter here, but in original program) and it always tries to replace the object with the lower index and drop the other one, so that N_active keeps working.

So after creating the new object, assigning it a new hash, I set the lower index of the two particles to it (a check of sim.particles shows that it has been assigned correctly) and delete the other one.
But weirdly it seems like afterwards (outside of the collision_resolve) it seems like the value goes back to what it was before:
https://gist.github.com/Findus23/c228493ea62bdd2034b13515a265f4b3#file-zz_output-txt-L18-L24
This allows a few years later another object to collide with 2486797600.

As said before, I think it is most likely that I do something wrong, but the fact that it also happens in this simplified version of the script mean that it probably more conceptual.

[hannorein]

Hm. I'm not very confident about the hash implementation for this use case. Is there an easy way to test if the problem occurs when you don't uses hashes? (e.g. by using a global variable instead)

[Findus23]

But what would I use instead to uniquely identify the particles (as indices might change if other objects collide)?

And I think the script never uses the hash for any logic, but just for the output.

[hannorein]

Fair enough.

Another potential non-rebound issue: Could your random number generator generate the same number twice? If you call it in your initialization routine and then later in the collision routine, you might not access the same seed.

[Findus23]

It definitely can. But isn't the chance that this happens in a simulation with at most 500 collisions 1/(2**32/500) which I think is 1/16 of the chance to be hit by a lightning in a year?

But this slightly bad feeling caused me to replace the random hash with incrementing integers in my main program (of course it is possible that I make the same mistake there with passing the object that also holds the current maximum)

[hannorein]

If for some reason the same seed gets used twice, then you get the same hash and the probability is more like 1/500.

[Findus23]

That makes sense and I agree that using a random hash isn't the best idea.

But the script also outputs all random hashes that are generated during the collisions and it never outputs the duplicate 2486797600.

[hannorein]

OK. That's not it then.

Is this an unexpected side-effect due to #494? If so, commenting out this line should fix it: https://github.com/hannorein/rebound/blob/main/src/integrator_mercurius.c#L366

[Findus23]

That might be possible. I upgraded to 3.16.0 about the same time this started to be weird (but I also did other changes at the same time, so I can't isolate it).

I quickly ran the linked script with 3.14.0 and I get other collisions, but at least no duplicates (https://gist.github.com/Findus23/c228493ea62bdd2034b13515a265f4b3#file-zz_output_3-14-0-txt)

[hannorein]

I think that makes sense. So maybe #494 should only be an option and turned off by default.

[Findus23]

I was wondering why this would affect my code as my testparticles do have a mass, but I just noticed that I lost the sim.testparticle_type = 1 when modifying the code recently.

[hannorein]

You're probably not the only one who might run into something related to this. So I'll make this feature optional and turn it off by default. Better to be slow and predictable, than fast and hard to debug.

[Findus23]

But I'm still not sure if there isn't a bug in #494 as @katvolk mentioned

if the test particles have mass, or the encounters involve two planets, then the planet evolution from the close-encounter routine is kept

and all particles in this scenario have a mass, so it should not have changed anything

[Rmelikyan]

Not related to the solution, but I do notice that you are not combining the two particles mass when you are merging the new planet.

I remember noticing that the conservation of momentum equation also scales the particles mass so it won't be exactly equal to p1.m + p2.m.

This may be something you want to check in on for the accuracy of your simulations!

[katvolk]

re #494, I think the check for test particle mass did depend on sim.testparticle_type = 1 being set (i.e., the check is whether there are close encounters between "active" particles or test particles of type 1, not actually explicitly a check for mass!

[Findus23]

Thanks @Rmelikyan,

You are indeed right. When I copied the simple example I forgot about this line (but it shouldn't change this issue).

In my real simulation I consider the mass and things are a bit more complex as my main work is on water (and stone) loss during collisions. If you are interested or have feedback about it, you can find it here: https://git.lw1.at/lw1/rebound-collisions/-/blob/0d2167cfdfd53997877025c3e8b69a864d6f191f/merge.py#L161-184

@katvolk Indeed when I set sim.testparticle_type = 1 (as it should be), I don't see any duplicates in a quick test run.

[Rmelikyan]

@Findus23 Yeah! Sorry to distract from the main point.

One thing I am finding while looking deeper at your gist example is this: if you search for the particle hash in index spot 12 (# 2486797600) in zz_output.txt you see that it is the target of your 3rd collision for that simulation.

In the next printout of the particle indices and hashes, you see that at index 12, the hash has changed, as expected, to the child hash.

BUT, the very next print out (after an integration step I think), the hash has reverted back! Is this a clue?

[Findus23]

@Rmelikyan No problem. There are definitely still many larger issues in my code. I haven't worked on something like that before and it will still take a few months until I have something that resembles a Master's thesis.

This is what I meant with

But weirdly it seems like afterwards (outside of the collision_resolve) it seems like the value goes back to what it was before:
https://gist.github.com/Findus23/c228493ea62bdd2034b13515a265f4b3#file-zz_output-txt-L18-L24

Unfortunately, the gist cuts off that part in the default view. But when printing the particles inside resolve() it is correct, but afterwards in the "normal" code, it reverted back.

[Rmelikyan]

Ok Good! I'm on the same page then 😂.

Maybe add another hash printout before the integrate call so we can be certain where the bug is. It could be occurring as the collision function returns or within the integrator call.

And my own question for you: Do you need to set the merged particle with a new hash? Can you not just record that particle xxx with hash yyy has encountered a collision? It seems by replacing the hash you might as well delete both p1 and p2 from the sim if not for the indexing drawbacks....

[Findus23]

@Rmelikyan

Isn't just before the integrate call exaclty the same as the one after the integrate call as the for loops back around directly after it?

for i in range(0, 500, 50):
    # <- here
    sim.integrate(i)
    print(i, sim.N)
    print([(p.index, p.hash.value) for p in sim.particles[:20]])
    # <- here
print("remaining", sim.N)

In theory one could keep hashes after collisions (e.g. the hash of the more massive body). But I think this only introduces more confusion as they aren't the same body any more afterwards (and e.g. have another water fraction). It also makes backtracing from the final body and plotting its mass over time easier if you only have to go up the collision tree.
I could in theory return 0 and add a completely new particle, but I feel safer that mercurius will handle it correctly (see #475 (comment)) this way (and as there are never more than two remnants after a collision in my case it doesn't matter)

[Rmelikyan]

You're right! my bad.

Ultimately i think it comes down to how you look at it. For me I probably wouldn't change the hash and would find a clear and sensical way to log the history of that particles collisions etc.

[hannorein]

Wow. What an active discussion! I love it!

Everything seems clear to me now. It is related to #494. The new hashes of the massive particles get overwritten after the close encounter because the massive particles should not have been affected by the collisions with massless particles. The possible solutions are:

• Add the missing sim.testparticle_type = 1 statement.
• Don't use hashes
• Make the modification to mercurius so massless test-particle encounters don't affect planet evolution #494 changes optional

[Findus23]

• Add the missing sim.testparticle_type = 1 statement.

That's the solution I will now be using as it is anyway what I wanted to do from the beginning