Mercury crashes after a collision with the central particle

[jilkoval]

It seems that Mercury crashes after a collision with the central particle occurs. I think this is becuase Mercury (always) removes the particle that hits the central body (mercury6_2.for, L228) and this is not comunicated to amuse. A conflict than appears in the evolve_model of interface.py (L1394). A minimal example follows.

from amuse.units import units
from amuse.datamodel import Particles
from amuse.ext.solarsystem import new_solar_system
from amuse.community.mercury.interface import Mercury

solar_system = new_solar_system()

stone = Particles(1)
stone.position = (3.,0.,0.) | units.AU
stone.velocity = (-10.,0.,0.) | units.kms
stone.mass = 0.|units.kg
stone.name = "stone"
stone.raidus = 100.|units.m

mercury = Mercury(redirection="file", redirect_file="m_redirection.out")
mercury.particles.add_particles(solar_system)
mercury.particles.add_particles(stone)
mercury.commit_particles()
mercury.parameters.info_file = "m_info_file.out"
mercury.parameters.timestep = (1.|units.day)

mercury.evolve_model(1.|units.yr)
mercury.stop()

The output info file 'm_info_file.out' has to exist, otherwise the example stalls. I think this is related to the issue #95 (and the earlier #91 where the solution requires calling commit_parameters() to set up the output files).

[jilkoval]

Actually, the particles get removed from the code here (amuse_helpers.for, L951, in routine mal_hcon, which is used for the default integrator=10). After commenting out the line, the particle is not removed and the code gets to an infinite loop (since the particle is identified to hit the central again and again).

[ipelupessy]

I think the best solution would be to include collision detection - there was some work on this but needs to be restarted...

[jilkoval]

A possible workaround is to set the radius of the central body (rcen in the mercury/src files) -- that is the attribute value of central_particle.radius in the amuse script -- to 0 (actually also some small value might work, as long as mce_cent is not able to identify the collision through its interpolation). I think that setting central_particle.radius=0 should not have effect on anything else in the code (it plays a role for the spin angular momentum of the central body, but this is probably not concerned if ones does not want to include collisions at athe first place).

[ipelupessy]

workaround availalble