src/cell_black_holes.h

@@ -55,9 +55,9 @@ struct cell_black_holes {
     struct task *density_ghost;
 
     /*! The ghost tasks related to BH swallowing */
-    struct task *swallow_ghost_0;
     struct task *swallow_ghost_1;
     struct task *swallow_ghost_2;
+    struct task *swallow_ghost_3;
 
     /*! Linked list of the tasks computing this cell's BH density. */
     struct link *density;

src/cell_convert_part.c

@@ -26,6 +26,7 @@
 #include "cell.h"
 
 /* Local headers. */
+#include "active.h"
 #include "engine.h"
 #include "hydro.h"
 #include "sink_properties.h"
@@ -220,6 +221,7 @@ struct spart *cell_add_spart(struct engine *e, struct cell *const c) {
 
     /* Update the spart->gpart links (shift by 1) */
     for (size_t i = 0; i < n_copy; ++i) {
+
 #ifdef SWIFT_DEBUG_CHECKS
       if (c->stars.parts[i + 1].gpart == NULL) {
         error("Incorrectly linked spart!");
@@ -264,6 +266,9 @@ struct spart *cell_add_spart(struct engine *e, struct cell *const c) {
   sp->ti_drift = e->ti_current;
 #endif
 
+  /* Give the new particle the correct depth */
+  cell_set_spart_h_depth(sp, c);
+
   /* Register that we used one of the free slots. */
   const size_t one = 1;
   atomic_sub(&e->s->nr_extra_sparts, one);
@@ -288,7 +293,7 @@ struct spart *cell_add_spart(struct engine *e, struct cell *const c) {
 struct sink *cell_add_sink(struct engine *e, struct cell *const c) {
   /* Perform some basic consitency checks */
   if (c->nodeID != engine_rank) error("Adding sink on a foreign node");
-  if (c->grav.ti_old_part != e->ti_current) error("Undrifted cell!");
+  if (c->sinks.ti_old_part != e->ti_current) error("Undrifted cell!");
   if (c->split) error("Addition of sink performed above the leaf level");
 
   /* Progeny number at each level */
@@ -353,6 +358,11 @@ struct sink *cell_add_sink(struct engine *e, struct cell *const c) {
 
     /* Update the sink->gpart links (shift by 1) */
     for (size_t i = 0; i < n_copy; ++i) {
+
+      // TODO: Matthieu figure out whether this is strictly needed
+      /* Skip inhibited (swallowed) sink particles */
+      if (sink_is_inhibited(&c->sinks.parts[i + 1], e)) continue;
+
 #ifdef SWIFT_DEBUG_CHECKS
       if (c->sinks.parts[i + 1].gpart == NULL) {
         error("Incorrectly linked sink!");
@@ -397,6 +407,9 @@ struct sink *cell_add_sink(struct engine *e, struct cell *const c) {
   sp->ti_drift = e->ti_current;
 #endif
 
+  /* Give the new particle the correct depth */
+  cell_set_sink_h_depth(sp, c);
+
   /* Register that we used one of the free slots. */
   const size_t one = 1;
   atomic_sub(&e->s->nr_extra_sinks, one);
@@ -489,6 +502,10 @@ struct gpart *cell_add_gpart(struct engine *e, struct cell *c) {
     /* Update the gpart->spart links (shift by 1) */
     struct gpart *gparts = c->grav.parts;
     for (size_t i = 0; i < n_copy; ++i) {
+
+      /* Skip inhibited particles */
+      if (gpart_is_inhibited(&c->grav.parts[i + 1], e)) continue;
+
       if (gparts[i + 1].type == swift_type_gas) {
         s->parts[-gparts[i + 1].id_or_neg_offset].gpart++;
       } else if (gparts[i + 1].type == swift_type_stars) {
@@ -572,7 +589,7 @@ void cell_remove_part(const struct engine *e, struct cell *c, struct part *p,
   /* Mark the gpart as inhibited and stand-alone */
   if (p->gpart) {
     p->gpart->time_bin = time_bin_inhibited;
-    p->gpart->id_or_neg_offset = p->id;
+    p->gpart->id_or_neg_offset = 1;
     p->gpart->type = swift_type_dark_matter;
   }
 
@@ -645,7 +662,7 @@ void cell_remove_spart(const struct engine *e, struct cell *c,
   sp->time_bin = time_bin_inhibited;
   if (sp->gpart) {
     sp->gpart->time_bin = time_bin_inhibited;
-    sp->gpart->id_or_neg_offset = sp->id;
+    sp->gpart->id_or_neg_offset = 1;
     sp->gpart->type = swift_type_dark_matter;
   }
 
@@ -684,7 +701,7 @@ void cell_remove_bpart(const struct engine *e, struct cell *c,
   bp->time_bin = time_bin_inhibited;
   if (bp->gpart) {
     bp->gpart->time_bin = time_bin_inhibited;
-    bp->gpart->id_or_neg_offset = bp->id;
+    bp->gpart->id_or_neg_offset = 1;
     bp->gpart->type = swift_type_dark_matter;
   }
 
@@ -722,7 +739,7 @@ void cell_remove_sink(const struct engine *e, struct cell *c,
   sink->time_bin = time_bin_inhibited;
   if (sink->gpart) {
     sink->gpart->time_bin = time_bin_inhibited;
-    sink->gpart->id_or_neg_offset = sink->id;
+    sink->gpart->id_or_neg_offset = 1;
     sink->gpart->type = swift_type_dark_matter;
   }
 
@@ -905,6 +922,9 @@ struct spart *cell_convert_part_to_spart(struct engine *e, struct cell *c,
   /* Set a smoothing length */
   sp->h = p->h;
 
+  /* Give the new particle the correct depth */
+  cell_set_spart_h_depth(sp, c);
+
   /* Here comes the Sun! */
   return sp;
 }
@@ -983,6 +1003,9 @@ struct spart *cell_spawn_new_spart_from_part(struct engine *e, struct cell *c,
   /* Set a smoothing length */
   sp->h = p->h;
 
+  /* Give the new particle the correct depth */
+  cell_set_spart_h_depth(sp, c);
+
   /* Here comes the Sun! */
   return sp;
 }
@@ -1051,10 +1074,9 @@ struct sink *cell_convert_part_to_sink(struct engine *e, struct cell *c,
 #ifdef SWIFT_DEBUG_CHECKS
   sp->ti_kick = gp->ti_kick;
   gp->ti_drift = sp->ti_drift;
-#endif
 
-  /* Set a smoothing length */
-  sp->r_cut = e->sink_properties->cut_off_radius;
+  message("A new sink (%lld) is born !", sp->id);
+#endif
 
   /* Here comes the Sink! */
   return sp;
@@ -1130,7 +1152,10 @@ struct spart *cell_spawn_new_spart_from_sink(struct engine *e, struct cell *c,
 #endif
 
   /* Set a smoothing length */
-  sp->h = s->r_cut;
+  sp->h = s->h;
+
+  /* Give the new particle the correct depth */
+  cell_set_spart_h_depth(sp, c);
 
   /* Here comes the Sun! */
   return sp;

src/cell_drift.c

@@ -27,6 +27,7 @@
 
 /* Local headers. */
 #include "active.h"
+#include "adaptive_softening.h"
 #include "drift.h"
 #include "feedback.h"
 #include "gravity.h"
@@ -135,6 +136,22 @@ void cell_set_ti_old_bpart(struct cell *c, const integertime_t ti) {
   }
 }
 
+/**
+ * @brief Recursively set the sinks' ti_old_part to the current time.
+ *
+ * @param c The cell to update.
+ * @param ti The current integer time.
+ */
+void cell_set_ti_old_sink(struct cell *c, const integertime_t ti) {
+
+  c->sinks.ti_old_part = ti;
+  if (c->split) {
+    for (int k = 0; k < 8; k++) {
+      if (c->progeny[k] != NULL) cell_set_ti_old_sink(c->progeny[k], ti);
+    }
+  }
+}
+
 /**
  * @brief Recursively drifts the #part in a cell hierarchy.
  *
@@ -315,6 +332,9 @@ void cell_drift_part(struct cell *c, const struct engine *e, int force,
       p->h = min(p->h, hydro_h_max);
       p->h = max(p->h, hydro_h_min);
 
+      /* Set the appropriate depth level for this particle */
+      cell_set_part_h_depth(p, c);
+
       /* Compute (square of) motion since last cell construction */
       const float dx2 = xp->x_diff[0] * xp->x_diff[0] +
                         xp->x_diff[1] * xp->x_diff[1] +
@@ -340,6 +360,7 @@ void cell_drift_part(struct cell *c, const struct engine *e, int force,
       /* Get ready for a density calculation */
       if (part_is_active(p, e)) {
         hydro_init_part(p, &e->s->hs);
+        adaptive_softening_init_part(p);
         mhd_init_part(p);
         black_holes_init_potential(&p->black_holes_data);
         chemistry_init_part(p, e->chemistry);
@@ -347,7 +368,7 @@ void cell_drift_part(struct cell *c, const struct engine *e, int force,
         tracers_after_init(p, xp, e->internal_units, e->physical_constants,
                            with_cosmology, e->cosmology, e->hydro_properties,
                            e->cooling_func, e->time);
-        sink_init_part(p);
+        sink_init_part(p, e->sink_properties);
         rt_init_part(p);
 
         /* Update the maximal active smoothing length in the cell */
@@ -567,6 +588,7 @@ void cell_drift_spart(struct cell *c, const struct engine *e, int force,
   const int periodic = e->s->periodic;
   const double dim[3] = {e->s->dim[0], e->s->dim[1], e->s->dim[2]};
   const int with_cosmology = (e->policy & engine_policy_cosmology);
+  const int with_rt = (e->policy & engine_policy_rt);
   const float stars_h_max = e->hydro_properties->h_max;
   const float stars_h_min = e->hydro_properties->h_min;
   const integertime_t ti_old_spart = c->stars.ti_old_part;
@@ -707,6 +729,9 @@ void cell_drift_spart(struct cell *c, const struct engine *e, int force,
       sp->h = min(sp->h, stars_h_max);
       sp->h = max(sp->h, stars_h_min);
 
+      /* Set the appropriate depth level for this particle */
+      cell_set_spart_h_depth(sp, c);
+
       /* Compute (square of) motion since last cell construction */
       const float dx2 = sp->x_diff[0] * sp->x_diff[0] +
                         sp->x_diff[1] * sp->x_diff[1] +
@@ -729,7 +754,8 @@ void cell_drift_spart(struct cell *c, const struct engine *e, int force,
         rt_init_spart(sp);
 
         /* Update the maximal active smoothing length in the cell */
-        cell_h_max_active = max(cell_h_max_active, sp->h);
+        if (feedback_is_active(sp, e) || with_rt)
+          cell_h_max_active = max(cell_h_max_active, sp->h);
       }
     }
 
@@ -911,6 +937,9 @@ void cell_drift_bpart(struct cell *c, const struct engine *e, int force,
       bp->h = min(bp->h, black_holes_h_max);
       bp->h = max(bp->h, black_holes_h_min);
 
+      /* Set the appropriate depth level for this particle */
+      cell_set_bpart_h_depth(bp, c);
+
       /* Compute (square of) motion since last cell construction */
       const float dx2 = bp->x_diff[0] * bp->x_diff[0] +
                         bp->x_diff[1] * bp->x_diff[1] +
@@ -968,13 +997,15 @@ void cell_drift_sink(struct cell *c, const struct engine *e, int force) {
   const int periodic = e->s->periodic;
   const double dim[3] = {e->s->dim[0], e->s->dim[1], e->s->dim[2]};
   const int with_cosmology = (e->policy & engine_policy_cosmology);
+  const float sinks_h_max = e->hydro_properties->h_max;
+  const float sinks_h_min = e->hydro_properties->h_min;
   const integertime_t ti_old_sink = c->sinks.ti_old_part;
   const integertime_t ti_current = e->ti_current;
   struct sink *const sinks = c->sinks.parts;
 
   float dx_max = 0.f, dx2_max = 0.f;
-  float cell_r_max = 0.f;
-  float cell_r_max_active = 0.f;
+  float cell_h_max = 0.f;
+  float cell_h_max_active = 0.f;
 
   /* Drift irrespective of cell flags? */
   force = (force || cell_get_flag(c, cell_flag_do_sink_drift));
@@ -994,7 +1025,7 @@ void cell_drift_sink(struct cell *c, const struct engine *e, int force) {
     cell_clear_flag(c, cell_flag_do_sink_drift | cell_flag_do_sink_sub_drift);
 
     /* Update the time of the last drift */
-    c->sinks.ti_old_part = ti_current;
+    cell_set_ti_old_sink(c, ti_current);
 
     return;
   }
@@ -1014,14 +1045,14 @@ void cell_drift_sink(struct cell *c, const struct engine *e, int force) {
 
         /* Update */
         dx_max = max(dx_max, cp->sinks.dx_max_part);
-        cell_r_max = max(cell_r_max, cp->sinks.r_cut_max);
-        cell_r_max_active = max(cell_r_max_active, cp->sinks.r_cut_max_active);
+        cell_h_max = max(cell_h_max, cp->sinks.h_max);
+        cell_h_max_active = max(cell_h_max_active, cp->sinks.h_max_active);
       }
     }
 
     /* Store the values */
-    c->sinks.r_cut_max = cell_r_max;
-    c->sinks.r_cut_max_active = cell_r_max_active;
+    c->sinks.h_max = cell_h_max;
+    c->sinks.h_max_active = cell_h_max_active;
     c->sinks.dx_max_part = dx_max;
 
     /* Update the time of the last drift */
@@ -1091,7 +1122,12 @@ void cell_drift_sink(struct cell *c, const struct engine *e, int force) {
         }
       }
 
-      /* sp->h does not need to be limited. */
+      /* Limit h to within the allowed range */
+      sink->h = min(sink->h, sinks_h_max);
+      sink->h = max(sink->h, sinks_h_min);
+
+      /* Set the appropriate depth level for this particle */
+      cell_set_sink_h_depth(sink, c);
 
       /* Compute (square of) motion since last cell construction */
       const float dx2 = sink->x_diff[0] * sink->x_diff[0] +
@@ -1100,7 +1136,7 @@ void cell_drift_sink(struct cell *c, const struct engine *e, int force) {
       dx2_max = max(dx2_max, dx2);
 
       /* Maximal smoothing length */
-      cell_r_max = max(cell_r_max, sink->r_cut);
+      cell_h_max = max(cell_h_max, sink->h);
 
       /* Mark the particle has not being swallowed */
       sink_mark_sink_as_not_swallowed(&sink->merger_data);
@@ -1109,7 +1145,7 @@ void cell_drift_sink(struct cell *c, const struct engine *e, int force) {
       if (sink_is_active(sink, e)) {
         sink_init_sink(sink);
 
-        cell_r_max_active = max(cell_r_max_active, sink->r_cut);
+        cell_h_max_active = max(cell_h_max_active, sink->h);
       }
     }
 
@@ -1117,8 +1153,8 @@ void cell_drift_sink(struct cell *c, const struct engine *e, int force) {
     dx_max = sqrtf(dx2_max);
 
     /* Store the values */
-    c->sinks.r_cut_max = cell_r_max;
-    c->sinks.r_cut_max_active = cell_r_max_active;
+    c->sinks.h_max = cell_h_max;
+    c->sinks.h_max_active = cell_h_max_active;
     c->sinks.dx_max_part = dx_max;
 
     /* Update the time of the last drift */

src/cell_grav.h

@@ -25,6 +25,7 @@
 #include <config.h>
 
 /* Local includes. */
+#include "gravity.h"
 #include "lock.h"
 #include "timeline.h"
 
@@ -34,11 +35,29 @@
 
 struct cell_grav {
 
-  /*! Pointer to the #gpart data. */
-  struct gpart *parts;
+  union {
 
-  /*! Pointer to the #spart data at rebuild time. */
-  struct gpart *parts_rebuild;
+    /*! Pointer to the #gpart data. */
+    struct gpart *parts;
+
+    /*! or #gpart_foreign data. */
+    struct gpart_foreign *parts_foreign;
+
+    /*! or #gpart_fof_foreign data. */
+    struct gpart_fof_foreign *parts_fof_foreign;
+  };
+
+  union {
+
+    /*! Pointer to the #gpart data at rebuild time. */
+    struct gpart *parts_rebuild;
+
+    /*! Pointer to the #gpart_foreign data at rebuild time. */
+    struct gpart_foreign *parts_foreign_rebuild;
+
+    /*! Pointer to the #gpart_fof_foreign data at rebuild time. */
+    struct gpart_fof_foreign *parts_fof_foreign_rebuild;
+  };
 
   /*! This cell's multipole. */
   struct gravity_tensors *multipole;

src/cell_grid.c

+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2024 Matthieu Schaller (schaller@strw.leidenuniv.nl)
+ *                             Yolan Uyttenhove (Yolan.Uyttenhove@UGent.be)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+
+/* Config parameters. */
+#include <config.h>
+
+/* Corresponding header */
+#include "cell_grid.h"
+
+/* Local headers */
+#include "engine.h"
+#include "error.h"
+#include "space_getsid.h"
+
+/**
+ * @brief Recursively free grid memory for cell.
+ *
+ * @param c The #cell.
+ */
+void cell_free_grid_rec(struct cell *c) {
+
+#ifndef MOVING_MESH
+  /* Nothing to do as we have no tessellations */
+#else
+#ifdef SWIFT_DEBUG_CHECKS
+  if (c->grid.construction_level != c && c->grid.voronoi != NULL)
+    error("Grid allocated, but not on grid construction level!");
+#endif
+  if (c->grid.construction_level == NULL) {
+    for (int k = 0; k < 8; k++)
+      if (c->progeny[k] != NULL) cell_free_grid_rec(c->progeny[k]);
+
+  } else if (c->grid.construction_level == c) {
+    cell_free_grid(c);
+  } else if (c->grid.construction_level != c) {
+    error("Somehow ended up below grid construction level!");
+  }
+#endif
+}
+
+/**
+ * @brief Updates the grid.self_completeness flag on this cell and its
+ * sub-cells.
+ *
+ * This cell satisfies the local completeness criterion for the Voronoi grid.
+ *
+ * A cell is defined as locally complete if, when we would split that cell in
+ * thirds along each dimension (i.e. in 27 smaller cells), every small cube
+ * would contain at least one particle.
+ *
+ * If a given cell and its direct neighbours on the same level in the AMR tree
+ * are self-complete, the Voronoi grid of that cell can completely be
+ * constructed using only particles of this cell and its direct neighbours,
+ * i.e. by doing a normal SWIFT neighbour loop.
+ *
+ * @param c The #cell to be checked
+ * @param force Whether to forcefully recompute the completeness if it was not
+ * invalidated.
+ * */
+void cell_grid_update_self_completeness(struct cell *c, int force) {
+  if (c == NULL) return;
+  if (!force && (c->grid.self_completeness != grid_invalidated_completeness))
+    return;
+
+  if (c->split) {
+    int all_complete = 1;
+
+    /* recurse */
+    for (int i = 0; all_complete && i < 8; i++) {
+      if (c->progeny[i] != NULL) {
+        cell_grid_update_self_completeness(c->progeny[i], force);
+        /* As long as all progeny is complete, this cell can safely be split for
+         * the grid construction (when not considering neighbouring cells) */
+        all_complete &=
+            (c->progeny[i]->grid.self_completeness == grid_complete);
+      }
+    }
+
+    /* If all sub-cells are complete, this cell is also complete. */
+    if (all_complete) {
+      c->grid.self_completeness = grid_complete;
+      /* We set complete to true for now */
+      c->grid.complete = 1;
+      /* We are done here */
+      return;
+    }
+  }
+
+  /* If this cell is not split, or not all subcells are complete, we need to
+   * check if this cell is complete by looping over all the particles. */
+
+  /* criterion = 0b111_111_111_111_111_111_111_111_111*/
+#ifdef HYDRO_DIMENSION_1D
+  const int criterion = (1 << 3) - 1;
+#elif defined(HYDRO_DIMENSION_2D)
+  const int criterion = (1 << 9) - 1;
+#elif defined(HYDRO_DIMENSION_3D)
+  const int criterion = (1 << 27) - 1;
+#else
+#error "Unknown hydro dimension"
+#endif
+  int flags = 0;
+  for (int i = 0; flags != criterion && i < c->hydro.count; i++) {
+    struct part *p = &c->hydro.parts[i];
+    int x_bin = (int)(3. * (p->x[0] - c->loc[0]) / c->width[0]);
+    int y_bin = (int)(3. * (p->x[1] - c->loc[1]) / c->width[1]);
+    int z_bin = (int)(3. * (p->x[2] - c->loc[2]) / c->width[2]);
+    if (x_bin >= 0 && x_bin < 3 && y_bin >= 0 && y_bin < 3 && z_bin >= 0 &&
+        z_bin < 3) {
+      flags |= 1 << (x_bin + 3 * y_bin + 9 * z_bin);
+    }
+  }
+
+  /* Set completeness flags accordingly */
+  if (flags == criterion) {
+    c->grid.self_completeness = grid_complete;
+    c->grid.complete = 1;
+  } else {
+    c->grid.self_completeness = grid_incomplete;
+    c->grid.complete = 0;
+  }
+}
+
+/**
+ * @brief (mapper) Updates the grid.self_completeness flag on this cell and its
+ * sub-cells.
+ *
+ * This function recomputes the self_completeness flag for all cells containing
+ * particles.
+ */
+void cell_grid_set_self_completeness_mapper(void *map_data, int num_elements,
+                                            void *extra_data) {
+  /* Extract the engine pointer. */
+  struct engine *e = (struct engine *)extra_data;
+
+  struct space *s = e->s;
+  const int nodeID = e->nodeID;
+  struct cell *cells = s->cells_top;
+
+  /* Loop through the elements, which are just byte offsets from NULL. */
+  for (int ind = 0; ind < num_elements; ind++) {
+    /* Get the cell index. */
+    const int cid = (size_t)(map_data) + ind;
+
+    /* Get the cell */
+    struct cell *c = &cells[cid];
+
+    /* A top level cell can be empty in 1D and 2D simulations, just skip it */
+    if (c->hydro.count == 0) {
+      continue;
+#ifdef SWIFT_DEBUG_CHECKS
+      if (hydro_dimension == 3)
+        error("Found empty top-level cell while running in 3D!");
+#endif
+    }
+    if (c->nodeID != nodeID) continue;
+
+    /* Set the splittable attribute for the moving mesh */
+    cell_grid_update_self_completeness(c, /*force*/ 1);
+  }
+}
+
+/**
+ * @brief Updates the grid.completeness flag for the given pair of cells and all
+ * recursive pairs of sub-cells.
+ *
+ * If one of the cells of the pair is not self-complete, or requires a rebuild,
+ * we mark the other cell in the pair as incomplete (it cannot construct its
+ * Voronoi grid on that level).
+ *
+ *
+ *
+ * @param ci The first #cell of the pair to be checked.
+ * @param cj The second #cell of the pair to be checked. If NULL, only check
+ * pairs of subcells from ci.
+ * @param sid The sort_id (direction) of the pair.
+ * @param e The #engine.
+ * */
+void cell_grid_set_pair_completeness(struct cell *restrict ci,
+                                     struct cell *restrict cj, int sid,
+                                     const struct engine *e) {
+
+  int ci_local = ci->nodeID == e->nodeID;
+  int cj_local = cj != NULL ? cj->nodeID == e->nodeID : 0;
+  /* Anything to do here? */
+  if (!ci_local && !cj_local) return;
+
+  /* Self or pair? */
+  if (cj == NULL) {
+    /* Self: Here we just need to recurse to hit all the pairs of sub-cells */
+    if (ci->split) {
+      /* recurse */
+      for (int k = 0; k < 7; k++) {
+        if (ci->progeny[k] != NULL) {
+          /* Self: Recurse for pairs of sub-cells of this sub-cell */
+          cell_grid_set_pair_completeness(ci->progeny[k], NULL, 0, e);
+
+          /* Recurse for pairs of sub-cells */
+          for (int l = k + 1; l < 8; l++) {
+            if (ci->progeny[l] != NULL) {
+              /* Get sid for pair */
+              int sid_sub = sub_sid_flag[k][l];
+              /* Pair: Recurse for pairs of sub-cells of this pair of sub-cells
+               */
+              cell_grid_set_pair_completeness(ci->progeny[k], ci->progeny[l],
+                                              sid_sub, e);
+            }
+          }
+        }
+      }
+    }
+  } else {
+    /* pair: Here we need to recurse further AND check whether one of the
+     * neighbouring cells invalidates the completeness of the other. */
+    struct cell_split_pair pairs = cell_split_pairs[sid];
+    if (ci->split && cj->split) {
+      /* recurse */
+      for (int i = 0; i < pairs.count; i++) {
+        struct cell *ci_sub = ci->progeny[pairs.pairs[i].pid];
+        struct cell *cj_sub = cj->progeny[pairs.pairs[i].pjd];
+        if (ci_sub == NULL || cj_sub == NULL) continue;
+        double shift[3];
+        int sid_sub =
+            space_getsid_and_swap_cells(e->s, &ci_sub, &cj_sub, shift);
+#ifdef SWIFT_DEBUG_CHECKS
+        assert(sid_sub == pairs.pairs[i].sid);
+#endif
+        cell_grid_set_pair_completeness(ci_sub, cj_sub, sid_sub, e);
+      }
+    } else if (!ci->split && cj->split) {
+      /* Set the completeness for the sub-cells of cj for this sid to 0 (they
+       * have no neighbouring cell on the same level for this SID) */
+      for (int i = 0; i < pairs.count; i++) {
+        int l = pairs.pairs[i].pjd;
+        if (cj->progeny[l] != NULL) cj->progeny[l]->grid.complete = 0;
+      }
+    } else if (!cj->split && ci->split) {
+      /* Set the completeness for the sub-cells of ci for this sid to 0 (they
+       * have no neighbouring cell on the same level for this SID) */
+      for (int i = 0; i < pairs.count; i++) {
+        int k = pairs.pairs[i].pid;
+        if (ci->progeny[k] != NULL) ci->progeny[k]->grid.complete = 0;
+      }
+    }
+
+    /* Update these cells' completeness flags (i.e. check whether the
+     * neighbouring cell invalidates completeness)
+     * We need to use atomics here, since multiple threads may change this at
+     * the same time. */
+    if (ci_local) {
+      atomic_and(&ci->grid.complete,
+                 !cell_grid_pair_invalidates_completeness(ci, cj));
+    }
+    if (cj_local) {
+      atomic_and(&cj->grid.complete,
+                 !cell_grid_pair_invalidates_completeness(cj, ci));
+    }
+  }
+}
+
+/**
+ * @brief (mapper) Sets the grid.completeness flag for all cells, by looping
+ * aggregating the self_completeness flags of all neighbours of each cell.
+ *
+ * A cell is considered complete if it and all its neighbours are
+ * self_complete. The Voronoi grid may be constructed at any level in the AMR
+ * tree as long as the cell at that level is complete.
+ * */
+void cell_set_grid_completeness_mapper(void *map_data, int num_elements,
+                                       void *extra_data) {
+  /* Extract the engine pointer. */
+  struct engine *e = (struct engine *)extra_data;
+  const int periodic = e->s->periodic;
+
+  struct space *s = e->s;
+  const int nodeID = e->nodeID;
+  const int *cdim = s->cdim;
+  struct cell *cells = s->cells_top;
+
+  /* Loop through the elements, which are just byte offsets from NULL, to set
+   * the neighbour flags. */
+  for (int ind = 0; ind < num_elements; ind++) {
+    /* Get the cell index. */
+    const int cid = (size_t)(map_data) + ind;
+
+    /* Integer indices of the cell in the top-level grid */
+    const int i = cid / (cdim[1] * cdim[2]);
+    const int j = (cid / cdim[2]) % cdim[1];
+    const int k = cid % cdim[2];
+
+    /* Get the cell */
+    struct cell *ci = &cells[cid];
+
+    /* Anything to do here? */
+    if (ci->hydro.count == 0) continue;
+
+    const int ci_local = ci->nodeID == nodeID;
+
+    /* Update completeness for all the pairs of sub cells of this cell */
+    if (ci_local) cell_grid_set_pair_completeness(ci, NULL, 0, e);
+
+    /* Now loop over all the neighbours of this cell to also update the
+     * completeness for pairs with this cell and all pairs of sub-cells */
+    for (int ii = -1; ii < 2; ii++) {
+      int iii = i + ii;
+      if (!periodic && (iii < 0 || iii >= cdim[0])) continue;
+      iii = (iii + cdim[0]) % cdim[0];
+      for (int jj = -1; jj < 2; jj++) {
+        int jjj = j + jj;
+        if (!periodic && (jjj < 0 || jjj >= cdim[1])) continue;
+        jjj = (jjj + cdim[1]) % cdim[1];
+        for (int kk = -1; kk < 2; kk++) {
+          int kkk = k + kk;
+          if (!periodic && (kkk < 0 || kkk >= cdim[2])) continue;
+          kkk = (kkk + cdim[2]) % cdim[2];
+
+          /* Get the neighbouring cell */
+          const int cjd = cell_getid(cdim, iii, jjj, kkk);
+          struct cell *cj = &cells[cjd];
+
+          /* Treat pairs only once. */
+          const int cj_local = cj->nodeID == nodeID;
+          if (cid >= cjd || cj->hydro.count == 0 || (!ci_local && !cj_local))
+            continue;
+
+          /* Update the completeness flag for this pair of cells and all pair of
+           * sub-cells */
+          int sid = (kk + 1) + 3 * ((jj + 1) + 3 * (ii + 1));
+          const int flip = runner_flip[sid];
+          sid = sortlistID[sid];
+          if (flip) {
+            cell_grid_set_pair_completeness(cj, ci, sid, e);
+          } else {
+            cell_grid_set_pair_completeness(ci, cj, sid, e);
+          }
+        }
+      }
+    } /* Now loop over all the neighbours of this cell */
+  } /* Loop through the elements, which are just byte offsets from NULL. */
+}
+
+/**
+ * @brief Select a suitable construction level for the Voronoi grid.
+ *
+ * The Voronoi grid is constructed at the lowest level for which the cell
+ * has more than #space_grid_split_threshold hydro particles *and* is marked
+ * as complete for the Voronoi construction.
+ *
+ * Cells below the construction level store a pointer to its higher level cell
+ * at the construction level.
+ *
+ * @param c The #cell
+ * @param construction_level NULL, if we are yet to encounter the suitable
+ * construction level, or a pointer to the parent-cell of c at the construction
+ * level.
+ */
+void cell_set_grid_construction_level(struct cell *c,
+                                      struct cell *construction_level) {
+
+  /* Above construction level? */
+  if (construction_level == NULL) {
+    /* Check if we can split this cell (i.e. all sub-cells are complete) */
+    int splittable = c->split && c->hydro.count > space_grid_split_threshold;
+    if (!c->grid.complete) {
+      /* Are we on the top level? */
+      if (c->top == c) {
+        warning("Found incomplete top level cell!");
+        splittable = 0;
+      } else {
+        error("Found incomplete cell above construction level!");
+      }
+    }
+    for (int k = 0; splittable && k < 8; k++) {
+      if (c->progeny[k] != NULL) splittable &= c->progeny[k]->grid.complete;
+    }
+
+    if (!splittable) {
+      /* This is the first time we encounter an unsplittable cell, meaning that
+       * it has too few particles to be split further or one of its progenitors
+       * is not complete. I.e. we have arrived at the construction level! */
+      construction_level = c;
+    }
+  }
+
+  /* Set the construction level of this cell */
+  c->grid.construction_level = construction_level;
+
+  /* Recurse. */
+  if (c->split)
+    for (int k = 0; k < 8; k++) {
+      if (c->progeny[k] != NULL)
+        cell_set_grid_construction_level(c->progeny[k], construction_level);
+    }
+}
+
+void cell_set_grid_construction_level_mapper(void *map_data, int num_elements,
+                                             void *extra_data) {
+  /* Extract the engine pointer. */
+  struct engine *e = (struct engine *)extra_data;
+
+  struct space *s = e->s;
+  const int nodeID = e->nodeID;
+  struct cell *cells = s->cells_top;
+
+  /* Loop through the elements, which are just byte offsets from NULL, to set
+   * the neighbour flags. */
+  for (int ind = 0; ind < num_elements; ind++) {
+
+    /* Get the cell index. */
+    const int cid = (size_t)(map_data) + ind;
+
+    /* Get the cell */
+    struct cell *ci = &cells[cid];
+
+    /* Anything to do here? */
+    if (ci->hydro.count == 0) continue;
+    const int ci_local = ci->nodeID == nodeID;
+
+    /* This cell's completeness flags are now set all the way down the cell
+     * hierarchy. We can now set the construction level. */
+    if (ci_local) {
+      cell_set_grid_construction_level(ci, NULL);
+    }
+  }
+}

src/cell_grid.h

+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2024 Matthieu Schaller (schaller@strw.leidenuniv.nl)
+ *                             Yolan Uyttenhove (Yolan.Uyttenhove@UGent.be)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+
+#ifndef SWIFTSIM_CELL_GRID_H
+#define SWIFTSIM_CELL_GRID_H
+
+/* Config parameters. */
+#include <config.h>
+
+/* Includes. */
+#include <stddef.h>
+
+/* Local includes */
+#include "const.h"
+#include "shadowswift/voronoi.h"
+#include "timeline.h"
+
+/*! @brief Enum indicating the completeness for the Voronoi mesh of this cell.
+ *
+ * A cell is considered complete when it and its neighbours on the same level in
+ * the AMR have at least one particle in every 1/27th cube of the cell (obtained
+ * by dividing cells in three along all axes).
+ *
+ * The Voronoi grid can safely be constructed on any level where the cell is
+ * complete. */
+enum grid_completeness {
+  grid_invalidated_completeness = 0,
+  grid_complete,
+  grid_incomplete,
+};
+
+struct cell_grid {
+  /*! Pointer to parent where the grid is constructed. */
+  struct cell *construction_level;
+
+  /*! Pointer to shallowest parent of this cell used in any pair construction
+   * task. Can be above the construction level of this cell. We need to drift at
+   * this level. */
+  struct cell *super;
+
+  /*! Whether this cell is complete (at least one particle in all 27 sub-cells
+   * if this cell is divided in thirds along each axis). */
+  enum grid_completeness self_completeness;
+
+  /*! Whether this cell is itself complete and has directly neighbouring cell
+   * on the same level in all directions which are also complete. */
+  int complete;
+
+#ifdef WITH_MPI
+  /*! Flags indicating whether we should send the faces for the corresponding
+   * SIDs over MPI */
+  int send_flags;
+#endif
+
+  /*! Pointer to the voronoi struct of this cell (if any) */
+  struct voronoi *voronoi;
+
+  /*! Pointer to this cells construction task. */
+  struct task *construction;
+
+  /*! Linked list of this cells outgoing construction synchronization tasks
+   * (i.e. for cells that need this cell for their construction task) */
+  struct link *sync_out;
+
+  /*! Linked list of this cells incoming construction synchronization tasks
+   * (i.e. cells needed for this cell's construction task) */
+  struct link *sync_in;
+
+  /*! Time of last construction */
+  integertime_t ti_old;
+};
+
+struct pcell_faces {
+  size_t counts[27];
+
+  struct voronoi_pair faces[];
+};
+
+/*! @brief Enum used to indicate whether a cell is above, below or on the
+ * construction level. Only used in the packed cell representation */
+enum grid_construction_level {
+  grid_above_construction_level,
+  grid_on_construction_level,
+  grid_below_construction_level
+};
+
+#endif  // SWIFTSIM_CELL_GRID_H

src/cell_hydro.h

@@ -106,6 +106,9 @@ struct cell_hydro {
     /*! Last (integer) time the cell's part were drifted forward in time. */
     integertime_t ti_old_part;
 
+    /*! Spin-lock for the case where we do an extra sort of the cell. */
+    swift_lock_type extra_sort_lock;
+
     /*! Max smoothing length of active particles in this cell. */
     float h_max_active;
 

src/cell_pack.c

@@ -25,6 +25,9 @@
 /* This object's header. */
 #include "cell.h"
 
+/* Local headers */
+#include "gravity.h"
+
 /**
  * @brief Pack the data of the given cell and all it's sub-cells.
  *
@@ -44,7 +47,7 @@ int cell_pack(struct cell *restrict c, struct pcell *restrict pc,
   pc->hydro.h_max = c->hydro.h_max;
   pc->stars.h_max = c->stars.h_max;
   pc->black_holes.h_max = c->black_holes.h_max;
-  pc->sinks.r_cut_max = c->sinks.r_cut_max;
+  pc->sinks.h_max = c->sinks.h_max;
 
   pc->hydro.ti_end_min = c->hydro.ti_end_min;
   pc->grav.ti_end_min = c->grav.ti_end_min;
@@ -68,6 +71,8 @@ int cell_pack(struct cell *restrict c, struct pcell *restrict pc,
   pc->black_holes.count = c->black_holes.count;
   pc->maxdepth = c->maxdepth;
 
+  pc->grid.self_completeness = c->grid.self_completeness;
+
   /* Copy the Multipole related information */
   if (with_gravity) {
     const struct gravity_tensors *mp = c->grav.multipole;
@@ -140,6 +145,47 @@ int cell_pack_tags(const struct cell *c, int *tags) {
 #endif
 }
 
+/**
+ * @brief Pack the extra grid info of the given cell and all it's sub-cells.
+ *
+ * @param c The #cell.
+ * @param info Pointer to an array of packed extra grid info (construction
+ * levels).
+ *
+ * @return The number of packed tags.
+ */
+int cell_pack_grid_extra(const struct cell *c,
+                         enum grid_construction_level *info) {
+#ifdef WITH_MPI
+
+  /* Start by packing the construction level of the current cell. */
+  if (c->grid.construction_level == NULL) {
+    info[0] = grid_above_construction_level;
+  } else if (c->grid.construction_level == c) {
+    info[0] = grid_on_construction_level;
+  } else {
+    info[0] = grid_below_construction_level;
+  }
+
+  /* Fill in the progeny, depth-first recursion. */
+  int count = 1;
+  for (int k = 0; k < 8; k++)
+    if (c->progeny[k] != NULL)
+      count += cell_pack_grid_extra(c->progeny[k], &info[count]);
+
+#ifdef SWIFT_DEBUG_CHECKS
+  if (c->mpi.pcell_size != count) error("Inconsistent info and pcell count!");
+#endif  // SWIFT_DEBUG_CHECKS
+
+  /* Return the number of packed tags used. */
+  return count;
+
+#else
+  error("SWIFT was not compiled with MPI support.");
+  return 0;
+#endif
+}
+
 void cell_pack_part_swallow(const struct cell *c,
                             struct black_holes_part_data *data) {
 
@@ -203,7 +249,7 @@ int cell_unpack(struct pcell *restrict pc, struct cell *restrict c,
   c->hydro.h_max = pc->hydro.h_max;
   c->stars.h_max = pc->stars.h_max;
   c->black_holes.h_max = pc->black_holes.h_max;
-  c->sinks.r_cut_max = pc->sinks.r_cut_max;
+  c->sinks.h_max = pc->sinks.h_max;
 
   c->hydro.ti_end_min = pc->hydro.ti_end_min;
   c->grav.ti_end_min = pc->grav.ti_end_min;
@@ -227,6 +273,8 @@ int cell_unpack(struct pcell *restrict pc, struct cell *restrict c,
   c->black_holes.count = pc->black_holes.count;
   c->maxdepth = pc->maxdepth;
 
+  c->grid.self_completeness = pc->grid.self_completeness;
+
 #ifdef SWIFT_DEBUG_CHECKS
   c->cellID = pc->cellID;
 #endif
@@ -258,6 +306,7 @@ int cell_unpack(struct pcell *restrict pc, struct cell *restrict c,
       temp->hydro.count = 0;
       temp->grav.count = 0;
       temp->stars.count = 0;
+      temp->sinks.count = 0;
       temp->loc[0] = c->loc[0];
       temp->loc[1] = c->loc[1];
       temp->loc[2] = c->loc[2];
@@ -265,6 +314,8 @@ int cell_unpack(struct pcell *restrict pc, struct cell *restrict c,
       temp->width[1] = c->width[1] / 2;
       temp->width[2] = c->width[2] / 2;
       temp->dmin = c->dmin / 2;
+      temp->h_min_allowed = temp->dmin * 0.5 * (1. / kernel_gamma);
+      temp->h_max_allowed = temp->dmin * (1. / kernel_gamma);
       if (k & 4) temp->loc[0] += temp->width[0];
       if (k & 2) temp->loc[1] += temp->width[1];
       if (k & 1) temp->loc[2] += temp->width[2];
@@ -274,6 +325,7 @@ int cell_unpack(struct pcell *restrict pc, struct cell *restrict c,
       temp->hydro.dx_max_sort = 0.f;
       temp->stars.dx_max_part = 0.f;
       temp->stars.dx_max_sort = 0.f;
+      temp->sinks.dx_max_part = 0.f;
       temp->black_holes.dx_max_part = 0.f;
       temp->nodeID = c->nodeID;
       temp->parent = c;
@@ -329,11 +381,67 @@ int cell_unpack_tags(const int *tags, struct cell *restrict c) {
 #endif
 }
 
+/**
+ * @brief Unpack the extra grid info of a given cell and its sub-cells.
+ *
+ * @param tags An array of extra grid info (construction levels).
+ * @param c The #cell in which to unpack the tags.
+ *
+ * @return The number of tags created.
+ */
+int cell_unpack_grid_extra(const enum grid_construction_level *info,
+                           struct cell *c, struct cell *construction_level) {
+#ifdef WITH_MPI
+
+  /* Unpack the current pcell. */
+  if (info[0] == grid_on_construction_level) {
+    construction_level = c;
+  } else if (info[0] == grid_above_construction_level) {
+#ifdef SWIFT_DEBUG_CHECKS
+    if (construction_level != NULL)
+      error(
+          "Above construction level, but construction level cell pointer is "
+          "not NULL!");
+#endif
+  } else {
+#ifdef SWIFT_DEBUG_CHECKS
+    if (info[0] != grid_below_construction_level)
+      error("Invalid construction level!");
+    if (construction_level == NULL || construction_level == c)
+      error("Invalid construction level cell pointer!");
+#endif
+  }
+
+  c->grid.construction_level = construction_level;
+
+  /* Number of new cells created. */
+  int count = 1;
+
+  /* Fill the progeny recursively, depth-first. */
+  for (int k = 0; k < 8; k++)
+    if (c->progeny[k] != NULL) {
+      count += cell_unpack_grid_extra(&info[count], c->progeny[k],
+                                      construction_level);
+    }
+
+#ifdef SWIFT_DEBUG_CHECKS
+  if (c->mpi.pcell_size != count) error("Inconsistent info and pcell count!");
+#endif  // SWIFT_DEBUG_CHECKS
+
+  /* Return the total number of unpacked tags. */
+  return count;
+
+#else
+  error("SWIFT was not compiled with MPI support.");
+  return 0;
+#endif
+}
+
 /**
  * @brief Pack the cell information about time-step sizes and displacements
  * of a cell hierarchy.
  *
- * @param c The #cells to pack.
+ * @param c The #cell's to pack.
  * @param pcells the packed cell structures to pack into.
  *
  * @return The number of cells that were packed.
@@ -356,6 +464,9 @@ int cell_pack_end_step(const struct cell *c, struct pcell_step *pcells) {
   pcells[0].black_holes.ti_end_min = c->black_holes.ti_end_min;
   pcells[0].black_holes.dx_max_part = c->black_holes.dx_max_part;
 
+  pcells[0].sinks.ti_end_min = c->sinks.ti_end_min;
+  pcells[0].sinks.dx_max_part = c->sinks.dx_max_part;
+
   /* Fill in the progeny, depth-first recursion. */
   int count = 1;
   for (int k = 0; k < 8; k++)
@@ -376,7 +487,7 @@ int cell_pack_end_step(const struct cell *c, struct pcell_step *pcells) {
  * @brief Unpack the cell information about time-step sizes and displacements
  * of a cell hierarchy.
  *
- * @param c The #cells to unpack into.
+ * @param c The #cell's to unpack into.
  * @param pcells the packed cell structures to unpack from.
  *
  * @return The number of cells that were packed.
@@ -400,6 +511,9 @@ int cell_unpack_end_step(struct cell *c, const struct pcell_step *pcells) {
   c->black_holes.ti_end_min = pcells[0].black_holes.ti_end_min;
   c->black_holes.dx_max_part = pcells[0].black_holes.dx_max_part;
 
+  c->sinks.ti_end_min = pcells[0].sinks.ti_end_min;
+  c->sinks.dx_max_part = pcells[0].sinks.dx_max_part;
+
   /* Fill in the progeny, depth-first recursion. */
   int count = 1;
   for (int k = 0; k < 8; k++)
@@ -460,6 +574,54 @@ void cell_unpack_timebin(struct cell *const c, timebin_t *const t) {
 #endif
 }
 
+/**
+ * @brief Pack the gpart information of the given cell.
+ *
+ * b needs to be aligned on SWIFT_CACHE_ALIGNMENT.
+ *
+ * @param c The #cell.
+ * @param b (output) The array of #gpart_foreign we pack into.
+ */
+void cell_pack_gpart(const struct cell *const c,
+                     struct gpart_foreign *const b) {
+
+#ifdef WITH_MPI
+
+  swift_declare_aligned_ptr(struct gpart_foreign, b_align, b,
+                            SWIFT_CACHE_ALIGNMENT);
+
+  for (int i = 0; i < c->grav.count; ++i)
+    gravity_foreign_copy(&b_align[i], &c->grav.parts[i]);
+
+#else
+  error("SWIFT was not compiled with MPI support.");
+#endif
+}
+
+/**
+ * @brief Pack the gpart FOF information of the given cell.
+ *
+ * b needs to be aligned on SWIFT_CACHE_ALIGNMENT.
+ *
+ * @param c The #cell.
+ * @param b (output) The array of #gpart_fof_foreign we pack into.
+ */
+void cell_pack_fof_gpart(const struct cell *const c,
+                         struct gpart_fof_foreign *const b) {
+
+#ifdef WITH_MPI
+
+  swift_declare_aligned_ptr(struct gpart_fof_foreign, b_align, b,
+                            SWIFT_CACHE_ALIGNMENT);
+
+  for (int i = 0; i < c->grav.count; ++i)
+    gravity_foreign_fof_copy(&b_align[i], &c->grav.parts[i]);
+
+#else
+  error("SWIFT was not compiled with MPI support.");
+#endif
+}
+
 /**
  * @brief Pack the multipole information of the given cell and all it's
  * sub-cells.
@@ -532,39 +694,108 @@ int cell_unpack_multipoles(struct cell *restrict c,
  *
  * @return The number of packed cells.
  */
-int cell_pack_sf_counts(struct cell *restrict c,
-                        struct pcell_sf *restrict pcells) {
+int cell_pack_sf_counts(struct cell *c, struct pcell_sf_stars *pcells) {
 
 #ifdef WITH_MPI
 
   /* Pack this cell's data. */
-  pcells[0].stars.delta_from_rebuild = c->stars.parts - c->stars.parts_rebuild;
-  pcells[0].stars.count = c->stars.count;
-  pcells[0].stars.dx_max_part = c->stars.dx_max_part;
-
-  /* Pack this cell's data. */
-  pcells[0].grav.delta_from_rebuild = c->grav.parts - c->grav.parts_rebuild;
-  pcells[0].grav.count = c->grav.count;
+  pcells[0].delta_from_rebuild = c->stars.parts - c->stars.parts_rebuild;
+  pcells[0].count = c->stars.count;
+  pcells[0].dx_max_part = c->stars.dx_max_part;
 
 #ifdef SWIFT_DEBUG_CHECKS
   /* Stars */
   if (c->stars.parts_rebuild == NULL)
     error("Star particles array at rebuild is NULL! c->depth=%d", c->depth);
 
-  if (pcells[0].stars.delta_from_rebuild < 0)
+  if (pcells[0].delta_from_rebuild < 0)
     error("Stars part pointer moved in the wrong direction!");
 
-  if (pcells[0].stars.delta_from_rebuild > 0 && c->depth == 0)
+  if (pcells[0].delta_from_rebuild > 0 && c->depth == 0)
     error("Shifting the top-level pointer is not allowed!");
+#endif
+
+  /* Fill in the progeny, depth-first recursion. */
+  int count = 1;
+  for (int k = 0; k < 8; k++)
+    if (c->progeny[k] != NULL) {
+      count += cell_pack_sf_counts(c->progeny[k], &pcells[count]);
+    }
+
+  /* Return the number of packed values. */
+  return count;
+
+#else
+  error("SWIFT was not compiled with MPI support.");
+  return 0;
+#endif
+}
+
+/**
+ * @brief Unpack the counts for star formation of a given cell and its
+ * sub-cells.
+ *
+ * @param c The #cell
+ * @param pcells The multipole information to unpack
+ *
+ * @return The number of cells created.
+ */
+int cell_unpack_sf_counts(struct cell *c, struct pcell_sf_stars *pcells) {
+
+#ifdef WITH_MPI
 
+#ifdef SWIFT_DEBUG_CHECKS
+  if (c->stars.parts_rebuild == NULL)
+    error("Star particles array at rebuild is NULL!");
+#endif
+
+  /* Unpack this cell's data. */
+  c->stars.count = pcells[0].count;
+  c->stars.parts = c->stars.parts_rebuild + pcells[0].delta_from_rebuild;
+  c->stars.dx_max_part = pcells[0].dx_max_part;
+
+  /* Fill in the progeny, depth-first recursion. */
+  int count = 1;
+  for (int k = 0; k < 8; k++)
+    if (c->progeny[k] != NULL) {
+      count += cell_unpack_sf_counts(c->progeny[k], &pcells[count]);
+    }
+
+  /* Return the number of packed values. */
+  return count;
+
+#else
+  error("SWIFT was not compiled with MPI support.");
+  return 0;
+#endif
+}
+
+/**
+ * @brief Pack the counts for star formation of the given cell and all it's
+ * sub-cells.
+ *
+ * @param c The #cell.
+ * @param pcells (output) The multipole information we pack into
+ *
+ * @return The number of packed cells.
+ */
+int cell_pack_grav_counts(struct cell *c, struct pcell_sf_grav *pcells) {
+
+#ifdef WITH_MPI
+
+  /* Pack this cell's data. */
+  pcells[0].delta_from_rebuild = c->grav.parts - c->grav.parts_rebuild;
+  pcells[0].count = c->grav.count;
+
+#ifdef SWIFT_DEBUG_CHECKS
   /* Grav */
   if (c->grav.parts_rebuild == NULL)
     error("Grav. particles array at rebuild is NULL! c->depth=%d", c->depth);
 
-  if (pcells[0].grav.delta_from_rebuild < 0)
+  if (pcells[0].delta_from_rebuild < 0)
     error("Grav part pointer moved in the wrong direction!");
 
-  if (pcells[0].grav.delta_from_rebuild > 0 && c->depth == 0)
+  if (pcells[0].delta_from_rebuild > 0 && c->depth == 0)
     error("Shifting the top-level pointer is not allowed!");
 #endif
 
@@ -572,7 +803,7 @@ int cell_pack_sf_counts(struct cell *restrict c,
   int count = 1;
   for (int k = 0; k < 8; k++)
     if (c->progeny[k] != NULL) {
-      count += cell_pack_sf_counts(c->progeny[k], &pcells[count]);
+      count += cell_pack_grav_counts(c->progeny[k], &pcells[count]);
     }
 
   /* Return the number of packed values. */
@@ -593,31 +824,24 @@ int cell_pack_sf_counts(struct cell *restrict c,
  *
  * @return The number of cells created.
  */
-int cell_unpack_sf_counts(struct cell *restrict c,
-                          struct pcell_sf *restrict pcells) {
+int cell_unpack_grav_counts(struct cell *c, struct pcell_sf_grav *pcells) {
 
 #ifdef WITH_MPI
 
 #ifdef SWIFT_DEBUG_CHECKS
   if (c->stars.parts_rebuild == NULL)
     error("Star particles array at rebuild is NULL!");
-  if (c->grav.parts_rebuild == NULL)
-    error("Grav particles array at rebuild is NULL!");
 #endif
 
   /* Unpack this cell's data. */
-  c->stars.count = pcells[0].stars.count;
-  c->stars.parts = c->stars.parts_rebuild + pcells[0].stars.delta_from_rebuild;
-  c->stars.dx_max_part = pcells[0].stars.dx_max_part;
-
-  c->grav.count = pcells[0].grav.count;
-  c->grav.parts = c->grav.parts_rebuild + pcells[0].grav.delta_from_rebuild;
+  c->grav.count = pcells[0].count;
+  c->grav.parts = c->grav.parts_rebuild + pcells[0].delta_from_rebuild;
 
   /* Fill in the progeny, depth-first recursion. */
   int count = 1;
   for (int k = 0; k < 8; k++)
     if (c->progeny[k] != NULL) {
-      count += cell_unpack_sf_counts(c->progeny[k], &pcells[count]);
+      count += cell_unpack_grav_counts(c->progeny[k], &pcells[count]);
     }
 
   /* Return the number of packed values. */

src/cell_sinks.h

@@ -42,6 +42,9 @@ struct cell_sinks {
     /*! Pointer to the #sink data. */
     struct sink *parts;
 
+    /*! Pointer to the #spart data at rebuild time. */
+    struct sink *parts_rebuild;
+
     /*! Linked list of the tasks computing this cell's sink swallow. */
     struct link *swallow;
 
@@ -57,6 +60,12 @@ struct cell_sinks {
     /*! Implicit tasks marking the entry of the sink block of tasks */
     struct task *sink_in;
 
+    /*! The sink ghost task itself */
+    struct task *density_ghost;
+
+    /*! Linked list of the tasks computing this cell's sink density. */
+    struct link *density;
+
     /*! Implicit tasks marking the end of sink swallow */
     struct task *sink_ghost1;
 
@@ -82,11 +91,12 @@ struct cell_sinks {
     /*! Nr of #sink this cell can hold after addition of new one. */
     int count_total;
 
-    /*! Max cut off radius of active particles in this cell. */
-    float r_cut_max_active;
+    /*! Max smoothing length of active particles in this cell.
+     */
+    float h_max_active;
 
-    /*! Values of r_cut_max before the drifts, used for sub-cell tasks. */
-    float r_cut_max_old;
+    /*! Values of h_max before the drifts, used for sub-cell tasks. */
+    float h_max_old;
 
     /*! Maximum part movement in this cell since last construction. */
     float dx_max_part;
@@ -108,8 +118,8 @@ struct cell_sinks {
   /*! Spin lock for various uses (#sink case). */
   swift_lock_type lock;
 
-  /*! Max cut off radius in this cell. */
-  float r_cut_max;
+  /*! Max smoothing length in this cell. */
+  float h_max;
 
   /*! Number of #sink updated in this cell. */
   int updated;

src/cell_split.c

@@ -384,6 +384,7 @@ void cell_split(struct cell *c, const ptrdiff_t parts_offset,
     c->progeny[k]->sinks.count = bucket_count[k];
     c->progeny[k]->sinks.count_total = c->progeny[k]->sinks.count;
     c->progeny[k]->sinks.parts = &c->sinks.parts[bucket_offset[k]];
+    c->progeny[k]->sinks.parts_rebuild = c->progeny[k]->sinks.parts;
   }
 
   /* Finally, do the same song and dance for the gparts. */
@@ -636,7 +637,8 @@ void cell_reorder_extra_sinks(struct cell *c, const ptrdiff_t sinks_offset) {
  * @param sinks The global array of #sink (for re-linking).
  */
 void cell_reorder_extra_gparts(struct cell *c, struct part *parts,
-                               struct spart *sparts, struct sink *sinks) {
+                               struct spart *sparts, struct sink *sinks,
+                               struct bpart *bparts) {
   struct gpart *gparts = c->grav.parts;
   const int count_real = c->grav.count;
 
@@ -668,6 +670,8 @@ void cell_reorder_extra_gparts(struct cell *c, struct part *parts,
         sinks[-gparts[i].id_or_neg_offset].gpart = &gparts[i];
       } else if (gparts[i].type == swift_type_stars) {
         sparts[-gparts[i].id_or_neg_offset].gpart = &gparts[i];
+      } else if (gparts[i].type == swift_type_black_hole) {
+        bparts[-gparts[i].id_or_neg_offset].gpart = &gparts[i];
       }
     }
   }

src/cell_unskip.c

@@ -847,7 +847,7 @@ void cell_activate_subcell_hydro_tasks(struct cell *ci, struct cell *cj,
 
     /* Get the orientation of the pair. */
     double shift[3];
-    const int sid = space_getsid(s->space, &ci, &cj, shift);
+    const int sid = space_getsid_and_swap_cells(s->space, &ci, &cj, shift);
 
     /* recurse? */
     if (cell_can_recurse_in_pair_hydro_task(ci) &&
@@ -959,7 +959,7 @@ void cell_activate_subcell_stars_tasks(struct cell *ci, struct cell *cj,
 
     /* Get the orientation of the pair. */
     double shift[3];
-    const int sid = space_getsid(s->space, &ci, &cj, shift);
+    const int sid = space_getsid_and_swap_cells(s->space, &ci, &cj, shift);
 
     const int ci_active = cell_need_activating_stars(ci, e, with_star_formation,
                                                      with_star_formation_sink);
@@ -970,8 +970,8 @@ void cell_activate_subcell_stars_tasks(struct cell *ci, struct cell *cj,
     if (!ci_active && !cj_active) return;
 
     /* recurse? */
-    if (cell_can_recurse_in_pair_stars_task(ci, cj) &&
-        cell_can_recurse_in_pair_stars_task(cj, ci)) {
+    if (cell_can_recurse_in_pair_stars_task(ci) &&
+        cell_can_recurse_in_pair_stars_task(cj)) {
 
       const struct cell_split_pair *csp = &cell_split_pairs[sid];
       for (int k = 0; k < csp->count; k++) {
@@ -1087,14 +1087,16 @@ void cell_activate_subcell_black_holes_tasks(struct cell *ci, struct cell *cj,
 
   /* Otherwise, pair interation */
   else {
-    /* Should we even bother? */
-    if (!cell_is_active_black_holes(ci, e) &&
-        !cell_is_active_black_holes(cj, e))
-      return;
 
     /* Get the orientation of the pair. */
     double shift[3];
-    const int sid = space_getsid(s->space, &ci, &cj, shift);
+    const int sid = space_getsid_and_swap_cells(s->space, &ci, &cj, shift);
+
+    const int ci_active = cell_is_active_black_holes(ci, e);
+    const int cj_active = cell_is_active_black_holes(cj, e);
+
+    /* Should we even bother? */
+    if (!ci_active && !cj_active) return;
 
     /* recurse? */
     if (cell_can_recurse_in_pair_black_holes_task(ci, cj) &&
@@ -1110,19 +1112,24 @@ void cell_activate_subcell_black_holes_tasks(struct cell *ci, struct cell *cj,
     }
 
     /* Otherwise, activate the drifts. */
-    else if (cell_is_active_black_holes(ci, e) ||
-             cell_is_active_black_holes(cj, e)) {
+    else if (ci_active || cj_active) {
+
+      /* Note we need to drift *both* BH cells to deal with BH<->BH swallows
+       * But we only need to drift the gas cell if the *other* cell has an
+       * active BH */
 
       /* Activate the drifts if the cells are local. */
       if (ci->nodeID == engine_rank) cell_activate_drift_bpart(ci, s);
-      if (cj->nodeID == engine_rank) cell_activate_drift_part(cj, s);
-      if (cj->nodeID == engine_rank && with_timestep_sync)
+      if (cj->nodeID == engine_rank && ci_active)
+        cell_activate_drift_part(cj, s);
+      if (cj->nodeID == engine_rank && ci_active && with_timestep_sync)
         cell_activate_sync_part(cj, s);
 
       /* Activate the drifts if the cells are local. */
-      if (ci->nodeID == engine_rank) cell_activate_drift_part(ci, s);
+      if (ci->nodeID == engine_rank && cj_active)
+        cell_activate_drift_part(ci, s);
       if (cj->nodeID == engine_rank) cell_activate_drift_bpart(cj, s);
-      if (ci->nodeID == engine_rank && with_timestep_sync)
+      if (ci->nodeID == engine_rank && cj_active && with_timestep_sync)
         cell_activate_sync_part(ci, s);
     }
   } /* Otherwise, pair interation */
@@ -1144,13 +1151,13 @@ void cell_activate_subcell_sinks_tasks(struct cell *ci, struct cell *cj,
 
   /* Store the current dx_max and h_max values. */
   ci->sinks.dx_max_part_old = ci->sinks.dx_max_part;
-  ci->sinks.r_cut_max_old = ci->sinks.r_cut_max;
+  ci->sinks.h_max_old = ci->sinks.h_max;
   ci->hydro.dx_max_part_old = ci->hydro.dx_max_part;
   ci->hydro.h_max_old = ci->hydro.h_max;
 
   if (cj != NULL) {
     cj->sinks.dx_max_part_old = cj->sinks.dx_max_part;
-    cj->sinks.r_cut_max_old = cj->sinks.r_cut_max;
+    cj->sinks.h_max_old = cj->sinks.h_max;
     cj->hydro.dx_max_part_old = cj->hydro.dx_max_part;
     cj->hydro.h_max_old = cj->hydro.h_max;
   }
@@ -1190,7 +1197,7 @@ void cell_activate_subcell_sinks_tasks(struct cell *ci, struct cell *cj,
   else {
     /* Get the orientation of the pair. */
     double shift[3];
-    const int sid = space_getsid(s->space, &ci, &cj, shift);
+    const int sid = space_getsid_and_swap_cells(s->space, &ci, &cj, shift);
 
     const int ci_active =
         cell_is_active_sinks(ci, e) || cell_is_active_hydro(ci, e);
@@ -1579,7 +1586,7 @@ void cell_activate_subcell_rt_tasks(struct cell *ci, struct cell *cj,
 
     /* Get the orientation of the pair. */
     double shift[3];
-    const int sid = space_getsid(s->space, &ci, &cj, shift);
+    const int sid = space_getsid_and_swap_cells(s->space, &ci, &cj, shift);
 
     /* recurse? */
     if (cell_can_recurse_in_pair_hydro_task(ci) &&
@@ -1625,11 +1632,9 @@ int cell_unskip_hydro_tasks(struct cell *c, struct scheduler *s) {
   const int with_feedback = e->policy & engine_policy_feedback;
   const int with_timestep_limiter =
       (e->policy & engine_policy_timestep_limiter);
-  const int with_sinks = e->policy & engine_policy_sinks;
 
 #ifdef WITH_MPI
   const int with_star_formation = e->policy & engine_policy_star_formation;
-  if (with_sinks) error("Cannot use sink tasks and MPI");
 #endif
   int rebuild = 0;
 
@@ -1653,20 +1658,18 @@ int cell_unskip_hydro_tasks(struct cell *c, struct scheduler *s) {
         (cj_active && cj_nodeID == nodeID)) {
       scheduler_activate(s, t);
 
-      /* Activate hydro drift */
+      /* Store current values of dx_max and h_max. */
       if (t->type == task_type_self) {
+        cell_activate_subcell_hydro_tasks(ci, NULL, s, with_timestep_limiter);
+
         if (ci_nodeID == nodeID) cell_activate_drift_part(ci, s);
         if (ci_nodeID == nodeID && with_timestep_limiter)
           cell_activate_limiter(ci, s);
       }
 
-      /* Set the correct sorting flags and activate hydro drifts */
+      /* Store current values of dx_max and h_max. */
       else if (t->type == task_type_pair) {
-        /* Store some values. */
-        atomic_or(&ci->hydro.requires_sorts, 1 << t->flags);
-        atomic_or(&cj->hydro.requires_sorts, 1 << t->flags);
-        ci->hydro.dx_max_sort_old = ci->hydro.dx_max_sort;
-        cj->hydro.dx_max_sort_old = cj->hydro.dx_max_sort;
+        cell_activate_subcell_hydro_tasks(ci, cj, s, with_timestep_limiter);
 
         /* Activate the drift tasks. */
         if (ci_nodeID == nodeID) cell_activate_drift_part(ci, s);
@@ -1677,25 +1680,11 @@ int cell_unskip_hydro_tasks(struct cell *c, struct scheduler *s) {
           cell_activate_limiter(ci, s);
         if (cj_nodeID == nodeID && with_timestep_limiter)
           cell_activate_limiter(cj, s);
-
-        /* Check the sorts and activate them if needed. */
-        cell_activate_hydro_sorts(ci, t->flags, s);
-        cell_activate_hydro_sorts(cj, t->flags, s);
-      }
-
-      /* Store current values of dx_max and h_max. */
-      else if (t->type == task_type_sub_self) {
-        cell_activate_subcell_hydro_tasks(ci, NULL, s, with_timestep_limiter);
-      }
-
-      /* Store current values of dx_max and h_max. */
-      else if (t->type == task_type_sub_pair) {
-        cell_activate_subcell_hydro_tasks(ci, cj, s, with_timestep_limiter);
       }
     }
 
     /* Only interested in pair interactions as of here. */
-    if (t->type == task_type_pair || t->type == task_type_sub_pair) {
+    if (t->type == task_type_pair) {
       /* Check whether there was too much particle motion, i.e. the
          cell neighbour conditions were violated. */
       if (cell_need_rebuild_for_hydro_pair(ci, cj)) rebuild = 1;
@@ -1930,9 +1919,6 @@ int cell_unskip_hydro_tasks(struct cell *c, struct scheduler *s) {
       cell_activate_star_formation_tasks(c->top, s, with_feedback);
       cell_activate_super_spart_drifts(c->top, s);
     }
-    if (with_sinks && c->top->sinks.star_formation_sink != NULL) {
-      cell_activate_star_formation_sink_tasks(c->top, s, with_feedback);
-    }
   }
   /* Additionally unskip force interactions between inactive local cell and
    * active remote cell. (The cell unskip will only be called for active cells,
@@ -1941,7 +1927,8 @@ int cell_unskip_hydro_tasks(struct cell *c, struct scheduler *s) {
 #if defined(MPI_SYMMETRIC_FORCE_INTERACTION) && defined(WITH_MPI)
     for (struct link *l = c->hydro.force; l != NULL; l = l->next) {
       struct task *t = l->t;
-      if (t->type != task_type_pair && t->type != task_type_sub_pair) continue;
+
+      if (t->type != task_type_pair) continue;
 
       struct cell *ci = l->t->ci;
       struct cell *cj = l->t->cj;
@@ -1956,30 +1943,8 @@ int cell_unskip_hydro_tasks(struct cell *c, struct scheduler *s) {
            ci_nodeID != nodeID)) {
         scheduler_activate(s, l->t);
 
-        if (t->type == task_type_pair) {
-          /* Store some values. */
-          atomic_or(&ci->hydro.requires_sorts, 1 << t->flags);
-          atomic_or(&cj->hydro.requires_sorts, 1 << t->flags);
-          ci->hydro.dx_max_sort_old = ci->hydro.dx_max_sort;
-          cj->hydro.dx_max_sort_old = cj->hydro.dx_max_sort;
-
-          /* Activate the drift tasks. */
-          if (ci_nodeID == nodeID) cell_activate_drift_part(ci, s);
-          if (cj_nodeID == nodeID) cell_activate_drift_part(cj, s);
-
-          /* Activate the limiter tasks. */
-          if (ci_nodeID == nodeID && with_timestep_limiter)
-            cell_activate_limiter(ci, s);
-          if (cj_nodeID == nodeID && with_timestep_limiter)
-            cell_activate_limiter(cj, s);
-
-          /* Check the sorts and activate them if needed. */
-          cell_activate_hydro_sorts(ci, t->flags, s);
-          cell_activate_hydro_sorts(cj, t->flags, s);
-        }
-
         /* Store current values of dx_max and h_max. */
-        else if (t->type == task_type_sub_pair) {
+        if (t->type == task_type_pair) {
           cell_activate_subcell_hydro_tasks(ci, cj, s, with_timestep_limiter);
         }
       }
@@ -2004,6 +1969,10 @@ int cell_unskip_gravity_tasks(struct cell *c, struct scheduler *s) {
   const int nodeID = e->nodeID;
   int rebuild = 0;
 
+#ifdef WITH_MPI
+  const int with_star_formation = e->policy & engine_policy_star_formation;
+#endif
+
   /* Un-skip the gravity tasks involved with this cell. */
   for (struct link *l = c->grav.grav; l != NULL; l = l->next) {
     struct task *t = l->t;
@@ -2050,6 +2019,9 @@ int cell_unskip_gravity_tasks(struct cell *c, struct scheduler *s) {
         /* Is the foreign cell active and will need stuff from us? */
         if (ci_active) {
 
+          scheduler_activate_pack(s, cj->mpi.pack, task_subtype_gpart,
+                                  ci_nodeID);
+
           scheduler_activate_send(s, cj->mpi.send, task_subtype_gpart,
                                   ci_nodeID);
 
@@ -2059,6 +2031,17 @@ int cell_unskip_gravity_tasks(struct cell *c, struct scheduler *s) {
           cell_activate_drift_gpart(cj, s);
         }
 
+        /* Propagating new star counts? */
+        if (with_star_formation) {
+          if (ci_active && ci->hydro.count > 0) {
+            scheduler_activate_recv(s, ci->mpi.recv, task_subtype_grav_counts);
+          }
+          if (cj_active && cj->hydro.count > 0) {
+            scheduler_activate_send(s, cj->mpi.send, task_subtype_grav_counts,
+                                    ci_nodeID);
+          }
+        }
+
       } else if (cj_nodeID != nodeID) {
         /* If the local cell is active, receive data from the foreign cell. */
         if (ci_active)
@@ -2067,6 +2050,9 @@ int cell_unskip_gravity_tasks(struct cell *c, struct scheduler *s) {
         /* Is the foreign cell active and will need stuff from us? */
         if (cj_active) {
 
+          scheduler_activate_pack(s, ci->mpi.pack, task_subtype_gpart,
+                                  cj_nodeID);
+
           scheduler_activate_send(s, ci->mpi.send, task_subtype_gpart,
                                   cj_nodeID);
 
@@ -2075,6 +2061,17 @@ int cell_unskip_gravity_tasks(struct cell *c, struct scheduler *s) {
              itself. */
           cell_activate_drift_gpart(ci, s);
         }
+
+        /* Propagating new star counts? */
+        if (with_star_formation) {
+          if (cj_active && cj->hydro.count > 0) {
+            scheduler_activate_recv(s, cj->mpi.recv, task_subtype_grav_counts);
+          }
+          if (ci_active && ci->hydro.count > 0) {
+            scheduler_activate_send(s, ci->mpi.send, task_subtype_grav_counts,
+                                    cj_nodeID);
+          }
+        }
       }
 #endif
     }
@@ -2177,82 +2174,40 @@ int cell_unskip_stars_tasks(struct cell *c, struct scheduler *s,
         (cj != NULL) && cell_need_activating_stars(cj, e, with_star_formation,
                                                    with_star_formation_sink);
 
-    /* Activate the drifts */
-    if (t->type == task_type_self && ci_active) {
-      cell_activate_drift_spart(ci, s);
-      cell_activate_drift_part(ci, s);
-      if (with_timestep_sync) cell_activate_sync_part(ci, s);
-    }
-
     /* Only activate tasks that involve a local active cell. */
     if ((ci_active || cj_active) &&
         (ci_nodeID == nodeID || cj_nodeID == nodeID)) {
       scheduler_activate(s, t);
 
-      if (t->type == task_type_pair) {
-        /* Activate stars_in for each cell that is part of
-         * a pair as to not miss any dependencies */
-        if (ci_nodeID == nodeID)
-          scheduler_activate(s, ci->hydro.super->stars.stars_in);
-        if (cj_nodeID == nodeID)
-          scheduler_activate(s, cj->hydro.super->stars.stars_in);
-
-        /* Do ci */
-        if (ci_active) {
-          /* stars for ci */
-          atomic_or(&ci->stars.requires_sorts, 1 << t->flags);
-          ci->stars.dx_max_sort_old = ci->stars.dx_max_sort;
-
-          /* hydro for cj */
-          atomic_or(&cj->hydro.requires_sorts, 1 << t->flags);
-          cj->hydro.dx_max_sort_old = cj->hydro.dx_max_sort;
-
-          /* Activate the drift tasks. */
-          if (ci_nodeID == nodeID) cell_activate_drift_spart(ci, s);
-          if (cj_nodeID == nodeID) cell_activate_drift_part(cj, s);
-          if (cj_nodeID == nodeID && with_timestep_sync)
-            cell_activate_sync_part(cj, s);
-
-          /* Check the sorts and activate them if needed. */
-          cell_activate_stars_sorts(ci, t->flags, s);
-          cell_activate_hydro_sorts(cj, t->flags, s);
-        }
-
-        /* Do cj */
-        if (cj_active) {
-          /* hydro for ci */
-          atomic_or(&ci->hydro.requires_sorts, 1 << t->flags);
-          ci->hydro.dx_max_sort_old = ci->hydro.dx_max_sort;
-
-          /* stars for cj */
-          atomic_or(&cj->stars.requires_sorts, 1 << t->flags);
-          cj->stars.dx_max_sort_old = cj->stars.dx_max_sort;
-
-          /* Activate the drift tasks. */
-          if (cj_nodeID == nodeID) cell_activate_drift_spart(cj, s);
-          if (ci_nodeID == nodeID) cell_activate_drift_part(ci, s);
-          if (ci_nodeID == nodeID && with_timestep_sync)
-            cell_activate_sync_part(ci, s);
-
-          /* Check the sorts and activate them if needed. */
-          cell_activate_hydro_sorts(ci, t->flags, s);
-          cell_activate_stars_sorts(cj, t->flags, s);
-        }
-      }
-
-      else if (t->type == task_type_sub_self) {
+      if (t->type == task_type_self) {
         cell_activate_subcell_stars_tasks(ci, NULL, s, with_star_formation,
                                           with_star_formation_sink,
                                           with_timestep_sync);
+
+        cell_activate_drift_spart(ci, s);
+        cell_activate_drift_part(ci, s);
+        if (with_timestep_sync) cell_activate_sync_part(ci, s);
       }
 
-      else if (t->type == task_type_sub_pair) {
+      else if (t->type == task_type_pair) {
         cell_activate_subcell_stars_tasks(ci, cj, s, with_star_formation,
                                           with_star_formation_sink,
                                           with_timestep_sync);
 
+        /* Activate the drift tasks. */
+        if (ci_nodeID == nodeID) cell_activate_drift_spart(ci, s);
+        if (cj_nodeID == nodeID) cell_activate_drift_part(cj, s);
+        if (cj_nodeID == nodeID && with_timestep_sync)
+          cell_activate_sync_part(cj, s);
+
+        /* Activate the drift tasks. */
+        if (cj_nodeID == nodeID) cell_activate_drift_spart(cj, s);
+        if (ci_nodeID == nodeID) cell_activate_drift_part(ci, s);
+        if (ci_nodeID == nodeID && with_timestep_sync)
+          cell_activate_sync_part(ci, s);
+
         /* Activate stars_in for each cell that is part of
-         * a sub_pair task as to not miss any dependencies */
+         * a pair task as to not miss any dependencies */
         if (ci_nodeID == nodeID)
           scheduler_activate(s, ci->hydro.super->stars.stars_in);
         if (cj_nodeID == nodeID)
@@ -2261,7 +2216,7 @@ int cell_unskip_stars_tasks(struct cell *c, struct scheduler *s,
     }
 
     /* Only interested in pair interactions as of here. */
-    if (t->type == task_type_pair || t->type == task_type_sub_pair) {
+    if (t->type == task_type_pair) {
       /* Check whether there was too much particle motion, i.e. the
          cell neighbour conditions were violated. */
       if (cell_need_rebuild_for_stars_pair(ci, cj)) rebuild = 1;
@@ -2374,11 +2329,7 @@ int cell_unskip_stars_tasks(struct cell *c, struct scheduler *s,
       scheduler_activate(s, t);
     }
 
-    else if (t->type == task_type_sub_self && ci_active) {
-      scheduler_activate(s, t);
-    }
-
-    else if (t->type == task_type_pair || t->type == task_type_sub_pair) {
+    else if (t->type == task_type_pair) {
       /* We only want to activate the task if the cell is active and is
          going to update some gas on the *local* node */
       if ((ci_nodeID == nodeID && cj_nodeID == nodeID) &&
@@ -2438,11 +2389,7 @@ int cell_unskip_stars_tasks(struct cell *c, struct scheduler *s,
       scheduler_activate(s, t);
     }
 
-    else if (t->type == task_type_sub_self && ci_active) {
-      scheduler_activate(s, t);
-    }
-
-    else if (t->type == task_type_pair || t->type == task_type_sub_pair) {
+    else if (t->type == task_type_pair) {
       /* We only want to activate the task if the cell is active and is
          going to update some gas on the *local* node */
       if ((ci_nodeID == nodeID && cj_nodeID == nodeID) &&
@@ -2484,15 +2431,11 @@ int cell_unskip_stars_tasks(struct cell *c, struct scheduler *s,
       scheduler_activate(s, t);
     }
 
-    else if (t->type == task_type_sub_self && ci_active) {
-      scheduler_activate(s, t);
-    }
-
-    else if (t->type == task_type_pair || t->type == task_type_sub_pair) {
+    else if (t->type == task_type_pair) {
 
       if (ci_active || cj_active) {
         /* Activate stars_out for each cell that is part of
-         * a pair/sub_pair task as to not miss any dependencies */
+         * a pair task as to not miss any dependencies */
         if (ci_nodeID == nodeID)
           scheduler_activate(s, ci->hydro.super->stars.stars_out);
         if (cj_nodeID == nodeID)
@@ -2563,7 +2506,8 @@ int cell_unskip_black_holes_tasks(struct cell *c, struct scheduler *s) {
   const int nodeID = e->nodeID;
   int rebuild = 0;
 
-  if (c->black_holes.drift != NULL && cell_is_active_black_holes(c, e)) {
+  if (c->black_holes.drift != NULL && c->black_holes.count > 0 &&
+      cell_is_active_black_holes(c, e)) {
     cell_activate_drift_bpart(c, s);
   }
 
@@ -2572,8 +2516,11 @@ int cell_unskip_black_holes_tasks(struct cell *c, struct scheduler *s) {
     struct task *t = l->t;
     struct cell *ci = t->ci;
     struct cell *cj = t->cj;
-    const int ci_active = cell_is_active_black_holes(ci, e);
-    const int cj_active = (cj != NULL) ? cell_is_active_black_holes(cj, e) : 0;
+    const int ci_active =
+        ci->black_holes.count > 0 && cell_is_active_black_holes(ci, e);
+    const int cj_active = (cj != NULL) ? (cj->black_holes.count > 0 &&
+                                          cell_is_active_black_holes(cj, e))
+                                       : 0;
 #ifdef WITH_MPI
     const int ci_nodeID = ci->nodeID;
     const int cj_nodeID = (cj != NULL) ? cj->nodeID : -1;
@@ -2588,117 +2535,129 @@ int cell_unskip_black_holes_tasks(struct cell *c, struct scheduler *s) {
 
       scheduler_activate(s, t);
 
-      /* Activate the drifts */
-      if (t->type == task_type_self) {
-        cell_activate_drift_part(ci, s);
-        cell_activate_drift_bpart(ci, s);
-        if (with_timestep_sync) cell_activate_sync_part(ci, s);
-      }
-
-      /* Activate the drifts */
-      else if (t->type == task_type_pair) {
-
-        /* Activate the drift tasks. */
-        if (ci_nodeID == nodeID) cell_activate_drift_bpart(ci, s);
-        if (ci_nodeID == nodeID) cell_activate_drift_part(ci, s);
-
-        if (cj_nodeID == nodeID) cell_activate_drift_part(cj, s);
-        if (cj_nodeID == nodeID) cell_activate_drift_bpart(cj, s);
-      }
-
       /* Store current values of dx_max and h_max. */
-      else if (t->type == task_type_sub_self) {
+      if (t->type == task_type_self) {
         cell_activate_subcell_black_holes_tasks(ci, NULL, s,
                                                 with_timestep_sync);
       }
 
       /* Store current values of dx_max and h_max. */
-      else if (t->type == task_type_sub_pair) {
+      else if (t->type == task_type_pair) {
         cell_activate_subcell_black_holes_tasks(ci, cj, s, with_timestep_sync);
+
+        /* Activate bh_in for each cell that is part of
+         * a pair task as to not miss any dependencies */
+        if (ci_nodeID == nodeID)
+          scheduler_activate(s, ci->hydro.super->black_holes.black_holes_in);
+        if (cj_nodeID == nodeID)
+          scheduler_activate(s, cj->hydro.super->black_holes.black_holes_in);
       }
     }
 
     /* Only interested in pair interactions as of here. */
-    if (t->type == task_type_pair || t->type == task_type_sub_pair) {
-
-      /* Activate bh_in for each cell that is part of
-       * a pair task as to not miss any dependencies */
-      if (ci_nodeID == nodeID)
-        scheduler_activate(s, ci->hydro.super->black_holes.black_holes_in);
-      if (cj_nodeID == nodeID)
-        scheduler_activate(s, cj->hydro.super->black_holes.black_holes_in);
+    if (t->type == task_type_pair) {
 
       /* Check whether there was too much particle motion, i.e. the
          cell neighbour conditions were violated. */
       if (cell_need_rebuild_for_black_holes_pair(ci, cj)) rebuild = 1;
       if (cell_need_rebuild_for_black_holes_pair(cj, ci)) rebuild = 1;
 
-      scheduler_activate(s, ci->hydro.super->black_holes.swallow_ghost_0);
-      scheduler_activate(s, cj->hydro.super->black_holes.swallow_ghost_0);
+      if (ci->hydro.super->black_holes.count > 0 && ci_active)
+        scheduler_activate(s, ci->hydro.super->black_holes.swallow_ghost_1);
+      if (cj->hydro.super->black_holes.count > 0 && cj_active)
+        scheduler_activate(s, cj->hydro.super->black_holes.swallow_ghost_1);
 
 #ifdef WITH_MPI
       /* Activate the send/recv tasks. */
       if (ci_nodeID != nodeID) {
 
-        /* Receive the foreign parts to compute BH accretion rates and do the
-         * swallowing */
-        scheduler_activate_recv(s, ci->mpi.recv, task_subtype_rho);
-        scheduler_activate_recv(s, ci->mpi.recv, task_subtype_part_swallow);
+        if (ci_active || cj_active) {
+          /* We must exchange the foreign BHs no matter the activity status */
+          scheduler_activate_recv(s, ci->mpi.recv, task_subtype_bpart_rho);
+          scheduler_activate_send(s, cj->mpi.send, task_subtype_bpart_rho,
+                                  ci_nodeID);
 
-        /* Send the local BHs to tag the particles to swallow and to do feedback
-         */
-        scheduler_activate_send(s, cj->mpi.send, task_subtype_bpart_rho,
-                                ci_nodeID);
-        scheduler_activate_send(s, cj->mpi.send, task_subtype_bpart_feedback,
-                                ci_nodeID);
+          /* Drift before you send */
+          if (cj->black_holes.count > 0) cell_activate_drift_bpart(cj, s);
+        }
 
-        /* Drift before you send */
-        cell_activate_drift_bpart(cj, s);
+        if (cj_active) {
 
-        /* Receive the foreign BHs to tag particles to swallow and for feedback
-         */
-        scheduler_activate_recv(s, ci->mpi.recv, task_subtype_bpart_rho);
-        scheduler_activate_recv(s, ci->mpi.recv, task_subtype_bpart_feedback);
+          /* Receive the foreign parts to compute BH accretion rates and do the
+           * swallowing */
+          scheduler_activate_recv(s, ci->mpi.recv, task_subtype_rho);
+          scheduler_activate_recv(s, ci->mpi.recv, task_subtype_part_swallow);
+          scheduler_activate_recv(s, ci->mpi.recv, task_subtype_bpart_merger);
 
-        /* Send the local part information */
-        scheduler_activate_send(s, cj->mpi.send, task_subtype_rho, ci_nodeID);
-        scheduler_activate_send(s, cj->mpi.send, task_subtype_part_swallow,
-                                ci_nodeID);
+          /* Send the local BHs to do feedback */
+          scheduler_activate_send(s, cj->mpi.send, task_subtype_bpart_feedback,
+                                  ci_nodeID);
 
-        /* Drift the cell which will be sent; note that not all sent
-           particles will be drifted, only those that are needed. */
-        cell_activate_drift_part(cj, s);
+          /* Drift before you send */
+          cell_activate_drift_bpart(cj, s);
+        }
+
+        if (ci_active) {
+
+          /* Receive the foreign BHs for feedback */
+          scheduler_activate_recv(s, ci->mpi.recv, task_subtype_bpart_feedback);
+
+          /* Send the local part information */
+          scheduler_activate_send(s, cj->mpi.send, task_subtype_rho, ci_nodeID);
+          scheduler_activate_send(s, cj->mpi.send, task_subtype_part_swallow,
+                                  ci_nodeID);
+          scheduler_activate_send(s, cj->mpi.send, task_subtype_bpart_merger,
+                                  ci_nodeID);
+
+          /* Drift the cell which will be sent; note that not all sent
+             particles will be drifted, only those that are needed. */
+          if (cj->hydro.count > 0) cell_activate_drift_part(cj, s);
+        }
 
       } else if (cj_nodeID != nodeID) {
 
-        /* Receive the foreign parts to compute BH accretion rates and do the
-         * swallowing */
-        scheduler_activate_recv(s, cj->mpi.recv, task_subtype_rho);
-        scheduler_activate_recv(s, cj->mpi.recv, task_subtype_part_swallow);
+        if (ci_active || cj_active) {
+          /* We must exchange the foreign BHs no matter the activity status */
+          scheduler_activate_recv(s, cj->mpi.recv, task_subtype_bpart_rho);
+          scheduler_activate_send(s, ci->mpi.send, task_subtype_bpart_rho,
+                                  cj_nodeID);
 
-        /* Send the local BHs to tag the particles to swallow and to do feedback
-         */
-        scheduler_activate_send(s, ci->mpi.send, task_subtype_bpart_rho,
-                                cj_nodeID);
-        scheduler_activate_send(s, ci->mpi.send, task_subtype_bpart_feedback,
-                                cj_nodeID);
+          /* Drift before you send */
+          if (ci->black_holes.count > 0) cell_activate_drift_bpart(ci, s);
+        }
 
-        /* Drift before you send */
-        cell_activate_drift_bpart(ci, s);
+        if (ci_active) {
 
-        /* Receive the foreign BHs to tag particles to swallow and for feedback
-         */
-        scheduler_activate_recv(s, cj->mpi.recv, task_subtype_bpart_rho);
-        scheduler_activate_recv(s, cj->mpi.recv, task_subtype_bpart_feedback);
+          /* Receive the foreign parts to compute BH accretion rates and do the
+           * swallowing */
+          scheduler_activate_recv(s, cj->mpi.recv, task_subtype_rho);
+          scheduler_activate_recv(s, cj->mpi.recv, task_subtype_part_swallow);
+          scheduler_activate_recv(s, cj->mpi.recv, task_subtype_bpart_merger);
 
-        /* Send the local part information */
-        scheduler_activate_send(s, ci->mpi.send, task_subtype_rho, cj_nodeID);
-        scheduler_activate_send(s, ci->mpi.send, task_subtype_part_swallow,
-                                cj_nodeID);
+          /* Send the local BHs to do feedback */
+          scheduler_activate_send(s, ci->mpi.send, task_subtype_bpart_feedback,
+                                  cj_nodeID);
 
-        /* Drift the cell which will be sent; note that not all sent
-           particles will be drifted, only those that are needed. */
-        cell_activate_drift_part(ci, s);
+          /* Drift before you send */
+          cell_activate_drift_bpart(ci, s);
+        }
+
+        if (cj_active) {
+
+          /* Receive the foreign BHs for feedback */
+          scheduler_activate_recv(s, cj->mpi.recv, task_subtype_bpart_feedback);
+
+          /* Send the local part information */
+          scheduler_activate_send(s, ci->mpi.send, task_subtype_rho, cj_nodeID);
+          scheduler_activate_send(s, ci->mpi.send, task_subtype_part_swallow,
+                                  cj_nodeID);
+          scheduler_activate_send(s, ci->mpi.send, task_subtype_bpart_merger,
+                                  cj_nodeID);
+
+          /* Drift the cell which will be sent; note that not all sent
+             particles will be drifted, only those that are needed. */
+          if (ci->hydro.count > 0) cell_activate_drift_part(ci, s);
+        }
       }
 #endif
     }
@@ -2794,9 +2753,9 @@ int cell_unskip_black_holes_tasks(struct cell *c, struct scheduler *s) {
 
       scheduler_activate(s, t);
 
-      if (t->type == task_type_pair || t->type == task_type_sub_pair) {
+      if (t->type == task_type_pair) {
         /* Activate bh_out for each cell that is part of
-         * a pair/sub_pair task as to not miss any dependencies */
+         * a pair task as to not miss any dependencies */
         if (ci_nodeID == nodeID)
           scheduler_activate(s, ci->hydro.super->black_holes.black_holes_out);
         if (cj_nodeID == nodeID)
@@ -2806,22 +2765,27 @@ 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)) {
-
-    /* If the cell doesn't have any pair/sub_pair type tasks,
-     * then we haven't unskipped all the implicit tasks yet. */
+  if (cell_is_active_black_holes(c, e)) {
     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.swallow_ghost_3 != NULL)
+      scheduler_activate(s, c->black_holes.swallow_ghost_3);
+  }
+  if (c->nodeID == nodeID && cell_is_active_black_holes(c, e)) {
     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)
       scheduler_activate(s, c->black_holes.black_holes_out);
+  }
+  if (c->nodeID == nodeID && c->black_holes.count > 0 &&
+      cell_is_active_black_holes(c, e)) {
+
+    /* If the cell doesn't have any pair type tasks,
+     * then we haven't unskipped all the implicit tasks yet. */
     if (c->kick1 != NULL) scheduler_activate(s, c->kick1);
     if (c->kick2 != NULL) scheduler_activate(s, c->kick2);
     if (c->timestep != NULL) scheduler_activate(s, c->timestep);
@@ -2848,6 +2812,7 @@ int cell_unskip_sinks_tasks(struct cell *c, struct scheduler *s) {
 
   struct engine *e = s->space->e;
   const int with_timestep_sync = (e->policy & engine_policy_timestep_sync);
+  const int with_feedback = e->policy & engine_policy_feedback;
   const int nodeID = e->nodeID;
   int rebuild = 0;
 
@@ -2856,11 +2821,12 @@ int cell_unskip_sinks_tasks(struct cell *c, struct scheduler *s) {
       cell_activate_drift_sink(c, s);
     }
 
-  /* Un-skip the star formation tasks involved with this cell. */
-  for (struct link *l = c->sinks.swallow; l != NULL; l = l->next) {
+  /* Un-skip the density tasks involved with this cell. */
+  for (struct link *l = c->sinks.density; l != NULL; l = l->next) {
     struct task *t = l->t;
     struct cell *ci = t->ci;
     struct cell *cj = t->cj;
+
 #ifdef WITH_MPI
     const int ci_nodeID = ci->nodeID;
     const int cj_nodeID = (cj != NULL) ? cj->nodeID : -1;
@@ -2871,97 +2837,49 @@ int cell_unskip_sinks_tasks(struct cell *c, struct scheduler *s) {
 
     const int ci_active =
         cell_is_active_sinks(ci, e) || cell_is_active_hydro(ci, e);
-
     const int cj_active = (cj != NULL) && (cell_is_active_sinks(cj, e) ||
                                            cell_is_active_hydro(cj, e));
 
-    /* Activate the drifts */
-    if (t->type == task_type_self && ci_active) {
-      cell_activate_drift_sink(ci, s);
-      cell_activate_drift_part(ci, s);
-      cell_activate_sink_formation_tasks(ci->top, s);
-      if (with_timestep_sync) cell_activate_sync_part(ci, s);
-    }
-
     /* Only activate tasks that involve a local active cell. */
     if ((ci_active || cj_active) &&
         (ci_nodeID == nodeID || cj_nodeID == nodeID)) {
-      scheduler_activate(s, t);
-
-      if (t->type == task_type_pair) {
-        /* For the mergers */
-        if (cj_nodeID == nodeID) {
-          cell_activate_drift_sink(cj, s);
-          cell_activate_sink_formation_tasks(cj->top, s);
-        }
-        if (ci_nodeID == nodeID) {
-          cell_activate_drift_sink(ci, s);
-          if (ci->top != cj->top) {
-            cell_activate_sink_formation_tasks(ci->top, s);
-          }
-        }
-
-        /* Do ci */
-        if (ci_active) {
-          /* hydro for cj */
-          atomic_or(&cj->hydro.requires_sorts, 1 << t->flags);
-          cj->hydro.dx_max_sort_old = cj->hydro.dx_max_sort;
-
-          /* Activate the drift tasks. */
-          if (cj_nodeID == nodeID) cell_activate_drift_part(cj, s);
-          if (cj_nodeID == nodeID && with_timestep_sync)
-            cell_activate_sync_part(cj, s);
-
-          /* Check the sorts and activate them if needed. */
-          cell_activate_hydro_sorts(cj, t->flags, s);
-        }
 
-        /* Do cj */
-        if (cj_active) {
-          /* hydro for ci */
-          atomic_or(&ci->hydro.requires_sorts, 1 << t->flags);
-          ci->hydro.dx_max_sort_old = ci->hydro.dx_max_sort;
-
-          /* Activate the drift tasks. */
-          if (ci_nodeID == nodeID) cell_activate_drift_part(ci, s);
-          if (ci_nodeID == nodeID && with_timestep_sync)
-            cell_activate_sync_part(ci, s);
-
-          /* Check the sorts and activate them if needed. */
-          cell_activate_hydro_sorts(ci, t->flags, s);
-        }
-      }
+      scheduler_activate(s, t);
 
-      else if (t->type == task_type_sub_self) {
+      /* Store current values of dx_max and h_max. */
+      if (t->type == task_type_self) {
         cell_activate_subcell_sinks_tasks(ci, NULL, s, with_timestep_sync);
       }
 
-      else if (t->type == task_type_sub_pair) {
+      /* Store current values of dx_max and h_max. */
+      else if (t->type == task_type_pair) {
         cell_activate_subcell_sinks_tasks(ci, cj, s, with_timestep_sync);
       }
     }
 
     /* Only interested in pair interactions as of here. */
-    if (t->type == task_type_pair || t->type == task_type_sub_pair) {
-      /* Check whether there was too much particle motion, i.e. the
-         cell neighbour conditions were violated. */
-      if (cell_need_rebuild_for_sinks_pair(ci, cj)) rebuild = 1;
-      if (cell_need_rebuild_for_sinks_pair(cj, ci)) rebuild = 1;
+    if (t->type == task_type_pair) {
 
-      /* Activate all sink_in tasks for each cell involved
-       * in pair/sub_pair type tasks */
+      /* Activate sink_in for each cell that is part of
+       * a pair task as to not miss any dependencies */
       if (ci_nodeID == nodeID)
         scheduler_activate(s, ci->hydro.super->sinks.sink_in);
       if (cj_nodeID == nodeID)
         scheduler_activate(s, cj->hydro.super->sinks.sink_in);
 
-#ifdef WITH_MPI
+      /* Check whether there was too much particle motion, i.e. the
+         cell neighbour conditions were violated. */
+      if (cell_need_rebuild_for_sinks_pair(ci, cj)) rebuild = 1;
+      if (cell_need_rebuild_for_sinks_pair(cj, ci)) rebuild = 1;
+
+#if defined(WITH_MPI) && !defined(SWIFT_DEBUG_CHECKS)
       error("TODO");
 #endif
     }
   }
 
-  for (struct link *l = c->sinks.do_sink_swallow; l != NULL; l = l->next) {
+  /* Un-skip the swallow tasks involved with this cell. */
+  for (struct link *l = c->sinks.swallow; l != NULL; l = l->next) {
     struct task *t = l->t;
     struct cell *ci = t->ci;
     struct cell *cj = t->cj;
@@ -2975,7 +2893,32 @@ int cell_unskip_sinks_tasks(struct cell *c, struct scheduler *s) {
 
     const int ci_active =
         cell_is_active_sinks(ci, e) || cell_is_active_hydro(ci, e);
+    const int cj_active = (cj != NULL) && (cell_is_active_sinks(cj, e) ||
+                                           cell_is_active_hydro(cj, e));
+
+    /* Only activate tasks that involve a local active cell. */
+    if ((ci_active || cj_active) &&
+        (ci_nodeID == nodeID || cj_nodeID == nodeID)) {
 
+      scheduler_activate(s, t);
+    }
+  }
+
+  /* Un-skip the do_sink_swallow tasks involved with this cell. */
+  for (struct link *l = c->sinks.do_sink_swallow; l != NULL; l = l->next) {
+    struct task *t = l->t;
+    struct cell *ci = t->ci;
+    struct cell *cj = t->cj;
+#ifdef WITH_MPI
+    const int ci_nodeID = ci->nodeID;
+    const int cj_nodeID = (cj != NULL) ? cj->nodeID : -1;
+#else
+    const int ci_nodeID = nodeID;
+    const int cj_nodeID = nodeID;
+#endif
+
+    const int ci_active =
+        cell_is_active_sinks(ci, e) || cell_is_active_hydro(ci, e);
     const int cj_active = (cj != NULL) && (cell_is_active_sinks(cj, e) ||
                                            cell_is_active_hydro(cj, e));
 
@@ -2986,6 +2929,7 @@ int cell_unskip_sinks_tasks(struct cell *c, struct scheduler *s) {
     }
   }
 
+  /* Un-skip the do_gas_swallow tasks involved with this cell. */
   for (struct link *l = c->sinks.do_gas_swallow; l != NULL; l = l->next) {
     struct task *t = l->t;
     struct cell *ci = t->ci;
@@ -3000,7 +2944,6 @@ int cell_unskip_sinks_tasks(struct cell *c, struct scheduler *s) {
 
     const int ci_active =
         cell_is_active_sinks(ci, e) || cell_is_active_hydro(ci, e);
-
     const int cj_active = (cj != NULL) && (cell_is_active_sinks(cj, e) ||
                                            cell_is_active_hydro(cj, e));
 
@@ -3009,9 +2952,9 @@ int cell_unskip_sinks_tasks(struct cell *c, struct scheduler *s) {
         (ci_nodeID == nodeID || cj_nodeID == nodeID)) {
       scheduler_activate(s, t);
 
-      if (t->type == task_type_pair || t->type == task_type_sub_pair) {
+      if (t->type == task_type_pair) {
         /* Activate sinks_out for each cell that is part of
-         * a pair/sub_pair task as to not miss any dependencies */
+         * a pair/pair task as to not miss any dependencies */
         if (ci_nodeID == nodeID)
           scheduler_activate(s, ci->hydro.super->sinks.sink_out);
         if (cj_nodeID == nodeID)
@@ -3025,11 +2968,21 @@ int cell_unskip_sinks_tasks(struct cell *c, struct scheduler *s) {
       (cell_is_active_sinks(c, e) || cell_is_active_hydro(c, e))) {
 
     if (c->sinks.sink_in != NULL) scheduler_activate(s, c->sinks.sink_in);
+    if (c->top->sinks.sink_formation != NULL) {
+      cell_activate_sink_formation_tasks(c->top, s);
+      cell_activate_super_sink_drifts(c->top, s);
+    }
+    if (c->sinks.density_ghost != NULL)
+      scheduler_activate(s, c->sinks.density_ghost);
     if (c->sinks.sink_ghost1 != NULL)
       scheduler_activate(s, c->sinks.sink_ghost1);
     if (c->sinks.sink_ghost2 != NULL)
       scheduler_activate(s, c->sinks.sink_ghost2);
     if (c->sinks.sink_out != NULL) scheduler_activate(s, c->sinks.sink_out);
+    if (c->top->sinks.star_formation_sink != NULL) {
+      cell_activate_star_formation_sink_tasks(c->top, s, with_feedback);
+      cell_activate_super_sink_drifts(c->top, s);
+    }
     if (c->kick1 != NULL) scheduler_activate(s, c->kick1);
     if (c->kick2 != NULL) scheduler_activate(s, c->kick2);
     if (c->timestep != NULL) scheduler_activate(s, c->timestep);
@@ -3039,7 +2992,6 @@ int cell_unskip_sinks_tasks(struct cell *c, struct scheduler *s) {
     if (c->csds != NULL) scheduler_activate(s, c->csds);
 #endif
   }
-
   return rebuild;
 }
 
@@ -3088,32 +3040,21 @@ int cell_unskip_rt_tasks(struct cell *c, struct scheduler *s,
       scheduler_activate(s, t);
 
       if (!sub_cycle) {
-        /* Activate sorts only during main/normal steps. */
-        if (t->type == task_type_pair) {
-          atomic_or(&ci->hydro.requires_sorts, 1 << t->flags);
-          atomic_or(&cj->hydro.requires_sorts, 1 << t->flags);
-          ci->hydro.dx_max_sort_old = ci->hydro.dx_max_sort;
-          cj->hydro.dx_max_sort_old = cj->hydro.dx_max_sort;
-
-          /* Check the sorts and activate them if needed. */
-          cell_activate_rt_sorts(ci, t->flags, s);
-          cell_activate_rt_sorts(cj, t->flags, s);
-        }
 
         /* Store current values of dx_max and h_max. */
-        else if (t->type == task_type_sub_self) {
+        if (t->type == task_type_self) {
           cell_activate_subcell_rt_tasks(ci, NULL, s, sub_cycle);
         }
 
         /* Store current values of dx_max and h_max. */
-        else if (t->type == task_type_sub_pair) {
+        else if (t->type == task_type_pair) {
           cell_activate_subcell_rt_tasks(ci, cj, s, sub_cycle);
         }
       }
     }
 
     /* Only interested in pair interactions as of here. */
-    if (t->type == task_type_pair || t->type == task_type_sub_pair) {
+    if (t->type == task_type_pair) {
 
 #ifdef WITH_MPI
 
@@ -3137,7 +3078,7 @@ int cell_unskip_rt_tasks(struct cell *c, struct scheduler *s,
             scheduler_activate_recv(s, ci->mpi.recv, task_subtype_rt_transport);
           }
         } else if (ci_active) {
-#ifdef MPI_SYMMETRIC_FORCE_INTERACTION
+#ifdef MPI_SYMMETRIC_FORCE_INTERACTION_RT
           /* If the local cell is inactive and the remote cell is active, we
            * still need to receive stuff to be able to do the force interaction
            * on this node as well.
@@ -3161,7 +3102,7 @@ int cell_unskip_rt_tasks(struct cell *c, struct scheduler *s,
                                     ci_nodeID);
           }
         } else if (cj_active) {
-#ifdef MPI_SYMMETRIC_FORCE_INTERACTION
+#ifdef MPI_SYMMETRIC_FORCE_INTERACTION_RT
           /* If the foreign cell is inactive, but the local cell is active,
            * we still need to send stuff to be able to do the force interaction
            * on both nodes.
@@ -3190,7 +3131,7 @@ int cell_unskip_rt_tasks(struct cell *c, struct scheduler *s,
             scheduler_activate_recv(s, cj->mpi.recv, task_subtype_rt_transport);
           }
         } else if (cj_active) {
-#ifdef MPI_SYMMETRIC_FORCE_INTERACTION
+#ifdef MPI_SYMMETRIC_FORCE_INTERACTION_RT
           /* If the local cell is inactive and the remote cell is active, we
            * still need to receive stuff to be able to do the force interaction
            * on this node as well.
@@ -3214,7 +3155,7 @@ int cell_unskip_rt_tasks(struct cell *c, struct scheduler *s,
                                     cj_nodeID);
           }
         } else if (ci_active) {
-#ifdef MPI_SYMMETRIC_FORCE_INTERACTION
+#ifdef MPI_SYMMETRIC_FORCE_INTERACTION_RT
           /* If the foreign cell is inactive, but the local cell is active,
            * we still need to send stuff to be able to do the force interaction
            * on both nodes
@@ -3252,10 +3193,10 @@ int cell_unskip_rt_tasks(struct cell *c, struct scheduler *s,
         (cj_active && cj_nodeID == nodeID)) {
       scheduler_activate(s, t);
 
-      if (t->type == task_type_pair || t->type == task_type_sub_pair) {
+      if (t->type == task_type_pair) {
 
         /* Activate transport_out for each cell that is part of
-         * a pair/sub_pair task as to not miss any dependencies */
+         * a pair task as to not miss any dependencies */
         if (ci_nodeID == nodeID)
           scheduler_activate(s, ci->hydro.super->rt.rt_transport_out);
         if (cj_nodeID == nodeID)
@@ -3275,14 +3216,13 @@ int cell_unskip_rt_tasks(struct cell *c, struct scheduler *s,
       if (c->rt.rt_tchem != NULL) scheduler_activate(s, c->rt.rt_tchem);
       if (c->rt.rt_out != NULL) scheduler_activate(s, c->rt.rt_out);
     } else {
-#if defined(MPI_SYMMETRIC_FORCE_INTERACTION) && defined(WITH_MPI)
+#if defined(MPI_SYMMETRIC_FORCE_INTERACTION_RT) && defined(WITH_MPI)
       /* Additionally unskip force interactions between inactive local cell and
        * active remote cell. (The cell unskip will only be called for active
        * cells, so, we have to do this now, from the active remote cell). */
       for (struct link *l = c->rt.rt_transport; l != NULL; l = l->next) {
         struct task *t = l->t;
-        if (t->type != task_type_pair && t->type != task_type_sub_pair)
-          continue;
+        if (t->type != task_type_pair) continue;
 
         struct cell *ci = l->t->ci;
         struct cell *cj = l->t->cj;
@@ -3298,20 +3238,9 @@ int cell_unskip_rt_tasks(struct cell *c, struct scheduler *s,
           scheduler_activate(s, l->t);
 
           if (!sub_cycle) {
-            /* Activate sorts only during main/normal steps. */
-            if (t->type == task_type_pair) {
-              atomic_or(&ci->hydro.requires_sorts, 1 << t->flags);
-              atomic_or(&cj->hydro.requires_sorts, 1 << t->flags);
-              ci->hydro.dx_max_sort_old = ci->hydro.dx_max_sort;
-              cj->hydro.dx_max_sort_old = cj->hydro.dx_max_sort;
-
-              /* Check the sorts and activate them if needed. */
-              cell_activate_rt_sorts(ci, t->flags, s);
-              cell_activate_rt_sorts(cj, t->flags, s);
-            }
 
             /* Store current values of dx_max and h_max. */
-            else if (t->type == task_type_sub_pair) {
+            if (t->type == task_type_pair) {
               cell_activate_subcell_rt_tasks(ci, cj, s, sub_cycle);
             }
           }

src/chemistry/AGORA/chemistry.h

@@ -294,10 +294,28 @@ __attribute__((always_inline)) INLINE static void chemistry_first_init_spart(
     const struct chemistry_global_data* data, struct spart* restrict sp) {
 
   for (int i = 0; i < AGORA_CHEMISTRY_ELEMENT_COUNT; i++) {
-    sp->chemistry_data.metal_mass_fraction[i] = data->initial_metallicities[i];
+    /* Bug fix (26.07.2024): Check that the initial me metallicities are non
+       negative. */
+    if (data->initial_metallicities[i] >= 0) {
+      /* Use the value from the parameter file */
+      sp->chemistry_data.metal_mass_fraction[i] =
+          data->initial_metallicities[i];
+    }
+    /* else : Use the value from the IC. We are reading the metal mass
+       fraction. So do not overwrite the metallicities */
   }
 }
 
+/**
+ * @brief Sets the chemistry properties of the sink particles to a valid start
+ * state.
+ *
+ * @param data The global chemistry information.
+ * @param sink Pointer to the sink particle data.
+ */
+__attribute__((always_inline)) INLINE static void chemistry_first_init_sink(
+    const struct chemistry_global_data* data, struct sink* restrict sink) {}
+
 /**
  * @brief Initialise the chemistry properties of a black hole with
  * the chemistry properties of the gas it is born from.

src/chemistry/AGORA/chemistry_iact.h

@@ -103,7 +103,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_chemistry(
 }
 
 /**
- * @brief do metal diffusion computation in the <FORCE LOOP>
+ * @brief do metal diffusion computation in the force loop
  * (symmetric version)
  *
  * @param r2 Comoving square distance between the two particles.
@@ -122,13 +122,13 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_chemistry(
  *
  */
 __attribute__((always_inline)) INLINE static void runner_iact_diffusion(
-    float r2, const float *dx, float hi, float hj, struct part *restrict pi,
-    struct part *restrict pj, const float a, const float H,
-    const float time_base, const integertime_t t_current,
+    const float r2, const float dx[3], const float hi, const float hj,
+    struct part *restrict pi, struct part *restrict pj, const float a,
+    const float H, const float time_base, const integertime_t t_current,
     const struct cosmology *cosmo, const int with_cosmology) {}
 
 /**
- * @brief do metal diffusion computation in the <FORCE LOOP>
+ * @brief do metal diffusion computation in the force loop
  * (nonsymmetric version)
  *
  * @param r2 Comoving square distance between the two particles.
@@ -147,9 +147,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_diffusion(
  *
  */
 __attribute__((always_inline)) INLINE static void runner_iact_nonsym_diffusion(
-    float r2, const float *dx, float hi, float hj, struct part *restrict pi,
-    struct part *restrict pj, const float a, const float H,
-    const float time_base, const integertime_t t_current,
+    const float r2, const float dx[3], const float hi, const float hj,
+    struct part *restrict pi, struct part *restrict pj, const float a,
+    const float H, const float time_base, const integertime_t t_current,
     const struct cosmology *cosmo, const int with_cosmology) {}
 
 #endif /* SWIFT_AGORA_CHEMISTRY_IACT_H */

src/chemistry/AGORA/chemistry_io.h

@@ -99,6 +99,19 @@ INLINE static int chemistry_write_sparticles(const struct spart* sparts,
   return 1;
 }
 
+/**
+ * @brief Specifies which sink fields to write to a dataset
+ *
+ * @param sinks The #sink array.
+ * @param list The list of i/o properties to write.
+ *
+ * @return Returns the number of fields to write.
+ */
+INLINE static int chemistry_write_sinkparticles(const struct sink* sinks,
+                                                struct io_props* list) {
+  return 0;
+}
+
 /**
  * @brief Specifies which bparticle fields to write to a dataset
  *

src/chemistry/EAGLE/chemistry.h

@@ -209,6 +209,16 @@ __attribute__((always_inline)) INLINE static void chemistry_first_init_spart(
   }
 }
 
+/**
+ * @brief Sets the chemistry properties of the sink particles to a valid start
+ * state.
+ *
+ * @param data The global chemistry information.
+ * @param sink Pointer to the sink particle data.
+ */
+__attribute__((always_inline)) INLINE static void chemistry_first_init_sink(
+    const struct chemistry_global_data* data, struct sink* restrict sink) {}
+
 /**
  * @brief Initialises the chemistry properties.
  *

src/chemistry/EAGLE/chemistry_iact.h

@@ -133,7 +133,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_chemistry(
 }
 
 /**
- * @brief do metal diffusion computation in the <FORCE LOOP>
+ * @brief do metal diffusion computation in the force loop
  * (symmetric version)
  *
  * @param r2 Comoving square distance between the two particles.
@@ -158,7 +158,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_diffusion(
     const struct cosmology *cosmo, const int with_cosmology) {}
 
 /**
- * @brief do metal diffusion computation in the <FORCE LOOP>
+ * @brief do metal diffusion computation in the force loop
  * (nonsymmetric version)
  *
  * @param r2 Comoving square distance between the two particles.

src/chemistry/EAGLE/chemistry_io.h

@@ -214,6 +214,19 @@ INLINE static int chemistry_write_sparticles(const struct spart* sparts,
   return 12;
 }
 
+/**
+ * @brief Specifies which sink fields to write to a dataset
+ *
+ * @param sinks The #sink array.
+ * @param list The list of i/o properties to write.
+ *
+ * @return Returns the number of fields to write.
+ */
+INLINE static int chemistry_write_sinkparticles(const struct sink* sinks,
+                                                struct io_props* list) {
+  return 0;
+}
+
 /**
  * @brief Specifies which black hole particle fields to write to a dataset
  *

src/chemistry/GEAR/chemistry.h

@@ -48,15 +48,52 @@
 INLINE static void chemistry_copy_star_formation_properties(
     struct part* p, const struct xpart* xp, struct spart* sp) {
 
+  /* gas mass after update */
   float mass = hydro_get_mass(p);
 
   /* Store the chemistry struct in the star particle */
-  for (int i = 0; i < GEAR_CHEMISTRY_ELEMENT_COUNT; i++) {
-    sp->chemistry_data.metal_mass_fraction[i] =
-        p->chemistry_data.smoothed_metal_mass_fraction[i];
+  for (int k = 0; k < GEAR_CHEMISTRY_ELEMENT_COUNT; k++) {
+    sp->chemistry_data.metal_mass_fraction[k] =
+        p->chemistry_data.smoothed_metal_mass_fraction[k];
 
     /* Remove the metals taken by the star. */
-    p->chemistry_data.metal_mass[i] *= mass / (mass + sp->mass);
+    p->chemistry_data.metal_mass[k] *= mass / (mass + sp->mass);
+  }
+}
+
+/**
+ * @brief Copies the chemistry properties of the sink particle over to the
+ * stellar particle.
+ *
+ * @param sink the sink particle with its properties.
+ * @param sp the new star particles.
+ */
+INLINE static void chemistry_copy_sink_properties_to_star(struct sink* sink,
+                                                          struct spart* sp) {
+
+  /* Store the chemistry struct in the star particle */
+  for (int k = 0; k < GEAR_CHEMISTRY_ELEMENT_COUNT; k++) {
+    sp->chemistry_data.metal_mass_fraction[k] =
+        sink->chemistry_data.metal_mass_fraction[k];
+  }
+}
+
+/**
+ * @brief Copies the chemistry properties of the gas particle over to the
+ * sink particle.
+ *
+ * @param p the gas particles.
+ * @param xp the additional properties of the gas particles.
+ * @param sink the new created star particle with its properties.
+ */
+INLINE static void chemistry_copy_sink_properties(const struct part* p,
+                                                  const struct xpart* xp,
+                                                  struct sink* sink) {
+
+  /* Store the chemistry struct in the star particle */
+  for (int i = 0; i < GEAR_CHEMISTRY_ELEMENT_COUNT; i++) {
+    sink->chemistry_data.metal_mass_fraction[i] =
+        p->chemistry_data.smoothed_metal_mass_fraction[i];
   }
 }
 
@@ -68,8 +105,9 @@ INLINE static void chemistry_copy_star_formation_properties(
  */
 static INLINE void chemistry_print_backend(
     const struct chemistry_global_data* data) {
-
-  message("Chemistry function is 'Gear'.");
+  if (engine_rank == 0) {
+    message("Chemistry function is 'Gear'.");
+  }
 }
 
 /**
@@ -246,10 +284,12 @@ static INLINE void chemistry_init_backend(struct swift_params* parameter_file,
   const float initial_metallicity = parser_get_param_float(
       parameter_file, "GEARChemistry:initial_metallicity");
 
-  if (initial_metallicity < 0) {
-    message("Setting the initial metallicity from the snapshot.");
-  } else {
-    message("Setting the initial metallicity from the parameter file.");
+  if (engine_rank == 0) {
+    if (initial_metallicity < 0) {
+      message("Setting the initial metallicity from the snapshot.");
+    } else {
+      message("Setting the initial metallicity from the parameter file.");
+    }
   }
 
   /* Set the initial metallicities */
@@ -426,7 +466,37 @@ __attribute__((always_inline)) INLINE static void chemistry_first_init_spart(
     const struct chemistry_global_data* data, struct spart* restrict sp) {
 
   for (int i = 0; i < GEAR_CHEMISTRY_ELEMENT_COUNT; i++) {
-    sp->chemistry_data.metal_mass_fraction[i] = data->initial_metallicities[i];
+    /* Bug fix (26.07.2024): Check that the initial me metallicities are non
+       negative. */
+    if (data->initial_metallicities[i] >= 0) {
+      /* Use the value from the parameter file */
+      sp->chemistry_data.metal_mass_fraction[i] =
+          data->initial_metallicities[i];
+    }
+    /* else : Use the value from the IC. We are reading the metal mass
+     fraction. So do not overwrite the metallicities */
+  }
+}
+
+/* Add chemistry first init sink ? */
+
+/**
+ * @brief Sets the chemistry properties of the sink particles to a valid start
+ * state.
+ *
+ * @param data The global chemistry information.
+ * @param sink Pointer to the sink particle data.
+ */
+__attribute__((always_inline)) INLINE static void chemistry_first_init_sink(
+    const struct chemistry_global_data* data, struct sink* restrict sink) {
+
+  for (int i = 0; i < GEAR_CHEMISTRY_ELEMENT_COUNT; i++) {
+    /* Use the value from the parameter file */
+    if (data->initial_metallicities[i] >= 0) {
+      sink->chemistry_data.metal_mass_fraction[i] =
+          data->initial_metallicities[i];
+    }
+    /* else : read the metallicities from the ICs. */
   }
 }
 
@@ -435,13 +505,17 @@ __attribute__((always_inline)) INLINE static void chemistry_first_init_spart(
  *
  * @param si_data The black hole data to add to.
  * @param sj_data The gas data to use.
- * @param gas_mass The mass of the gas particle.
+ * @param mi_old The mass of the #sink i before accreting the #part p.
  */
 __attribute__((always_inline)) INLINE static void chemistry_add_sink_to_sink(
-    struct chemistry_sink_data* si_data,
-    const struct chemistry_sink_data* sj_data) {
+    struct sink* si, const struct sink* sj, const double mi_old) {
+
+  for (int k = 0; k < GEAR_CHEMISTRY_ELEMENT_COUNT; k++) {
+    double mk = si->chemistry_data.metal_mass_fraction[k] * mi_old +
+                sj->chemistry_data.metal_mass_fraction[k] * sj->mass;
 
-  // To be implemented.
+    si->chemistry_data.metal_mass_fraction[k] = mk / si->mass;
+  }
 }
 
 /**
@@ -449,13 +523,20 @@ __attribute__((always_inline)) INLINE static void chemistry_add_sink_to_sink(
  *
  * @param sp_data The sink data to add to.
  * @param p_data The gas data to use.
- * @param gas_mass The mass of the gas particle.
+ * @param ms_old The mass of the #sink before accreting the #part p.
  */
 __attribute__((always_inline)) INLINE static void chemistry_add_part_to_sink(
-    struct chemistry_sink_data* sp_data,
-    const struct chemistry_part_data* p_data, const double gas_mass) {
+    struct sink* s, const struct part* p, const double ms_old) {
+
+  /* gas mass */
+  const float mass = hydro_get_mass(p);
 
-  // To be implemented.
+  for (int k = 0; k < GEAR_CHEMISTRY_ELEMENT_COUNT; k++) {
+    double mk = s->chemistry_data.metal_mass_fraction[k] * ms_old +
+                p->chemistry_data.smoothed_metal_mass_fraction[k] * mass;
+
+    s->chemistry_data.metal_mass_fraction[k] = mk / s->mass;
+  }
 }
 
 /**
@@ -564,6 +645,20 @@ chemistry_get_star_total_iron_mass_fraction_for_feedback(
   return sp->chemistry_data.metal_mass_fraction[0];
 }
 
+/**
+ * @brief Returns the total iron mass fraction of the
+ * sink particle to be used in feedback/enrichment related routines.
+ * We assume iron to be stored at index 0.
+ *
+ * @param sp Pointer to the particle data.
+ */
+__attribute__((always_inline)) INLINE static double
+chemistry_get_sink_total_iron_mass_fraction_for_feedback(
+    const struct sink* restrict sink) {
+
+  return sink->chemistry_data.metal_mass_fraction[0];
+}
+
 /**
  * @brief Returns the abundances (metal mass fraction) of the
  * star particle to be used in feedback/enrichment related routines.

src/chemistry/GEAR/chemistry_iact.h

@@ -108,7 +108,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_chemistry(
 }
 
 /**
- * @brief do metal diffusion computation in the <FORCE LOOP>
+ * @brief do metal diffusion computation in the force loop
  * (symmetric version)
  *
  * @param r2 Comoving square distance between the two particles.
@@ -133,7 +133,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_diffusion(
     const struct cosmology *cosmo, const int with_cosmology) {}
 
 /**
- * @brief do metal diffusion computation in the <FORCE LOOP>
+ * @brief do metal diffusion computation in the force loop
  * (nonsymmetric version)
  *
  * @param r2 Comoving square distance between the two particles.

src/chemistry/GEAR/chemistry_io.h

@@ -107,6 +107,26 @@ INLINE static int chemistry_write_sparticles(const struct spart* sparts,
   return 1;
 }
 
+/**
+ * @brief Specifies which sink fields to write to a dataset
+ *
+ * @param sinks The #sink array.
+ * @param list The list of i/o properties to write.
+ *
+ * @return Returns the number of fields to write.
+ */
+INLINE static int chemistry_write_sinkparticles(const struct sink* sinks,
+                                                struct io_props* list) {
+
+  /* List what we want to write */
+  list[0] = io_make_output_field(
+      "MetalMassFractions", DOUBLE, GEAR_CHEMISTRY_ELEMENT_COUNT,
+      UNIT_CONV_NO_UNITS, 0.f, sinks, chemistry_data.metal_mass_fraction,
+      "Mass fraction of each element");
+
+  return 1;
+}
+
 /**
  * @brief Specifies which black hole particle fields to write to a dataset
  *