Issues with MERCURIUS and collision handling #475
Comments
|
Maybe I'm missing something, but I think what you want to do is implement this as a collision_resolve function. When the exception is thrown, the simulation stops immediately. That's nice for a quick test to see if particles collide, but not really useful if you want to continue with the integration. You'd need to piece together the particle data from some internal arrays, which could be messy. In short, use a collision_resolve function pointer and implement your logic there. |
|
Apart from the things I mentioned above (and a few older misunderstandings of mine #418) doing all the collision handling in Python worked really well for me (the part that is calculating the water fractions and then removing and adding the particles afterwards is of course missing in above example). And halting the simulation on collisions, getting some water retention estimations via some external methods and then continuing with the new particles is exactly what I would want to do. With the estimation code already being in Python I assumed it would be most straight-forward (and avoid additional errors in the implementation) to keep the merging code in Python. |
|
That's not a problem! The collision_resolve function pointer can point to a python function! |
|
Okay, I totally misunderstood the "pointer". So instead of def merge_particles(sim: Simulation):
...
try:
sim.integrate(t)
except Collision:
merge_particles(sim)you suggest instead to directly run def merge_particles(sim: POINTER_REB_SIM, collision: reb_collision)-> int:
...
sim.collision_resolve = merge_particles
sim.integrate(t)? If I understand it correctly, this means I no longer need the hack to find the two particles as collision.p1 and collision.p2 contain their indices, but I can no longer access the original Python And even the collision.p1 don't seem helpful as I can no longer access |
|
You want to access the sim pointer that is given to your function. Because it's a pointer, you need to use the syntax You can change the two particles involved in the collisions. Depending on the return value of your merge_particle function, you can keep both, one, the other, or neither of the particles. This way, MERCURIUS will make all the necessary changes to the various internal arrays so that the simulation can continue. (But make sure to run a few quick tests) If you want to do something more complicated, for example adding a whole swarm of new particles when a collisions occurs, you need to dig deep into the code and it can be complicated. If that's what you want to do, it might be best to chat over zoom to brain storm what the best way to implement this would be. |
|
Many thanks, I am starting to understand. So it is more of a matter of def merge_particles(sim_p: POINTER_REB_SIM, collision: reb_collision) -> int:
sim: Simulation = sim_p.contents
p1: Particle = sim.particles[collision.p1],
p2: Particle = sim.particles[collision.p2]
... # modify p1 and p2 to fit the collision outcome (but I have no way to access some_var for this)
return 0 # or 1 or 2 to delete one of them
some_var = SomeObject()
sim.collision_resolve = merge_particles
sim.integrate(t)In theory this should be enough for my use case (there will never be more than two remaining particles). I'll have to think a bit more about if this might work (one issue I still have is that I can't access variables from outside the integration) and do some tests if everything works as expected inside the But at least this solves the third issue and potentially also another unrelated one (having to remove and reorder all particles after a collision in case a semi-active particle collides with an active one and N_active changes therefore) |
|
You should be able to access variables: |
|
Okay, my issue was that this wouldn't work if from merge import really_merge_particles
def merge_particles(sim_p: POINTER_REB_SIM, collision: reb_collision) -> int:
return really_merge_particles(sim_p, collision, some_var)
some_var = SomeObject()
sim.collision_resolve = merge_particles
sim.integrate(t)One last question before I'll try out all of this in a few days: I assume rebound doesn't use the hash anywhere internally, but rather the pointers to p1 and p2, so it should be fine if their hash changes during the collision, right? (I'll check just in case) |
|
Yes, you can change the hash at any time without consequences. I'll close this issue for now, but feel free to reopen it if you run into other issues! |
Hello,
This time I have stumbled across a bit more complex issue, so this one will be a bit longer (I still tried to break it down to a minimal example).
A bit of background: I'm currently working on my Master's Thesis and my main plan would be to reproduce https://ui.adsabs.harvard.edu/abs/2020A%26A...634A..76B/abstract, but with a bit simplified collision handling (letting rebound do most of the work) and interpolation of a collision catalogue for water retention fractions instead of dedicated SPH simulations.
So the core part would be simulating late planetary formation and tracking water fractions over collisions. As this takes a long time, but also has a lot of collisions MERCURIUS seems like the obvious solution.
So this is the basis of the program that runs the simulation and on an encounter calculates collision speeds and angle
uncollapse for data.py
Sorry for the rather verbose file, but I wanted to keep it reproducable independent on any code that reads input files.
The program should run for a few seconds and output something like this:
Now I have a few issues with the program:
sim.ri_mercurius.modealways returns 0 (so whfast) even in situations where I assume that mercurius has already switched to IAS15. Also sim.dt stays the same (and shows the initial whfast dt even during collision), but I assume this is intentional.p.lastcollision == sim.tas shown in CloseEncounters.ipynb doesn't seem to work, butp.lastcollision >= sim.t - sim.dtseems to reliably find the two colliding bodies.Am I missing some obvious limitation or am I doing something wrong? I tried increasing
sim.ri_mercurius.hillfacbut even with a higher value (e.g. 30), but simulation takes a lot longer (which makes sense as the particles are quite close to each others; weirdly the program gets even quite a bit slower than just using IAS15), but I don't think it changes the distance when the simulation halts.The text was updated successfully, but these errors were encountered: