collision detection with bridge

[jilkoval]

I am having troubles with collision_detection stopping condition while using bridge. After detecting a collision, I remove one of the colliding particles and evolve the system further, but the same collision is detected again. It seems to me that .collision_detection.is_set() is not set to False for the evolve_model after the first collision.

I might be missing something basic here however. Any comments welcome! Please see example code below (I tested several gravity codes, the same problem appears for all).

from amuse.community.rebound.interface import Rebound
from amuse.community.huayno.interface import Huayno
from amuse.community.ph4.interface import ph4
from amuse.units import units, nbody_system
from amuse.datamodel import Particles
from amuse.couple import bridge
from amuse.ext.galactic_potentials import MWpotentialBovy2015

collinders = Particles(3)
collinders.mass = [5.0, 0.0, 0.0] | units.MSun
collinders.position = [(500.,0.,0.), (501.,0.,0.), (502.,0.,0.)] | units.AU
collinders.radius = [0.3, 0.3, 0.4] | units.AU
collinders.velocity = [(0.,0.,0.), (0.,0.,0.), (0.,0.,0.)] | units.kms

converter = nbody_system.nbody_to_si(1|units.MSun,1|units.AU)

#gravity = Rebound(converter)
#gravity.parameters.integrator = "whfast"
#gravity.parameters.timestep = 0.05|units.day

#gravity = Huayno(converter)
#gravity.parameters.timestep_parameter = 0.0001

gravity = ph4(converter)
gravity.parameters.timestep = 0.05|units.day

gravity.stopping_conditions.collision_detection.enable()
gravity.particles.add_particles(collinders)
gravity.commit_particles()

gravity_with_bridge = bridge.Bridge()
external_potential = MWpotentialBovy2015()
gravity_with_bridge.add_system(gravity, (external_potential,), do_sync=True)
gravity_with_bridge.timestep = 0.1|units.day

t_end = 60.|units.day
dt = 5.|units.day

t = 0.|units.day
while t < t_end:
    #gravity.evolve_model(t)
    print " >> evolving >", t, "n_particles =", len(gravity.particles)
    gravity_with_bridge.evolve_model(t)
    #print gravity.stopping_conditions.collision_detection.is_set()
    if gravity.stopping_conditions.collision_detection.is_set():
        #print gravity.stopping_conditions.collision_detection.is_set()
        t_collision = gravity_with_bridge.model_time
        n_particles_in_collision = 2*len(gravity.stopping_conditions.collision_detection.particles(0))
        print " xx collision at", t_collision, "n_particles_in_collisions =", n_particles_in_collision 
        if n_particles_in_collision == 0:
            print " !! only one particle in collision > break"
            break
        else:
            print gravity.stopping_conditions.collision_detection.particles(0)
            print gravity.stopping_conditions.collision_detection.particles(1)
            gravity.particles.remove_particles(gravity.stopping_conditions.collision_detection.particles(0))
            #gravity_with_bridge.particles.remove_particles(gravity.stopping_conditions.collision_detection.particles(0))
            #print gravity.particles
        continue
    t += dt

gravity.stop()
gravity_with_bridge.stop()

[jilkoval]

The gravity code needs to be synchronized with the bridge model_time after the collision: gravity.evolve_model(t_collision_bridge)

The complete fixed while loop:

while t < t_end:
    #gravity.evolve_model(t)
    print " >> evolving >", t.in_(units.day), "n_particles =", len(gravity.particles)
    gravity_with_bridge.evolve_model(t)
    #print gravity.stopping_conditions.collision_detection.is_set()
    if gravity.stopping_conditions.collision_detection.is_set():
        t_collision_bridge = gravity_with_bridge.model_time
        t_collision_gravity = gravity.model_time
        n_particles_in_collision = 2*len(gravity.stopping_conditions.collision_detection.particles(0))
        print " xx collision at", t_collision_gravity, "n_particles_in_collisions =", n_particles_in_collision 
        print "    t_gravity=", gravity.model_time.in_(units.day), "t_bridge=", gravity_with_bridge.model_time.in_(units.day)
        if n_particles_in_collision == 0:
            print " !! only one particle in collision > break"
            break
        else:
            print gravity.stopping_conditions.collision_detection.particles(0)
            print gravity.stopping_conditions.collision_detection.particles(1)
            gravity.particles.remove_particles(gravity.stopping_conditions.collision_detection.particles(0))
            #gravity_with_bridge.particles.remove_particles(gravity.stopping_conditions.collision_detection.particles(0))
            #print gravity.particles
        if (t_collision_gravity > t_collision_bridge):
            t=t_collision_gravity
        else:
            gravity.evolve_model(t_collision_bridge)
        continue
    t += dt