Version that correctly flags the BHs to swallow and removes them before dying...

Version that correctly flags the BHs to swallow and removes them before dying in the time-step calculation function.

src/black_holes/EAGLE/black_holes_iact.h

@@ -230,6 +230,8 @@ runner_iact_nonsym_bh_bh_swallow(const float r2, const float *dx,
      * smoothing length */
     if (v2_pec < G * M / (kernel_gamma * h)) {
 
+      // MATTHIEU: Also add distance check: factor * Plummer softening
+
       /* This particle is swallowed by the BH with the largest ID of all the
        * candidates wanting to swallow it */
       if (bj->merger_data.swallow_id < bi->id) {

src/cell.c

@@ -1988,7 +1988,8 @@ void cell_clean_links(struct cell *c, void *data) {
   c->stars.feedback = NULL;
   c->black_holes.density = NULL;
   c->black_holes.swallow = NULL;
-  c->black_holes.do_swallow = NULL;
+  c->black_holes.do_gas_swallow = NULL;
+  c->black_holes.do_bh_swallow = NULL;
   c->black_holes.feedback = NULL;
 }
 
@@ -4001,7 +4002,7 @@ int cell_unskip_black_holes_tasks(struct cell *c, struct scheduler *s) {
   }
 
   /* Un-skip the swallow tasks involved with this cell. */
-  for (struct link *l = c->black_holes.do_swallow; l != NULL; l = l->next) {
+  for (struct link *l = c->black_holes.do_gas_swallow; l != NULL; l = l->next) {
     struct task *t = l->t;
     struct cell *ci = t->ci;
     struct cell *cj = t->cj;
@@ -4049,23 +4050,14 @@ int cell_unskip_black_holes_tasks(struct cell *c, struct scheduler *s) {
   /* Unskip all the other task types. */
   if (c->nodeID == nodeID && cell_is_active_black_holes(c, e)) {
 
-    /* for (struct link *l = c->black_holes.swallow; l != NULL; l = l->next)
-     */
-    /*   scheduler_activate(s, l->t); */
-    /* for (struct link *l = c->black_holes.do_swallow; l != NULL; l =
-     * l->next)
-     */
-    /*   scheduler_activate(s, l->t); */
-    /* for (struct link *l = c->black_holes.feedback; l != NULL; l = l->next)
-     */
-    /*   scheduler_activate(s, l->t); */
-
     if (c->black_holes.density_ghost != NULL)
       scheduler_activate(s, c->black_holes.density_ghost);
     if (c->black_holes.swallow_ghost[0] != NULL)
       scheduler_activate(s, c->black_holes.swallow_ghost[0]);
     if (c->black_holes.swallow_ghost[1] != NULL)
       scheduler_activate(s, c->black_holes.swallow_ghost[1]);
+    if (c->black_holes.swallow_ghost[2] != NULL)
+      scheduler_activate(s, c->black_holes.swallow_ghost[2]);
     if (c->black_holes.black_holes_in != NULL)
       scheduler_activate(s, c->black_holes.black_holes_in);
     if (c->black_holes.black_holes_out != NULL)

src/cell.h

@@ -651,7 +651,7 @@ struct cell {
     struct task *density_ghost;
 
     /*! The star ghost task itself */
-    struct task *swallow_ghost[2];
+    struct task *swallow_ghost[3];
 
     /*! Linked list of the tasks computing this cell's BH density. */
     struct link *density;
@@ -661,7 +661,10 @@ struct cell {
     struct link *swallow;
 
     /*! Linked list of the tasks processing the particles to swallow */
-    struct link *do_swallow;
+    struct link *do_gas_swallow;
+
+    /*! Linked list of the tasks processing the particles to swallow */
+    struct link *do_bh_swallow;
 
     /*! Linked list of the tasks computing this cell's BH feedback. */
     struct link *feedback;

src/engine.c

@@ -3327,11 +3327,13 @@ void engine_skip_force_and_kick(struct engine *e) {
         t->type == task_type_extra_ghost ||
         t->type == task_type_bh_swallow_ghost1 ||
         t->type == task_type_bh_swallow_ghost2 ||
+        t->type == task_type_bh_swallow_ghost3 ||
         t->subtype == task_subtype_gradient ||
         t->subtype == task_subtype_stars_feedback ||
         t->subtype == task_subtype_bh_feedback ||
         t->subtype == task_subtype_bh_swallow ||
-        t->subtype == task_subtype_do_swallow ||
+        t->subtype == task_subtype_do_gas_swallow ||
+        t->subtype == task_subtype_do_bh_swallow ||
         t->subtype == task_subtype_bpart_rho ||
         t->subtype == task_subtype_part_swallow ||
         t->subtype == task_subtype_bpart_swallow ||

src/engine_maketasks.c

@@ -672,11 +672,10 @@ void engine_addtasks_recv_black_holes(struct engine *e, struct cell *c,
 
     for (struct link *l = c->black_holes.swallow; l != NULL; l = l->next) {
       scheduler_addunlock(s, t_rho, l->t);
-      // scheduler_addunlock(s, l->t, t_swallow);
       scheduler_addunlock(s, l->t, t_gas_swallow);
     }
-    for (struct link *l = c->black_holes.do_swallow; l != NULL; l = l->next) {
-      // scheduler_addunlock(s, t_swallow, l->t);
+    for (struct link *l = c->black_holes.do_gas_swallow; l != NULL;
+         l = l->next) {
       scheduler_addunlock(s, t_gas_swallow, l->t);
       scheduler_addunlock(s, l->t, t_feedback);
     }
@@ -1093,8 +1092,12 @@ void engine_make_hierarchical_tasks_hydro(struct engine *e, struct cell *c) {
         c->black_holes.density_ghost = scheduler_addtask(
             s, task_type_bh_density_ghost, task_subtype_none, 0, 0, c, NULL);
 
-        c->black_holes.swallow_ghost[1] = scheduler_addtask(
-            s, task_type_bh_swallow_ghost2, task_subtype_none, 0, 0, c, NULL);
+        c->black_holes.swallow_ghost[1] =
+            scheduler_addtask(s, task_type_bh_swallow_ghost2, task_subtype_none,
+                              0, /* implicit =*/1, c, NULL);
+
+        c->black_holes.swallow_ghost[2] = scheduler_addtask(
+            s, task_type_bh_swallow_ghost3, task_subtype_none, 0, 0, c, NULL);
 
         scheduler_addunlock(s, c->super->kick2, c->black_holes.black_holes_in);
         scheduler_addunlock(s, c->black_holes.black_holes_out,
@@ -1718,7 +1721,8 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
   struct task *t_star_feedback = NULL;
   struct task *t_bh_density = NULL;
   struct task *t_bh_swallow = NULL;
-  struct task *t_do_swallow = NULL;
+  struct task *t_do_gas_swallow = NULL;
+  struct task *t_do_bh_swallow = NULL;
   struct task *t_bh_feedback = NULL;
 
   for (int ind = 0; ind < num_elements; ind++) {
@@ -1774,8 +1778,12 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
             sched, task_type_self, task_subtype_bh_density, flags, 0, ci, NULL);
         t_bh_swallow = scheduler_addtask(
             sched, task_type_self, task_subtype_bh_swallow, flags, 0, ci, NULL);
-        t_do_swallow = scheduler_addtask(
-            sched, task_type_self, task_subtype_do_swallow, flags, 0, ci, NULL);
+        t_do_gas_swallow =
+            scheduler_addtask(sched, task_type_self,
+                              task_subtype_do_gas_swallow, flags, 0, ci, NULL);
+        t_do_bh_swallow =
+            scheduler_addtask(sched, task_type_self, task_subtype_do_bh_swallow,
+                              flags, 0, ci, NULL);
         t_bh_feedback =
             scheduler_addtask(sched, task_type_self, task_subtype_bh_feedback,
                               flags, 0, ci, NULL);
@@ -1793,7 +1801,8 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
       if (with_black_holes && bcount_i > 0) {
         engine_addlink(e, &ci->black_holes.density, t_bh_density);
         engine_addlink(e, &ci->black_holes.swallow, t_bh_swallow);
-        engine_addlink(e, &ci->black_holes.do_swallow, t_do_swallow);
+        engine_addlink(e, &ci->black_holes.do_gas_swallow, t_do_gas_swallow);
+        engine_addlink(e, &ci->black_holes.do_bh_swallow, t_do_bh_swallow);
         engine_addlink(e, &ci->black_holes.feedback, t_bh_feedback);
       }
 
@@ -1851,13 +1860,20 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
         scheduler_addunlock(sched, t_bh_swallow,
                             ci->hydro.super->black_holes.swallow_ghost[0]);
 
-        scheduler_addunlock(
-            sched, ci->hydro.super->black_holes.swallow_ghost[0], t_do_swallow);
-        scheduler_addunlock(sched, t_do_swallow,
+        scheduler_addunlock(sched,
+                            ci->hydro.super->black_holes.swallow_ghost[0],
+                            t_do_gas_swallow);
+        scheduler_addunlock(sched, t_do_gas_swallow,
                             ci->hydro.super->black_holes.swallow_ghost[1]);
 
         scheduler_addunlock(sched,
                             ci->hydro.super->black_holes.swallow_ghost[1],
+                            t_do_bh_swallow);
+        scheduler_addunlock(sched, t_do_bh_swallow,
+                            ci->hydro.super->black_holes.swallow_ghost[2]);
+
+        scheduler_addunlock(sched,
+                            ci->hydro.super->black_holes.swallow_ghost[2],
                             t_bh_feedback);
         scheduler_addunlock(sched, t_bh_feedback,
                             ci->hydro.super->black_holes.black_holes_out);
@@ -1916,8 +1932,12 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
             sched, task_type_pair, task_subtype_bh_density, flags, 0, ci, cj);
         t_bh_swallow = scheduler_addtask(
             sched, task_type_pair, task_subtype_bh_swallow, flags, 0, ci, cj);
-        t_do_swallow = scheduler_addtask(
-            sched, task_type_pair, task_subtype_do_swallow, flags, 0, ci, cj);
+        t_do_gas_swallow =
+            scheduler_addtask(sched, task_type_pair,
+                              task_subtype_do_gas_swallow, flags, 0, ci, cj);
+        t_do_bh_swallow =
+            scheduler_addtask(sched, task_type_pair, task_subtype_do_bh_swallow,
+                              flags, 0, ci, cj);
         t_bh_feedback = scheduler_addtask(
             sched, task_type_pair, task_subtype_bh_feedback, flags, 0, ci, cj);
       }
@@ -1939,8 +1959,10 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
         engine_addlink(e, &cj->black_holes.density, t_bh_density);
         engine_addlink(e, &ci->black_holes.swallow, t_bh_swallow);
         engine_addlink(e, &cj->black_holes.swallow, t_bh_swallow);
-        engine_addlink(e, &ci->black_holes.do_swallow, t_do_swallow);
-        engine_addlink(e, &cj->black_holes.do_swallow, t_do_swallow);
+        engine_addlink(e, &ci->black_holes.do_gas_swallow, t_do_gas_swallow);
+        engine_addlink(e, &cj->black_holes.do_gas_swallow, t_do_gas_swallow);
+        engine_addlink(e, &ci->black_holes.do_bh_swallow, t_do_bh_swallow);
+        engine_addlink(e, &cj->black_holes.do_bh_swallow, t_do_bh_swallow);
         engine_addlink(e, &ci->black_holes.feedback, t_bh_feedback);
         engine_addlink(e, &cj->black_holes.feedback, t_bh_feedback);
       }
@@ -2030,12 +2052,18 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
 
           scheduler_addunlock(sched,
                               ci->hydro.super->black_holes.swallow_ghost[0],
-                              t_do_swallow);
-          scheduler_addunlock(sched, t_do_swallow,
+                              t_do_gas_swallow);
+          scheduler_addunlock(sched, t_do_gas_swallow,
                               ci->hydro.super->black_holes.swallow_ghost[1]);
 
           scheduler_addunlock(sched,
                               ci->hydro.super->black_holes.swallow_ghost[1],
+                              t_do_bh_swallow);
+          scheduler_addunlock(sched, t_do_bh_swallow,
+                              ci->hydro.super->black_holes.swallow_ghost[2]);
+
+          scheduler_addunlock(sched,
+                              ci->hydro.super->black_holes.swallow_ghost[2],
                               t_bh_feedback);
           scheduler_addunlock(sched, t_bh_feedback,
                               ci->hydro.super->black_holes.black_holes_out);
@@ -2102,12 +2130,18 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
 
             scheduler_addunlock(sched,
                                 cj->hydro.super->black_holes.swallow_ghost[0],
-                                t_do_swallow);
-            scheduler_addunlock(sched, t_do_swallow,
+                                t_do_gas_swallow);
+            scheduler_addunlock(sched, t_do_gas_swallow,
                                 cj->hydro.super->black_holes.swallow_ghost[1]);
 
             scheduler_addunlock(sched,
                                 cj->hydro.super->black_holes.swallow_ghost[1],
+                                t_do_bh_swallow);
+            scheduler_addunlock(sched, t_do_bh_swallow,
+                                cj->hydro.super->black_holes.swallow_ghost[2]);
+
+            scheduler_addunlock(sched,
+                                cj->hydro.super->black_holes.swallow_ghost[2],
                                 t_bh_feedback);
             scheduler_addunlock(sched, t_bh_feedback,
                                 cj->hydro.super->black_holes.black_holes_out);
@@ -2172,9 +2206,13 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
             scheduler_addtask(sched, task_type_sub_self,
                               task_subtype_bh_swallow, flags, 0, ci, NULL);
 
-        t_do_swallow =
+        t_do_gas_swallow =
+            scheduler_addtask(sched, task_type_sub_self,
+                              task_subtype_do_gas_swallow, flags, 0, ci, NULL);
+
+        t_do_bh_swallow =
             scheduler_addtask(sched, task_type_sub_self,
-                              task_subtype_do_swallow, flags, 0, ci, NULL);
+                              task_subtype_do_bh_swallow, flags, 0, ci, NULL);
 
         t_bh_feedback =
             scheduler_addtask(sched, task_type_sub_self,
@@ -2193,7 +2231,8 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
       if (with_black_holes && bcount_i > 0) {
         engine_addlink(e, &ci->black_holes.density, t_bh_density);
         engine_addlink(e, &ci->black_holes.swallow, t_bh_swallow);
-        engine_addlink(e, &ci->black_holes.do_swallow, t_do_swallow);
+        engine_addlink(e, &ci->black_holes.do_gas_swallow, t_do_gas_swallow);
+        engine_addlink(e, &ci->black_holes.do_bh_swallow, t_do_bh_swallow);
         engine_addlink(e, &ci->black_holes.feedback, t_bh_feedback);
       }
 
@@ -2257,13 +2296,20 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
         scheduler_addunlock(sched, t_bh_swallow,
                             ci->hydro.super->black_holes.swallow_ghost[0]);
 
-        scheduler_addunlock(
-            sched, ci->hydro.super->black_holes.swallow_ghost[0], t_do_swallow);
-        scheduler_addunlock(sched, t_do_swallow,
+        scheduler_addunlock(sched,
+                            ci->hydro.super->black_holes.swallow_ghost[0],
+                            t_do_gas_swallow);
+        scheduler_addunlock(sched, t_do_gas_swallow,
                             ci->hydro.super->black_holes.swallow_ghost[1]);
 
         scheduler_addunlock(sched,
                             ci->hydro.super->black_holes.swallow_ghost[1],
+                            t_do_bh_swallow);
+        scheduler_addunlock(sched, t_do_bh_swallow,
+                            ci->hydro.super->black_holes.swallow_ghost[2]);
+
+        scheduler_addunlock(sched,
+                            ci->hydro.super->black_holes.swallow_ghost[2],
                             t_bh_feedback);
         scheduler_addunlock(sched, t_bh_feedback,
                             ci->hydro.super->black_holes.black_holes_out);
@@ -2326,9 +2372,12 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
         t_bh_swallow =
             scheduler_addtask(sched, task_type_sub_pair,
                               task_subtype_bh_swallow, flags, 0, ci, cj);
-        t_do_swallow =
+        t_do_gas_swallow =
             scheduler_addtask(sched, task_type_sub_pair,
-                              task_subtype_do_swallow, flags, 0, ci, cj);
+                              task_subtype_do_gas_swallow, flags, 0, ci, cj);
+        t_do_bh_swallow =
+            scheduler_addtask(sched, task_type_sub_pair,
+                              task_subtype_do_bh_swallow, flags, 0, ci, cj);
         t_bh_feedback =
             scheduler_addtask(sched, task_type_sub_pair,
                               task_subtype_bh_feedback, flags, 0, ci, cj);
@@ -2351,8 +2400,10 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
         engine_addlink(e, &cj->black_holes.density, t_bh_density);
         engine_addlink(e, &ci->black_holes.swallow, t_bh_swallow);
         engine_addlink(e, &cj->black_holes.swallow, t_bh_swallow);
-        engine_addlink(e, &ci->black_holes.do_swallow, t_do_swallow);
-        engine_addlink(e, &cj->black_holes.do_swallow, t_do_swallow);
+        engine_addlink(e, &ci->black_holes.do_gas_swallow, t_do_gas_swallow);
+        engine_addlink(e, &cj->black_holes.do_gas_swallow, t_do_gas_swallow);
+        engine_addlink(e, &ci->black_holes.do_bh_swallow, t_do_bh_swallow);
+        engine_addlink(e, &cj->black_holes.do_bh_swallow, t_do_bh_swallow);
         engine_addlink(e, &ci->black_holes.feedback, t_bh_feedback);
         engine_addlink(e, &cj->black_holes.feedback, t_bh_feedback);
       }
@@ -2441,12 +2492,18 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
 
           scheduler_addunlock(sched,
                               ci->hydro.super->black_holes.swallow_ghost[0],
-                              t_do_swallow);
-          scheduler_addunlock(sched, t_do_swallow,
+                              t_do_gas_swallow);
+          scheduler_addunlock(sched, t_do_gas_swallow,
                               ci->hydro.super->black_holes.swallow_ghost[1]);
 
           scheduler_addunlock(sched,
                               ci->hydro.super->black_holes.swallow_ghost[1],
+                              t_do_bh_swallow);
+          scheduler_addunlock(sched, t_do_bh_swallow,
+                              ci->hydro.super->black_holes.swallow_ghost[2]);
+
+          scheduler_addunlock(sched,
+                              ci->hydro.super->black_holes.swallow_ghost[2],
                               t_bh_feedback);
           scheduler_addunlock(sched, t_bh_feedback,
                               ci->hydro.super->black_holes.black_holes_out);
@@ -2515,12 +2572,18 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
 
             scheduler_addunlock(sched,
                                 cj->hydro.super->black_holes.swallow_ghost[0],
-                                t_do_swallow);
-            scheduler_addunlock(sched, t_do_swallow,
+                                t_do_gas_swallow);
+            scheduler_addunlock(sched, t_do_gas_swallow,
                                 cj->hydro.super->black_holes.swallow_ghost[1]);
 
             scheduler_addunlock(sched,
                                 cj->hydro.super->black_holes.swallow_ghost[1],
+                                t_do_bh_swallow);
+            scheduler_addunlock(sched, t_do_bh_swallow,
+                                cj->hydro.super->black_holes.swallow_ghost[2]);
+
+            scheduler_addunlock(sched,
+                                cj->hydro.super->black_holes.swallow_ghost[2],
                                 t_bh_feedback);
             scheduler_addunlock(sched, t_bh_feedback,
                                 cj->hydro.super->black_holes.black_holes_out);

src/engine_marktasks.c

@@ -207,14 +207,26 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
       }
 
       else if (t_type == task_type_self &&
-               t_subtype == task_subtype_do_swallow) {
+               t_subtype == task_subtype_do_gas_swallow) {
         if (ci_active_black_holes) {
           scheduler_activate(s, t);
         }
       }
 
       else if (t_type == task_type_sub_self &&
-               t_subtype == task_subtype_do_swallow) {
+               t_subtype == task_subtype_do_gas_swallow) {
+        if (ci_active_black_holes) scheduler_activate(s, t);
+      }
+
+      else if (t_type == task_type_self &&
+               t_subtype == task_subtype_do_bh_swallow) {
+        if (ci_active_black_holes) {
+          scheduler_activate(s, t);
+        }
+      }
+
+      else if (t_type == task_type_sub_self &&
+               t_subtype == task_subtype_do_bh_swallow) {
         if (ci_active_black_holes) scheduler_activate(s, t);
       }
 
@@ -404,7 +416,7 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
       /* Black_Holes density */
       else if ((t_subtype == task_subtype_bh_density ||
                 t_subtype == task_subtype_bh_swallow ||
-                t_subtype == task_subtype_do_swallow ||
+                t_subtype == task_subtype_do_gas_swallow ||
                 t_subtype == task_subtype_bh_feedback) &&
                (ci_active_black_holes || cj_active_black_holes) &&
                (ci_nodeID == nodeID || cj_nodeID == nodeID)) {
@@ -892,7 +904,8 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
     /* Black hole ghost tasks ? */
     else if (t_type == task_type_bh_density_ghost ||
              t_type == task_type_bh_swallow_ghost1 ||
-             t_type == task_type_bh_swallow_ghost2) {
+             t_type == task_type_bh_swallow_ghost2 ||
+             t_type == task_type_bh_swallow_ghost3) {
       if (cell_is_active_black_holes(t->ci, e)) scheduler_activate(s, t);
     }
 

src/runner.c

@@ -3198,20 +3198,23 @@ void runner_do_timestep(struct runner *r, struct cell *c, int timer) {
         } else { /* gpart is inactive */
 
           /* Count the number of inhibited particles */
-          if (gpart_is_inhibited(gp, e)) g_inhibited++;
+          if (gpart_is_inhibited(gp, e)) {
+            g_inhibited++;
+          } else {
 
-          const integertime_t ti_end =
-              get_integer_time_end(ti_current, gp->time_bin);
+            const integertime_t ti_end =
+                get_integer_time_end(ti_current, gp->time_bin);
 
-          /* What is the next sync-point ? */
-          ti_gravity_end_min = min(ti_end, ti_gravity_end_min);
-          ti_gravity_end_max = max(ti_end, ti_gravity_end_max);
+            /* What is the next sync-point ? */
+            ti_gravity_end_min = min(ti_end, ti_gravity_end_min);
+            ti_gravity_end_max = max(ti_end, ti_gravity_end_max);
 
-          const integertime_t ti_beg =
-              get_integer_time_begin(ti_current + 1, gp->time_bin);
+            const integertime_t ti_beg =
+                get_integer_time_begin(ti_current + 1, gp->time_bin);
 
-          /* What is the next starting point for this cell ? */
-          ti_gravity_beg_max = max(ti_beg, ti_gravity_beg_max);
+            /* What is the next starting point for this cell ? */
+            ti_gravity_beg_max = max(ti_beg, ti_gravity_beg_max);
+          }
         }
       }
     }
@@ -3319,22 +3322,26 @@ void runner_do_timestep(struct runner *r, struct cell *c, int timer) {
       } else {
 
         /* Count the number of inhibited particles */
-        if (bpart_is_inhibited(bp, e)) ++b_inhibited;
+        if (bpart_is_inhibited(bp, e)) {
+          message("BH inhibited!");
+          ++b_inhibited;
+        } else {
 
-        const integertime_t ti_end =
-            get_integer_time_end(ti_current, bp->time_bin);
+          const integertime_t ti_end =
+              get_integer_time_end(ti_current, bp->time_bin);
 
-        const integertime_t ti_beg =
-            get_integer_time_begin(ti_current + 1, bp->time_bin);
+          const integertime_t ti_beg =
+              get_integer_time_begin(ti_current + 1, bp->time_bin);
 
-        ti_black_holes_end_min = min(ti_end, ti_black_holes_end_min);
-        ti_black_holes_end_max = max(ti_end, ti_black_holes_end_max);
-        ti_gravity_end_min = min(ti_end, ti_gravity_end_min);
-        ti_gravity_end_max = max(ti_end, ti_gravity_end_max);
+          ti_black_holes_end_min = min(ti_end, ti_black_holes_end_min);
+          ti_black_holes_end_max = max(ti_end, ti_black_holes_end_max);
+          ti_gravity_end_min = min(ti_end, ti_gravity_end_min);
+          ti_gravity_end_max = max(ti_end, ti_gravity_end_max);
 
-        /* What is the next starting point for this cell ? */
-        ti_black_holes_beg_max = max(ti_beg, ti_black_holes_beg_max);
-        ti_gravity_beg_max = max(ti_beg, ti_gravity_beg_max);
+          /* What is the next starting point for this cell ? */
+          ti_black_holes_beg_max = max(ti_beg, ti_black_holes_beg_max);
+          ti_gravity_beg_max = max(ti_beg, ti_gravity_beg_max);
+        }
       }
     }
 
@@ -3731,7 +3738,7 @@ void runner_do_end_grav_force(struct runner *r, struct cell *c, int timer) {
  * @param c The #cell.
  * @param timer Are we timing this?
  */
-void runner_do_swallow(struct runner *r, struct cell *c, int timer) {
+void runner_do_gas_swallow(struct runner *r, struct cell *c, int timer) {
 
   struct engine *e = r->e;
   struct space *s = e->s;
@@ -3759,7 +3766,7 @@ void runner_do_swallow(struct runner *r, struct cell *c, int timer) {
       if (c->progeny[k] != NULL) {
         struct cell *restrict cp = c->progeny[k];
 
-        runner_do_swallow(r, cp, 0);
+        runner_do_gas_swallow(r, cp, 0);
       }
     }
   } else {
@@ -3883,7 +3890,7 @@ void runner_do_swallow(struct runner *r, struct cell *c, int timer) {
  * @param c The #cell.
  * @param timer Are we timing this?
  */
-void runner_do_swallow_self(struct runner *r, struct cell *c, int timer) {
+void runner_do_gas_swallow_self(struct runner *r, struct cell *c, int timer) {
 
 #ifdef SWIFT_DEBUG_CHECKS
   if (c->nodeID != r->e->nodeID) error("Running self task on foreign node");
@@ -3891,7 +3898,7 @@ void runner_do_swallow_self(struct runner *r, struct cell *c, int timer) {
     error("Running self task on inactive cell");
 #endif
 
-  runner_do_swallow(r, c, timer);
+  runner_do_gas_swallow(r, c, timer);
 }
 
 /**
@@ -3902,8 +3909,8 @@ void runner_do_swallow_self(struct runner *r, struct cell *c, int timer) {
  * @param cj Second #cell.
  * @param timer Are we timing this?
  */
-void runner_do_swallow_pair(struct runner *r, struct cell *ci, struct cell *cj,
-                            int timer) {
+void runner_do_gas_swallow_pair(struct runner *r, struct cell *ci,
+                                struct cell *cj, int timer) {
 
   const struct engine *e = r->e;
 
@@ -3916,8 +3923,213 @@ void runner_do_swallow_pair(struct runner *r, struct cell *ci, struct cell *cj,
 
   /* Run the swallowing loop only in the cell that is the neighbour of the
    * active BH */
-  if (cell_is_active_black_holes(cj, e)) runner_do_swallow(r, ci, timer);
-  if (cell_is_active_black_holes(ci, e)) runner_do_swallow(r, cj, timer);
+  if (cell_is_active_black_holes(cj, e)) runner_do_gas_swallow(r, ci, timer);
+  if (cell_is_active_black_holes(ci, e)) runner_do_gas_swallow(r, cj, timer);
+}
+
+/**
+ * @brief Process all the BH particles in a cell that have been flagged for
+ * swallowing by a black hole.
+ *
+ * This is done by recursing down to the leaf-level and skipping the sub-cells
+ * that have not been drifted as they would not have any particles with
+ * swallowing flag. We then loop over the particles with a flag and look into
+ * the space-wide list of black holes for the particle with the corresponding
+ * ID. If found, the BH swallows the BH particle and the BH particle is
+ * removed. If the cell is local, we may be looking for a foreign BH, in which
+ * case, we do not update the BH (that will be done on its node) but just remove
+ * the BH particle.
+ *
+ * @param r The thread #runner.
+ * @param c The #cell.
+ * @param timer Are we timing this?
+ */
+void runner_do_bh_swallow(struct runner *r, struct cell *c, int timer) {
+
+  struct engine *e = r->e;
+  struct space *s = e->s;
+  struct bpart *bparts = s->bparts;
+  const size_t nr_bpart = s->nr_bparts;
+#ifdef WITH_MPI
+  struct bpart *bparts_foreign = s->bparts_foreign;
+  const size_t nr_bparts_foreign = s->nr_bparts_foreign;
+#endif
+
+  struct bpart *cell_bparts = c->black_holes.parts;
+
+  /* Early abort?
+   * (We only want cells for which we drifted the gas as these are
+   * the only ones that could have gas particles that have been flagged
+   * for swallowing) */
+  if (c->black_holes.count == 0 ||
+      c->black_holes.ti_old_part != e->ti_current) {
+    return;
+  }
+
+  /* Loop over the progeny ? */
+  if (c->split) {
+    for (int k = 0; k < 8; k++) {
+      if (c->progeny[k] != NULL) {
+        struct cell *restrict cp = c->progeny[k];
+
+        runner_do_bh_swallow(r, cp, 0);
+      }
+    }
+  } else {
+
+    /* Loop over all the gas particles in the cell
+     * Note that the cell (and hence the parts) may be local or foreign. */
+    const size_t nr_cell_bparts = c->black_holes.count;
+    for (size_t k = 0; k < nr_cell_bparts; k++) {
+
+      /* Get a handle on the part. */
+      struct bpart *const cell_bp = &cell_bparts[k];
+
+      /* Ignore inhibited particles (they have already been removed!) */
+      if (bpart_is_inhibited(cell_bp, e)) continue;
+
+      /* Get the ID of the black holes that will swallow this part */
+      const long long swallow_id =
+          black_holes_get_bpart_swallow_id(&cell_bp->merger_data);
+
+      /* Has this particle been flagged for swallowing? */
+      if (swallow_id >= 0) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+        if (cell_bp->