hannorein / rebound

💫 An open-source multi-purpose N-body code.
https://rebound.readthedocs.io/
GNU General Public License v3.0
861 stars 221 forks source link

lastcolision time #358

Closed MilevaE closed 6 years ago

MilevaE commented 6 years ago

Hi,

I'm having problems with the lastcolision function. I have defined an integration like that:

end = 1000 #final time
Noutputs = 1000 #outputs 
times = np.linspace(0.,end,Noutputs,endpoint=False)
totalN = np.zeros(Noutputs)
for i, timex in enumerate(times):
    try:
        sim.integrate(timex)
        totalN[i] = sim.N_active
    except rebound.Collision:
        collided = []
        for p in sim.particles:
            if p.lastcollision == sim.t:
                collided.append(p.index)
            #print sim.t
        water_update(sim)

So, what I'm doing here is to define my own collision resolve rutine, after a collision exception is caught by the exception. My problem is that depending if I use IAS15 or Mercurius, lastcollision is not working. In both the exception is caught, so I guess that the collision happens and both integrators "see" it, but using Mercurius, lastcollision is 0 instead the time when the collision happened.

Any idea why I could be having this problem? suggestions to solve it?

I would prefer use Mercurius because is much faster than IAS15..

best, Mil

hannorein commented 6 years ago

Thanks for catching that. I just pushed an update which should fix the issue. Feel free to reopen the issue if that's not the case!

MilevaE commented 6 years ago

Hi Hanno,

thanks a lot for your super-fast answer. I upgraded and works partially. Now seems that both values sim.t and lastcollision have different accuracy, so they don't satisfy the condition if p.lastcollision == sim.t:

here what I obtain:

sim.t = 695.726
lastcollision = 695.725297086

So, the collision is caught, but I still can not identify which particles were involved.

Thanks a lot for your help

hannorein commented 6 years ago

That's expected as Mercurius takes many small timesteps during close encounters. You'd need to add some extra logic in your collision routine to say if the collision happened within the last timestep.

MilevaE commented 6 years ago

ok! thanks a lot!