Question on removing particles in 'additional forces' function

[ryanwhite1]

Hi all,

I'm having problems with removing particles from a simulation within an additional forces function. What I'm trying to do is remove a particle if its eccentricity is larger than some critical value (0.8 in my case). I won't attach my exact code here (unless requested) as my function is quite long and the code overall relies on data from external scripts, but here is a rough outline of what I'm doing:

void disk_forces(struct reb_simulation* r){
    // this is my additional forces function
    const double G = r->G;
    int N = r->N;
    struct reb_particle* const particles = r->particles;
    struct reb_particle com = particles[0]; // calculate migration forces with respect to center of mass;

    for (int i = 1; i < N; i++){
        struct reb_particle* p = &(particles[i]); // get the particle
        int p_hash = p->hash;
        // get orbital properties of the particle
        const double dx = p->x-com.x;
        const double dy = p->y-com.y;
        const double dz = p->z-com.z;
        const double dvx = p->vx-com.vx;
        const double dvy = p->vy-com.vy;
        const double dvz = p->vz-com.vz;
        const double mass = p->m;
        double radius = sqrt(dx*dx + dy*dy + dz*dz);
        const double vr = (dx*dvx + dy*dvy + dz*dvz);   // dot product of v and r

        // ~80 other lines of code changing the particle's acceleration

        double mu = G*(com.m + mass);
        double vel = sqrt(dvx*dvx + dvy*dvy + dvz*dvz);

        // calculate eccentricity
        double ex = 1. / mu * ((vel*vel - mu / radius) * dx - vr * dvx);
        double ey = 1. / mu * ((vel*vel - mu / radius) * dy - vr * dvy);
        double ez = 1. / mu * ((vel*vel - mu / radius) * dz - vr * dvz);
        double e = sqrt(ex*ex + ey*ey + ez*ez);
        if (e > 0.8){               // if eccentricity >0.8, want to eject the particle by its hash
            printf("\nParticle ejected!, %.4e, %.4e\n", e, mass);
            reb_simulation_remove_particle_by_hash(r, p_hash, 1);
        }
}

The problem is that if a particle at index i is ejected, all particles from index i onwards are removed from the simulation and I'm not sure why.
I'm not very confident with C programming hence why I'm hesitant to label this a bug yet instead of my bad programming.
I could probably think of a work around to avoid ejecting particles in the additional forces function, but I'd like to know if there is an obvious error in what I'm doing! I'm happy to give more info as needed :)

Environment

• REBOUND Version: 4.0.2
• API interface: C
• Operating System: WSL2 Ubuntu

Physics
I'm simulating migration of massive particles within a custom accretion disk.
I've non-dimensionalised all of my units so that the characteristic length and timescales are ~1.
I'm simulating on the order of 10^1 particles, and for ~10^5 time units.
I'm also using the IAS15 integrator.

[ryanwhite1]

In writing this I'd actually thought of a solution! The solution was to check for ejections in a separate function called on in the heartbeat function instead of checking in the additional forces function.
I feel a bit silly, but I'll close this for posterity!

[hannorein]

Your solution is absolutely correct. You can't remove particles in the force functions. You can only do it once per timestep in the heartbeat function. Congrats to finding the solution yourself! 😉