src/cooling/grackle/cooling.h

@@ -105,6 +105,16 @@ INLINE static float cooling_get_subgrid_density(const struct part* p,
   return -1.f;
 }
 
+/**
+ * @brief Record the time when cooling was switched off for a particle.
+ *
+ * @param p #part data.
+ * @param xp Pointer to the #xpart data.
+ * @param time The time when the cooling was switched off.
+ */
+void cooling_set_part_time_cooling_off(struct part* p, struct xpart* xp,
+                                       const double time);
+
 float cooling_get_radiated_energy(const struct xpart* restrict xp);
 void cooling_print_backend(const struct cooling_function_data* cooling);
 

src/cooling/grackle/cooling_io.h

@@ -25,6 +25,7 @@
 #include "cooling_struct.h"
 #include "io_properties.h"
 #include "physical_constants.h"
+#include "units.h"
 
 #ifdef HAVE_HDF5
 
@@ -140,7 +141,7 @@ __attribute__((always_inline)) INLINE static int cooling_write_particles(
  */
 __attribute__((always_inline)) INLINE static void cooling_read_parameters(
     struct swift_params* parameter_file, struct cooling_function_data* cooling,
-    const struct phys_const* phys_const) {
+    const struct phys_const* phys_const, const struct unit_system* us) {
 
   parser_get_param_string(parameter_file, "GrackleCooling:cloudy_table",
                           cooling->cloudy_table);
@@ -207,6 +208,36 @@ __attribute__((always_inline)) INLINE static void cooling_read_parameters(
   cooling->HydrogenFractionByMass = parser_get_opt_param_double(
       parameter_file, "GrackleCooling:HydrogenFractionByMass", 0.76);
 
+  cooling->initial_nHII_to_nH_ratio = parser_get_opt_param_double(
+      parameter_file, "GrackleCooling:initial_nHII_to_nH_ratio", -1);
+
+  cooling->initial_nHeI_to_nH_ratio = parser_get_opt_param_double(
+      parameter_file, "GrackleCooling:initial_nHeI_to_nH_ratio", -1);
+
+  cooling->initial_nHeII_to_nH_ratio = parser_get_opt_param_double(
+      parameter_file, "GrackleCooling:initial_nHeII_to_nH_ratio", -1);
+
+  cooling->initial_nHeIII_to_nH_ratio = parser_get_opt_param_double(
+      parameter_file, "GrackleCooling:initial_nHeIII_to_nH_ratio", -1);
+
+  cooling->initial_nDI_to_nH_ratio = parser_get_opt_param_double(
+      parameter_file, "GrackleCooling:initial_nDI_to_nH_ratio", -1);
+
+  cooling->initial_nDII_to_nH_ratio = parser_get_opt_param_double(
+      parameter_file, "GrackleCooling:initial_nDII_to_nH_ratio", -1);
+
+  cooling->initial_nHM_to_nH_ratio = parser_get_opt_param_double(
+      parameter_file, "GrackleCooling:initial_nHM_to_nH_ratio", -1);
+
+  cooling->initial_nH2I_to_nH_ratio = parser_get_opt_param_double(
+      parameter_file, "GrackleCooling:initial_nH2I_to_nH_ratio", -1);
+
+  cooling->initial_nH2II_to_nH_ratio = parser_get_opt_param_double(
+      parameter_file, "GrackleCooling:initial_nH2II_to_nH_ratio", -1);
+
+  cooling->initial_nHDI_to_nH_ratio = parser_get_opt_param_double(
+      parameter_file, "GrackleCooling:initial_nHDI_to_nH_ratio", -1);
+
   /* Self shielding */
   cooling->self_shielding_method = parser_get_opt_param_int(
       parameter_file, "GrackleCooling:self_shielding_method", 0);
@@ -216,9 +247,6 @@ __attribute__((always_inline)) INLINE static void cooling_read_parameters(
         parameter_file, "GrackleCooling:self_shielding_threshold_atom_per_cm3");
   }
 
-  cooling->HydrogenFractionByMass = parser_get_opt_param_double(
-      parameter_file, "GrackleCooling:HydrogenFractionByMass", 0.76);
-
   /* Initial step convergence */
   cooling->max_step = parser_get_opt_param_int(
       parameter_file, "GrackleCooling:max_steps", 10000);
@@ -230,6 +258,16 @@ __attribute__((always_inline)) INLINE static void cooling_read_parameters(
   cooling->thermal_time = parser_get_param_double(
       parameter_file, "GrackleCooling:thermal_time_myr");
   cooling->thermal_time *= phys_const->const_year * 1e6;
+
+  /* Maximal cooling density (in acc) */
+  cooling->cooling_density_max = parser_get_param_double(
+      parameter_file, "GrackleCooling:maximal_density_Hpcm3");
+
+  /* Convert from acc to intenal units */
+  const double m_p_cgs = phys_const->const_proton_mass *
+                         units_cgs_conversion_factor(us, UNIT_CONV_MASS);
+  cooling->cooling_density_max *=
+      m_p_cgs / units_cgs_conversion_factor(us, UNIT_CONV_DENSITY);
 }
 
 #endif /* SWIFT_COOLING_GRACKLE_IO_H */

src/cooling/grackle/cooling_properties.h

@@ -103,6 +103,36 @@ struct cooling_function_data {
   /*! Hydrogen fraction by mass */
   float HydrogenFractionByMass;
 
+  /*! initial nHII to nH ratio (number density ratio) */
+  float initial_nHII_to_nH_ratio;
+
+  /*! initial nHeI to nH ratio (number density ratio) */
+  float initial_nHeI_to_nH_ratio;
+
+  /*! initial nHeII to nH ratio (number density ratio) */
+  float initial_nHeII_to_nH_ratio;
+
+  /*! initial nHeIII to nH ratio (number density ratio) */
+  float initial_nHeIII_to_nH_ratio;
+
+  /*! initial nDI to nH ratio (number density ratio) */
+  float initial_nDI_to_nH_ratio;
+
+  /*! initial nDII to nH ratio (number density ratio) */
+  float initial_nDII_to_nH_ratio;
+
+  /*! initial nHM to nH ratio (number density ratio) */
+  float initial_nHM_to_nH_ratio;
+
+  /*! initial nH2I to nH ratio (number density ratio) */
+  float initial_nH2I_to_nH_ratio;
+
+  /*! initial nH2II to nH ratio (number density ratio) */
+  float initial_nH2II_to_nH_ratio;
+
+  /*! initial nHDI to nH ratio (number density ratio) */
+  float initial_nHDI_to_nH_ratio;
+
   /*! Self shielding method (1 -> 3 for grackle's ones, 0 for none and -1 for
    * GEAR) */
   int self_shielding_method;
@@ -121,6 +151,9 @@ struct cooling_function_data {
 
   /*! Duration for switching off cooling after an event (e.g. supernovae) */
   double thermal_time;
+
+  /*! Maximal allowed density for cooling (in internal units). */
+  double cooling_density_max;
 };
 
 #endif /* SWIFT_COOLING_PROPERTIES_GRACKLE_H */

src/cooling/none/cooling.h

@@ -277,6 +277,19 @@ __attribute__((always_inline)) INLINE static float cooling_get_radiated_energy(
   return 0.f;
 }
 
+/**
+ * @brief Record the time when cooling was switched off for a particle.
+ *
+ * In none model, this function does nothing.
+ *
+ * @param p #part data.
+ * @param xp Pointer to the #xpart data.
+ * @param time The time when the cooling was switched off.
+ */
+INLINE static void cooling_set_part_time_cooling_off(struct part* p,
+                                                     struct xpart* xp,
+                                                     const double time) {}
+
 /**
  * @brief Split the coolong content of a particle into n pieces
  *

src/cosmology.c

@@ -813,7 +813,9 @@ void cosmology_init_tables(struct cosmology *c) {
 
     /* Choose value of expansion factor for check */
     const double dloga = (c->log_a_end - c->log_a_begin) / (n - 1);
-    const double a = exp(c->log_a_begin + dloga * i);
+    double a = exp(c->log_a_begin + dloga * i);
+    a = fmax(a, c->a_begin);
+    a = fmin(a, c->a_end);
 
     /* Verify that converting expansion factor to time and back recovers the
      * original value */
@@ -830,11 +832,10 @@ void cosmology_init_tables(struct cosmology *c) {
     if (frac_error > max_error_distance) max_error_distance = frac_error;
   }
 
-  message("Max fractional error in a to age of universe round trip = %16.8e\n",
+  message("Max fractional error in a to age of universe round trip = %16.8e",
           max_error_time);
-  message(
-      "Max fractional error in a to comoving distance round trip = %16.8e\n",
-      max_error_distance);
+  message("Max fractional error in a to comoving distance round trip = %16.8e",
+          max_error_distance);
 
 #endif /* SWIFT_DEBUG_CHECKS */
 

src/cosmology.h

@@ -113,7 +113,7 @@ struct cosmology {
   /*! Scale-factor at the previous time-step */
   double a_old;
 
-  /*! Redshit at the previous time-step */
+  /*! Redshift at the previous time-step */
   double z_old;
 
   /*------------------------------------------------------------------ */

src/csds.c

@@ -1232,7 +1232,7 @@ void csds_struct_restore(struct csds_writer *log, FILE *stream) {
 
   /* Restore the pointers */
   for (int i = 0; i < swift_type_count; i++) {
-    if (log->field_pointers == NULL) continue;
+    if (log->field_pointers[i] == NULL) continue;
 
     log->field_pointers[i] =
         log->list_fields + (log->field_pointers[i] - old_list_fields);

src/csds_io.h

@@ -165,12 +165,10 @@ struct csds_field {
  * @param conversion_func The conversion function.
  * @param field_size The size of the field to write.
  */
-#define csds_define_field_from_function_hydro(csds_field, field_name,     \
-                                              conversion_func, size)      \
-  {                                                                       \
-    csds_define_field_from_function_general(csds_field, field_name,       \
-                                            conversion_func, size, hydro) \
-  }
+#define csds_define_field_from_function_hydro(csds_field, field_name, \
+                                              conversion_func, size)  \
+  {csds_define_field_from_function_general(csds_field, field_name,    \
+                                           conversion_func, size, hydro)}
 
 /**
  * @brief Define a field from a function for stars.
@@ -180,12 +178,10 @@ struct csds_field {
  * @param conversion_func The conversion function.
  * @param field_size The size of the field to write.
  */
-#define csds_define_field_from_function_stars(csds_field, field_name,     \
-                                              conversion_func, size)      \
-  {                                                                       \
-    csds_define_field_from_function_general(csds_field, field_name,       \
-                                            conversion_func, size, stars) \
-  }
+#define csds_define_field_from_function_stars(csds_field, field_name, \
+                                              conversion_func, size)  \
+  {csds_define_field_from_function_general(csds_field, field_name,    \
+                                           conversion_func, size, stars)}
 
 /**
  * @brief Define a field from a function for gravity.
@@ -195,12 +191,10 @@ struct csds_field {
  * @param conversion_func The conversion function.
  * @param field_size The size of the field to write.
  */
-#define csds_define_field_from_function_gravity(csds_field, field_name,  \
-                                                conversion_func, size)   \
-  {                                                                      \
-    csds_define_field_from_function_general(csds_field, field_name,      \
-                                            conversion_func, size, grav) \
-  }
+#define csds_define_field_from_function_gravity(csds_field, field_name, \
+                                                conversion_func, size)  \
+  {csds_define_field_from_function_general(csds_field, field_name,      \
+                                           conversion_func, size, grav)}
 
 void csds_write_description(struct csds_writer *log, struct engine *e);
 

src/debug.c

@@ -70,10 +70,12 @@
 #include "./hydro/Phantom/hydro_debug.h"
 #elif defined(GIZMO_MFV_SPH) || defined(GIZMO_MFM_SPH)
 #include "./hydro/Gizmo/hydro_debug.h"
-#elif defined(SHADOWFAX_SPH)
+#elif defined(SHADOWSWIFT)
 #include "./hydro/Shadowswift/hydro_debug.h"
 #elif defined(PLANETARY_SPH)
 #include "./hydro/Planetary/hydro_debug.h"
+#elif defined(REMIX_SPH)
+#include "./hydro/REMIX/hydro_debug.h"
 #elif defined(SPHENIX_SPH)
 #include "./hydro/SPHENIX/hydro_debug.h"
 #elif defined(GASOLINE_SPH)
@@ -242,8 +244,8 @@ int checkSpacehmax(struct space *s) {
   float cell_sinks_h_max = 0.0f;
   for (int k = 0; k < s->nr_cells; k++) {
     if (s->cells_top[k].nodeID == s->e->nodeID &&
-        s->cells_top[k].sinks.r_cut_max > cell_sinks_h_max) {
-      cell_sinks_h_max = s->cells_top[k].sinks.r_cut_max;
+        s->cells_top[k].sinks.h_max > cell_sinks_h_max) {
+      cell_sinks_h_max = s->cells_top[k].sinks.h_max;
     }
   }
 
@@ -266,8 +268,8 @@ int checkSpacehmax(struct space *s) {
   /* Now all the sinks. */
   float sink_h_max = 0.0f;
   for (size_t k = 0; k < s->nr_sinks; k++) {
-    if (s->sinks[k].r_cut > sink_h_max) {
-      sink_h_max = s->sinks[k].r_cut;
+    if (s->sinks[k].h > sink_h_max) {
+      sink_h_max = s->sinks[k].h;
     }
   }
 
@@ -315,17 +317,17 @@ int checkSpacehmax(struct space *s) {
   /* sink */
   for (int k = 0; k < s->nr_cells; k++) {
     if (s->cells_top[k].nodeID == s->e->nodeID) {
-      if (s->cells_top[k].sinks.r_cut_max > sink_h_max) {
+      if (s->cells_top[k].sinks.h_max > sink_h_max) {
         message("cell %d is inconsistent (%f > %f)", k,
-                s->cells_top[k].sinks.r_cut_max, sink_h_max);
+                s->cells_top[k].sinks.h_max, sink_h_max);
       }
     }
   }
 
   for (size_t k = 0; k < s->nr_sinks; k++) {
-    if (s->sinks[k].r_cut > cell_sinks_h_max) {
+    if (s->sinks[k].h > cell_sinks_h_max) {
       message("spart %lld is inconsistent (%f > %f)", s->sinks[k].id,
-              s->sinks[k].r_cut, cell_sinks_h_max);
+              s->sinks[k].h, cell_sinks_h_max);
     }
   }
 
@@ -436,7 +438,7 @@ int checkCellhdxmax(const struct cell *c, int *depth) {
                       sp->x_diff[1] * sp->x_diff[1] +
                       sp->x_diff[2] * sp->x_diff[2];
 
-    sinks_h_max = max(sinks_h_max, sp->r_cut);
+    sinks_h_max = max(sinks_h_max, sp->h);
     sinks_dx_max = max(sinks_dx_max, sqrtf(dx2));
   }
 
@@ -476,9 +478,9 @@ int checkCellhdxmax(const struct cell *c, int *depth) {
     result = 0;
   }
 
-  if (c->sinks.r_cut_max != sinks_h_max) {
+  if (c->sinks.h_max != sinks_h_max) {
     message("%d Inconsistent sinks_h_max: cell %f != parts %f", *depth,
-            c->sinks.r_cut_max, sinks_h_max);
+            c->sinks.h_max, sinks_h_max);
     message("location: %f %f %f", c->loc[0], c->loc[1], c->loc[2]);
     result = 0;
   }

src/distributed_io.c

@@ -56,13 +56,17 @@
 #include "sink_io.h"
 #include "star_formation_io.h"
 #include "stars_io.h"
+#include "swift_lustre_api.h"
 #include "tools.h"
 #include "units.h"
 #include "version.h"
 #include "xmf.h"
 
 /* Are we timing the i/o? */
-//#define IO_SPEED_MEASUREMENT
+// #define IO_SPEED_MEASUREMENT
+
+/* Max number of entries that can be written for a given particle type */
+static const int io_max_size_output_list = 100;
 
 /**
  * @brief Writes a data array in given HDF5 group.
@@ -236,6 +240,9 @@ void write_distributed_array(
   io_write_attribute_f(h_data, "a-scale exponent", props.scale_factor_exponent);
   io_write_attribute_s(h_data, "Expression for physical CGS units", buffer);
   io_write_attribute_s(h_data, "Lossy compression filter", comp_buffer);
+  io_write_attribute_b(h_data, "Value stored as physical", props.is_physical);
+  io_write_attribute_b(h_data, "Property can be converted to comoving",
+                       props.is_convertible_to_comoving);
 
   /* Write the actual number this conversion factor corresponds to */
   const double factor =
@@ -401,6 +408,9 @@ void write_array_virtual(struct engine* e, hid_t grp, const char* fileName_base,
   io_write_attribute_f(h_data, "a-scale exponent", props.scale_factor_exponent);
   io_write_attribute_s(h_data, "Expression for physical CGS units", buffer);
   io_write_attribute_s(h_data, "Lossy compression filter", comp_buffer);
+  io_write_attribute_b(h_data, "Value stored as physical", props.is_physical);
+  io_write_attribute_b(h_data, "Property can be converted to comoving",
+                       props.is_convertible_to_comoving);
 
   /* Write the actual number this conversion factor corresponds to */
   const double factor =
@@ -499,8 +509,14 @@ void write_virtual_file(struct engine* e, const char* fileName_base,
    * specific output */
   xmf_write_outputheader(xmfFile, fileName, e->time);
 
+  /* Set the minimal API version to avoid issues with advanced features */
+  hid_t h_props = H5Pcreate(H5P_FILE_ACCESS);
+  herr_t err = H5Pset_libver_bounds(h_props, HDF5_LOWEST_FILE_FORMAT_VERSION,
+                                    HDF5_HIGHEST_FILE_FORMAT_VERSION);
+  if (err < 0) error("Error setting the hdf5 API version");
+
   /* Open HDF5 file with the chosen parameters */
-  hid_t h_file = H5Fcreate(fileName, H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
+  hid_t h_file = H5Fcreate(fileName, H5F_ACC_TRUNC, H5P_DEFAULT, h_props);
   if (h_file < 0) error("Error while opening file '%s'.", fileName);
 
   /* Open header to write simulation properties */
@@ -642,8 +658,8 @@ void write_virtual_file(struct engine* e, const char* fileName_base,
     io_write_attribute_ll(h_grp, "TotalNumberOfParticles", N_total[ptype]);
 
     int num_fields = 0;
-    struct io_props list[100];
-    bzero(list, 100 * sizeof(struct io_props));
+    struct io_props list[io_max_size_output_list];
+    bzero(list, io_max_size_output_list * sizeof(struct io_props));
 
     /* Write particle fields from the particle structure */
     switch (ptype) {
@@ -684,6 +700,15 @@ void write_virtual_file(struct engine* e, const char* fileName_base,
         error("Particle Type %d not yet supported. Aborting", ptype);
     }
 
+    /* Verify we are not going to crash when writing below */
+    if (num_fields >= io_max_size_output_list)
+      error("Too many fields to write for particle type %d", ptype);
+    for (int i = 0; i < num_fields; ++i) {
+      if (!list[i].is_used) error("List of field contains an empty entry!");
+      if (!list[i].dimension)
+        error("Dimension of field '%s' is <= 1!", list[i].name);
+    }
+
     /* Did the user specify a non-standard default for the entire particle
      * type? */
     const enum lossy_compression_schemes compression_level_current_default =
@@ -723,8 +748,9 @@ void write_virtual_file(struct engine* e, const char* fileName_base,
   /* Write LXMF file descriptor */
   xmf_write_outputfooter(xmfFile, e->snapshot_output_count, e->time);
 
-  /* Close the file for now */
+  /* Close the file */
   H5Fclose(h_file);
+  H5Pclose(h_props);
 
 #else
   error(
@@ -943,7 +969,7 @@ void write_output_distributed(struct engine* e,
   }
 
   /* Compute offset in the file and total number of particles */
-  const long long N[swift_type_count] = {
+  long long N[swift_type_count] = {
       Ngas_written,   Ndm_written,         Ndm_background, Nsinks_written,
       Nstars_written, Nblackholes_written, Ndm_neutrino};
 
@@ -967,26 +993,59 @@ void write_output_distributed(struct engine* e,
   };
 
   /* Use a single Lustre stripe with a rank-based OST offset? */
-  if (e->snapshot_lustre_OST_count != 0) {
-
-    /* Use a random offset to avoid placing things in the same OSTs. We do
-     * this to keep the use of OSTs balanced, much like using -1 for the
-     * stripe. */
-    int offset = rand() % e->snapshot_lustre_OST_count;
-    MPI_Bcast(&offset, 1, MPI_INT, 0, MPI_COMM_WORLD);
-
-    char string[1200];
-    sprintf(string, "lfs setstripe -c 1 -i %d %s",
-            ((e->nodeID + offset) % e->snapshot_lustre_OST_count), fileName);
-    const int result = system(string);
-    if (result != 0) {
-      message("lfs setstripe command returned error code %d", result);
+  if (e->snapshot_lustre_OST_checks != 0) {
+
+    /* Gather information about the current state of the OSTs. */
+    struct swift_ost_store ost_infos;
+
+    /* Select good OSTs sorted by free space. */
+    if (e->nodeID == 0) {
+      swift_ost_select(&ost_infos, fileName, e->snapshot_lustre_OST_free,
+                       e->snapshot_lustre_OST_test, e->verbose);
+    }
+
+    /* Distribute the OST information. */
+    MPI_Bcast(&ost_infos, sizeof(struct swift_ost_store), MPI_BYTE, 0,
+              MPI_COMM_WORLD);
+
+    /* Need to make space for the OSTs and copy those locally. If the count is
+     * zero this is probably not a lustre mount. */
+    if (ost_infos.count > 0) {
+      if (e->nodeID != 0) swift_ost_store_alloc(&ost_infos, ost_infos.size);
+      MPI_Bcast(ost_infos.infos, sizeof(struct swift_ost_info) * ost_infos.size,
+                MPI_BYTE, 0, MPI_COMM_WORLD);
+
+      /* We now know how many OSTs are available, each rank should attempt to
+       * use a different one, but overtime we should try not to use the same
+       * ones. Culling will order things by free space so we should get some
+       * reordering of those if we do this process each time. */
+      int dummy = e->nodeID;
+      int offset = swift_ost_next(&ost_infos, &dummy, 1);
+
+      /* And create the file with a stripe of 1 on the OST. */
+      const int result = swift_create_striped_file(fileName, offset, 1, &dummy);
+      if (result != 0) message("failed to set stripe of snapshot");
+
+      /* Finished with this. */
+      swift_ost_store_free(&ost_infos);
+    } else if (e->nodeID == 0) {
+      swift_ost_store_free(&ost_infos);
+
+      /* Don't try this again until next launch. */
+      e->snapshot_lustre_OST_checks = 0;
+      message("Disabling further lustre OST checks");
     }
   }
 
+  /* Set the minimal API version to avoid issues with advanced features */
+  hid_t h_props = H5Pcreate(H5P_FILE_ACCESS);
+  herr_t err = H5Pset_libver_bounds(h_props, HDF5_LOWEST_FILE_FORMAT_VERSION,
+                                    HDF5_HIGHEST_FILE_FORMAT_VERSION);
+  if (err < 0) error("Error setting the hdf5 API version");
+
   /* Open file */
   /* message("Opening file '%s'.", fileName); */
-  h_file = H5Fcreate(fileName, H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
+  h_file = H5Fcreate(fileName, H5F_ACC_TRUNC, H5P_DEFAULT, h_props);
   if (h_file < 0) error("Error while opening file '%s'.", fileName);
 
   /* Open header to write simulation properties */
@@ -1147,8 +1206,8 @@ void write_output_distributed(struct engine* e,
     io_write_attribute_ll(h_grp, "TotalNumberOfParticles", N_total[ptype]);
 
     int num_fields = 0;
-    struct io_props list[100];
-    bzero(list, 100 * sizeof(struct io_props));
+    struct io_props list[io_max_size_output_list];
+    bzero(list, io_max_size_output_list * sizeof(struct io_props));
     size_t Nparticles = 0;
 
     struct part* parts_written = NULL;
@@ -1416,6 +1475,15 @@ void write_output_distributed(struct engine* e,
         error("Particle Type %d not yet supported. Aborting", ptype);
     }
 
+    /* Verify we are not going to crash when writing below */
+    if (num_fields >= io_max_size_output_list)
+      error("Too many fields to write for particle type %d", ptype);
+    for (int i = 0; i < num_fields; ++i) {
+      if (!list[i].is_used) error("List of field contains an empty entry!");
+      if (!list[i].dimension)
+        error("Dimension of field '%s' is <= 1!", list[i].name);
+    }
+
     /* Did the user specify a non-standard default for the entire particle
      * type? */
     const enum lossy_compression_schemes compression_level_current_default =
@@ -1463,6 +1531,7 @@ void write_output_distributed(struct engine* e,
 
   /* Close file */
   H5Fclose(h_file);
+  H5Pclose(h_props);
 
 #if H5_VERSION_GE(1, 10, 0)
 
@@ -1485,8 +1554,13 @@ void write_output_distributed(struct engine* e,
     char fileName_virtual[1030];
     sprintf(fileName_virtual, "%s.hdf5", fileName_base);
 
+    h_props = H5Pcreate(H5P_FILE_ACCESS);
+    err = H5Pset_libver_bounds(h_props, HDF5_LOWEST_FILE_FORMAT_VERSION,
+                               HDF5_HIGHEST_FILE_FORMAT_VERSION);
+    if (err < 0) error("Error setting the hdf5 API version");
+
     /* Open the snapshot on rank 0 */
-    h_file_cells = H5Fopen(fileName_virtual, H5F_ACC_RDWR, H5P_DEFAULT);
+    h_file_cells = H5Fopen(fileName_virtual, H5F_ACC_RDWR, h_props);
     if (h_file_cells < 0)
       error("Error while opening file '%s' on rank %d.", fileName_virtual,
             mpi_rank);
@@ -1514,6 +1588,7 @@ void write_output_distributed(struct engine* e,
   if (mpi_rank == 0) {
     H5Gclose(h_grp_cells);
     H5Fclose(h_file_cells);
+    H5Pclose(h_props);
   }
 
 #endif

src/drift.h

@@ -23,6 +23,7 @@
 #include <config.h>
 
 /* Local headers. */
+#include "adaptive_softening.h"
 #include "black_holes.h"
 #include "const.h"
 #include "debug.h"
@@ -205,6 +206,7 @@ __attribute__((always_inline)) INLINE static void drift_part(
   mhd_predict_extra(p, xp, dt_drift, dt_therm, cosmo, hydro_props,
                     entropy_floor);
   rt_predict_extra(p, xp, dt_drift);
+  if (p->gpart) gravity_update_softening(p->gpart, p, e->gravity_properties);
 
   /* Compute offsets since last cell construction */
   for (int k = 0; k < 3; k++) {
@@ -390,6 +392,25 @@ __attribute__((always_inline)) INLINE static void drift_sink(
   sink->ti_drift = ti_current;
 #endif
 
+#ifdef SWIFT_FIXED_BOUNDARY_PARTICLES
+
+  /* Get the ID of the gpart */
+  const long long id = sink->id;
+
+  /* Cancel the velocity of the particles */
+  if (id < SWIFT_FIXED_BOUNDARY_PARTICLES) {
+
+    /* Don't move! */
+    sink->v[0] = 0.f;
+    sink->v[1] = 0.f;
+    sink->v[2] = 0.f;
+  }
+#endif
+
+#ifdef WITH_LIGHTCONE
+  error("Lightcone treatment of sinks needs implementing");
+#endif
+
   /* Drift... */
   sink->x[0] += sink->v[0] * dt_drift;
   sink->x[1] += sink->v[1] * dt_drift;

src/engine.c

@@ -134,7 +134,9 @@ const char *engine_policy_names[] = {"none",
                                      "line of sight",
                                      "sink",
                                      "rt",
-                                     "power spectra"};
+                                     "power spectra",
+                                     "moving mesh",
+                                     "moving mesh hydro"};
 
 const int engine_default_snapshot_subsample[swift_type_count] = {0};
 
@@ -322,7 +324,7 @@ void engine_repartition_trigger(struct engine *e) {
           e->usertime_last_step, e->systime_last_step, (double)resident,
           e->local_deadtime / (e->nr_threads * e->wallclock_time)};
       double timemems[e->nr_nodes * 4];
-      MPI_Gather(&timemem, 4, MPI_DOUBLE, timemems, 4, MPI_DOUBLE, 0,
+      MPI_Gather(timemem, 4, MPI_DOUBLE, timemems, 4, MPI_DOUBLE, 0,
                  MPI_COMM_WORLD);
       if (e->nodeID == 0) {
 
@@ -491,6 +493,34 @@ void engine_exchange_cells(struct engine *e) {
 #endif
 }
 
+/**
+ * @brief Exchange extra information for the grid construction with other nodes.
+ *
+ * @param e The #engine.
+ */
+void engine_exchange_grid_extra(struct engine *e) {
+
+#ifdef WITH_MPI
+
+#ifdef SWIFT_DEBUG_CHECKS
+  if (!(e->policy & engine_policy_grid))
+    error("Not running with grid, but trying to exchange grid information!");
+#endif
+
+  const ticks tic = getticks();
+
+  /* Exchange the grid info with neighbouring ranks. */
+  proxy_grid_extra_exchange(e->proxies, e->nr_proxies, e->s);
+
+  if (e->verbose)
+    message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
+            clocks_getunit());
+
+#else
+  error("SWIFT was not compiled with MPI support.");
+#endif
+}
+
 /**
  * @brief Exchanges the top-level multipoles between all the nodes
  * such that every node has a multipole for each top-level cell.
@@ -747,16 +777,18 @@ void engine_allocate_foreign_particles(struct engine *e, const int fof) {
 #ifdef WITH_MPI
 
   const int nr_proxies = e->nr_proxies;
-  const int with_hydro = e->policy & engine_policy_hydro;
+  const int with_hydro =
+      e->policy & (engine_policy_hydro | engine_policy_grid_hydro);
   const int with_stars = e->policy & engine_policy_stars;
   const int with_black_holes = e->policy & engine_policy_black_holes;
+  const int with_sinks = e->policy & engine_policy_sinks;
   struct space *s = e->s;
   ticks tic = getticks();
 
   /* Count the number of particles we need to import and re-allocate
      the buffer if needed. */
   size_t count_parts_in = 0, count_gparts_in = 0, count_sparts_in = 0,
-         count_bparts_in = 0;
+         count_bparts_in = 0, count_sinks_in = 0;
   for (int k = 0; k < nr_proxies; k++) {
     for (int j = 0; j < e->proxies[k].nr_cells_in; j++) {
 
@@ -775,6 +807,11 @@ void engine_allocate_foreign_particles(struct engine *e, const int fof) {
 
       /* For black holes, we just use the numbers in the top-level cells */
       count_bparts_in += e->proxies[k].cells_in[j]->black_holes.count;
+
+      /* For sinks, we just use the numbers in the top-level cells + some
+         extra space */
+      count_sinks_in +=
+          e->proxies[k].cells_in[j]->sinks.count + space_extra_sinks;
     }
   }
 
@@ -801,16 +838,27 @@ void engine_allocate_foreign_particles(struct engine *e, const int fof) {
 
   /* Allocate space for the foreign particles we will receive */
   size_t old_size_gparts_foreign = s->size_gparts_foreign;
-  if (count_gparts_in > s->size_gparts_foreign) {
+  if (!fof && count_gparts_in > s->size_gparts_foreign) {
     if (s->gparts_foreign != NULL)
       swift_free("gparts_foreign", s->gparts_foreign);
     s->size_gparts_foreign = engine_foreign_alloc_margin * count_gparts_in;
-    if (swift_memalign("gparts_foreign", (void **)&s->gparts_foreign,
-                       gpart_align,
-                       sizeof(struct gpart) * s->size_gparts_foreign) != 0)
+    if (swift_memalign(
+            "gparts_foreign", (void **)&s->gparts_foreign, gpart_align,
+            sizeof(struct gpart_foreign) * s->size_gparts_foreign) != 0)
       error("Failed to allocate foreign gpart data.");
   }
 
+  /* Allocate space for the foreign FOF particles we will receive */
+  if (fof && count_gparts_in > s->size_gparts_foreign) {
+    if (s->gparts_foreign != NULL)
+      swift_free("gparts_fof_foreign", s->gparts_fof_foreign);
+    s->size_gparts_foreign = engine_foreign_alloc_margin * count_gparts_in;
+    if (swift_memalign(
+            "gparts_fof_foreign", (void **)&s->gparts_fof_foreign, gpart_align,
+            sizeof(struct gpart_fof_foreign) * s->size_gparts_foreign) != 0)
+      error("Failed to allocate foreign FOF gpart data.");
+  }
+
   /* Allocate space for the foreign particles we will receive */
   size_t old_size_sparts_foreign = s->size_sparts_foreign;
   if (!fof && count_sparts_in > s->size_sparts_foreign) {
@@ -821,11 +869,15 @@ void engine_allocate_foreign_particles(struct engine *e, const int fof) {
                        spart_align,
                        sizeof(struct spart) * s->size_sparts_foreign) != 0)
       error("Failed to allocate foreign spart data.");
+
+#ifdef SWIFT_DEBUG_CHECKS
     bzero(s->sparts_foreign, s->size_sparts_foreign * sizeof(struct spart));
+
     for (size_t i = 0; i < s->size_sparts_foreign; ++i) {
       s->sparts_foreign[i].time_bin = time_bin_not_created;
       s->sparts_foreign[i].id = -43;
     }
+#endif
   }
 
   /* Allocate space for the foreign particles we will receive */
@@ -840,44 +892,80 @@ void engine_allocate_foreign_particles(struct engine *e, const int fof) {
       error("Failed to allocate foreign bpart data.");
   }
 
+  /* Allocate space for the foreign particles we will receive */
+  size_t old_size_sinks_foreign = s->size_sinks_foreign;
+  if (!fof && count_sinks_in > s->size_sinks_foreign) {
+    if (s->sinks_foreign != NULL) swift_free("sinks_foreign", s->sinks_foreign);
+    s->size_sinks_foreign = engine_foreign_alloc_margin * count_sinks_in;
+    if (swift_memalign("sinks_foreign", (void **)&s->sinks_foreign, sink_align,
+                       sizeof(struct sink) * s->size_sinks_foreign) != 0)
+      error("Failed to allocate foreign sink data.");
+
+#ifdef SWIFT_DEBUG_CHECKS
+    bzero(s->sinks_foreign, s->size_sinks_foreign * sizeof(struct sink));
+
+    /* Note: If you ever see a sink particle with id = -666, the following
+       lines is the ones that sets the ID to this value. */
+    for (size_t i = 0; i < s->size_sinks_foreign; ++i) {
+      s->sinks_foreign[i].time_bin = time_bin_not_created;
+      s->sinks_foreign[i].id = -666;
+    }
+#endif
+  }
+
   if (e->verbose) {
     message(
-        "Allocating %zd/%zd/%zd/%zd foreign part/gpart/spart/bpart "
-        "(%zd/%zd/%zd/%zd MB)",
+        "Allocating %zd/%zd/%zd/%zd/%zd foreign part/gpart/spart/bpart/sink "
+        "(%zd/%zd/%zd/%zd/%zd MB)",
         s->size_parts_foreign, s->size_gparts_foreign, s->size_sparts_foreign,
-        s->size_bparts_foreign,
+        s->size_bparts_foreign, s->size_sinks_foreign,
         s->size_parts_foreign * sizeof(struct part) / (1024 * 1024),
-        s->size_gparts_foreign * sizeof(struct gpart) / (1024 * 1024),
+        s->size_gparts_foreign * sizeof(struct gpart_foreign) / (1024 * 1024),
         s->size_sparts_foreign * sizeof(struct spart) / (1024 * 1024),
-        s->size_bparts_foreign * sizeof(struct bpart) / (1024 * 1024));
+        s->size_bparts_foreign * sizeof(struct bpart) / (1024 * 1024),
+        s->size_sinks_foreign * sizeof(struct sink) / (1024 * 1024));
 
     if ((s->size_parts_foreign - old_size_parts_foreign) > 0 ||
         (s->size_gparts_foreign - old_size_gparts_foreign) > 0 ||
         (s->size_sparts_foreign - old_size_sparts_foreign) > 0 ||
-        (s->size_bparts_foreign - old_size_bparts_foreign) > 0) {
+        (s->size_bparts_foreign - old_size_bparts_foreign) > 0 ||
+        (s->size_sinks_foreign - old_size_sinks_foreign) > 0) {
       message(
-          "Re-allocations %zd/%zd/%zd/%zd part/gpart/spart/bpart "
-          "(%zd/%zd/%zd/%zd MB)",
+          "Re-allocations %zd/%zd/%zd/%zd/%zd part/gpart/spart/bpart/sink "
+          "(%zd/%zd/%zd/%zd/%zd MB)",
           (s->size_parts_foreign - old_size_parts_foreign),
           (s->size_gparts_foreign - old_size_gparts_foreign),
           (s->size_sparts_foreign - old_size_sparts_foreign),
           (s->size_bparts_foreign - old_size_bparts_foreign),
+          (s->size_sinks_foreign - old_size_sinks_foreign),
           (s->size_parts_foreign - old_size_parts_foreign) *
               sizeof(struct part) / (1024 * 1024),
-          (s->size_gparts_foreign - old_size_gparts_foreign) *
-              sizeof(struct gpart) / (1024 * 1024),
+          fof ? (s->size_gparts_foreign - old_size_gparts_foreign) *
+                    sizeof(struct gpart_fof_foreign) / (1024 * 1024)
+              : (s->size_gparts_foreign - old_size_gparts_foreign) *
+                    sizeof(struct gpart_foreign) / (1024 * 1024),
           (s->size_sparts_foreign - old_size_sparts_foreign) *
               sizeof(struct spart) / (1024 * 1024),
           (s->size_bparts_foreign - old_size_bparts_foreign) *
-              sizeof(struct bpart) / (1024 * 1024));
+              sizeof(struct bpart) / (1024 * 1024),
+          (s->size_sinks_foreign - old_size_sinks_foreign) *
+              sizeof(struct sink) / (1024 * 1024));
     }
   }
 
+  if (e->verbose)
+    message("Allocating and zeroing arrays took %.3f %s.",
+            clocks_from_ticks(getticks() - tic), clocks_getunit());
+
+  tic = getticks();
+
   /* Unpack the cells and link to the particle data. */
   struct part *parts = s->parts_foreign;
-  struct gpart *gparts = s->gparts_foreign;
+  struct gpart_foreign *gparts_foreign = s->gparts_foreign;
+  struct gpart_fof_foreign *gparts_fof_foreign = s->gparts_fof_foreign;
   struct spart *sparts = s->sparts_foreign;
   struct bpart *bparts = s->bparts_foreign;
+  struct sink *sinks = s->sinks_foreign;
   for (int k = 0; k < nr_proxies; k++) {
     for (int j = 0; j < e->proxies[k].nr_cells_in; j++) {
 
@@ -888,11 +976,18 @@ void engine_allocate_foreign_particles(struct engine *e, const int fof) {
         parts = &parts[count_parts];
       }
 
-      if (e->proxies[k].cells_in_type[j] & proxy_cell_type_gravity) {
+      if (fof && e->proxies[k].cells_in_type[j] & proxy_cell_type_gravity) {
+
+        const size_t count_gparts = cell_link_foreign_fof_gparts(
+            e->proxies[k].cells_in[j], gparts_fof_foreign);
+        gparts_fof_foreign = &gparts_fof_foreign[count_gparts];
+      }
+
+      if (!fof && e->proxies[k].cells_in_type[j] & proxy_cell_type_gravity) {
 
         const size_t count_gparts =
-            cell_link_foreign_gparts(e->proxies[k].cells_in[j], gparts);
-        gparts = &gparts[count_gparts];
+            cell_link_foreign_gparts(e->proxies[k].cells_in[j], gparts_foreign);
+        gparts_foreign = &gparts_foreign[count_gparts];
       }
 
       if (!fof && with_stars) {
@@ -909,14 +1004,26 @@ void engine_allocate_foreign_particles(struct engine *e, const int fof) {
         cell_link_bparts(e->proxies[k].cells_in[j], bparts);
         bparts = &bparts[e->proxies[k].cells_in[j]->black_holes.count];
       }
+
+      if (!fof && with_sinks) {
+
+        /* For sinks, we just use the numbers in the top-level cells */
+        cell_link_sinks(e->proxies[k].cells_in[j], sinks);
+        sinks =
+            &sinks[e->proxies[k].cells_in[j]->sinks.count + space_extra_sinks];
+      }
     }
   }
 
   /* Update the counters */
   s->nr_parts_foreign = parts - s->parts_foreign;
-  s->nr_gparts_foreign = gparts - s->gparts_foreign;
+  if (fof)
+    s->nr_gparts_foreign = gparts_fof_foreign - s->gparts_fof_foreign;
+  else
+    s->nr_gparts_foreign = gparts_foreign - s->gparts_foreign;
   s->nr_sparts_foreign = sparts - s->sparts_foreign;
   s->nr_bparts_foreign = bparts - s->bparts_foreign;
+  s->nr_sinks_foreign = sinks - s->sinks_foreign;
 
   if (e->verbose)
     message("Recursively linking foreign arrays took %.3f %s.",
@@ -1016,6 +1123,19 @@ void engine_print_task_counts(const struct engine *e) {
   message("nr_sparts = %zu.", e->s->nr_sparts);
   message("nr_bparts = %zu.", e->s->nr_bparts);
 
+#if defined(SWIFT_DEBUG_CHECKS) && defined(WITH_MPI)
+  if (e->verbose == 2) {
+    /* check that the global number of sends matches the global number of
+       recvs */
+    int global_counts[2] = {counts[task_type_send], counts[task_type_recv]};
+    MPI_Allreduce(MPI_IN_PLACE, global_counts, 2, MPI_INT, MPI_SUM,
+                  MPI_COMM_WORLD);
+    if (global_counts[0] != global_counts[1])
+      error("Missing communications (%i sends, %i recvs)!", global_counts[0],
+            global_counts[1]);
+  }
+#endif
+
   if (e->verbose)
     message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
             clocks_getunit());
@@ -1104,8 +1224,15 @@ int engine_estimate_nr_tasks(const struct engine *e) {
 #endif
   }
   if (e->policy & engine_policy_sinks) {
-    /* 1 drift, 2 kicks, 1 time-step, 1 sink formation */
-    n1 += 5;
+    /* 1 drift, 2 kicks, 1 time-step, 1 sink formation     | 5
+       density: 1 self + 13 pairs                          | 14
+       swallow: 1 self + 13 pairs                          | 14
+       do_gas_swallow: 1 self + 13 pairs                   | 14
+       do_sink_swallow: 1 self + 13 pairs                  | 14
+       ghosts: density_ghost, sink_ghost_1, sink_ghost_2   | 3
+       implicit: sink_in,  sink_out                        | 2 */
+    n1 += 66;
+    n2 += 3;
     if (e->policy & engine_policy_stars) {
       /* 1 star formation */
       n1 += 1;
@@ -1130,6 +1257,35 @@ int engine_estimate_nr_tasks(const struct engine *e) {
     n1 += 39;
 #ifdef WITH_MPI
     n1 += 2;
+#endif
+  }
+  if (e->policy & engine_policy_grid) {
+    /* Grid construction: 1 self + 26 (asymmetric) pairs + 1 ghost + 1 sort */
+    n1 += 29;
+    n2 += 3;
+#ifdef WITH_MPI
+    n1 += 3;
+#endif
+  }
+  if (e->policy & engine_policy_grid_hydro) {
+    /* slope estimate: 1 self + 13 pairs (on average)       |   14
+     * others: 1 ghosts, 2 kicks, 1 drift, 1 timestep       | +  7
+     * Total:                                               =   21 */
+    n1 += 21;
+    n2 += 2;
+#ifdef EXTRA_HYDRO_LOOP
+    /* slope limiter: 1 self + 13 pairs                     | + 14
+     * flux: 1 self + 13 pairs                              | + 14
+     * others: 2 ghost.                                     | +  2
+     * Total:                                               =   30  */
+    n1 += 30;
+    n2 += 3;
+#endif
+#ifdef WITH_MPI
+    n1 += 1;
+#ifdef EXTRA_HYDRO_LOOP
+    n1 += 1;
+#endif
 #endif
   }
 
@@ -1145,7 +1301,8 @@ int engine_estimate_nr_tasks(const struct engine *e) {
     struct cell *c = &e->s->cells_top[k];
 
     /* Any cells with particles will have tasks (local & foreign). */
-    int nparts = c->hydro.count + c->grav.count + c->stars.count;
+    int nparts = c->hydro.count + c->grav.count + c->stars.count +
+                 c->black_holes.count + c->sinks.count;
     if (nparts > 0) {
       ntop++;
       ncells++;
@@ -1281,7 +1438,7 @@ void engine_rebuild(struct engine *e, const int repartitioned,
                 MPI_COMM_WORLD);
   MPI_Allreduce(MPI_IN_PLACE, &e->s->max_softening, 1, MPI_FLOAT, MPI_MAX,
                 MPI_COMM_WORLD);
-  MPI_Allreduce(MPI_IN_PLACE, &e->s->max_mpole_power,
+  MPI_Allreduce(MPI_IN_PLACE, e->s->max_mpole_power,
                 SELF_GRAVITY_MULTIPOLE_ORDER + 1, MPI_FLOAT, MPI_MAX,
                 MPI_COMM_WORLD);
 #endif
@@ -1307,6 +1464,12 @@ void engine_rebuild(struct engine *e, const int repartitioned,
   /* Initial cleaning up session ? */
   if (clean_smoothing_length_values) space_sanitize(e->s);
 
+  /* Set the initial completeness flag for the moving mesh (before exchange) */
+  if (e->policy & engine_policy_grid) {
+    cell_grid_set_self_completeness_mapper(e->s->cells_top, e->s->nr_cells,
+                                           NULL);
+  }
+
 /* If in parallel, exchange the cell structure, top-level and neighbouring
  * multipoles. To achieve this, free the foreign particle buffers first. */
 #ifdef WITH_MPI
@@ -1330,8 +1493,27 @@ void engine_rebuild(struct engine *e, const int repartitioned,
     if (counter != e->total_nr_gparts)
       error("Total particles in multipoles inconsistent with engine");
   }
+  if (e->policy & engine_policy_grid) {
+    for (int i = 0; i < e->s->nr_cells; i++) {
+      const struct cell *ci = &e->s->cells_top[i];
+      if (ci->hydro.count > 0 && !(ci->grid.self_completeness == grid_complete))
+        error("Encountered incomplete top level cell!");
+    }
+  }
 #endif
 
+  /* Set the grid construction level, is needed before splitting the tasks */
+  if (e->policy & engine_policy_grid) {
+    /* Set the completeness and construction level */
+    threadpool_map(&e->threadpool, cell_set_grid_completeness_mapper, NULL,
+                   e->s->nr_cells, 1, threadpool_auto_chunk_size, e);
+    threadpool_map(&e->threadpool, cell_set_grid_construction_level_mapper,
+                   NULL, e->s->nr_cells, 1, threadpool_auto_chunk_size, e);
+#ifdef WITH_MPI
+    engine_exchange_grid_extra(e);
+#endif
+  }
+
   /* Re-build the tasks. */
   engine_maketasks(e);
 
@@ -1362,9 +1544,9 @@ void engine_rebuild(struct engine *e, const int repartitioned,
   space_check_unskip_flags(e->s);
 #endif
 
-  /* Run through the tasks and mark as skip or not. */
-  if (engine_marktasks(e))
-    error("engine_marktasks failed after space_rebuild.");
+  /* Run through the cells, and their tasks to mark as unskipped. */
+  engine_unskip(e);
+  if (e->forcerebuild) error("engine_unskip faled after a rebuild!");
 
   /* Print the status of the system */
   if (e->verbose) engine_print_task_counts(e);
@@ -1643,6 +1825,7 @@ void engine_skip_force_and_kick(struct engine *e) {
         t->subtype == task_subtype_part_prep1 ||
         t->subtype == task_subtype_spart_prep2 ||
         t->subtype == task_subtype_sf_counts ||
+        t->subtype == task_subtype_grav_counts ||
         t->subtype == task_subtype_rt_gradient ||
         t->subtype == task_subtype_rt_transport)
       t->skip = 1;
@@ -1771,6 +1954,47 @@ void engine_get_max_ids(struct engine *e) {
 #endif
 }
 
+/**
+ * @brief Gather the information about the top-level cells whose time-step has
+ * changed and activate the communications required to synchonize the
+ * time-steps.
+ *
+ * @param e The #engine.
+ */
+void engine_synchronize_times(struct engine *e) {
+
+#ifdef WITH_MPI
+
+  const ticks tic_start = getticks();
+
+  /* Collect which top-level cells have been updated */
+  MPI_Allreduce(MPI_IN_PLACE, e->s->cells_top_updated, e->s->nr_cells, MPI_CHAR,
+                MPI_SUM, MPI_COMM_WORLD);
+
+  /* Activate tend communications involving the cells that have changed. */
+  for (int i = 0; i < e->s->nr_cells; ++i) {
+
+    if (e->s->cells_top_updated[i]) {
+
+      struct cell *c = &e->s->cells_top[i];
+      scheduler_activate_all_subtype(&e->sched, c->mpi.send, task_subtype_tend);
+      scheduler_activate_all_subtype(&e->sched, c->mpi.recv, task_subtype_tend);
+    }
+  }
+
+  if (e->verbose)
+    message("Gathering and activating tend took %.3f %s.",
+            clocks_from_ticks(getticks() - tic_start), clocks_getunit());
+
+  TIMER_TIC;
+  engine_launch(e, "tend");
+  TIMER_TOC(timer_runners);
+
+#else
+  error("SWIFT was not compiled with MPI support.");
+#endif
+}
+
 /**
  * @brief Run the radiative transfer sub-cycles outside the
  * regular time-steps.
@@ -1820,13 +2044,14 @@ void engine_run_rt_sub_cycles(struct engine *e) {
   /* Get some time variables for printouts. Don't update the ones in the
    * engine like in the regular step, or the outputs in the regular steps
    * will be wrong. */
-  /* think cosmology one day: needs adapting here */
-  /* Also needs adapting further below - we print out current values of a
-   * and z. They need to be updated in the engine. */
-  if (e->policy & engine_policy_cosmology)
-    error("Can't run RT subcycling with cosmology yet");
-  const double dt_subcycle = rt_step_size * e->time_base;
-  double time = e->ti_current_subcycle * e->time_base + e->time_begin;
+  double dt_subcycle;
+  if (e->policy & engine_policy_cosmology) {
+    dt_subcycle =
+        cosmology_get_delta_time(e->cosmology, e->ti_current, e->ti_rt_end_min);
+  } else {
+    dt_subcycle = rt_step_size * e->time_base;
+  }
+  double time = e->time;
 
   /* Keep track and accumulate the deadtime over all sub-cycles. */
   /* We need to manually put this back in the engine struct when
@@ -1886,10 +2111,20 @@ void engine_run_rt_sub_cycles(struct engine *e) {
     e->max_active_bin_subcycle = get_max_active_bin(e->ti_current_subcycle);
     e->min_active_bin_subcycle =
         get_min_active_bin(e->ti_current_subcycle, ti_subcycle_old);
-    /* think cosmology one day: needs adapting here */
-    if (e->policy & engine_policy_cosmology)
-      error("Can't run RT subcycling with cosmology yet");
-    time = e->ti_current_subcycle * e->time_base + e->time_begin;
+
+    /* Update rt properties */
+    rt_props_update(e->rt_props, e->internal_units, e->cosmology);
+
+    if (e->policy & engine_policy_cosmology) {
+      double time_old = time;
+      cosmology_update(
+          e->cosmology, e->physical_constants,
+          e->ti_current_subcycle);  // Update cosmological parameters
+      time = e->cosmology->time;    // Grab new cosmology time
+      dt_subcycle = time - time_old;
+    } else {
+      time = e->ti_current_subcycle * e->time_base + e->time_begin;
+    }
 
     /* Do the actual work now. */
     engine_unskip_rt_sub_cycle(e);
@@ -1991,6 +2226,10 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs,
   if (e->nodeID == 0) message("Setting particles to a valid state...");
   engine_first_init_particles(e);
 
+  /* Initialise the particle splitting mechanism */
+  if (e->hydro_properties->particle_splitting)
+    engine_init_split_gas_particles(e);
+
   if (e->nodeID == 0)
     message("Computing initial gas densities and approximate gravity.");
 
@@ -2028,6 +2267,9 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs,
     cooling_update(e->physical_constants, e->cosmology, e->pressure_floor_props,
                    e->cooling_func, e->s, e->time);
 
+  if (e->policy & engine_policy_rt)
+    rt_props_update(e->rt_props, e->internal_units, e->cosmology);
+
 #ifdef WITH_CSDS
   if (e->policy & engine_policy_csds) {
     /* Mark the first time step in the particle csds file. */
@@ -2045,6 +2287,9 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs,
   }
 #endif
 
+  /* Zero the list of cells that have had their time-step updated */
+  bzero(e->s->cells_top_updated, e->s->nr_cells * sizeof(char));
+
   /* Now, launch the calculation */
   TIMER_TIC;
   engine_launch(e, "tasks");
@@ -2069,7 +2314,7 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs,
 
     /* Correct what we did (e.g. in PE-SPH, need to recompute rho_bar) */
     if (hydro_need_extra_init_loop) {
-      engine_marktasks(e);
+      engine_unskip(e);
       engine_skip_force_and_kick(e);
       engine_launch(e, "tasks");
     }
@@ -2133,11 +2378,19 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs,
   scheduler_write_cell_dependencies(&e->sched, e->verbose, e->step);
   if (e->nodeID == 0) scheduler_write_task_level(&e->sched, e->step);
 
+  /* Zero the list of cells that have had their time-step updated */
+  bzero(e->s->cells_top_updated, e->s->nr_cells * sizeof(char));
+
   /* Run the 0th time-step */
   TIMER_TIC2;
   engine_launch(e, "tasks");
   TIMER_TOC2(timer_runners);
 
+  /* When running over MPI, synchronize top-level cells */
+#ifdef WITH_MPI
+  engine_synchronize_times(e);
+#endif
+
 #ifdef SWIFT_HYDRO_DENSITY_CHECKS
   /* Run the brute-force hydro calculation for some parts */
   if (e->policy & engine_policy_hydro)
@@ -2157,6 +2410,15 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs,
     stars_exact_density_check(e->s, e, /*rel_tol=*/1e-3);
 #endif
 
+#ifdef SWIFT_SINK_DENSITY_CHECKS
+  /* Run the brute-force sink calculation for some sinks */
+  if (e->policy & engine_policy_sinks) sink_exact_density_compute(e->s, e);
+
+  /* Check the accuracy of the sink calculation */
+  if (e->policy & engine_policy_sinks)
+    sink_exact_density_check(e->s, e, /*rel_tol=*/1e-3);
+#endif
+
 #ifdef SWIFT_GRAVITY_FORCE_CHECKS
   /* Check the accuracy of the gravity calculation */
   if (e->policy & engine_policy_self_gravity)
@@ -2269,12 +2531,12 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs,
     for (int i = 0; i < s->nr_cells; i++) {
       struct cell *c = &s->cells_top[i];
       if (c->nodeID == engine_rank && c->sinks.count > 0) {
-        float sink_h_max = c->sinks.parts[0].r_cut;
+        float sink_h_max = c->sinks.parts[0].h;
         for (int k = 1; k < c->sinks.count; k++) {
-          if (c->sinks.parts[k].r_cut > sink_h_max)
-            sink_h_max = c->sinks.parts[k].r_cut;
+          if (c->sinks.parts[k].h > sink_h_max)
+            sink_h_max = c->sinks.parts[k].h;
         }
-        c->sinks.r_cut_max = max(sink_h_max, c->sinks.r_cut_max);
+        c->sinks.h_max = max(sink_h_max, c->sinks.h_max);
       }
     }
   }
@@ -2461,6 +2723,10 @@ int engine_step(struct engine *e) {
     hydro_props_update(e->hydro_properties, e->gravity_properties,
                        e->cosmology);
 
+  /* Update the rt properties */
+  if (e->policy & engine_policy_rt)
+    rt_props_update(e->rt_props, e->internal_units, e->cosmology);
+
   /* Check for any snapshot triggers */
   engine_io_check_snapshot_triggers(e);
 
@@ -2658,11 +2924,19 @@ int engine_step(struct engine *e) {
      want to lose the data from the tasks) */
   space_reset_ghost_histograms(e->s);
 
+  /* Zero the list of cells that have had their time-step updated */
+  bzero(e->s->cells_top_updated, e->s->nr_cells * sizeof(char));
+
   /* Start all the tasks. */
   TIMER_TIC;
   engine_launch(e, "tasks");
   TIMER_TOC(timer_runners);
 
+  /* When running over MPI, synchronize top-level cells */
+#ifdef WITH_MPI
+  engine_synchronize_times(e);
+#endif
+
   /* Now record the CPU times used by the tasks. */
 #ifdef WITH_MPI
   double end_usertime = 0.0;
@@ -2691,6 +2965,15 @@ int engine_step(struct engine *e) {
     stars_exact_density_check(e->s, e, /*rel_tol=*/1e-2);
 #endif
 
+#ifdef SWIFT_SINK_DENSITY_CHECKS
+  /* Run the brute-force sink calculation for some sinks */
+  if (e->policy & engine_policy_sinks) sink_exact_density_compute(e->s, e);
+
+  /* Check the accuracy of the sink calculation */
+  if (e->policy & engine_policy_sinks)
+    sink_exact_density_check(e->s, e, /*rel_tol=*/1e-2);
+#endif
+
 #ifdef SWIFT_GRAVITY_FORCE_CHECKS
   /* Check if we want to run force checks this timestep. */
   if (e->policy & engine_policy_self_gravity) {
@@ -2888,6 +3171,7 @@ void engine_split(struct engine *e, struct partition *initial_partition) {
   if (e->verbose)
     message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
             clocks_getunit());
+
 #else
   error("SWIFT was not compiled with MPI support.");
 #endif
@@ -2999,8 +3283,9 @@ void engine_pin(void) {
   threadpool_set_affinity_mask(entry_affinity);
 
   int pin;
-  for (pin = 0; pin < CPU_SETSIZE && !CPU_ISSET(pin, entry_affinity); ++pin)
-    ;
+  for (pin = 0; pin < CPU_SETSIZE && !CPU_ISSET(pin, entry_affinity); ++pin) {
+    /* Nothing to do here */
+  }
 
   cpu_set_t affinity;
   CPU_ZERO(&affinity);
@@ -3217,8 +3502,12 @@ void engine_init(
       parser_get_opt_param_int(params, "Snapshots:compression", 0);
   e->snapshot_distributed =
       parser_get_opt_param_int(params, "Snapshots:distributed", 0);
-  e->snapshot_lustre_OST_count =
-      parser_get_opt_param_int(params, "Snapshots:lustre_OST_count", 0);
+  e->snapshot_lustre_OST_checks =
+      parser_get_opt_param_int(params, "Snapshots:lustre_OST_checks", 0);
+  e->snapshot_lustre_OST_free =
+      parser_get_opt_param_int(params, "Snapshots:lustre_OST_free", 0);
+  e->snapshot_lustre_OST_test =
+      parser_get_opt_param_int(params, "Snapshots:lustre_OST_test", 0);
   e->snapshot_invoke_stf =
       parser_get_opt_param_int(params, "Snapshots:invoke_stf", 0);
   e->snapshot_invoke_fof =
@@ -3587,6 +3876,9 @@ void engine_recompute_displacement_constraint(struct engine *e) {
 
     /* Apply the dimensionless factor */
     e->dt_max_RMS_displacement = dt * e->max_RMS_displacement_factor;
+    if (e->dt_max_RMS_displacement == 0.f) {
+      error("Setting dt_max_RMS_displacement to 0!");
+    }
 
     if (e->verbose)
       message("max_dt_RMS_displacement = %e", e->dt_max_RMS_displacement);
@@ -3943,7 +4235,7 @@ void engine_struct_restore(struct engine *e, FILE *stream) {
   struct rt_props *rt_properties =
       (struct rt_props *)malloc(sizeof(struct rt_props));
   rt_struct_restore(rt_properties, stream, e->physical_constants,
-                    e->internal_units);
+                    e->internal_units, cosmo);
   e->rt_props = rt_properties;
 
   struct black_holes_props *black_holes_properties =

src/engine.h

@@ -89,8 +89,10 @@ enum engine_policy {
   engine_policy_sinks = (1 << 25),
   engine_policy_rt = (1 << 26),
   engine_policy_power_spectra = (1 << 27),
+  engine_policy_grid = (1 << 28),
+  engine_policy_grid_hydro = (1 << 29),
 };
-#define engine_maxpolicy 28
+#define engine_maxpolicy 30
 extern const char *engine_policy_names[engine_maxpolicy + 1];
 
 /**
@@ -348,7 +350,9 @@ struct engine {
   float snapshot_subsample_fraction[swift_type_count];
   int snapshot_run_on_dump;
   int snapshot_distributed;
-  int snapshot_lustre_OST_count;
+  int snapshot_lustre_OST_checks;
+  int snapshot_lustre_OST_free;
+  int snapshot_lustre_OST_test;
   int snapshot_compression;
   int snapshot_invoke_stf;
   int snapshot_invoke_fof;
@@ -587,8 +591,16 @@ struct engine {
   /* Whether to dump restart files after the last step. */
   int restart_onexit;
 
-  /* Number of Lustre OSTs on the system to use as rank-based striping offset */
-  int restart_lustre_OST_count;
+  /* Perform OST checks and assign each restart file a stripe on the basis of
+   * most free space first. */
+  int restart_lustre_OST_checks;
+
+  /* Free space that an OST should have to be used, -1 makes this
+   * the rss size. In MiB so we can use an int and human sized. */
+  int restart_lustre_OST_free;
+
+  /* Whether to check is OSTs are writable, if not then they are not used. */
+  int restart_lustre_OST_test;
 
   /* Do we free the foreign data before writing restart files? */
   int free_foreign_when_dumping_restart;
@@ -682,7 +694,7 @@ struct engine {
 /* Function prototypes, engine.c. */
 void engine_addlink(struct engine *e, struct link **l, struct task *t);
 void engine_barrier(struct engine *e);
-void engine_compute_next_snapshot_time(struct engine *e);
+void engine_compute_next_snapshot_time(struct engine *e, const int restart);
 void engine_compute_next_stf_time(struct engine *e);
 void engine_compute_next_fof_time(struct engine *e);
 void engine_compute_next_statistics_time(struct engine *e);
@@ -747,7 +759,8 @@ void engine_exchange_strays(struct engine *e, const size_t offset_parts,
                             size_t *Ngpart, const size_t offset_sparts,
                             const int *ind_spart, size_t *Nspart,
                             const size_t offset_bparts, const int *ind_bpart,
-                            size_t *Nbpart);
+                            size_t *Nbpart, const size_t offset_sinks,
+                            const int *ind_sink, size_t *Nsink);
 void engine_rebuild(struct engine *e, int redistributed, int clean_h_values);
 void engine_repartition(struct engine *e);
 void engine_repartition_trigger(struct engine *e);
@@ -772,11 +785,9 @@ void engine_maketasks(struct engine *e);
 /* Function prototypes, engine_maketasks.c. */
 void engine_make_fof_tasks(struct engine *e);
 
-/* Function prototypes, engine_marktasks.c. */
-int engine_marktasks(struct engine *e);
-
 /* Function prototypes, engine_split_particles.c. */
 void engine_split_gas_particles(struct engine *e);
+void engine_init_split_gas_particles(struct engine *e);
 
 #ifdef HAVE_SETAFFINITY
 cpu_set_t *engine_entry_affinity(void);

src/engine_config.c

@@ -37,6 +37,7 @@
 
 /* Local headers. */
 #include "fof.h"
+#include "line_of_sight.h"
 #include "mpiuse.h"
 #include "part.h"
 #include "pressure_floor.h"
@@ -311,8 +312,9 @@ void engine_config(int restart, int fof, struct engine *e,
     int skip = 0;
     for (int k = 0; k < nr_affinity_cores; k++) {
       int c;
-      for (c = skip; c < CPU_SETSIZE && !CPU_ISSET(c, entry_affinity); ++c)
-        ;
+      for (c = skip; c < CPU_SETSIZE && !CPU_ISSET(c, entry_affinity); ++c) {
+        /* Nothing to do here */
+      }
       cpuid[k] = c;
       skip = c + 1;
     }
@@ -407,8 +409,8 @@ void engine_config(int restart, int fof, struct engine *e,
     /* Make sure the corresponding policy is set and make space for the proxies
      */
     e->policy |= engine_policy_mpi;
-    if ((e->proxies = (struct proxy *)calloc(sizeof(struct proxy),
-                                             engine_maxproxies)) == NULL)
+    if ((e->proxies = (struct proxy *)calloc(engine_maxproxies,
+                                             sizeof(struct proxy))) == NULL)
       error("Failed to allocate memory for proxies.");
     e->nr_proxies = 0;
 
@@ -787,14 +789,16 @@ void engine_config(int restart, int fof, struct engine *e,
 #endif
 
     /* Find the time of the first snapshot output */
-    engine_compute_next_snapshot_time(e);
+    engine_compute_next_snapshot_time(e, restart);
 
     /* Find the time of the first statistics output */
     engine_compute_next_statistics_time(e);
 
-    /* Find the time of the first line of sight output */
+    /* Find the time of the first line of sight output
+     * and verify the outputs */
     if (e->policy & engine_policy_line_of_sight) {
       engine_compute_next_los_time(e);
+      los_io_output_check(e);
     }
 
     /* Find the time of the first stf output */
@@ -827,9 +831,13 @@ void engine_config(int restart, int fof, struct engine *e,
      * on restart. */
     e->restart_onexit = parser_get_opt_param_int(params, "Restarts:onexit", 0);
 
-    /* Read the number of Lustre OSTs to distribute the restart files over */
-    e->restart_lustre_OST_count =
-        parser_get_opt_param_int(params, "Restarts:lustre_OST_count", 0);
+    /* Lustre OST options. Disabled by default. */
+    e->restart_lustre_OST_checks =
+        parser_get_opt_param_int(params, "Restarts:lustre_OST_checks", 0);
+    e->restart_lustre_OST_free =
+        parser_get_opt_param_int(params, "Restarts:lustre_OST_free", 0);
+    e->restart_lustre_OST_test =
+        parser_get_opt_param_int(params, "Restarts:lustre_OST_test", 0);
 
     /* Hours between restart dumps. Can be changed on restart. */
     float dhours =
@@ -946,6 +954,8 @@ void engine_config(int restart, int fof, struct engine *e,
         params, "Scheduler:cell_extra_gparts", space_extra_gparts_default);
     space_extra_bparts = parser_get_opt_param_int(
         params, "Scheduler:cell_extra_bparts", space_extra_bparts_default);
+    space_extra_sinks = parser_get_opt_param_int(
+        params, "Scheduler:cell_extra_sinks", space_extra_sinks_default);
 
     /* Do we want any spare particles for on the fly creation?
        This condition should be the same than in space.c */

src/engine_fof.c

@@ -48,10 +48,13 @@ void engine_activate_gpart_comms(struct engine *e) {
 
     struct task *t = &tasks[k];
 
-    if ((t->type == task_type_send) && (t->subtype == task_subtype_gpart)) {
+    if ((t->type == task_type_send) && (t->subtype == task_subtype_fof)) {
+      scheduler_activate(s, t);
+    } else if ((t->type == task_type_pack) &&
+               (t->subtype == task_subtype_fof)) {
       scheduler_activate(s, t);
     } else if ((t->type == task_type_recv) &&
-               (t->subtype == task_subtype_gpart)) {
+               (t->subtype == task_subtype_fof)) {
       scheduler_activate(s, t);
     } else {
       t->skip = 1;
@@ -179,26 +182,34 @@ void engine_fof(struct engine *e, const int dump_results,
 
   /* Compute the local<->foreign group links (nothing to do without MPI)*/
   fof_search_foreign_cells(e->fof_properties, e->s);
+
+  /* Link the foreign fragments and finalise global group list (nothing to do
+   * without MPI) */
+  fof_link_foreign_fragments(e->fof_properties, e->s);
 #endif
 
-  /* Compute the attachable->linkable links */
-  fof_link_attachable_particles(e->fof_properties, e->s);
+  /* Count the total number of (linkable) groups,
+   * assign unique group IDs and set the IDs in the particles */
+  fof_assign_group_ids(e->fof_properties, e->s);
 
 #ifdef WITH_MPI
 
-  /* Free the foreign particles */
-  space_free_foreign_parts(e->s, /*clear pointers=*/1);
+  /* Activate the tasks exchanging all the required gparts */
+  engine_activate_gpart_comms(e);
+
+  /* Perform send and receive tasks. */
+  engine_launch(e, "fof comms");
 
 #endif
 
-  /* Finish the operations attaching the attachables to their groups */
-  fof_finalise_attachables(e->fof_properties, e->s);
+  /* Compute the attachable->linkable links */
+  fof_link_attachable_particles(e->fof_properties, e->s);
 
 #ifdef WITH_MPI
 
-  /* Link the foreign fragments and finalise global group list (nothing to do
-   * without MPI) */
-  fof_link_foreign_fragments(e->fof_properties, e->s);
+  /* Free the foreign particles */
+  space_free_foreign_parts(e->s, /*clear pointers=*/1);
+
 #endif
 
   /* Compute group properties and act on the results
@@ -216,6 +227,13 @@ void engine_fof(struct engine *e, const int dump_results,
   /* ... and find the next FOF time */
   if (seed_black_holes) engine_compute_next_fof_time(e);
 
+  /* Free everything we allocated */
+  fof_free_arrays(e->fof_properties);
+
+#ifdef WITH_MPI
+  MPI_Barrier(MPI_COMM_WORLD);
+#endif
+
   /* Restore the foreign buffers as they were*/
   if (foreign_buffers_allocated) {
 #ifdef WITH_MPI

src/engine_io.c

@@ -565,7 +565,7 @@ void engine_io(struct engine *e) {
 #endif
 
         /* ... and find the next output time */
-        engine_compute_next_snapshot_time(e);
+        engine_compute_next_snapshot_time(e, /*restart=*/0);
         break;
 
       case output_statistics:
@@ -694,17 +694,77 @@ void engine_io(struct engine *e) {
   e->time = time;
 }
 
+/**
+ * @brief Set the value of the recording trigger windows based
+ * on the user's desires and the time to the next snapshot.
+ *
+ * @param e The #engine.
+ */
+void engine_set_and_verify_snapshot_triggers(struct engine *e) {
+
+  integertime_t ti_next_snap = e->ti_next_snapshot;
+  if (ti_next_snap == -1) ti_next_snap = max_nr_timesteps;
+
+  /* Time until the next snapshot */
+  double time_to_next_snap;
+  if (e->policy & engine_policy_cosmology) {
+    time_to_next_snap =
+        cosmology_get_delta_time(e->cosmology, e->ti_current, ti_next_snap);
+  } else {
+    time_to_next_snap = (ti_next_snap - e->ti_current) * e->time_base;
+  }
+
+  /* Do we need to reduce any of the recording trigger times?
+   * Or can we set them with the user's desired range? */
+  for (int k = 0; k < num_snapshot_triggers_part; ++k) {
+    if (e->snapshot_recording_triggers_desired_part[k] > 0) {
+      if (e->snapshot_recording_triggers_desired_part[k] > time_to_next_snap) {
+        e->snapshot_recording_triggers_part[k] = time_to_next_snap;
+      } else {
+        e->snapshot_recording_triggers_part[k] =
+            e->snapshot_recording_triggers_desired_part[k];
+      }
+    }
+  }
+  for (int k = 0; k < num_snapshot_triggers_spart; ++k) {
+    if (e->snapshot_recording_triggers_desired_spart[k] > 0) {
+      if (e->snapshot_recording_triggers_desired_spart[k] > time_to_next_snap) {
+        e->snapshot_recording_triggers_spart[k] = time_to_next_snap;
+      } else {
+        e->snapshot_recording_triggers_spart[k] =
+            e->snapshot_recording_triggers_desired_spart[k];
+      }
+    }
+  }
+  for (int k = 0; k < num_snapshot_triggers_bpart; ++k) {
+    if (e->snapshot_recording_triggers_desired_bpart[k] > 0) {
+      if (e->snapshot_recording_triggers_desired_bpart[k] > time_to_next_snap) {
+        e->snapshot_recording_triggers_bpart[k] = time_to_next_snap;
+      } else {
+        e->snapshot_recording_triggers_bpart[k] =
+            e->snapshot_recording_triggers_desired_bpart[k];
+      }
+    }
+  }
+}
+
 /**
  * @brief Computes the next time (on the time line) for a dump
  *
  * @param e The #engine.
+ * @param restart Are we calling this upon a restart event?
  */
-void engine_compute_next_snapshot_time(struct engine *e) {
+void engine_compute_next_snapshot_time(struct engine *e, const int restart) {
 
   /* Do output_list file case */
   if (e->output_list_snapshots) {
     output_list_read_next_time(e->output_list_snapshots, e, "snapshots",
                                &e->ti_next_snapshot);
+
+    /* Unless we are restarting, check the allowed recording trigger time */
+    if (!restart) engine_set_and_verify_snapshot_triggers(e);
+
+    /* All done in the list case */
     return;
   }
 
@@ -762,49 +822,9 @@ void engine_compute_next_snapshot_time(struct engine *e) {
         message("Next snapshot time set to t=%e.", next_snapshot_time);
     }
 
-    /* Time until the next snapshot */
-    double time_to_next_snap;
-    if (e->policy & engine_policy_cosmology) {
-      time_to_next_snap = cosmology_get_delta_time(e->cosmology, e->ti_current,
-                                                   e->ti_next_snapshot);
-    } else {
-      time_to_next_snap = (e->ti_next_snapshot - e->ti_current) * e->time_base;
-    }
-
-    /* Do we need to reduce any of the recording trigger times? */
-    for (int k = 0; k < num_snapshot_triggers_part; ++k) {
-      if (e->snapshot_recording_triggers_desired_part[k] > 0) {
-        if (e->snapshot_recording_triggers_desired_part[k] >
-            time_to_next_snap) {
-          e->snapshot_recording_triggers_part[k] = time_to_next_snap;
-        } else {
-          e->snapshot_recording_triggers_part[k] =
-              e->snapshot_recording_triggers_desired_part[k];
-        }
-      }
-    }
-    for (int k = 0; k < num_snapshot_triggers_spart; ++k) {
-      if (e->snapshot_recording_triggers_desired_spart[k] > 0) {
-        if (e->snapshot_recording_triggers_desired_spart[k] >
-            time_to_next_snap) {
-          e->snapshot_recording_triggers_spart[k] = time_to_next_snap;
-        } else {
-          e->snapshot_recording_triggers_spart[k] =
-              e->snapshot_recording_triggers_desired_spart[k];
-        }
-      }
-    }
-    for (int k = 0; k < num_snapshot_triggers_bpart; ++k) {
-      if (e->snapshot_recording_triggers_desired_bpart[k] > 0) {
-        if (e->snapshot_recording_triggers_desired_bpart[k] >
-            time_to_next_snap) {
-          e->snapshot_recording_triggers_bpart[k] = time_to_next_snap;
-        } else {
-          e->snapshot_recording_triggers_bpart[k] =
-              e->snapshot_recording_triggers_desired_bpart[k];
-        }
-      }
-    }
+    /* Unless we are restarting, set the recording triggers accordingly for the
+     * next output */
+    if (!restart) engine_set_and_verify_snapshot_triggers(e);
   }
 }
 
@@ -1169,7 +1189,7 @@ void engine_init_output_lists(struct engine *e, struct swift_params *params,
     if (e->output_list_snapshots->select_output_on)
       output_list_check_selection(e->output_list_snapshots, output_options);
 
-    engine_compute_next_snapshot_time(e);
+    engine_compute_next_snapshot_time(e, /*restart=*/0);
 
     if (e->policy & engine_policy_cosmology)
       e->a_first_snapshot =
@@ -1265,12 +1285,15 @@ void engine_io_check_snapshot_triggers(struct engine *e) {
   const int with_cosmology = (e->policy & engine_policy_cosmology);
 
   /* Time until the next snapshot */
+  integertime_t ti_next_snap = e->ti_next_snapshot;
+  if (ti_next_snap == -1) ti_next_snap = max_nr_timesteps;
+
   double time_to_next_snap;
   if (e->policy & engine_policy_cosmology) {
-    time_to_next_snap = cosmology_get_delta_time(e->cosmology, e->ti_current,
-                                                 e->ti_next_snapshot);
+    time_to_next_snap =
+        cosmology_get_delta_time(e->cosmology, e->ti_current, ti_next_snap);
   } else {
-    time_to_next_snap = (e->ti_next_snapshot - e->ti_current) * e->time_base;
+    time_to_next_snap = (ti_next_snap - e->ti_current) * e->time_base;
   }
 
   /* Should any not yet switched on trigger be activated? (part version) */
@@ -1305,10 +1328,10 @@ void engine_io_check_snapshot_triggers(struct engine *e) {
         /* Time from the start of the particle's step to the snapshot */
         double total_time;
         if (with_cosmology) {
-          total_time = cosmology_get_delta_time(e->cosmology, ti_begin,
-                                                e->ti_next_snapshot);
+          total_time =
+              cosmology_get_delta_time(e->cosmology, ti_begin, ti_next_snap);
         } else {
-          total_time = (e->ti_next_snapshot - ti_begin) * e->time_base;
+          total_time = (ti_next_snap - ti_begin) * e->time_base;
         }
 
         /* Time to deduct = time since the start of the step - trigger time */
@@ -1366,10 +1389,10 @@ void engine_io_check_snapshot_triggers(struct engine *e) {
         /* Time from the start of the particle's step to the snapshot */
         double total_time;
         if (with_cosmology) {
-          total_time = cosmology_get_delta_time(e->cosmology, ti_begin,
-                                                e->ti_next_snapshot);
+          total_time =
+              cosmology_get_delta_time(e->cosmology, ti_begin, ti_next_snap);
         } else {
-          total_time = (e->ti_next_snapshot - ti_begin) * e->time_base;
+          total_time = (ti_next_snap - ti_begin) * e->time_base;
         }
 
         /* Time to deduct = time since the start of the step - trigger time */
@@ -1426,10 +1449,10 @@ void engine_io_check_snapshot_triggers(struct engine *e) {
         /* Time from the start of the particle's step to the snapshot */
         double total_time;
         if (with_cosmology) {
-          total_time = cosmology_get_delta_time(e->cosmology, ti_begin,
-                                                e->ti_next_snapshot);
+          total_time =
+              cosmology_get_delta_time(e->cosmology, ti_begin, ti_next_snap);
         } else {
-          total_time = (e->ti_next_snapshot - ti_begin) * e->time_base;
+          total_time = (ti_next_snap - ti_begin) * e->time_base;
         }
 
         /* Time to deduct = time since the start of the step - trigger time */

src/engine_maketasks.c

@@ -43,6 +43,7 @@
 #include "engine.h"
 
 /* Local headers. */
+#include "adaptive_softening.h"
 #include "atomic.h"
 #include "cell.h"
 #include "clocks.h"
@@ -66,10 +67,18 @@ extern int engine_max_parts_per_cooling;
  * @param e The #engine.
  * @param ci The sending #cell.
  * @param cj Dummy cell containing the nodeID of the receiving node.
+ * @param t_pack_grav The grav packing #task, if it has already been created.
  * @param t_grav The send_grav #task, if it has already been created.
+ * @param t_pack_fof The fof packing #task, if it has already been created.
+ * @param t_fof The send_fof #task, if it has already been created.
+ * @param with_fof Are we running with FOF?
  */
 void engine_addtasks_send_gravity(struct engine *e, struct cell *ci,
-                                  struct cell *cj, struct task *t_grav) {
+                                  struct cell *cj, struct task *t_grav_counts,
+                                  struct task *t_grav, struct task *t_pack_grav,
+                                  struct task *t_fof, struct task *t_pack_fof,
+                                  const int with_fof,
+                                  const int with_star_formation) {
 
 #ifdef WITH_MPI
   struct link *l = NULL;
@@ -79,6 +88,19 @@ void engine_addtasks_send_gravity(struct engine *e, struct cell *ci,
   /* Early abort (are we below the level where tasks are)? */
   if (!cell_get_flag(ci, cell_flag_has_tasks)) return;
 
+  if (t_grav_counts == NULL && with_star_formation && ci->hydro.count > 0) {
+#ifdef SWIFT_DEBUG_CHECKS
+    if (ci->depth != 0)
+      error(
+          "Attaching a grav_count task at a non-top level c->depth=%d "
+          "c->count=%d",
+          ci->depth, ci->hydro.count);
+#endif
+    t_grav_counts = scheduler_addtask(
+        s, task_type_send, task_subtype_grav_counts, ci->mpi.tag, 0, ci, cj);
+    scheduler_addunlock(s, ci->hydro.star_formation, t_grav_counts);
+  }
+
   /* Check if any of the gravity tasks are for the target node. */
   for (l = ci->grav.grav; l != NULL; l = l->next)
     if (l->t->ci->nodeID == nodeID ||
@@ -97,22 +119,55 @@ void engine_addtasks_send_gravity(struct engine *e, struct cell *ci,
       t_grav = scheduler_addtask(s, task_type_send, task_subtype_gpart,
                                  ci->mpi.tag, 0, ci, cj);
 
+      t_pack_grav = scheduler_addtask(s, task_type_pack, task_subtype_gpart, 0,
+                                      0, ci, cj);
+
+      /* Pack before you send */
+      scheduler_addunlock(s, t_pack_grav, t_grav);
+
+      if (with_fof) {
+        t_fof = scheduler_addtask(s, task_type_send, task_subtype_fof,
+                                  ci->mpi.tag, 0, ci, cj);
+
+        t_pack_fof = scheduler_addtask(s, task_type_pack, task_subtype_fof, 0,
+                                       0, ci, cj);
+
+        /* Pack before you send */
+        scheduler_addunlock(s, t_pack_fof, t_fof);
+      }
+
       /* The sends should unlock the down pass. */
       scheduler_addunlock(s, t_grav, ci->grav.super->grav.down);
 
-      /* Drift before you send */
-      scheduler_addunlock(s, ci->grav.super->grav.drift, t_grav);
+      if (with_star_formation && ci->top->hydro.count > 0)
+        scheduler_addunlock(s, t_grav, ci->top->hydro.star_formation);
+
+      /* Drift before you pack + send */
+      scheduler_addunlock(s, ci->grav.super->grav.drift, t_pack_grav);
+
+      if (gravity_after_hydro_density)
+        scheduler_addunlock(s, ci->grav.super->grav.init_out, t_pack_grav);
     }
 
     /* Add them to the local cell. */
     engine_addlink(e, &ci->mpi.send, t_grav);
+    engine_addlink(e, &ci->mpi.pack, t_pack_grav);
+    if (with_fof) {
+      engine_addlink(e, &ci->mpi.send, t_fof);
+      engine_addlink(e, &ci->mpi.pack, t_pack_fof);
+    }
+    if (with_star_formation && ci->hydro.count > 0) {
+      engine_addlink(e, &ci->mpi.send, t_grav_counts);
+    }
   }
 
   /* Recurse? */
   if (ci->split)
     for (int k = 0; k < 8; k++)
       if (ci->progeny[k] != NULL)
-        engine_addtasks_send_gravity(e, ci->progeny[k], cj, t_grav);
+        engine_addtasks_send_gravity(e, ci->progeny[k], cj, t_grav_counts,
+                                     t_grav, t_pack_grav, t_fof, t_pack_fof,
+                                     with_fof, with_star_formation);
 
 #else
   error("SWIFT was not compiled with MPI support.");
@@ -338,14 +393,13 @@ void engine_addtasks_send_stars(struct engine *e, struct cell *ci,
                                 struct cell *cj, struct task *t_density,
                                 struct task *t_prep2, struct task *t_sf_counts,
                                 const int with_star_formation) {
-#ifdef SWIFT_DEBUG_CHECKS
+#ifdef WITH_MPI
+#if !defined(SWIFT_DEBUG_CHECKS)
   if (e->policy & engine_policy_sinks && e->policy & engine_policy_stars) {
-    error("TODO");
+    error("TODO: Star formation sink over MPI");
   }
 #endif
 
-#ifdef WITH_MPI
-
   struct link *l = NULL;
   struct scheduler *s = &e->sched;
   const int nodeID = cj->nodeID;
@@ -501,24 +555,24 @@ void engine_addtasks_send_black_holes(struct engine *e, struct cell *ci,
       scheduler_addunlock(s, t_feedback,
                           ci->hydro.super->black_holes.black_holes_out);
 
-      scheduler_addunlock(s, ci->hydro.super->black_holes.swallow_ghost_2,
+      scheduler_addunlock(s, ci->hydro.super->black_holes.swallow_ghost_3,
                           t_feedback);
 
       /* Ghost before you send */
       scheduler_addunlock(s, ci->hydro.super->black_holes.drift, t_rho);
       scheduler_addunlock(s, ci->hydro.super->black_holes.density_ghost, t_rho);
       scheduler_addunlock(s, t_rho,
-                          ci->hydro.super->black_holes.swallow_ghost_0);
+                          ci->hydro.super->black_holes.swallow_ghost_1);
 
-      scheduler_addunlock(s, ci->hydro.super->black_holes.swallow_ghost_0,
+      scheduler_addunlock(s, ci->hydro.super->black_holes.swallow_ghost_1,
                           t_bh_merger);
       scheduler_addunlock(s, t_bh_merger,
-                          ci->hydro.super->black_holes.swallow_ghost_2);
+                          ci->hydro.super->black_holes.swallow_ghost_3);
 
-      scheduler_addunlock(s, ci->hydro.super->black_holes.swallow_ghost_0,
+      scheduler_addunlock(s, ci->hydro.super->black_holes.swallow_ghost_1,
                           t_gas_swallow);
       scheduler_addunlock(s, t_gas_swallow,
-                          ci->hydro.super->black_holes.swallow_ghost_1);
+                          ci->hydro.super->black_holes.swallow_ghost_2);
     }
 
     engine_addlink(e, &ci->mpi.send, t_rho);
@@ -889,13 +943,13 @@ void engine_addtasks_recv_stars(struct engine *e, struct cell *c,
                                 struct task *t_sf_counts,
                                 struct task *const tend,
                                 const int with_star_formation) {
-#ifdef SWIFT_DEBUG_CHECKS
+#ifdef WITH_MPI
+#if !defined(SWIFT_DEBUG_CHECKS)
   if (e->policy & engine_policy_sinks && e->policy & engine_policy_stars) {
-    error("TODO");
+    error("TODO: Star formation sink over MPI");
   }
 #endif
 
-#ifdef WITH_MPI
   struct scheduler *s = &e->sched;
 
   /* Early abort (are we below the level where tasks are)? */
@@ -1105,11 +1159,15 @@ void engine_addtasks_recv_black_holes(struct engine *e, struct cell *c,
  * @param e The #engine.
  * @param c The foreign #cell.
  * @param t_grav The recv_gpart #task, if it has already been created.
+ * @param t_fof The recv_fof #task, if it has already been created.
  * @param tend The top-level time-step communication #task.
+ * @param with_fof Are we running with FOF?
  */
 void engine_addtasks_recv_gravity(struct engine *e, struct cell *c,
-                                  struct task *t_grav,
-                                  struct task *const tend) {
+                                  struct task *t_grav_counts,
+                                  struct task *t_grav, struct task *t_fof,
+                                  struct task *const tend, const int with_fof,
+                                  const int with_star_formation) {
 
 #ifdef WITH_MPI
   struct scheduler *s = &e->sched;
@@ -1117,6 +1175,19 @@ void engine_addtasks_recv_gravity(struct engine *e, struct cell *c,
   /* Early abort (are we below the level where tasks are)? */
   if (!cell_get_flag(c, cell_flag_has_tasks)) return;
 
+  if (t_grav_counts == NULL && with_star_formation && c->hydro.count > 0) {
+#ifdef SWIFT_DEBUG_CHECKS
+    if (c->depth != 0)
+      error(
+          "Attaching a grav_count task at a non-top level c->depth=%d "
+          "c->count=%d",
+          c->depth, c->hydro.count);
+#endif
+
+    t_grav_counts = scheduler_addtask(
+        s, task_type_recv, task_subtype_grav_counts, c->mpi.tag, 0, c, NULL);
+  }
+
   /* Have we reached a level where there are any gravity tasks ? */
   if (t_grav == NULL && c->grav.grav != NULL) {
 
@@ -1128,11 +1199,22 @@ void engine_addtasks_recv_gravity(struct engine *e, struct cell *c,
     /* Create the tasks. */
     t_grav = scheduler_addtask(s, task_type_recv, task_subtype_gpart,
                                c->mpi.tag, 0, c, NULL);
+
+    if (t_grav_counts != NULL) scheduler_addunlock(s, t_grav, t_grav_counts);
+
+    if (with_fof)
+      t_fof = scheduler_addtask(s, task_type_recv, task_subtype_fof, c->mpi.tag,
+                                0, c, NULL);
   }
 
   /* If we have tasks, link them. */
   if (t_grav != NULL) {
     engine_addlink(e, &c->mpi.recv, t_grav);
+    if (with_fof) engine_addlink(e, &c->mpi.recv, t_fof);
+
+    if (with_star_formation && c->hydro.count > 0) {
+      engine_addlink(e, &c->mpi.recv, t_grav_counts);
+    }
 
     for (struct link *l = c->grav.grav; l != NULL; l = l->next) {
       scheduler_addunlock(s, t_grav, l->t);
@@ -1144,7 +1226,9 @@ void engine_addtasks_recv_gravity(struct engine *e, struct cell *c,
   if (c->split)
     for (int k = 0; k < 8; k++)
       if (c->progeny[k] != NULL)
-        engine_addtasks_recv_gravity(e, c->progeny[k], t_grav, tend);
+        engine_addtasks_recv_gravity(e, c->progeny[k], t_grav_counts, t_grav,
+                                     t_fof, tend, with_fof,
+                                     with_star_formation);
 
 #else
   error("SWIFT was not compiled with MPI support.");
@@ -1272,6 +1356,11 @@ void engine_make_hierarchical_tasks_common(struct engine *e, struct cell *c) {
         scheduler_addunlock(s, c->top->sinks.star_formation_sink, c->timestep);
       }
 
+      /* Subgrid tasks: sinks formation */
+      if (with_sinks) {
+        scheduler_addunlock(s, c->kick2, c->top->sinks.sink_formation);
+      }
+
       /* Time-step limiter */
       if (with_timestep_limiter) {
 
@@ -1379,6 +1468,11 @@ void engine_make_hierarchical_tasks_gravity(struct engine *e, struct cell *c) {
         scheduler_addunlock(s, c->grav.long_range, c->grav.down);
         scheduler_addunlock(s, c->grav.down, c->grav.super->grav.end_force);
 
+        /* With adaptive softening, force the hydro density to complete first */
+        if (gravity_after_hydro_density && c->hydro.super == c) {
+          scheduler_addunlock(s, c->hydro.ghost_out, c->grav.init_out);
+        }
+
         /* Link in the implicit tasks */
         scheduler_addunlock(s, c->grav.init, c->grav.init_out);
         scheduler_addunlock(s, c->grav.drift, c->grav.drift_out);
@@ -1499,6 +1593,7 @@ void engine_make_hierarchical_tasks_hydro(struct engine *e, struct cell *c,
   const int with_star_formation_sink = (with_sinks && with_stars);
   const int with_black_holes = (e->policy & engine_policy_black_holes);
   const int with_rt = (e->policy & engine_policy_rt);
+  const int with_timestep_sync = (e->policy & engine_policy_timestep_sync);
 #ifdef WITH_CSDS
   const int with_csds = (e->policy & engine_policy_csds);
 #endif
@@ -1509,22 +1604,10 @@ void engine_make_hierarchical_tasks_hydro(struct engine *e, struct cell *c,
   if ((c->nodeID == e->nodeID) && (star_resort_cell == NULL) &&
       (c->depth == engine_star_resort_task_depth || c->hydro.super == c)) {
 
-    /* Star formation */
-    if (with_feedback && c->hydro.count > 0 && with_star_formation) {
-
-      /* Record this is the level where we re-sort */
-      star_resort_cell = c;
-
-      c->hydro.stars_resort = scheduler_addtask(
-          s, task_type_stars_resort, task_subtype_none, 0, 0, c, NULL);
-
-      scheduler_addunlock(s, c->top->hydro.star_formation,
-                          c->hydro.stars_resort);
-    }
-
-    /* Star formation from sinks */
-    if (with_feedback && with_star_formation_sink &&
-        (c->hydro.count > 0 || c->sinks.count > 0)) {
+    /* If star formation from gas or sinks has happened, we need to resort */
+    if (with_feedback && ((c->hydro.count > 0 && with_star_formation) ||
+                          ((c->hydro.count > 0 || c->sinks.count > 0) &&
+                           with_star_formation_sink))) {
 
       /* Record this is the level where we re-sort */
       star_resort_cell = c;
@@ -1532,8 +1615,20 @@ void engine_make_hierarchical_tasks_hydro(struct engine *e, struct cell *c,
       c->hydro.stars_resort = scheduler_addtask(
           s, task_type_stars_resort, task_subtype_none, 0, 0, c, NULL);
 
-      scheduler_addunlock(s, c->top->sinks.star_formation_sink,
-                          c->hydro.stars_resort);
+      /* Now add the relevant unlocks */
+      /* If we can make stars, we should wait until SF is done before resorting
+       */
+      if (with_star_formation && c->hydro.count > 0) {
+        scheduler_addunlock(s, c->top->hydro.star_formation,
+                            c->hydro.stars_resort);
+      }
+      /* If we can make sinks or spawn from existing ones, we should wait until
+      SF is done before resorting */
+      if (with_star_formation_sink &&
+          (c->hydro.count > 0 || c->sinks.count > 0)) {
+        scheduler_addunlock(s, c->top->sinks.star_formation_sink,
+                            c->hydro.stars_resort);
+      }
     }
   }
 
@@ -1550,7 +1645,7 @@ void engine_make_hierarchical_tasks_hydro(struct engine *e, struct cell *c,
     }
 
     if (with_black_holes) {
-      c->black_holes.swallow_ghost_0 =
+      c->black_holes.swallow_ghost_1 =
           scheduler_addtask(s, task_type_bh_swallow_ghost1, task_subtype_none,
                             0, /* implicit =*/1, c, NULL);
     }
@@ -1601,6 +1696,9 @@ void engine_make_hierarchical_tasks_hydro(struct engine *e, struct cell *c,
             scheduler_addtask(s, task_type_sink_in, task_subtype_none, 0,
                               /* implicit = */ 1, c, NULL);
 
+        c->sinks.density_ghost = scheduler_addtask(
+            s, task_type_sink_density_ghost, task_subtype_none, 0, 0, c, NULL);
+
         c->sinks.sink_ghost1 =
             scheduler_addtask(s, task_type_sink_ghost1, task_subtype_none, 0,
                               /* implicit = */ 1, c, NULL);
@@ -1616,6 +1714,9 @@ void engine_make_hierarchical_tasks_hydro(struct engine *e, struct cell *c,
         /* Link to the main tasks */
         scheduler_addunlock(s, c->super->kick2, c->sinks.sink_in);
         scheduler_addunlock(s, c->sinks.sink_out, c->super->timestep);
+        scheduler_addunlock(s, c->top->sinks.sink_formation, c->sinks.sink_in);
+        if (with_timestep_sync)
+          scheduler_addunlock(s, c->sinks.sink_out, c->super->timestep_sync);
 
         if (with_stars &&
             (c->top->hydro.count > 0 || c->top->sinks.count > 0)) {
@@ -1690,13 +1791,9 @@ void engine_make_hierarchical_tasks_hydro(struct engine *e, struct cell *c,
         scheduler_addunlock(s, c->stars.stars_out, c->super->timestep);
 
         /* Star formation*/
-        if (with_feedback && c->hydro.count > 0 && with_star_formation) {
-          scheduler_addunlock(s, star_resort_cell->hydro.stars_resort,
-                              c->stars.stars_in);
-        }
-        /* Star formation from sinks */
-        if (with_feedback && with_star_formation_sink &&
-            (c->hydro.count > 0 || c->sinks.count > 0)) {
+        if (with_feedback && ((c->hydro.count > 0 && with_star_formation) ||
+                              ((c->hydro.count > 0 || c->sinks.count > 0) &&
+                               with_star_formation_sink))) {
           scheduler_addunlock(s, star_resort_cell->hydro.stars_resort,
                               c->stars.stars_in);
         }
@@ -1785,11 +1882,11 @@ void engine_make_hierarchical_tasks_hydro(struct engine *e, struct cell *c,
         c->black_holes.density_ghost = scheduler_addtask(
             s, task_type_bh_density_ghost, task_subtype_none, 0, 0, c, NULL);
 
-        c->black_holes.swallow_ghost_1 =
+        c->black_holes.swallow_ghost_2 =
             scheduler_addtask(s, task_type_bh_swallow_ghost2, task_subtype_none,
                               0, /* implicit =*/1, c, NULL);
 
-        c->black_holes.swallow_ghost_2 = scheduler_addtask(
+        c->black_holes.swallow_ghost_3 = scheduler_addtask(
             s, task_type_bh_swallow_ghost3, task_subtype_none, 0, 0, c, NULL);
 
 #ifdef WITH_CSDS
@@ -1811,7 +1908,7 @@ void engine_make_hierarchical_tasks_hydro(struct engine *e, struct cell *c,
         /* Make sure we don't start swallowing gas particles before the stars
            have converged on their smoothing lengths. */
         scheduler_addunlock(s, c->stars.density_ghost,
-                            c->black_holes.swallow_ghost_0);
+                            c->black_holes.swallow_ghost_1);
       }
     }
   } else { /* We are above the super-cell so need to go deeper */
@@ -2098,8 +2195,9 @@ void engine_count_and_link_tasks_mapper(void *map_data, int num_elements,
 
       /* Link self tasks to cells. */
     } else if (t_type == task_type_self) {
+#ifdef SWIFT_DEBUG_CHECKS
       atomic_inc(&ci->nr_tasks);
-
+#endif
       if (t_subtype == task_subtype_density) {
         engine_addlink(e, &ci->hydro.density, t);
       } else if (t_subtype == task_subtype_grav) {
@@ -2110,38 +2208,10 @@ void engine_count_and_link_tasks_mapper(void *map_data, int num_elements,
 
       /* Link pair tasks to cells. */
     } else if (t_type == task_type_pair) {
-      atomic_inc(&ci->nr_tasks);
-      atomic_inc(&cj->nr_tasks);
-
-      if (t_subtype == task_subtype_density) {
-        engine_addlink(e, &ci->hydro.density, t);
-        engine_addlink(e, &cj->hydro.density, t);
-      } else if (t_subtype == task_subtype_grav) {
-        engine_addlink(e, &ci->grav.grav, t);
-        engine_addlink(e, &cj->grav.grav, t);
-      }
 #ifdef SWIFT_DEBUG_CHECKS
-      else if (t_subtype == task_subtype_external_grav) {
-        error("Found a pair/external-gravity task...");
-      }
-#endif
-
-      /* Link sub-self tasks to cells. */
-    } else if (t_type == task_type_sub_self) {
-      atomic_inc(&ci->nr_tasks);
-
-      if (t_subtype == task_subtype_density) {
-        engine_addlink(e, &ci->hydro.density, t);
-      } else if (t_subtype == task_subtype_grav) {
-        engine_addlink(e, &ci->grav.grav, t);
-      } else if (t_subtype == task_subtype_external_grav) {
-        engine_addlink(e, &ci->grav.grav, t);
-      }
-
-      /* Link sub-pair tasks to cells. */
-    } else if (t_type == task_type_sub_pair) {
       atomic_inc(&ci->nr_tasks);
       atomic_inc(&cj->nr_tasks);
+#endif
 
       if (t_subtype == task_subtype_density) {
         engine_addlink(e, &ci->hydro.density, t);
@@ -2170,24 +2240,28 @@ void engine_count_and_link_tasks_mapper(void *map_data, int num_elements,
 /**
  * @brief Creates all the task dependencies for the gravity
  *
- * @param e The #engine
+ * @param map_data The task array passed to this pool thread.
+ * @param num_elements The number of tasks in this pool thread.
+ * @param extra_data Pointer to the #engine.
  */
-void engine_link_gravity_tasks(struct engine *e) {
+void engine_link_gravity_tasks_mapper(void *map_data, int num_elements,
+                                      void *extra_data) {
 
+  struct task *tasks = (struct task *)map_data;
+  struct engine *e = (struct engine *)extra_data;
   struct scheduler *sched = &e->sched;
   const int nodeID = e->nodeID;
-  const int nr_tasks = sched->nr_tasks;
 
-  for (int k = 0; k < nr_tasks; k++) {
+  for (int k = 0; k < num_elements; k++) {
 
     /* Get a pointer to the task. */
-    struct task *t = &sched->tasks[k];
+    struct task *t = &tasks[k];
 
     if (t->type == task_type_none) continue;
 
     /* Get the cells we act on */
-    struct cell *ci = t->ci;
-    struct cell *cj = t->cj;
+    struct cell *restrict ci = t->ci;
+    struct cell *restrict cj = t->cj;
     const enum task_types t_type = t->type;
     const enum task_subtypes t_subtype = t->subtype;
 
@@ -2229,7 +2303,8 @@ void engine_link_gravity_tasks(struct engine *e) {
     }
 
     /* Self-interaction for external gravity ? */
-    if (t_type == task_type_self && t_subtype == task_subtype_external_grav) {
+    else if (t_type == task_type_self &&
+             t_subtype == task_subtype_external_grav) {
 
 #ifdef SWIFT_DEBUG_CHECKS
       if (ci_nodeID != nodeID) error("Non-local self task");
@@ -2264,7 +2339,7 @@ void engine_link_gravity_tasks(struct engine *e) {
     }
 
     /* Otherwise, sub-self interaction? */
-    else if (t_type == task_type_sub_self && t_subtype == task_subtype_grav) {
+    else if (t_type == task_type_self && t_subtype == task_subtype_grav) {
 
 #ifdef SWIFT_DEBUG_CHECKS
       if (ci_nodeID != nodeID) error("Non-local sub-self task");
@@ -2277,7 +2352,7 @@ void engine_link_gravity_tasks(struct engine *e) {
     }
 
     /* Sub-self-interaction for external gravity ? */
-    else if (t_type == task_type_sub_self &&
+    else if (t_type == task_type_self &&
              t_subtype == task_subtype_external_grav) {
 
 #ifdef SWIFT_DEBUG_CHECKS
@@ -2290,7 +2365,7 @@ void engine_link_gravity_tasks(struct engine *e) {
     }
 
     /* Otherwise, sub-pair interaction? */
-    else if (t_type == task_type_sub_pair && t_subtype == task_subtype_grav) {
+    else if (t_type == task_type_pair && t_subtype == task_subtype_grav) {
 
       if (ci_nodeID == nodeID) {
 
@@ -2425,6 +2500,7 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
   struct task *t_do_gas_swallow = NULL;
   struct task *t_do_bh_swallow = NULL;
   struct task *t_bh_feedback = NULL;
+  struct task *t_sink_density = NULL;
   struct task *t_sink_swallow = NULL;
   struct task *t_rt_gradient = NULL;
   struct task *t_rt_transport = NULL;
@@ -2457,19 +2533,20 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
       scheduler_addunlock(sched, ci->hydro.super->stars.drift, t);
     }
 
-    /* Self-interaction? */
+    /* Otherwise, self interaction? */
     else if (t_type == task_type_self && t_subtype == task_subtype_density) {
 
       const int bcount_i = ci->black_holes.count;
 
-      /* Make the self-density tasks depend on the drift only. */
+      /* Make all density tasks depend on the drift and sorts. */
       scheduler_addunlock(sched, ci->hydro.super->hydro.drift, t);
+      scheduler_addunlock(sched, ci->hydro.super->hydro.sorts, t);
 
-      /* Task for the second hydro loop, */
+      /* Start by constructing the task for the second hydro loop */
       t_force = scheduler_addtask(sched, task_type_self, task_subtype_force,
                                   flags, 0, ci, NULL);
 
-      /* the task for the time-step limiter */
+      /* and the task for the time-step limiter */
       if (with_timestep_limiter) {
         t_limiter = scheduler_addtask(sched, task_type_self,
                                       task_subtype_limiter, flags, 0, ci, NULL);
@@ -2496,6 +2573,9 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
 
       /* The sink tasks */
       if (with_sink) {
+        t_sink_density =
+            scheduler_addtask(sched, task_type_self, task_subtype_sink_density,
+                              flags, 0, ci, NULL);
         t_sink_swallow =
             scheduler_addtask(sched, task_type_self, task_subtype_sink_swallow,
                               flags, 0, ci, NULL);
@@ -2513,12 +2593,15 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
             sched, task_type_self, task_subtype_bh_density, flags, 0, ci, NULL);
         t_bh_swallow = scheduler_addtask(
             sched, task_type_self, task_subtype_bh_swallow, flags, 0, ci, NULL);
+
         t_do_gas_swallow =
             scheduler_addtask(sched, task_type_self,
                               task_subtype_do_gas_swallow, flags, 0, ci, NULL);
+
         t_do_bh_swallow =
             scheduler_addtask(sched, task_type_self, task_subtype_do_bh_swallow,
                               flags, 0, ci, NULL);
+
         t_bh_feedback =
             scheduler_addtask(sched, task_type_self, task_subtype_bh_feedback,
                               flags, 0, ci, NULL);
@@ -2533,7 +2616,7 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
                               flags, 0, ci, NULL);
       }
 
-      /* Link the tasks to the cells */
+      /* Add the link between the new loop and the cell */
       engine_addlink(e, &ci->hydro.force, t_force);
       if (with_timestep_limiter) {
         engine_addlink(e, &ci->hydro.limiter, t_limiter);
@@ -2547,6 +2630,7 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
 #endif
       }
       if (with_sink) {
+        engine_addlink(e, &ci->sinks.density, t_sink_density);
         engine_addlink(e, &ci->sinks.swallow, t_sink_swallow);
         engine_addlink(e, &ci->sinks.do_sink_swallow, t_sink_do_sink_swallow);
         engine_addlink(e, &ci->sinks.do_gas_swallow, t_sink_do_gas_swallow);
@@ -2565,20 +2649,22 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
 
 #ifdef EXTRA_HYDRO_LOOP
 
-      /* Same work for the additional hydro loop */
+      /* Start by constructing the task for the second and third hydro loop */
       t_gradient = scheduler_addtask(sched, task_type_self,
                                      task_subtype_gradient, flags, 0, ci, NULL);
 
-      /* Add the link between the new loops and the cell */
+      /* Add the link between the new loop and the cell */
       engine_addlink(e, &ci->hydro.gradient, t_gradient);
 
-      /* Now, build all the dependencies for the hydro */
+      /* Now, build all the dependencies for the hydro for the cells */
+      /* that are local and are not descendant of the same super_hydro-cells */
       engine_make_hydro_loops_dependencies(sched, t, t_gradient, t_force,
                                            t_limiter, ci, with_cooling,
                                            with_timestep_limiter);
 #else
 
-      /* Now, build all the dependencies for the hydro */
+      /* Now, build all the dependencies for the hydro for the cells */
+      /* that are local and are not descendant of the same super_hydro-cells */
       engine_make_hydro_loops_dependencies(sched, t, t_force, t_limiter, ci,
                                            with_cooling, with_timestep_limiter);
 #endif
@@ -2594,8 +2680,12 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
 
         scheduler_addunlock(sched, ci->hydro.super->stars.drift,
                             t_star_density);
+        scheduler_addunlock(sched, ci->hydro.super->stars.sorts,
+                            t_star_density);
         scheduler_addunlock(sched, ci->hydro.super->hydro.drift,
                             t_star_density);
+        scheduler_addunlock(sched, ci->hydro.super->hydro.sorts,
+                            t_star_density);
         scheduler_addunlock(sched, ci->hydro.super->stars.stars_in,
                             t_star_density);
         scheduler_addunlock(sched, t_star_density,
@@ -2623,20 +2713,21 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
                             ci->hydro.super->stars.stars_out);
       }
 
-      /* The sink's tasks. */
       if (with_sink) {
 
-        /* Do the sink_formation */
+        /* Sink density */
         scheduler_addunlock(sched, ci->hydro.super->sinks.drift,
-                            ci->top->sinks.sink_formation);
+                            t_sink_density);
         scheduler_addunlock(sched, ci->hydro.super->hydro.drift,
-                            ci->top->sinks.sink_formation);
+                            t_sink_density);
         scheduler_addunlock(sched, ci->hydro.super->sinks.sink_in,
-                            ci->top->sinks.sink_formation);
-        scheduler_addunlock(sched, ci->top->sinks.sink_formation,
-                            t_sink_swallow);
+                            t_sink_density);
+        scheduler_addunlock(sched, t_sink_density,
+                            ci->hydro.super->sinks.density_ghost);
 
         /* Do the sink_swallow */
+        scheduler_addunlock(sched, ci->hydro.super->sinks.density_ghost,
+                            t_sink_swallow);
         scheduler_addunlock(sched, t_sink_swallow,
                             ci->hydro.super->sinks.sink_ghost1);
 
@@ -2670,27 +2761,27 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
         scheduler_addunlock(sched, ci->hydro.super->black_holes.density_ghost,
                             t_bh_swallow);
         scheduler_addunlock(sched, t_bh_swallow,
-                            ci->hydro.super->black_holes.swallow_ghost_0);
+                            ci->hydro.super->black_holes.swallow_ghost_1);
 
-        scheduler_addunlock(sched, ci->hydro.super->black_holes.swallow_ghost_0,
+        scheduler_addunlock(sched, ci->hydro.super->black_holes.swallow_ghost_1,
                             t_do_gas_swallow);
         scheduler_addunlock(sched, t_do_gas_swallow,
-                            ci->hydro.super->black_holes.swallow_ghost_1);
+                            ci->hydro.super->black_holes.swallow_ghost_2);
 
-        scheduler_addunlock(sched, ci->hydro.super->black_holes.swallow_ghost_1,
+        scheduler_addunlock(sched, ci->hydro.super->black_holes.swallow_ghost_2,
                             t_do_bh_swallow);
         scheduler_addunlock(sched, t_do_bh_swallow,
-                            ci->hydro.super->black_holes.swallow_ghost_2);
+                            ci->hydro.super->black_holes.swallow_ghost_3);
 
-        scheduler_addunlock(sched, ci->hydro.super->black_holes.swallow_ghost_2,
+        scheduler_addunlock(sched, ci->hydro.super->black_holes.swallow_ghost_3,
                             t_bh_feedback);
         scheduler_addunlock(sched, t_bh_feedback,
                             ci->hydro.super->black_holes.black_holes_out);
       }
 
       if (with_timestep_limiter) {
-        scheduler_addunlock(sched, ci->super->timestep, t_limiter);
         scheduler_addunlock(sched, ci->hydro.super->hydro.drift, t_limiter);
+        scheduler_addunlock(sched, ci->super->timestep, t_limiter);
         scheduler_addunlock(sched, t_limiter, ci->super->kick1);
         scheduler_addunlock(sched, t_limiter, ci->super->timestep_limiter);
       }
@@ -2704,6 +2795,7 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
 
       if (with_rt) {
         scheduler_addunlock(sched, ci->hydro.super->hydro.drift, t_rt_gradient);
+        scheduler_addunlock(sched, ci->hydro.super->hydro.sorts, t_rt_gradient);
         scheduler_addunlock(sched, ci->hydro.super->rt.rt_ghost1,
                             t_rt_gradient);
         scheduler_addunlock(sched, t_rt_gradient,
@@ -2766,6 +2858,7 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
         t_star_feedback =
             scheduler_addtask(sched, task_type_pair,
                               task_subtype_stars_feedback, flags, 0, ci, cj);
+
 #ifdef EXTRA_STAR_LOOPS
         t_star_prep1 = scheduler_addtask(
             sched, task_type_pair, task_subtype_stars_prep1, flags, 0, ci, cj);
@@ -2776,6 +2869,8 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
 
       /* The sink tasks */
       if (with_sink) {
+        t_sink_density = scheduler_addtask(
+            sched, task_type_pair, task_subtype_sink_density, flags, 0, ci, cj);
         t_sink_swallow = scheduler_addtask(
             sched, task_type_pair, task_subtype_sink_swallow, flags, 0, ci, cj);
         t_sink_do_sink_swallow = scheduler_addtask(
@@ -2823,8 +2918,8 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
                                 t_rt_transport);
         }
         /* We need to ensure that a local inactive cell is sorted before
-         * the interaction in the transport loop. Local cells don't have an
-         * rt_sorts task. */
+         * the interaction in the transport loop. Local cells don't have
+         * an rt_sort task. */
         if (ci->hydro.super->hydro.sorts != NULL)
           scheduler_addunlock(sched, ci->hydro.super->hydro.sorts,
                               t_rt_transport);
@@ -2854,13 +2949,12 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
 #endif
       }
       if (with_sink) {
-        /* Formation */
+        engine_addlink(e, &ci->sinks.density, t_sink_density);
+        engine_addlink(e, &cj->sinks.density, t_sink_density);
         engine_addlink(e, &ci->sinks.swallow, t_sink_swallow);
         engine_addlink(e, &cj->sinks.swallow, t_sink_swallow);
-        /* Merger */
         engine_addlink(e, &ci->sinks.do_sink_swallow, t_sink_do_sink_swallow);
         engine_addlink(e, &cj->sinks.do_sink_swallow, t_sink_do_sink_swallow);
-        /* Accretion */
         engine_addlink(e, &ci->sinks.do_gas_swallow, t_sink_do_gas_swallow);
         engine_addlink(e, &cj->sinks.do_gas_swallow, t_sink_do_gas_swallow);
       }
@@ -2924,15 +3018,14 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
       if (with_feedback) {
         scheduler_addunlock(sched, ci->hydro.super->hydro.sorts,
                             t_star_density);
-
         if (ci->hydro.super != cj->hydro.super) {
           scheduler_addunlock(sched, cj->hydro.super->hydro.sorts,
                               t_star_density);
         }
       }
+
       if (with_rt) {
         scheduler_addunlock(sched, ci->hydro.super->hydro.sorts, t_rt_gradient);
-
         if (ci->hydro.super != cj->hydro.super) {
           scheduler_addunlock(sched, cj->hydro.super->hydro.sorts,
                               t_rt_gradient);
@@ -2948,10 +3041,10 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
             scheduler_addunlock(sched, ci->hydro.super->hydro.cooling_out,
                                 t_star_density);
 
-          scheduler_addunlock(sched, ci->hydro.super->stars.drift,
-                              t_star_density);
           scheduler_addunlock(sched, ci->hydro.super->stars.sorts,
                               t_star_density);
+          scheduler_addunlock(sched, ci->hydro.super->stars.drift,
+                              t_star_density);
           scheduler_addunlock(sched, ci->hydro.super->hydro.drift,
                               t_star_density);
           scheduler_addunlock(sched, ci->hydro.super->stars.stars_in,
@@ -2983,17 +3076,19 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
 
         if (with_sink) {
 
-          /* Do the sink_formation */
+          /* Sink density */
           scheduler_addunlock(sched, ci->hydro.super->sinks.drift,
-                              ci->top->sinks.sink_formation);
+                              t_sink_density);
           scheduler_addunlock(sched, ci->hydro.super->hydro.drift,
-                              ci->top->sinks.sink_formation);
+                              t_sink_density);
           scheduler_addunlock(sched, ci->hydro.super->sinks.sink_in,
-                              ci->top->sinks.sink_formation);
-          scheduler_addunlock(sched, ci->top->sinks.sink_formation,
-                              t_sink_swallow);
+                              t_sink_density);
+          scheduler_addunlock(sched, t_sink_density,
+                              ci->hydro.super->sinks.density_ghost);
 
           /* Do the sink_swallow */
+          scheduler_addunlock(sched, ci->hydro.super->sinks.density_ghost,
+                              t_sink_swallow);
           scheduler_addunlock(sched, t_sink_swallow,
                               ci->hydro.super->sinks.sink_ghost1);
 
@@ -3028,22 +3123,22 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
           scheduler_addunlock(sched, ci->hydro.super->black_holes.density_ghost,
                               t_bh_swallow);
           scheduler_addunlock(sched, t_bh_swallow,
-                              ci->hydro.super->black_holes.swallow_ghost_0);
+                              ci->hydro.super->black_holes.swallow_ghost_1);
 
           scheduler_addunlock(sched,
-                              ci->hydro.super->black_holes.swallow_ghost_0,
+                              ci->hydro.super->black_holes.swallow_ghost_1,
                               t_do_gas_swallow);
           scheduler_addunlock(sched, t_do_gas_swallow,
-                              ci->hydro.super->black_holes.swallow_ghost_1);
+                              ci->hydro.super->black_holes.swallow_ghost_2);
 
           scheduler_addunlock(sched,
-                              ci->hydro.super->black_holes.swallow_ghost_1,
+                              ci->hydro.super->black_holes.swallow_ghost_2,
                               t_do_bh_swallow);
           scheduler_addunlock(sched, t_do_bh_swallow,
-                              ci->hydro.super->black_holes.swallow_ghost_2);
+                              ci->hydro.super->black_holes.swallow_ghost_3);
 
           scheduler_addunlock(sched,
-                              ci->hydro.super->black_holes.swallow_ghost_2,
+                              ci->hydro.super->black_holes.swallow_ghost_3,
                               t_bh_feedback);
           scheduler_addunlock(sched, t_bh_feedback,
                               ci->hydro.super->black_holes.black_holes_out);
@@ -3076,8 +3171,8 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
           scheduler_addunlock(sched, t_rt_transport,
                               ci->hydro.super->rt.rt_transport_out);
         }
+      } else /* ci->nodeID != nodeID */ {
 
-      } else /*(ci->nodeID != nodeID) */ {
         if (with_feedback) {
 #ifdef EXTRA_STAR_LOOPS
           scheduler_addunlock(sched, ci->hydro.super->stars.sorts,
@@ -3086,10 +3181,10 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
           scheduler_addunlock(sched, ci->hydro.super->stars.sorts,
                               t_star_feedback);
         }
-
         if (with_black_holes && (bcount_i > 0 || bcount_j > 0)) {
+
           scheduler_addunlock(sched, t_bh_swallow,
-                              ci->hydro.super->black_holes.swallow_ghost_0);
+                              ci->hydro.super->black_holes.swallow_ghost_1);
         }
       }
 
@@ -3140,17 +3235,19 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
 
           if (with_sink) {
 
-            /* Do the sink_formation */
+            /* Sink density */
             scheduler_addunlock(sched, cj->hydro.super->sinks.drift,
-                                cj->top->sinks.sink_formation);
+                                t_sink_density);
             scheduler_addunlock(sched, cj->hydro.super->hydro.drift,
-                                cj->top->sinks.sink_formation);
+                                t_sink_density);
             scheduler_addunlock(sched, cj->hydro.super->sinks.sink_in,
-                                cj->top->sinks.sink_formation);
-            scheduler_addunlock(sched, cj->top->sinks.sink_formation,
-                                t_sink_swallow);
+                                t_sink_density);
+            scheduler_addunlock(sched, t_sink_density,
+                                cj->hydro.super->sinks.density_ghost);
 
             /* Do the sink_swallow */
+            scheduler_addunlock(sched, cj->hydro.super->sinks.density_ghost,
+                                t_sink_swallow);
             scheduler_addunlock(sched, t_sink_swallow,
                                 cj->hydro.super->sinks.sink_ghost1);
 
@@ -3187,27 +3284,26 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
                                 cj->hydro.super->black_holes.density_ghost,
                                 t_bh_swallow);
             scheduler_addunlock(sched, t_bh_swallow,
-                                cj->hydro.super->black_holes.swallow_ghost_0);
+                                cj->hydro.super->black_holes.swallow_ghost_1);
 
             scheduler_addunlock(sched,
-                                cj->hydro.super->black_holes.swallow_ghost_0,
+                                cj->hydro.super->black_holes.swallow_ghost_1,
                                 t_do_gas_swallow);
             scheduler_addunlock(sched, t_do_gas_swallow,
-                                cj->hydro.super->black_holes.swallow_ghost_1);
+                                cj->hydro.super->black_holes.swallow_ghost_2);
 
             scheduler_addunlock(sched,
-                                cj->hydro.super->black_holes.swallow_ghost_1,
+                                cj->hydro.super->black_holes.swallow_ghost_2,
                                 t_do_bh_swallow);
             scheduler_addunlock(sched, t_do_bh_swallow,
-                                cj->hydro.super->black_holes.swallow_ghost_2);
+                                cj->hydro.super->black_holes.swallow_ghost_3);
 
             scheduler_addunlock(sched,
-                                cj->hydro.super->black_holes.swallow_ghost_2,
+                                cj->hydro.super->black_holes.swallow_ghost_3,
                                 t_bh_feedback);
             scheduler_addunlock(sched, t_bh_feedback,
                                 cj->hydro.super->black_holes.black_holes_out);
           }
-
           if (with_rt) {
             scheduler_addunlock(sched, cj->hydro.super->hydro.drift,
                                 t_rt_gradient);
@@ -3243,8 +3339,7 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
             scheduler_addunlock(sched, t_bh_feedback, cj->super->timestep_sync);
           }
         }
-
-      } else /*(cj->nodeID != nodeID) */ {
+      } else /* cj->nodeID != nodeID */ {
         if (with_feedback) {
 #ifdef EXTRA_STAR_LOOPS
           scheduler_addunlock(sched, cj->hydro.super->stars.sorts,
@@ -3255,895 +3350,69 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
         }
 
         if (with_black_holes && (bcount_i > 0 || bcount_j > 0)) {
-
           scheduler_addunlock(sched, t_bh_swallow,
-                              cj->hydro.super->black_holes.swallow_ghost_0);
+                              cj->hydro.super->black_holes.swallow_ghost_1);
         }
       }
     }
+  }
+}
 
-    /* Otherwise, sub-self interaction? */
-    else if (t_type == task_type_sub_self &&
-             t_subtype == task_subtype_density) {
+/**
+ * @brief Constructs the top-level pair tasks for the first hydro loop over
+ * neighbours
+ *
+ * Here we construct all the tasks for all possible neighbouring non-empty
+ * local cells in the hierarchy. No dependencies are being added thus far.
+ * Additional loop over neighbours can later be added by simply duplicating
+ * all the tasks created by this function.
+ *
+ * @param map_data Offset of first two indices disguised as a pointer.
+ * @param num_elements Number of cells to traverse.
+ * @param extra_data The #engine.
+ */
+void engine_make_hydroloop_tasks_mapper(void *map_data, int num_elements,
+                                        void *extra_data) {
 
-      const int bcount_i = ci->black_holes.count;
+  /* Extract the engine pointer. */
+  struct engine *e = (struct engine *)extra_data;
+  const int periodic = e->s->periodic;
+  const int with_feedback = (e->policy & engine_policy_feedback);
+  const int with_stars = (e->policy & engine_policy_stars);
+  const int with_sinks = (e->policy & engine_policy_sinks);
+  const int with_black_holes = (e->policy & engine_policy_black_holes);
 
-      /* Make all density tasks depend on the drift and sorts. */
-      scheduler_addunlock(sched, ci->hydro.super->hydro.drift, t);
-      scheduler_addunlock(sched, ci->hydro.super->hydro.sorts, t);
+  struct space *s = e->s;
+  struct scheduler *sched = &e->sched;
+  const int nodeID = e->nodeID;
+  const int *cdim = s->cdim;
+  struct cell *cells = s->cells_top;
 
-      /* Start by constructing the task for the second hydro loop */
-      t_force = scheduler_addtask(sched, task_type_sub_self, task_subtype_force,
-                                  flags, 0, ci, NULL);
+  /* Loop through the elements, which are just byte offsets from NULL. */
+  for (int ind = 0; ind < num_elements; ind++) {
 
-      /* and the task for the time-step limiter */
-      if (with_timestep_limiter) {
-        t_limiter = scheduler_addtask(sched, task_type_sub_self,
-                                      task_subtype_limiter, flags, 0, ci, NULL);
-      }
+    /* Get the cell index. */
+    const int cid = (size_t)(map_data) + ind;
 
-      /* The stellar feedback tasks */
-      if (with_feedback) {
-        t_star_density =
-            scheduler_addtask(sched, task_type_sub_self,
-                              task_subtype_stars_density, flags, 0, ci, NULL);
-        t_star_feedback =
-            scheduler_addtask(sched, task_type_sub_self,
-                              task_subtype_stars_feedback, flags, 0, ci, NULL);
+    /* 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];
 
-#ifdef EXTRA_STAR_LOOPS
-        t_star_prep1 =
-            scheduler_addtask(sched, task_type_sub_self,
-                              task_subtype_stars_prep1, flags, 0, ci, NULL);
-        t_star_prep2 =
-            scheduler_addtask(sched, task_type_sub_self,
-                              task_subtype_stars_prep2, flags, 0, ci, NULL);
-#endif
-      }
+    /* Get the cell */
+    struct cell *ci = &cells[cid];
 
-      /* The sink tasks */
-      if (with_sink) {
-        t_sink_swallow =
-            scheduler_addtask(sched, task_type_sub_self,
-                              task_subtype_sink_swallow, flags, 0, ci, NULL);
-        t_sink_do_sink_swallow = scheduler_addtask(
-            sched, task_type_sub_self, task_subtype_sink_do_sink_swallow, flags,
-            0, ci, NULL);
-        t_sink_do_gas_swallow = scheduler_addtask(
-            sched, task_type_sub_self, task_subtype_sink_do_gas_swallow, flags,
-            0, ci, NULL);
-      }
+    /* Skip cells without hydro or star particles */
+    if ((ci->hydro.count == 0) && (!with_stars || ci->stars.count == 0) &&
+        (!with_sinks || ci->sinks.count == 0) &&
+        (!with_black_holes || ci->black_holes.count == 0))
+      continue;
 
-      /* The black hole feedback tasks */
-      if (with_black_holes && bcount_i > 0) {
-        t_bh_density =
-            scheduler_addtask(sched, task_type_sub_self,
-                              task_subtype_bh_density, flags, 0, ci, NULL);
-        t_bh_swallow =
-            scheduler_addtask(sched, task_type_sub_self,
-                              task_subtype_bh_swallow, flags, 0, ci, NULL);
-
-        t_do_gas_swallow =
-            scheduler_addtask(sched, task_type_sub_self,
-                              task_subtype_do_gas_swallow, flags, 0, ci, NULL);
-
-        t_do_bh_swallow =
-            scheduler_addtask(sched, task_type_sub_self,
-                              task_subtype_do_bh_swallow, flags, 0, ci, NULL);
-
-        t_bh_feedback =
-            scheduler_addtask(sched, task_type_sub_self,
-                              task_subtype_bh_feedback, flags, 0, ci, NULL);
-      }
-
-      if (with_rt) {
-        t_rt_gradient =
-            scheduler_addtask(sched, task_type_sub_self,
-                              task_subtype_rt_gradient, flags, 0, ci, NULL);
-        t_rt_transport =
-            scheduler_addtask(sched, task_type_sub_self,
-                              task_subtype_rt_transport, flags, 0, ci, NULL);
-      }
-
-      /* Add the link between the new loop and the cell */
-      engine_addlink(e, &ci->hydro.force, t_force);
-      if (with_timestep_limiter) {
-        engine_addlink(e, &ci->hydro.limiter, t_limiter);
-      }
-      if (with_feedback) {
-        engine_addlink(e, &ci->stars.density, t_star_density);
-        engine_addlink(e, &ci->stars.feedback, t_star_feedback);
-#ifdef EXTRA_STAR_LOOPS
-        engine_addlink(e, &ci->stars.prepare1, t_star_prep1);
-        engine_addlink(e, &ci->stars.prepare2, t_star_prep2);
-#endif
-      }
-      if (with_sink) {
-        engine_addlink(e, &ci->sinks.swallow, t_sink_swallow);
-        engine_addlink(e, &ci->sinks.do_sink_swallow, t_sink_do_sink_swallow);
-        engine_addlink(e, &ci->sinks.do_gas_swallow, t_sink_do_gas_swallow);
-      }
-      if (with_black_holes && bcount_i > 0) {
-        engine_addlink(e, &ci->black_holes.density, t_bh_density);
-        engine_addlink(e, &ci->black_holes.swallow, t_bh_swallow);
-        engine_addlink(e, &ci->black_holes.do_gas_swallow, t_do_gas_swallow);
-        engine_addlink(e, &ci->black_holes.do_bh_swallow, t_do_bh_swallow);
-        engine_addlink(e, &ci->black_holes.feedback, t_bh_feedback);
-      }
-      if (with_rt) {
-        engine_addlink(e, &ci->rt.rt_gradient, t_rt_gradient);
-        engine_addlink(e, &ci->rt.rt_transport, t_rt_transport);
-      }
-
-#ifdef EXTRA_HYDRO_LOOP
-
-      /* Start by constructing the task for the second and third hydro loop */
-      t_gradient = scheduler_addtask(sched, task_type_sub_self,
-                                     task_subtype_gradient, flags, 0, ci, NULL);
-
-      /* Add the link between the new loop and the cell */
-      engine_addlink(e, &ci->hydro.gradient, t_gradient);
-
-      /* Now, build all the dependencies for the hydro for the cells */
-      /* that are local and are not descendant of the same super_hydro-cells */
-      engine_make_hydro_loops_dependencies(sched, t, t_gradient, t_force,
-                                           t_limiter, ci, with_cooling,
-                                           with_timestep_limiter);
-#else
-
-      /* Now, build all the dependencies for the hydro for the cells */
-      /* that are local and are not descendant of the same super_hydro-cells */
-      engine_make_hydro_loops_dependencies(sched, t, t_force, t_limiter, ci,
-                                           with_cooling, with_timestep_limiter);
-#endif
-
-      /* Create the task dependencies */
-      scheduler_addunlock(sched, t_force, ci->hydro.super->hydro.end_force);
-
-      if (with_feedback) {
-
-        if (with_cooling)
-          scheduler_addunlock(sched, ci->hydro.super->hydro.cooling_out,
-                              t_star_density);
-
-        scheduler_addunlock(sched, ci->hydro.super->stars.drift,
-                            t_star_density);
-        scheduler_addunlock(sched, ci->hydro.super->stars.sorts,
-                            t_star_density);
-        scheduler_addunlock(sched, ci->hydro.super->hydro.drift,
-                            t_star_density);
-        scheduler_addunlock(sched, ci->hydro.super->hydro.sorts,
-                            t_star_density);
-        scheduler_addunlock(sched, ci->hydro.super->stars.stars_in,
-                            t_star_density);
-        scheduler_addunlock(sched, t_star_density,
-                            ci->hydro.super->stars.density_ghost);
-#ifdef EXTRA_STAR_LOOPS
-        scheduler_addunlock(sched, ci->hydro.super->stars.density_ghost,
-                            t_star_prep1);
-        scheduler_addunlock(sched, t_star_prep1,
-                            ci->hydro.super->stars.prep1_ghost);
-        scheduler_addunlock(sched, t_star_prep1,
-                            ci->hydro.super->hydro.prep1_ghost);
-        scheduler_addunlock(sched, ci->hydro.super->stars.prep1_ghost,
-                            t_star_prep2);
-        scheduler_addunlock(sched, ci->hydro.super->hydro.prep1_ghost,
-                            t_star_prep2);
-        scheduler_addunlock(sched, t_star_prep2,
-                            ci->hydro.super->stars.prep2_ghost);
-        scheduler_addunlock(sched, ci->hydro.super->stars.prep2_ghost,
-                            t_star_feedback);
-#else
-        scheduler_addunlock(sched, ci->hydro.super->stars.density_ghost,
-                            t_star_feedback);
-#endif
-        scheduler_addunlock(sched, t_star_feedback,
-                            ci->hydro.super->stars.stars_out);
-      }
-
-      if (with_sink) {
-
-        /* Do the sink_formation */
-        scheduler_addunlock(sched, ci->hydro.super->sinks.drift,
-                            ci->top->sinks.sink_formation);
-        scheduler_addunlock(sched, ci->hydro.super->hydro.drift,
-                            ci->top->sinks.sink_formation);
-        scheduler_addunlock(sched, ci->hydro.super->sinks.sink_in,
-                            ci->top->sinks.sink_formation);
-        scheduler_addunlock(sched, ci->top->sinks.sink_formation,
-                            t_sink_swallow);
-
-        /* Do the sink_swallow */
-        scheduler_addunlock(sched, t_sink_swallow,
-                            ci->hydro.super->sinks.sink_ghost1);
-
-        /* Do the sink_do_gas_swallow */
-        scheduler_addunlock(sched, ci->hydro.super->sinks.sink_ghost1,
-                            t_sink_do_gas_swallow);
-        scheduler_addunlock(sched, t_sink_do_gas_swallow,
-                            ci->hydro.super->sinks.sink_ghost2);
-
-        /* Do the sink_do_sink_swallow */
-        scheduler_addunlock(sched, ci->hydro.super->sinks.sink_ghost2,
-                            t_sink_do_sink_swallow);
-        scheduler_addunlock(sched, t_sink_do_sink_swallow,
-                            ci->hydro.super->sinks.sink_out);
-      }
-
-      if (with_black_holes && bcount_i > 0) {
-
-        if (with_cooling)
-          scheduler_addunlock(sched, ci->hydro.super->hydro.cooling_out,
-                              t_bh_density);
-
-        scheduler_addunlock(sched, ci->hydro.super->black_holes.drift,
-                            t_bh_density);
-        scheduler_addunlock(sched, ci->hydro.super->hydro.drift, t_bh_density);
-        scheduler_addunlock(sched, ci->hydro.super->black_holes.black_holes_in,
-                            t_bh_density);
-        scheduler_addunlock(sched, t_bh_density,
-                            ci->hydro.super->black_holes.density_ghost);
-
-        scheduler_addunlock(sched, ci->hydro.super->black_holes.density_ghost,
-                            t_bh_swallow);
-        scheduler_addunlock(sched, t_bh_swallow,
-                            ci->hydro.super->black_holes.swallow_ghost_0);
-
-        scheduler_addunlock(sched, ci->hydro.super->black_holes.swallow_ghost_0,
-                            t_do_gas_swallow);
-        scheduler_addunlock(sched, t_do_gas_swallow,
-                            ci->hydro.super->black_holes.swallow_ghost_1);
-
-        scheduler_addunlock(sched, ci->hydro.super->black_holes.swallow_ghost_1,
-                            t_do_bh_swallow);
-        scheduler_addunlock(sched, t_do_bh_swallow,
-                            ci->hydro.super->black_holes.swallow_ghost_2);
-
-        scheduler_addunlock(sched, ci->hydro.super->black_holes.swallow_ghost_2,
-                            t_bh_feedback);
-        scheduler_addunlock(sched, t_bh_feedback,
-                            ci->hydro.super->black_holes.black_holes_out);
-      }
-
-      if (with_timestep_limiter) {
-        scheduler_addunlock(sched, ci->hydro.super->hydro.drift, t_limiter);
-        scheduler_addunlock(sched, ci->super->timestep, t_limiter);
-        scheduler_addunlock(sched, t_limiter, ci->super->kick1);
-        scheduler_addunlock(sched, t_limiter, ci->super->timestep_limiter);
-      }
-
-      if (with_timestep_sync && with_feedback) {
-        scheduler_addunlock(sched, t_star_feedback, ci->super->timestep_sync);
-      }
-      if (with_timestep_sync && with_black_holes && bcount_i > 0) {
-        scheduler_addunlock(sched, t_bh_feedback, ci->super->timestep_sync);
-      }
-
-      if (with_rt) {
-        scheduler_addunlock(sched, ci->hydro.super->hydro.drift, t_rt_gradient);
-        scheduler_addunlock(sched, ci->hydro.super->hydro.sorts, t_rt_gradient);
-        scheduler_addunlock(sched, ci->hydro.super->rt.rt_ghost1,
-                            t_rt_gradient);
-        scheduler_addunlock(sched, t_rt_gradient,
-                            ci->hydro.super->rt.rt_ghost2);
-        scheduler_addunlock(sched, ci->hydro.super->rt.rt_ghost2,
-                            t_rt_transport);
-        scheduler_addunlock(sched, t_rt_transport,
-                            ci->hydro.super->rt.rt_transport_out);
-      }
-    }
-
-    /* Otherwise, sub-pair interaction? */
-    else if (t_type == task_type_sub_pair &&
-             t_subtype == task_subtype_density) {
-
-      const int bcount_i = ci->black_holes.count;
-      const int bcount_j = cj->black_holes.count;
-
-      /* Make all density tasks depend on the drift */
-      if (ci->nodeID == nodeID) {
-        scheduler_addunlock(sched, ci->hydro.super->hydro.drift, t);
-      }
-      if ((cj->nodeID == nodeID) && (ci->hydro.super != cj->hydro.super)) {
-        scheduler_addunlock(sched, cj->hydro.super->hydro.drift, t);
-      }
-
-      /* Make all density tasks depend on the sorts */
-      scheduler_addunlock(sched, ci->hydro.super->hydro.sorts, t);
-      if (ci->hydro.super != cj->hydro.super) {
-        scheduler_addunlock(sched, cj->hydro.super->hydro.sorts, t);
-      }
-
-      /* New task for the force */
-      t_force = scheduler_addtask(sched, task_type_sub_pair, task_subtype_force,
-                                  flags, 0, ci, cj);
-
-#ifdef MPI_SYMMETRIC_FORCE_INTERACTION
-      /* The order of operations for an inactive local cell interacting
-       * with an active foreign cell is not guaranteed because the density
-       * (and gradient) iact loops don't exist in that case. So we need
-       * an explicit dependency here to have sorted cells. */
-
-      /* Make all force tasks depend on the sorts */
-      scheduler_addunlock(sched, ci->hydro.super->hydro.sorts, t_force);
-      if (ci->hydro.super != cj->hydro.super) {
-        scheduler_addunlock(sched, cj->hydro.super->hydro.sorts, t_force);
-      }
-#endif
-
-      /* and the task for the time-step limiter */
-      if (with_timestep_limiter) {
-        t_limiter = scheduler_addtask(sched, task_type_sub_pair,
-                                      task_subtype_limiter, flags, 0, ci, cj);
-      }
-
-      /* The stellar feedback tasks */
-      if (with_feedback) {
-        t_star_density =
-            scheduler_addtask(sched, task_type_sub_pair,
-                              task_subtype_stars_density, flags, 0, ci, cj);
-        t_star_feedback =
-            scheduler_addtask(sched, task_type_sub_pair,
-                              task_subtype_stars_feedback, flags, 0, ci, cj);
-
-#ifdef EXTRA_STAR_LOOPS
-        t_star_prep1 =
-            scheduler_addtask(sched, task_type_sub_pair,
-                              task_subtype_stars_prep1, flags, 0, ci, cj);
-        t_star_prep2 =
-            scheduler_addtask(sched, task_type_sub_pair,
-                              task_subtype_stars_prep2, flags, 0, ci, cj);
-#endif
-      }
-
-      /* The sink tasks */
-      if (with_sink) {
-        t_sink_swallow =
-            scheduler_addtask(sched, task_type_sub_pair,
-                              task_subtype_sink_swallow, flags, 0, ci, cj);
-        t_sink_do_sink_swallow = scheduler_addtask(
-            sched, task_type_sub_pair, task_subtype_sink_do_sink_swallow, flags,
-            0, ci, cj);
-        t_sink_do_gas_swallow = scheduler_addtask(
-            sched, task_type_sub_pair, task_subtype_sink_do_gas_swallow, flags,
-            0, ci, cj);
-      }
-
-      /* The black hole feedback tasks */
-      if (with_black_holes && (bcount_i > 0 || bcount_j > 0)) {
-        t_bh_density =
-            scheduler_addtask(sched, task_type_sub_pair,
-                              task_subtype_bh_density, flags, 0, ci, cj);
-        t_bh_swallow =
-            scheduler_addtask(sched, task_type_sub_pair,
-                              task_subtype_bh_swallow, flags, 0, ci, cj);
-        t_do_gas_swallow =
-            scheduler_addtask(sched, task_type_sub_pair,
-                              task_subtype_do_gas_swallow, flags, 0, ci, cj);
-        t_do_bh_swallow =
-            scheduler_addtask(sched, task_type_sub_pair,
-                              task_subtype_do_bh_swallow, flags, 0, ci, cj);
-        t_bh_feedback =
-            scheduler_addtask(sched, task_type_sub_pair,
-                              task_subtype_bh_feedback, flags, 0, ci, cj);
-      }
-
-      if (with_rt) {
-        t_rt_gradient =
-            scheduler_addtask(sched, task_type_sub_pair,
-                              task_subtype_rt_gradient, flags, 0, ci, cj);
-        t_rt_transport =
-            scheduler_addtask(sched, task_type_sub_pair,
-                              task_subtype_rt_transport, flags, 0, ci, cj);
-#ifdef MPI_SYMMETRIC_FORCE_INTERACTION
-        /* The order of operations for an inactive local cell interacting
-         * with an active foreign cell is not guaranteed because the gradient
-         * iact loops don't exist in that case. So we need an explicit
-         * dependency here to have sorted cells. */
-
-        /* Make all force tasks depend on the sorts */
-        if (ci->hydro.super->rt.rt_sorts != NULL)
-          scheduler_addunlock(sched, ci->hydro.super->rt.rt_sorts,
-                              t_rt_transport);
-        if (ci->hydro.super != cj->hydro.super) {
-          if (cj->hydro.super->rt.rt_sorts != NULL)
-            scheduler_addunlock(sched, cj->hydro.super->rt.rt_sorts,
-                                t_rt_transport);
-        }
-        /* We need to ensure that a local inactive cell is sorted before
-         * the interaction in the transport loop. Local cells don't have
-         * an rt_sort task. */
-        if (ci->hydro.super->hydro.sorts != NULL)
-          scheduler_addunlock(sched, ci->hydro.super->hydro.sorts,
-                              t_rt_transport);
-        if ((ci->hydro.super != cj->hydro.super) &&
-            (cj->hydro.super->hydro.sorts != NULL))
-          scheduler_addunlock(sched, cj->hydro.super->hydro.sorts,
-                              t_rt_transport);
-#endif
-      }
-
-      engine_addlink(e, &ci->hydro.force, t_force);
-      engine_addlink(e, &cj->hydro.force, t_force);
-      if (with_timestep_limiter) {
-        engine_addlink(e, &ci->hydro.limiter, t_limiter);
-        engine_addlink(e, &cj->hydro.limiter, t_limiter);
-      }
-      if (with_feedback) {
-        engine_addlink(e, &ci->stars.density, t_star_density);
-        engine_addlink(e, &cj->stars.density, t_star_density);
-        engine_addlink(e, &ci->stars.feedback, t_star_feedback);
-        engine_addlink(e, &cj->stars.feedback, t_star_feedback);
-#ifdef EXTRA_STAR_LOOPS
-        engine_addlink(e, &ci->stars.prepare1, t_star_prep1);
-        engine_addlink(e, &cj->stars.prepare1, t_star_prep1);
-        engine_addlink(e, &ci->stars.prepare2, t_star_prep2);
-        engine_addlink(e, &cj->stars.prepare2, t_star_prep2);
-#endif
-      }
-      if (with_sink) {
-        /* Formation */
-        engine_addlink(e, &ci->sinks.swallow, t_sink_swallow);
-        engine_addlink(e, &cj->sinks.swallow, t_sink_swallow);
-
-        /* Merger */
-        engine_addlink(e, &ci->sinks.do_sink_swallow, t_sink_do_sink_swallow);
-        engine_addlink(e, &cj->sinks.do_sink_swallow, t_sink_do_sink_swallow);
-
-        /* Accretion */
-        engine_addlink(e, &ci->sinks.do_gas_swallow, t_sink_do_gas_swallow);
-        engine_addlink(e, &cj->sinks.do_gas_swallow, t_sink_do_gas_swallow);
-      }
-      if (with_black_holes && (bcount_i > 0 || bcount_j > 0)) {
-        engine_addlink(e, &ci->black_holes.density, t_bh_density);
-        engine_addlink(e, &cj->black_holes.density, t_bh_density);
-        engine_addlink(e, &ci->black_holes.swallow, t_bh_swallow);
-        engine_addlink(e, &cj->black_holes.swallow, t_bh_swallow);
-        engine_addlink(e, &ci->black_holes.do_gas_swallow, t_do_gas_swallow);
-        engine_addlink(e, &cj->black_holes.do_gas_swallow, t_do_gas_swallow);
-        engine_addlink(e, &ci->black_holes.do_bh_swallow, t_do_bh_swallow);
-        engine_addlink(e, &cj->black_holes.do_bh_swallow, t_do_bh_swallow);
-        engine_addlink(e, &ci->black_holes.feedback, t_bh_feedback);
-        engine_addlink(e, &cj->black_holes.feedback, t_bh_feedback);
-      }
-      if (with_rt) {
-        engine_addlink(e, &ci->rt.rt_gradient, t_rt_gradient);
-        engine_addlink(e, &cj->rt.rt_gradient, t_rt_gradient);
-        engine_addlink(e, &ci->rt.rt_transport, t_rt_transport);
-        engine_addlink(e, &cj->rt.rt_transport, t_rt_transport);
-      }
-
-#ifdef EXTRA_HYDRO_LOOP
-
-      /* Start by constructing the task for the second and third hydro loop */
-      t_gradient = scheduler_addtask(sched, task_type_sub_pair,
-                                     task_subtype_gradient, flags, 0, ci, cj);
-
-      /* Add the link between the new loop and both cells */
-      engine_addlink(e, &ci->hydro.gradient, t_gradient);
-      engine_addlink(e, &cj->hydro.gradient, t_gradient);
-
-      /* Now, build all the dependencies for the hydro for the cells */
-      /* that are local and are not descendant of the same super_hydro-cells */
-      if (ci->nodeID == nodeID) {
-        engine_make_hydro_loops_dependencies(sched, t, t_gradient, t_force,
-                                             t_limiter, ci, with_cooling,
-                                             with_timestep_limiter);
-      }
-      if ((cj->nodeID == nodeID) && (ci->hydro.super != cj->hydro.super)) {
-        engine_make_hydro_loops_dependencies(sched, t, t_gradient, t_force,
-                                             t_limiter, cj, with_cooling,
-                                             with_timestep_limiter);
-      }
-#else
-
-      /* Now, build all the dependencies for the hydro for the cells */
-      /* that are local and are not descendant of the same super_hydro-cells */
-      if (ci->nodeID == nodeID) {
-        engine_make_hydro_loops_dependencies(sched, t, t_force, t_limiter, ci,
-                                             with_cooling,
-                                             with_timestep_limiter);
-      }
-      if ((cj->nodeID == nodeID) && (ci->hydro.super != cj->hydro.super)) {
-        engine_make_hydro_loops_dependencies(sched, t, t_force, t_limiter, cj,
-                                             with_cooling,
-                                             with_timestep_limiter);
-      }
-#endif
-
-      if (with_feedback) {
-        scheduler_addunlock(sched, ci->hydro.super->hydro.sorts,
-                            t_star_density);
-        if (ci->hydro.super != cj->hydro.super) {
-          scheduler_addunlock(sched, cj->hydro.super->hydro.sorts,
-                              t_star_density);
-        }
-      }
-
-      if (with_rt) {
-        scheduler_addunlock(sched, ci->hydro.super->hydro.sorts, t_rt_gradient);
-        if (ci->hydro.super != cj->hydro.super) {
-          scheduler_addunlock(sched, cj->hydro.super->hydro.sorts,
-                              t_rt_gradient);
-        }
-      }
-
-      if (ci->nodeID == nodeID) {
-        scheduler_addunlock(sched, t_force, ci->hydro.super->hydro.end_force);
-
-        if (with_feedback) {
-
-          if (with_cooling)
-            scheduler_addunlock(sched, ci->hydro.super->hydro.cooling_out,
-                                t_star_density);
-
-          scheduler_addunlock(sched, ci->hydro.super->stars.sorts,
-                              t_star_density);
-          scheduler_addunlock(sched, ci->hydro.super->stars.drift,
-                              t_star_density);
-          scheduler_addunlock(sched, ci->hydro.super->hydro.drift,
-                              t_star_density);
-          scheduler_addunlock(sched, ci->hydro.super->stars.stars_in,
-                              t_star_density);
-          scheduler_addunlock(sched, t_star_density,
-                              ci->hydro.super->stars.density_ghost);
-#ifdef EXTRA_STAR_LOOPS
-          scheduler_addunlock(sched, ci->hydro.super->stars.density_ghost,
-                              t_star_prep1);
-          scheduler_addunlock(sched, t_star_prep1,
-                              ci->hydro.super->stars.prep1_ghost);
-          scheduler_addunlock(sched, t_star_prep1,
-                              ci->hydro.super->hydro.prep1_ghost);
-          scheduler_addunlock(sched, ci->hydro.super->stars.prep1_ghost,
-                              t_star_prep2);
-          scheduler_addunlock(sched, ci->hydro.super->hydro.prep1_ghost,
-                              t_star_prep2);
-          scheduler_addunlock(sched, t_star_prep2,
-                              ci->hydro.super->stars.prep2_ghost);
-          scheduler_addunlock(sched, ci->hydro.super->stars.prep2_ghost,
-                              t_star_feedback);
-#else
-          scheduler_addunlock(sched, ci->hydro.super->stars.density_ghost,
-                              t_star_feedback);
-#endif
-          scheduler_addunlock(sched, t_star_feedback,
-                              ci->hydro.super->stars.stars_out);
-        }
-
-        if (with_sink) {
-
-          /* Do the sink_formation */
-          scheduler_addunlock(sched, ci->hydro.super->sinks.drift,
-                              ci->top->sinks.sink_formation);
-          scheduler_addunlock(sched, ci->hydro.super->hydro.drift,
-                              ci->top->sinks.sink_formation);
-          scheduler_addunlock(sched, ci->hydro.super->sinks.sink_in,
-                              ci->top->sinks.sink_formation);
-          scheduler_addunlock(sched, ci->top->sinks.sink_formation,
-                              t_sink_swallow);
-
-          /* Do the sink_swallow */
-          scheduler_addunlock(sched, t_sink_swallow,
-                              ci->hydro.super->sinks.sink_ghost1);
-
-          /* Do the sink_do_gas_swallow */
-          scheduler_addunlock(sched, ci->hydro.super->sinks.sink_ghost1,
-                              t_sink_do_gas_swallow);
-          scheduler_addunlock(sched, t_sink_do_gas_swallow,
-                              ci->hydro.super->sinks.sink_ghost2);
-
-          /* Do the sink_do_sink_swallow */
-          scheduler_addunlock(sched, ci->hydro.super->sinks.sink_ghost2,
-                              t_sink_do_sink_swallow);
-          scheduler_addunlock(sched, t_sink_do_sink_swallow,
-                              ci->hydro.super->sinks.sink_out);
-        }
-
-        if (with_black_holes && (bcount_i > 0 || bcount_j > 0)) {
-
-          if (with_cooling)
-            scheduler_addunlock(sched, ci->hydro.super->hydro.cooling_out,
-                                t_bh_density);
-
-          scheduler_addunlock(sched, ci->hydro.super->black_holes.drift,
-                              t_bh_density);
-          scheduler_addunlock(sched, ci->hydro.super->hydro.drift,
-                              t_bh_density);
-          scheduler_addunlock(
-              sched, ci->hydro.super->black_holes.black_holes_in, t_bh_density);
-          scheduler_addunlock(sched, t_bh_density,
-                              ci->hydro.super->black_holes.density_ghost);
-
-          scheduler_addunlock(sched, ci->hydro.super->black_holes.density_ghost,
-                              t_bh_swallow);
-          scheduler_addunlock(sched, t_bh_swallow,
-                              ci->hydro.super->black_holes.swallow_ghost_0);
-
-          scheduler_addunlock(sched,
-                              ci->hydro.super->black_holes.swallow_ghost_0,
-                              t_do_gas_swallow);
-          scheduler_addunlock(sched, t_do_gas_swallow,
-                              ci->hydro.super->black_holes.swallow_ghost_1);
-
-          scheduler_addunlock(sched,
-                              ci->hydro.super->black_holes.swallow_ghost_1,
-                              t_do_bh_swallow);
-          scheduler_addunlock(sched, t_do_bh_swallow,
-                              ci->hydro.super->black_holes.swallow_ghost_2);
-
-          scheduler_addunlock(sched,
-                              ci->hydro.super->black_holes.swallow_ghost_2,
-                              t_bh_feedback);
-          scheduler_addunlock(sched, t_bh_feedback,
-                              ci->hydro.super->black_holes.black_holes_out);
-        }
-
-        if (with_timestep_limiter) {
-          scheduler_addunlock(sched, ci->hydro.super->hydro.drift, t_limiter);
-          scheduler_addunlock(sched, ci->super->timestep, t_limiter);
-          scheduler_addunlock(sched, t_limiter, ci->super->kick1);
-          scheduler_addunlock(sched, t_limiter, ci->super->timestep_limiter);
-        }
-
-        if (with_timestep_sync && with_feedback) {
-          scheduler_addunlock(sched, t_star_feedback, ci->super->timestep_sync);
-        }
-        if (with_timestep_sync && with_black_holes &&
-            (bcount_i > 0 || bcount_j > 0)) {
-          scheduler_addunlock(sched, t_bh_feedback, ci->super->timestep_sync);
-        }
-
-        if (with_rt) {
-          scheduler_addunlock(sched, ci->hydro.super->hydro.drift,
-                              t_rt_gradient);
-          scheduler_addunlock(sched, ci->hydro.super->rt.rt_ghost1,
-                              t_rt_gradient);
-          scheduler_addunlock(sched, t_rt_gradient,
-                              ci->hydro.super->rt.rt_ghost2);
-          scheduler_addunlock(sched, ci->hydro.super->rt.rt_ghost2,
-                              t_rt_transport);
-          scheduler_addunlock(sched, t_rt_transport,
-                              ci->hydro.super->rt.rt_transport_out);
-        }
-      } else /* ci->nodeID != nodeID */ {
-
-        if (with_feedback) {
-#ifdef EXTRA_STAR_LOOPS
-          scheduler_addunlock(sched, ci->hydro.super->stars.sorts,
-                              t_star_prep1);
-#endif
-          scheduler_addunlock(sched, ci->hydro.super->stars.sorts,
-                              t_star_feedback);
-        }
-        if (with_black_holes && (bcount_i > 0 || bcount_j > 0)) {
-
-          scheduler_addunlock(sched, t_bh_swallow,
-                              ci->hydro.super->black_holes.swallow_ghost_0);
-        }
-      }
-
-      if (cj->nodeID == nodeID) {
-
-        if (ci->hydro.super != cj->hydro.super) {
-
-          scheduler_addunlock(sched, t_force, cj->hydro.super->hydro.end_force);
-
-          if (with_feedback) {
-
-            if (with_cooling)
-              scheduler_addunlock(sched, cj->hydro.super->hydro.cooling_out,
-                                  t_star_density);
-
-            scheduler_addunlock(sched, cj->hydro.super->stars.sorts,
-                                t_star_density);
-            scheduler_addunlock(sched, cj->hydro.super->stars.drift,
-                                t_star_density);
-            scheduler_addunlock(sched, cj->hydro.super->hydro.drift,
-                                t_star_density);
-            scheduler_addunlock(sched, cj->hydro.super->stars.stars_in,
-                                t_star_density);
-            scheduler_addunlock(sched, t_star_density,
-                                cj->hydro.super->stars.density_ghost);
-#ifdef EXTRA_STAR_LOOPS
-            scheduler_addunlock(sched, cj->hydro.super->stars.density_ghost,
-                                t_star_prep1);
-            scheduler_addunlock(sched, t_star_prep1,
-                                cj->hydro.super->stars.prep1_ghost);
-            scheduler_addunlock(sched, t_star_prep1,
-                                cj->hydro.super->hydro.prep1_ghost);
-            scheduler_addunlock(sched, cj->hydro.super->stars.prep1_ghost,
-                                t_star_prep2);
-            scheduler_addunlock(sched, cj->hydro.super->hydro.prep1_ghost,
-                                t_star_prep2);
-            scheduler_addunlock(sched, t_star_prep2,
-                                cj->hydro.super->stars.prep2_ghost);
-            scheduler_addunlock(sched, cj->hydro.super->stars.prep2_ghost,
-                                t_star_feedback);
-#else
-            scheduler_addunlock(sched, cj->hydro.super->stars.density_ghost,
-                                t_star_feedback);
-#endif
-            scheduler_addunlock(sched, t_star_feedback,
-                                cj->hydro.super->stars.stars_out);
-          }
-          if (with_sink) {
-
-            /* Do the sink_formation */
-            scheduler_addunlock(sched, cj->hydro.super->sinks.drift,
-                                cj->top->sinks.sink_formation);
-            scheduler_addunlock(sched, cj->hydro.super->hydro.drift,
-                                cj->top->sinks.sink_formation);
-            scheduler_addunlock(sched, cj->hydro.super->sinks.sink_in,
-                                cj->top->sinks.sink_formation);
-            scheduler_addunlock(sched, cj->top->sinks.sink_formation,
-                                t_sink_swallow);
-
-            /* Do the sink_swallow */
-            scheduler_addunlock(sched, t_sink_swallow,
-                                cj->hydro.super->sinks.sink_ghost1);
-
-            /* Do the sink_do_gas_swallow */
-            scheduler_addunlock(sched, cj->hydro.super->sinks.sink_ghost1,
-                                t_sink_do_gas_swallow);
-            scheduler_addunlock(sched, t_sink_do_gas_swallow,
-                                cj->hydro.super->sinks.sink_ghost2);
-
-            /* Do the sink_do_sink_swallow */
-            scheduler_addunlock(sched, cj->hydro.super->sinks.sink_ghost2,
-                                t_sink_do_sink_swallow);
-            scheduler_addunlock(sched, t_sink_do_sink_swallow,
-                                cj->hydro.super->sinks.sink_out);
-          }
-
-          if (with_black_holes && (bcount_i > 0 || bcount_j > 0)) {
-
-            if (with_cooling)
-              scheduler_addunlock(sched, cj->hydro.super->hydro.cooling_out,
-                                  t_bh_density);
-
-            scheduler_addunlock(sched, cj->hydro.super->black_holes.drift,
-                                t_bh_density);
-            scheduler_addunlock(sched, cj->hydro.super->hydro.drift,
-                                t_bh_density);
-            scheduler_addunlock(sched,
-                                cj->hydro.super->black_holes.black_holes_in,
-                                t_bh_density);
-            scheduler_addunlock(sched, t_bh_density,
-                                cj->hydro.super->black_holes.density_ghost);
-
-            scheduler_addunlock(sched,
-                                cj->hydro.super->black_holes.density_ghost,
-                                t_bh_swallow);
-            scheduler_addunlock(sched, t_bh_swallow,
-                                cj->hydro.super->black_holes.swallow_ghost_0);
-
-            scheduler_addunlock(sched,
-                                cj->hydro.super->black_holes.swallow_ghost_0,
-                                t_do_gas_swallow);
-            scheduler_addunlock(sched, t_do_gas_swallow,
-                                cj->hydro.super->black_holes.swallow_ghost_1);
-
-            scheduler_addunlock(sched,
-                                cj->hydro.super->black_holes.swallow_ghost_1,
-                                t_do_bh_swallow);
-            scheduler_addunlock(sched, t_do_bh_swallow,
-                                cj->hydro.super->black_holes.swallow_ghost_2);
-
-            scheduler_addunlock(sched,
-                                cj->hydro.super->black_holes.swallow_ghost_2,
-                                t_bh_feedback);
-            scheduler_addunlock(sched, t_bh_feedback,
-                                cj->hydro.super->black_holes.black_holes_out);
-          }
-          if (with_rt) {
-            scheduler_addunlock(sched, cj->hydro.super->hydro.drift,
-                                t_rt_gradient);
-            scheduler_addunlock(sched, cj->hydro.super->rt.rt_ghost1,
-                                t_rt_gradient);
-            scheduler_addunlock(sched, t_rt_gradient,
-                                cj->hydro.super->rt.rt_ghost2);
-            scheduler_addunlock(sched, cj->hydro.super->rt.rt_ghost2,
-                                t_rt_transport);
-            scheduler_addunlock(sched, t_rt_transport,
-                                cj->hydro.super->rt.rt_transport_out);
-          }
-
-          if (with_timestep_limiter) {
-            scheduler_addunlock(sched, cj->hydro.super->hydro.drift, t_limiter);
-          }
-        }
-
-        if (ci->super != cj->super) {
-
-          if (with_timestep_limiter) {
-            scheduler_addunlock(sched, cj->super->timestep, t_limiter);
-            scheduler_addunlock(sched, t_limiter, cj->super->kick1);
-            scheduler_addunlock(sched, t_limiter, cj->super->timestep_limiter);
-          }
-
-          if (with_timestep_sync && with_feedback) {
-            scheduler_addunlock(sched, t_star_feedback,
-                                cj->super->timestep_sync);
-          }
-          if (with_timestep_sync && with_black_holes &&
-              (bcount_i > 0 || bcount_j > 0)) {
-            scheduler_addunlock(sched, t_bh_feedback, cj->super->timestep_sync);
-          }
-        }
-      } else /* cj->nodeID != nodeID */ {
-        if (with_feedback) {
-#ifdef EXTRA_STAR_LOOPS
-          scheduler_addunlock(sched, cj->hydro.super->stars.sorts,
-                              t_star_prep1);
-#endif
-          scheduler_addunlock(sched, cj->hydro.super->stars.sorts,
-                              t_star_feedback);
-        }
-
-        if (with_black_holes && (bcount_i > 0 || bcount_j > 0)) {
-          scheduler_addunlock(sched, t_bh_swallow,
-                              cj->hydro.super->black_holes.swallow_ghost_0);
-        }
-      }
-    }
-  }
-}
-
-/**
- * @brief Constructs the top-level pair tasks for the first hydro loop over
- * neighbours
- *
- * Here we construct all the tasks for all possible neighbouring non-empty
- * local cells in the hierarchy. No dependencies are being added thus far.
- * Additional loop over neighbours can later be added by simply duplicating
- * all the tasks created by this function.
- *
- * @param map_data Offset of first two indices disguised as a pointer.
- * @param num_elements Number of cells to traverse.
- * @param extra_data The #engine.
- */
-void engine_make_hydroloop_tasks_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;
-  const int with_feedback = (e->policy & engine_policy_feedback);
-  const int with_stars = (e->policy & engine_policy_stars);
-  const int with_sinks = (e->policy & engine_policy_sinks);
-  const int with_black_holes = (e->policy & engine_policy_black_holes);
-
-  struct space *s = e->s;
-  struct scheduler *sched = &e->sched;
-  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. */
-  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];
-
-    /* Skip cells without hydro or star particles */
-    if ((ci->hydro.count == 0) && (!with_stars || ci->stars.count == 0) &&
-        (!with_sinks || ci->sinks.count == 0) &&
-        (!with_black_holes || ci->black_holes.count == 0))
-      continue;
-
-    /* If the cell is local build a self-interaction */
-    if (ci->nodeID == nodeID) {
-      scheduler_addtask(sched, task_type_self, task_subtype_density, 0, 0, ci,
-                        NULL);
-    }
+    /* If the cell is local build a self-interaction */
+    if (ci->nodeID == nodeID) {
+      scheduler_addtask(sched, task_type_self, task_subtype_density, 0, 0, ci,
+                        NULL);
+    }
 
     /* Now loop over all the neighbours of this cell */
     for (int ii = -1; ii < 2; ii++) {
@@ -4291,11 +3560,12 @@ void engine_addtasks_send_mapper(void *map_data, int num_elements,
   const int with_feedback = (e->policy & engine_policy_feedback);
   const int with_sync = (e->policy & engine_policy_timestep_sync);
   const int with_rt = (e->policy & engine_policy_rt);
+  const int with_fof = (e->policy & engine_policy_fof);
   struct cell_type_pair *cell_type_pairs = (struct cell_type_pair *)map_data;
 
-#ifdef SWIFT_DEBUG_CHECKS
+#if defined(WITH_MPI) && !defined(SWIFT_DEBUG_CHECKS)
   if (e->policy & engine_policy_sinks) {
-    error("TODO");
+    error("TODO: Sink MPI tasks are not implemented yet!");
   }
 #endif
 
@@ -4349,7 +3619,10 @@ void engine_addtasks_send_mapper(void *map_data, int num_elements,
      * connection. */
     if ((e->policy & engine_policy_self_gravity) &&
         (type & proxy_cell_type_gravity))
-      engine_addtasks_send_gravity(e, ci, cj, /*t_grav=*/NULL);
+      engine_addtasks_send_gravity(e, ci, cj, /*t_grav_count=*/NULL,
+                                   /*t_grav=*/NULL, /*t_pack_grav=*/NULL,
+                                   /*t_fof=*/NULL, /*t_pack_fof=*/NULL,
+                                   with_fof, with_star_formation);
   }
 }
 
@@ -4362,12 +3635,13 @@ void engine_addtasks_recv_mapper(void *map_data, int num_elements,
   const int with_feedback = (e->policy & engine_policy_feedback);
   const int with_black_holes = (e->policy & engine_policy_black_holes);
   const int with_sync = (e->policy & engine_policy_timestep_sync);
+  const int with_fof = (e->policy & engine_policy_fof);
   const int with_rt = (e->policy & engine_policy_rt);
   struct cell_type_pair *cell_type_pairs = (struct cell_type_pair *)map_data;
 
-#ifdef SWIFT_DEBUG_CHECKS
+#if defined(WITH_MPI) && !defined(SWIFT_DEBUG_CHECKS)
   if (e->policy & engine_policy_sinks) {
-    error("TODO");
+    error("TODO: Sink MPI tasks are not implemented yet!");
   }
 #endif
 
@@ -4445,7 +3719,9 @@ void engine_addtasks_recv_mapper(void *map_data, int num_elements,
      * connection. */
     if ((e->policy & engine_policy_self_gravity) &&
         (type & proxy_cell_type_gravity))
-      engine_addtasks_recv_gravity(e, ci, /*t_grav=*/NULL, tend);
+      engine_addtasks_recv_gravity(e, ci, /*t_grav_count=*/NULL,
+                                   /*t_grav=*/NULL, /*t_fof=*/NULL, tend,
+                                   with_fof, with_star_formation);
   }
 }
 
@@ -4545,7 +3821,10 @@ void engine_make_fof_tasks(struct engine *e) {
   struct scheduler *sched = &e->sched;
   ticks tic = getticks();
 
-  if (e->restarting) error("Running FOF on a restart step!");
+  if (e->restarting) {
+    /* Re-set the scheduler. */
+    scheduler_reset(sched, engine_estimate_nr_tasks(e));
+  }
 
   /* Construct a FOF loop over neighbours */
   if (e->policy & engine_policy_fof)
@@ -4721,8 +4000,13 @@ void engine_maketasks(struct engine *e) {
   tic2 = getticks();
 
   /* Add the dependencies for the gravity stuff */
-  if (e->policy & (engine_policy_self_gravity | engine_policy_external_gravity))
-    engine_link_gravity_tasks(e);
+  if (e->policy &
+      (engine_policy_self_gravity | engine_policy_external_gravity)) {
+
+    threadpool_map(&e->threadpool, engine_link_gravity_tasks_mapper,
+                   e->sched.tasks, e->sched.nr_tasks, sizeof(struct task),
+                   threadpool_auto_chunk_size, e);
+  }
 
   if (e->verbose)
     message("Linking gravity tasks took %.3f %s.",
@@ -4840,7 +4124,7 @@ void engine_maketasks(struct engine *e) {
   tic2 = getticks();
 
   /* Set the unlocks per task. */
-  scheduler_set_unlocks(sched);
+  scheduler_set_unlocks(sched, &e->threadpool);
 
   if (e->verbose)
     message("Setting unlocks took %.3f %s.",

src/engine_marktasks.c

-/*******************************************************************************
- * This file is part of SWIFT.
- * Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
- *                    Matthieu Schaller (schaller@strw.leidenuniv.nl)
- *               2015 Peter W. Draper (p.w.draper@durham.ac.uk)
- *                    Angus Lepper (angus.lepper@ed.ac.uk)
- *               2016 John A. Regan (john.a.regan@durham.ac.uk)
- *                    Tom Theuns (tom.theuns@durham.ac.uk)
- *
- * 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>
-
-/* Some standard headers. */
-#include <stdlib.h>
-#include <unistd.h>
-
-/* MPI headers. */
-#ifdef WITH_MPI
-#include <mpi.h>
-#endif
-
-/* Load the profiler header, if needed. */
-#ifdef WITH_PROFILER
-#include <gperftools/profiler.h>
-#endif
-
-/* This object's header. */
-#include "engine.h"
-
-/* Local headers. */
-#include "active.h"
-#include "atomic.h"
-#include "cell.h"
-#include "clocks.h"
-#include "cycle.h"
-#include "debug.h"
-#include "error.h"
-#include "feedback.h"
-#include "proxy.h"
-#include "timers.h"
-
-/**
- * @brief Mark tasks to be un-skipped and set the sort flags accordingly.
- *        Threadpool mapper function.
- *
- * @param map_data pointer to the tasks
- * @param num_elements number of tasks
- * @param extra_data pointer to int that will define if a rebuild is needed.
- */
-void engine_marktasks_mapper(void *map_data, int num_elements,
-                             void *extra_data) {
-  /* Unpack the arguments. */
-  struct task *tasks = (struct task *)map_data;
-  size_t *rebuild_space = &((size_t *)extra_data)[1];
-  struct scheduler *s = (struct scheduler *)(((size_t *)extra_data)[2]);
-  struct engine *e = (struct engine *)((size_t *)extra_data)[0];
-  const int nodeID = e->nodeID;
-  const int with_timestep_limiter = e->policy & engine_policy_timestep_limiter;
-  const int with_timestep_sync = e->policy & engine_policy_timestep_sync;
-  const int with_sinks = e->policy & engine_policy_sinks;
-  const int with_stars = e->policy & engine_policy_stars;
-  const int with_star_formation = (e->policy & engine_policy_star_formation);
-  const int with_star_formation_sink = with_sinks && with_stars;
-  const int with_feedback = e->policy & engine_policy_feedback;
-
-  for (int ind = 0; ind < num_elements; ind++) {
-
-    /* Get basic task information */
-    struct task *t = &tasks[ind];
-    const enum task_types t_type = t->type;
-    const enum task_subtypes t_subtype = t->subtype;
-
-    /* Single-cell task? */
-    if (t_type == task_type_self || t_type == task_type_sub_self) {
-
-      /* Local pointer. */
-      struct cell *ci = t->ci;
-
-#ifdef SWIFT_DEBUG_CHECKS
-      if (ci->nodeID != nodeID) error("Non-local self task found");
-#endif
-
-      const int ci_active_hydro = cell_is_active_hydro(ci, e);
-      const int ci_active_gravity = cell_is_active_gravity(ci, e);
-      const int ci_active_black_holes = cell_is_active_black_holes(ci, e);
-      const int ci_active_sinks =
-          cell_is_active_sinks(ci, e) || ci_active_hydro;
-      const int ci_active_stars = cell_need_activating_stars(
-          ci, e, with_star_formation, with_star_formation_sink);
-      const int ci_active_rt = cell_is_rt_active(ci, e);
-
-      /* Activate the hydro drift */
-      if (t_type == task_type_self && t_subtype == task_subtype_density) {
-        if (ci_active_hydro) {
-          scheduler_activate(s, t);
-          cell_activate_drift_part(ci, s);
-          if (with_timestep_limiter) cell_activate_limiter(ci, s);
-        }
-      }
-
-      /* Store current values of dx_max and h_max. */
-      else if (t_type == task_type_sub_self &&
-               t_subtype == task_subtype_density) {
-        if (ci_active_hydro) {
-          scheduler_activate(s, t);
-          cell_activate_subcell_hydro_tasks(ci, NULL, s, with_timestep_limiter);
-          if (with_timestep_limiter) cell_activate_limiter(ci, s);
-        }
-      }
-
-      else if (t_type == task_type_self && t_subtype == task_subtype_force) {
-        if (ci_active_hydro) scheduler_activate(s, t);
-      }
-
-      else if (t_type == task_type_sub_self &&
-               t_subtype == task_subtype_force) {
-        if (ci_active_hydro) scheduler_activate(s, t);
-      }
-
-      else if (t->type == task_type_self &&
-               t->subtype == task_subtype_limiter) {
-        if (ci_active_hydro) scheduler_activate(s, t);
-      }
-
-      else if (t->type == task_type_sub_self &&
-               t->subtype == task_subtype_limiter) {
-        if (ci_active_hydro) scheduler_activate(s, t);
-      }
-
-      else if (t_type == task_type_self && t_subtype == task_subtype_gradient) {
-        if (ci_active_hydro) scheduler_activate(s, t);
-      }
-
-      else if (t_type == task_type_sub_self &&
-               t_subtype == task_subtype_gradient) {
-        if (ci_active_hydro) scheduler_activate(s, t);
-      }
-
-      /* Activate the star density */
-      else if (t_type == task_type_self &&
-               t_subtype == task_subtype_stars_density) {
-        if (ci_active_stars) {
-          scheduler_activate(s, t);
-          cell_activate_drift_part(ci, s);
-          cell_activate_drift_spart(ci, s);
-          if (with_timestep_sync) cell_activate_sync_part(ci, s);
-        }
-      }
-
-      /* Store current values of dx_max and h_max. */
-      else if (t_type == task_type_sub_self &&
-               t_subtype == task_subtype_stars_density) {
-        if (ci_active_stars) {
-          scheduler_activate(s, t);
-          cell_activate_subcell_stars_tasks(ci, NULL, s, with_star_formation,
-                                            with_star_formation_sink,
-                                            with_timestep_sync);
-        }
-      }
-
-      else if (t_type == task_type_self &&
-               t_subtype == task_subtype_stars_prep1) {
-        if (ci_active_stars) {
-          scheduler_activate(s, t);
-        }
-      }
-
-      else if (t_type == task_type_sub_self &&
-               t_subtype == task_subtype_stars_prep1) {
-        if (ci_active_stars) scheduler_activate(s, t);
-      }
-
-      else if (t_type == task_type_self &&
-               t_subtype == task_subtype_stars_prep2) {
-        if (ci_active_stars) {
-          scheduler_activate(s, t);
-        }
-      }
-
-      else if (t_type == task_type_sub_self &&
-               t_subtype == task_subtype_stars_prep2) {
-        if (ci_active_stars) scheduler_activate(s, t);
-      }
-
-      else if (t_type == task_type_self &&
-               t_subtype == task_subtype_stars_feedback) {
-        if (ci_active_stars) {
-          scheduler_activate(s, t);
-        }
-      }
-
-      else if (t_type == task_type_sub_self &&
-               t_subtype == task_subtype_stars_feedback) {
-        if (ci_active_stars) scheduler_activate(s, t);
-      }
-
-      /* Activate the sink formation */
-      else if (t_type == task_type_self &&
-               t_subtype == task_subtype_sink_swallow) {
-        if (ci_active_sinks) {
-          scheduler_activate(s, t);
-          cell_activate_drift_part(ci, s);
-          cell_activate_drift_sink(ci, s);
-          cell_activate_sink_formation_tasks(ci->top, s);
-          if (with_timestep_sync) cell_activate_sync_part(ci, s);
-        }
-      }
-
-      /* Store current values of dx_max and h_max. */
-      else if (t_type == task_type_sub_self &&
-               t_subtype == task_subtype_sink_swallow) {
-        if (ci_active_sinks) {
-          scheduler_activate(s, t);
-          cell_activate_subcell_sinks_tasks(ci, NULL, s, with_timestep_sync);
-        }
-      }
-
-      /* Activate the sink merger */
-      else if (t_type == task_type_self &&
-               t_subtype == task_subtype_sink_do_sink_swallow) {
-        if (ci_active_sinks) {
-          scheduler_activate(s, t);
-        }
-      }
-
-      else if (t_type == task_type_sub_self &&
-               t_subtype == task_subtype_sink_do_sink_swallow) {
-        if (ci_active_sinks) {
-          scheduler_activate(s, t);
-        }
-      }
-
-      /* Activate the sink accretion */
-      else if (t_type == task_type_self &&
-               t_subtype == task_subtype_sink_do_gas_swallow) {
-        if (ci_active_sinks) {
-          scheduler_activate(s, t);
-        }
-      }
-
-      else if (t_type == task_type_sub_self &&
-               t_subtype == task_subtype_sink_do_gas_swallow) {
-        if (ci_active_sinks) {
-          scheduler_activate(s, t);
-        }
-      }
-
-      /* Activate the black hole density */
-      else if (t_type == task_type_self &&
-               t_subtype == task_subtype_bh_density) {
-        if (ci_active_black_holes) {
-          scheduler_activate(s, t);
-          cell_activate_drift_part(ci, s);
-          cell_activate_drift_bpart(ci, s);
-          if (with_timestep_sync) cell_activate_sync_part(ci, s);
-        }
-      }
-
-      /* Store current values of dx_max and h_max. */
-      else if (t_type == task_type_sub_self &&
-               t_subtype == task_subtype_bh_density) {
-        if (ci_active_black_holes) {
-          scheduler_activate(s, t);
-          cell_activate_subcell_black_holes_tasks(ci, NULL, s,
-                                                  with_timestep_sync);
-        }
-      }
-
-      else if (t_type == task_type_self &&
-               t_subtype == task_subtype_bh_swallow) {
-        if (ci_active_black_holes) {
-          scheduler_activate(s, t);
-        }
-      }
-
-      else if (t_type == task_type_sub_self &&
-               t_subtype == task_subtype_bh_swallow) {
-        if (ci_active_black_holes) scheduler_activate(s, t);
-      }
-
-      else if (t_type == task_type_self &&
-               t_subtype == task_subtype_do_gas_swallow) {
-        if (ci_active_black_holes) {
-          scheduler_activate(s, t);
-        }
-      }
-
-      else if (t_type == task_type_sub_self &&
-               t_subtype == task_subtype_do_gas_swallow) {
-        if (ci_active_black_holes) scheduler_activate(s, t);
-      }
-
-      else if (t_type == task_type_self &&
-               t_subtype == task_subtype_do_bh_swallow) {
-        if (ci_active_black_holes) {
-          scheduler_activate(s, t);
-        }
-      }
-
-      else if (t_type == task_type_sub_self &&
-               t_subtype == task_subtype_do_bh_swallow) {
-        if (ci_active_black_holes) scheduler_activate(s, t);
-      }
-
-      else if (t_type == task_type_self &&
-               t_subtype == task_subtype_bh_feedback) {
-        if (ci_active_black_holes) {
-          scheduler_activate(s, t);
-        }
-      }
-
-      else if (t_type == task_type_sub_self &&
-               t_subtype == task_subtype_bh_feedback) {
-        if (ci_active_black_holes) scheduler_activate(s, t);
-      }
-
-      /* Activate the gravity drift */
-      else if (t_type == task_type_self && t_subtype == task_subtype_grav) {
-        if (ci_active_gravity) {
-          scheduler_activate(s, t);
-          cell_activate_subcell_grav_tasks(t->ci, NULL, s);
-        }
-      }
-
-      /* Activate the gravity drift */
-      else if (t_type == task_type_self &&
-               t_subtype == task_subtype_external_grav) {
-        if (ci_active_gravity) {
-          scheduler_activate(s, t);
-          cell_activate_drift_gpart(t->ci, s);
-        }
-      }
-
-      /* Activate RT tasks */
-      else if (t_type == task_type_self &&
-               t_subtype == task_subtype_rt_gradient) {
-        if (ci_active_rt) scheduler_activate(s, t);
-      }
-
-      else if (t_type == task_type_sub_self &&
-               t_subtype == task_subtype_rt_gradient) {
-        if (ci_active_rt) {
-          scheduler_activate(s, t);
-          cell_activate_subcell_rt_tasks(ci, NULL, s, /*sub_cycle=*/0);
-        }
-      }
-
-      else if (t_subtype == task_subtype_rt_transport) {
-        if (ci_active_rt) scheduler_activate(s, t);
-      }
-
-#ifdef SWIFT_DEBUG_CHECKS
-      else {
-        error("Invalid task type / sub-type encountered");
-      }
-#endif
-    }
-
-    /* Pair? */
-    else if (t_type == task_type_pair || t_type == task_type_sub_pair) {
-
-      /* Local pointers. */
-      struct cell *ci = t->ci;
-      struct cell *cj = t->cj;
-#ifdef WITH_MPI
-      const int ci_nodeID = ci->nodeID;
-      const int cj_nodeID = cj->nodeID;
-#else
-      const int ci_nodeID = nodeID;
-      const int cj_nodeID = nodeID;
-#endif
-      const int ci_active_hydro = cell_is_active_hydro(ci, e);
-      const int cj_active_hydro = cell_is_active_hydro(cj, e);
-
-      const int ci_active_gravity = cell_is_active_gravity(ci, e);
-      const int cj_active_gravity = cell_is_active_gravity(cj, e);
-
-      const int ci_active_black_holes = cell_is_active_black_holes(ci, e);
-      const int cj_active_black_holes = cell_is_active_black_holes(cj, e);
-
-      const int ci_active_sinks =
-          cell_is_active_sinks(ci, e) || ci_active_hydro;
-      const int cj_active_sinks =
-          cell_is_active_sinks(cj, e) || cj_active_hydro;
-
-      const int ci_active_stars = cell_need_activating_stars(
-          ci, e, with_star_formation, with_star_formation_sink);
-      const int cj_active_stars = cell_need_activating_stars(
-          cj, e, with_star_formation, with_star_formation_sink);
-
-      const int ci_active_rt = cell_is_rt_active(ci, e);
-      const int cj_active_rt = cell_is_rt_active(cj, e);
-
-      /* Only activate tasks that involve a local active cell. */
-      if ((t_subtype == task_subtype_density ||
-           t_subtype == task_subtype_gradient ||
-           t_subtype == task_subtype_limiter ||
-           t_subtype == task_subtype_force) &&
-          ((ci_active_hydro && ci_nodeID == nodeID) ||
-           (cj_active_hydro && cj_nodeID == nodeID))) {
-
-        scheduler_activate(s, t);
-
-        /* Set the correct sorting flags */
-        if (t_type == task_type_pair && t_subtype == task_subtype_density) {
-
-          /* 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 hydro drift tasks. */
-          if (ci_nodeID == nodeID) cell_activate_drift_part(ci, s);
-          if (cj_nodeID == nodeID) cell_activate_drift_part(cj, s);
-
-          /* And the limiter */
-          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 &&
-                 t_subtype == task_subtype_density) {
-          cell_activate_subcell_hydro_tasks(t->ci, t->cj, s,
-                                            with_timestep_limiter);
-        }
-      }
-
-      /* Stars density */
-      else if ((t_subtype == task_subtype_stars_density) &&
-               (ci_active_stars || cj_active_stars) &&
-               (ci_nodeID == nodeID || cj_nodeID == nodeID)) {
-
-        scheduler_activate(s, t);
-
-        /* Set the correct sorting flags */
-        if (t_type == task_type_pair) {
-
-          /* Add stars_in dependencies 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->stars.stars_in);
-          if (cj_nodeID == nodeID)
-            scheduler_activate(s, cj->hydro.super->stars.stars_in);
-
-          /* Do ci */
-          if (ci_active_stars) {
-
-            /* 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_hydro_sorts(cj, t->flags, s);
-            cell_activate_stars_sorts(ci, t->flags, s);
-          }
-
-          /* Do cj */
-          if (cj_active_stars) {
-
-            /* 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 (ci_nodeID == nodeID) cell_activate_drift_part(ci, s);
-            if (cj_nodeID == nodeID) cell_activate_drift_spart(cj, 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);
-          }
-        }
-
-        /* Store current values of dx_max and h_max. */
-        else if (t_type == task_type_sub_pair &&
-                 t_subtype == task_subtype_stars_density) {
-
-          /* Add stars_in dependencies for each cell that is part of
-           * a pair/sub_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)
-            scheduler_activate(s, cj->hydro.super->stars.stars_in);
-
-          cell_activate_subcell_stars_tasks(ci, cj, s, with_star_formation,
-                                            with_star_formation_sink,
-                                            with_timestep_sync);
-        }
-      }
-
-      /* Stars prep1 */
-      else if (t_subtype == task_subtype_stars_prep1) {
-
-        /* 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) &&
-            (ci_active_stars || cj_active_stars)) {
-
-          scheduler_activate(s, t);
-
-          /* If there are active sparts in ci, activate hydro ghost in cj */
-          if (ci_active_stars)
-            scheduler_activate(s, cj->hydro.super->hydro.prep1_ghost);
-          /* If there are active sparts in cj, activate hydro ghost in ci */
-          if (cj_active_stars)
-            scheduler_activate(s, ci->hydro.super->hydro.prep1_ghost);
-
-        } else if ((ci_nodeID == nodeID && cj_nodeID != nodeID) &&
-                   (cj_active_stars)) {
-
-          scheduler_activate(s, t);
-          /* If there are active sparts in cj, activate hydro ghost in ci */
-          scheduler_activate(s, ci->hydro.super->hydro.prep1_ghost);
-
-        } else if ((ci_nodeID != nodeID && cj_nodeID == nodeID) &&
-                   (ci_active_stars)) {
-
-          scheduler_activate(s, t);
-          /* If there are active sparts in ci, activate hydro ghost in cj */
-          scheduler_activate(s, cj->hydro.super->hydro.prep1_ghost);
-        }
-      }
-
-      /* Stars prep2 */
-      else if (t_subtype == task_subtype_stars_prep2) {
-
-        /* We only want to activate the task if the cell is active and is
-           going to update some sparts on the *local* node */
-        if ((ci_nodeID == nodeID && cj_nodeID == nodeID) &&
-            (ci_active_stars || cj_active_stars)) {
-
-          scheduler_activate(s, t);
-
-        } else if ((ci_nodeID == nodeID && cj_nodeID != nodeID) &&
-                   (ci_active_stars)) {
-
-          scheduler_activate(s, t);
-
-        } else if ((ci_nodeID != nodeID && cj_nodeID == nodeID) &&
-                   (cj_active_stars)) {
-
-          scheduler_activate(s, t);
-        }
-      }
-
-      /* Stars feedback */
-      else if (t_subtype == task_subtype_stars_feedback) {
-
-        /* 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) &&
-            (ci_active_stars || cj_active_stars)) {
-
-          scheduler_activate(s, t);
-
-        } else if ((ci_nodeID == nodeID && cj_nodeID != nodeID) &&
-                   (cj_active_stars)) {
-
-          scheduler_activate(s, t);
-
-        } else if ((ci_nodeID != nodeID && cj_nodeID == nodeID) &&
-                   (ci_active_stars)) {
-
-          scheduler_activate(s, t);
-        }
-
-        if (t->type == task_type_pair || t->type == task_type_sub_pair) {
-          /* Add stars_out dependencies for each cell that is part of
-           * a pair/sub_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)
-            scheduler_activate(s, cj->hydro.super->stars.stars_out);
-        }
-      }
-
-      /* Black_Holes density */
-      else if ((t_subtype == task_subtype_bh_density ||
-                t_subtype == task_subtype_bh_swallow ||
-                t_subtype == task_subtype_do_gas_swallow ||
-                t_subtype == task_subtype_do_bh_swallow ||
-                t_subtype == task_subtype_bh_feedback) &&
-               (ci_active_black_holes || cj_active_black_holes) &&
-               (ci_nodeID == nodeID || cj_nodeID == nodeID)) {
-
-        scheduler_activate(s, t);
-
-        /* Set the correct drifting flags */
-        if (t_type == task_type_pair && t_subtype == task_subtype_bh_density) {
-          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);
-
-          /* 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 || t_type == task_type_sub_pair) &&
-            t_subtype == task_subtype_bh_feedback) {
-          /* Add bh_out dependencies for each cell that is part of
-           * a pair/sub_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)
-            scheduler_activate(s, cj->hydro.super->black_holes.black_holes_out);
-        }
-
-        /* Store current values of dx_max and h_max. */
-        else if (t_type == task_type_sub_pair &&
-                 t_subtype == task_subtype_bh_density) {
-          cell_activate_subcell_black_holes_tasks(ci, cj, s,
-                                                  with_timestep_sync);
-          /* Activate bh_in for each cell that is part of
-           * a sub_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);
-        }
-      }
-
-      /* Gravity */
-      else if ((t_subtype == task_subtype_grav) &&
-               ((ci_active_gravity && ci_nodeID == nodeID) ||
-                (cj_active_gravity && cj_nodeID == nodeID))) {
-
-        scheduler_activate(s, t);
-
-        if (t_type == task_type_pair && t_subtype == task_subtype_grav) {
-          /* Activate the gravity drift */
-          cell_activate_subcell_grav_tasks(t->ci, t->cj, s);
-        }
-
-#ifdef SWIFT_DEBUG_CHECKS
-        else if (t_type == task_type_sub_pair &&
-                 t_subtype == task_subtype_grav) {
-          error("Invalid task sub-type encountered");
-        }
-#endif
-      }
-
-      /* Sink formation */
-      else if ((t_subtype == task_subtype_sink_swallow ||
-                t_subtype == task_subtype_sink_do_sink_swallow ||
-                t_subtype == task_subtype_sink_do_gas_swallow) &&
-               (ci_active_sinks || cj_active_sinks) &&
-               (ci_nodeID == nodeID || cj_nodeID == nodeID)) {
-
-        scheduler_activate(s, t);
-
-        /* Set the correct sorting flags */
-        if (t_type == task_type_pair &&
-            t_subtype == task_subtype_sink_swallow) {
-
-          /* Activate the sink drift for the sink merger */
-          if (ci_nodeID == nodeID) {
-            cell_activate_drift_sink(ci, s);
-            cell_activate_sink_formation_tasks(ci->top, s);
-            /* Activate all sink_in tasks for each cell involved
-             * in pair type tasks */
-            scheduler_activate(s, ci->hydro.super->sinks.sink_in);
-          }
-
-          if (cj_nodeID == nodeID) {
-            cell_activate_drift_sink(cj, s);
-            if (ci->top != cj->top) {
-              cell_activate_sink_formation_tasks(cj->top, s);
-            }
-            /* Activate all sink_in tasks for each cell involved
-             * in pair type tasks */
-            scheduler_activate(s, cj->hydro.super->sinks.sink_in);
-          }
-
-          /* Do ci */
-          if (ci_active_sinks) {
-
-            /* 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_sinks) {
-
-            /* 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. */
-            /* Activate the sink drift for the merger */
-            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);
-          }
-        }
-
-        else if (t_type == task_type_sub_pair &&
-                 t_subtype == task_subtype_sink_swallow) {
-          /* Activate all sink_in tasks for each cell involved
-           * in sub_pair type tasks */
-          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);
-
-          /* Store current values of dx_max and h_max. */
-          cell_activate_subcell_sinks_tasks(ci, cj, s, with_timestep_sync);
-        }
-
-        else if ((t_type == task_type_pair || t_type == task_type_sub_pair) &&
-                 t_subtype == task_subtype_sink_do_gas_swallow) {
-          /* Activate sinks_out for each cell that is part of
-           * a pair/sub_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)
-            scheduler_activate(s, cj->hydro.super->sinks.sink_out);
-        }
-      }
-
-      /* RT gradient and transport tasks */
-      else if (t_subtype == task_subtype_rt_gradient) {
-
-        /* 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 && ci_active_rt) ||
-            (cj_nodeID == nodeID && cj_active_rt)) {
-
-          scheduler_activate(s, t);
-
-          /* Set the correct sorting flags */
-          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;
-
-            /* 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) {
-            cell_activate_subcell_rt_tasks(ci, cj, s, /*sub_cycle=*/0);
-          }
-        }
-      }
-
-      else if (t_subtype == task_subtype_rt_transport) {
-        /* 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 && ci_active_rt) ||
-            (cj_nodeID == nodeID && cj_active_rt)) {
-
-          /* The gradient and transport task subtypes mirror the hydro tasks.
-           * Therefore all the (subcell) sorts and drifts should already have
-           * been activated properly in the hydro part of the activation. */
-          scheduler_activate(s, t);
-
-          if (t_type == task_type_pair || t_type == task_type_sub_pair) {
-            /* Activate transport_out for each cell that is part of
-             * a pair/sub_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)
-              scheduler_activate(s, cj->hydro.super->rt.rt_transport_out);
-          }
-        }
-      }
-
-      /* Pair tasks between inactive local cells and active remote cells. */
-      if ((ci_nodeID != nodeID && cj_nodeID == nodeID && ci_active_hydro &&
-           !cj_active_hydro) ||
-          (ci_nodeID == nodeID && cj_nodeID != nodeID && !ci_active_hydro &&
-           cj_active_hydro)) {
-
-#if defined(WITH_MPI) && defined(MPI_SYMMETRIC_FORCE_INTERACTION)
-        if (t_subtype == task_subtype_force) {
-
-          scheduler_activate(s, t);
-
-          /* Set the correct sorting flags */
-          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 hydro drift tasks. */
-            if (ci_nodeID == nodeID) cell_activate_drift_part(ci, s);
-            if (cj_nodeID == nodeID) cell_activate_drift_part(cj, s);
-
-            /* And the limiter */
-            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) {
-            cell_activate_subcell_hydro_tasks(t->ci, t->cj, s,
-                                              with_timestep_limiter);
-          }
-        }
-      }
-
-      /* Pair tasks between inactive local cells and active remote cells. */
-      if ((ci_nodeID != nodeID && cj_nodeID == nodeID && ci_active_rt &&
-           !cj_active_rt) ||
-          (ci_nodeID == nodeID && cj_nodeID != nodeID && !ci_active_rt &&
-           cj_active_rt)) {
-
-        if (t_subtype == task_subtype_rt_transport) {
-
-          scheduler_activate(s, t);
-
-          /* Set the correct sorting flags */
-          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;
-
-            /* 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) {
-            cell_activate_subcell_rt_tasks(ci, cj, s, /*sub_cycle=*/0);
-          }
-        }
-#endif
-      }
-
-      /* Only interested in density tasks as of here. */
-      if (t_subtype == task_subtype_density) {
-
-        /* Too much particle movement? */
-        if (cell_need_rebuild_for_hydro_pair(ci, cj)) *rebuild_space = 1;
-
-#ifdef WITH_MPI
-        /* Activate the send/recv tasks. */
-        if (ci_nodeID != nodeID) {
-
-          /* If the local cell is active, receive data from the foreign cell. */
-          if (cj_active_hydro) {
-            scheduler_activate_recv(s, ci->mpi.recv, task_subtype_xv);
-            if (ci_active_hydro) {
-              scheduler_activate_recv(s, ci->mpi.recv, task_subtype_rho);
-#ifdef EXTRA_HYDRO_LOOP
-              scheduler_activate_recv(s, ci->mpi.recv, task_subtype_gradient);
-#endif
-            }
-          }
-          /* 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. */
-          else if (ci_active_hydro) {
-#ifdef MPI_SYMMETRIC_FORCE_INTERACTION
-            /* NOTE: (yuyttenh, 09/2022) Since the particle communications send
-             * over whole particles currently, just activating the gradient
-             * send/recieve should be enough for now. The remote active
-             * particles are only needed for the sorts and the flux exchange on
-             * the node of the inactive cell, so sending over the xv and
-             * gradient suffices. If at any point the commutications change, we
-             * should probably also send over the rho separately. */
-            scheduler_activate_recv(s, ci->mpi.recv, task_subtype_xv);
-#ifndef EXTRA_HYDRO_LOOP
-            scheduler_activate_recv(s, ci->mpi.recv, task_subtype_rho);
-#else
-            scheduler_activate_recv(s, ci->mpi.recv, task_subtype_gradient);
-#endif
-#endif
-          }
-
-          /* If the foreign cell is active, we want its particles for the
-           * limiter */
-          if (ci_active_hydro && with_timestep_limiter) {
-            scheduler_activate_recv(s, ci->mpi.recv, task_subtype_limiter);
-            scheduler_activate_unpack(s, ci->mpi.unpack, task_subtype_limiter);
-          }
-
-          /* Is the foreign cell active and will need stuff from us? */
-          if (ci_active_hydro) {
-            struct link *l = scheduler_activate_send(
-                s, cj->mpi.send, task_subtype_xv, ci_nodeID);
-
-            /* Drift the cell which will be sent at the level at which it is
-               sent, i.e. drift the cell specified in the send task (l->t)
-               itself. */
-            cell_activate_drift_part(l->t->ci, s);
-
-            /* If the local cell is also active, more stuff will be needed. */
-            if (cj_active_hydro) {
-              scheduler_activate_send(s, cj->mpi.send, task_subtype_rho,
-                                      ci_nodeID);
-
-#ifdef EXTRA_HYDRO_LOOP
-              scheduler_activate_send(s, cj->mpi.send, task_subtype_gradient,
-                                      ci_nodeID);
-#endif
-            }
-          }
-          /* 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 */
-          else if (cj_active_hydro) {
-#ifdef MPI_SYMMETRIC_FORCE_INTERACTION
-            /* See NOTE on line 867 */
-            struct link *l = scheduler_activate_send(
-                s, cj->mpi.send, task_subtype_xv, ci_nodeID);
-            /* Drift the cell which will be sent at the level at which it is
-             * sent, i.e. drift the cell specified in the send task (l->t)
-             * itself. */
-            cell_activate_drift_part(l->t->ci, s);
-#ifndef EXTRA_HYDRO_LOOP
-            scheduler_activate_send(s, cj->mpi.send, task_subtype_rho,
-                                    ci_nodeID);
-#else
-            scheduler_activate_send(s, cj->mpi.send, task_subtype_gradient,
-                                    ci_nodeID);
-#endif
-#endif
-          }
-
-          /* If the local cell is active, send its particles for the limiting.
-           */
-          if (cj_active_hydro && with_timestep_limiter) {
-            scheduler_activate_send(s, cj->mpi.send, task_subtype_limiter,
-                                    ci_nodeID);
-            scheduler_activate_pack(s, cj->mpi.pack, task_subtype_limiter,
-                                    ci_nodeID);
-          }
-
-          /* Propagating new star counts? */
-          if (with_star_formation_sink) error("TODO");
-          if (with_star_formation && with_feedback) {
-            if (ci_active_hydro && ci->hydro.count > 0) {
-              scheduler_activate_recv(s, ci->mpi.recv, task_subtype_sf_counts);
-            }
-            if (cj_active_hydro && cj->hydro.count > 0) {
-              scheduler_activate_send(s, cj->mpi.send, task_subtype_sf_counts,
-                                      ci_nodeID);
-            }
-          }
-
-        } else if (cj_nodeID != nodeID) {
-
-          /* If the local cell is active, receive data from the foreign cell. */
-          if (ci_active_hydro) {
-
-            scheduler_activate_recv(s, cj->mpi.recv, task_subtype_xv);
-            if (cj_active_hydro) {
-              scheduler_activate_recv(s, cj->mpi.recv, task_subtype_rho);
-#ifdef EXTRA_HYDRO_LOOP
-              scheduler_activate_recv(s, cj->mpi.recv, task_subtype_gradient);
-#endif
-            }
-          }
-          /* 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. */
-          else if (cj_active_hydro) {
-#ifdef MPI_SYMMETRIC_FORCE_INTERACTION
-            /* See NOTE on line 867. */
-            scheduler_activate_recv(s, cj->mpi.recv, task_subtype_xv);
-#ifndef EXTRA_HYDRO_LOOP
-            scheduler_activate_recv(s, cj->mpi.recv, task_subtype_rho);
-#else
-            scheduler_activate_recv(s, cj->mpi.recv, task_subtype_gradient);
-#endif
-#endif
-          }
-
-          /* If the foreign cell is active, we want its particles for the
-           * limiter */
-          if (cj_active_hydro && with_timestep_limiter) {
-            scheduler_activate_recv(s, cj->mpi.recv, task_subtype_limiter);
-            scheduler_activate_unpack(s, cj->mpi.unpack, task_subtype_limiter);
-          }
-
-          /* Is the foreign cell active and will need stuff from us? */
-          if (cj_active_hydro) {
-
-            struct link *l = scheduler_activate_send(
-                s, ci->mpi.send, task_subtype_xv, cj_nodeID);
-
-            /* Drift the cell which will be sent at the level at which it is
-               sent, i.e. drift the cell specified in the send task (l->t)
-               itself. */
-            cell_activate_drift_part(l->t->ci, s);
-
-            /* If the local cell is also active, more stuff will be needed. */
-            if (ci_active_hydro) {
-
-              scheduler_activate_send(s, ci->mpi.send, task_subtype_rho,
-                                      cj_nodeID);
-
-#ifdef EXTRA_HYDRO_LOOP
-              scheduler_activate_send(s, ci->mpi.send, task_subtype_gradient,
-                                      cj_nodeID);
-#endif
-            }
-          }
-          /* 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 */
-          else if (ci_active_hydro) {
-#ifdef MPI_SYMMETRIC_FORCE_INTERACTION
-            /* See NOTE on line 867. */
-            struct link *l = scheduler_activate_send(
-                s, ci->mpi.send, task_subtype_xv, cj_nodeID);
-            /* Drift the cell which will be sent at the level at which it is
-             * sent, i.e. drift the cell specified in the send task (l->t)
-             * itself. */
-            cell_activate_drift_part(l->t->ci, s);
-#ifndef EXTRA_HYDRO_LOOP
-            scheduler_activate_send(s, ci->mpi.send, task_subtype_rho,
-                                    cj_nodeID);
-#else
-            scheduler_activate_send(s, ci->mpi.send, task_subtype_gradient,
-                                    cj_nodeID);
-#endif
-#endif
-          }
-
-          /* If the local cell is active, send its particles for the limiting.
-           */
-          if (ci_active_hydro && with_timestep_limiter) {
-            scheduler_activate_send(s, ci->mpi.send, task_subtype_limiter,
-                                    cj_nodeID);
-            scheduler_activate_pack(s, ci->mpi.pack, task_subtype_limiter,
-                                    cj_nodeID);
-          }
-
-          /* Propagating new star counts? */
-          if (with_star_formation_sink) error("TODO");
-          if (with_star_formation && with_feedback) {
-            if (cj_active_hydro && cj->hydro.count > 0) {
-              scheduler_activate_recv(s, cj->mpi.recv, task_subtype_sf_counts);
-            }
-            if (ci_active_hydro && ci->hydro.count > 0) {
-              scheduler_activate_send(s, ci->mpi.send, task_subtype_sf_counts,
-                                      cj_nodeID);
-            }
-          }
-        }
-#endif
-      }
-
-      /* Only interested in stars_density tasks as of here. */
-      else if (t->subtype == task_subtype_stars_density) {
-
-        /* Too much particle movement? */
-        if (cell_need_rebuild_for_stars_pair(ci, cj)) *rebuild_space = 1;
-        if (cell_need_rebuild_for_stars_pair(cj, ci)) *rebuild_space = 1;
-
-#ifdef WITH_MPI
-        /* Activate the send/recv tasks. */
-        if (ci_nodeID != nodeID) {
-
-          if (cj_active_stars) {
-            scheduler_activate_recv(s, ci->mpi.recv, task_subtype_xv);
-            scheduler_activate_recv(s, ci->mpi.recv, task_subtype_rho);
-#ifdef EXTRA_STAR_LOOPS
-            scheduler_activate_recv(s, ci->mpi.recv, task_subtype_part_prep1);
-#endif
-
-            /* If the local cell is active, more stuff will be needed. */
-            scheduler_activate_send(s, cj->mpi.send, task_subtype_spart_density,
-                                    ci_nodeID);
-#ifdef EXTRA_STAR_LOOPS
-            scheduler_activate_send(s, cj->mpi.send, task_subtype_spart_prep2,
-                                    ci_nodeID);
-#endif
-            cell_activate_drift_spart(cj, s);
-          }
-
-          if (ci_active_stars) {
-            scheduler_activate_recv(s, ci->mpi.recv,
-                                    task_subtype_spart_density);
-#ifdef EXTRA_STAR_LOOPS
-            scheduler_activate_recv(s, ci->mpi.recv, task_subtype_spart_prep2);
-#endif
-
-            /* Is the foreign cell active and will need stuff from us? */
-            scheduler_activate_send(s, cj->mpi.send, task_subtype_xv,
-                                    ci_nodeID);
-            scheduler_activate_send(s, cj->mpi.send, task_subtype_rho,
-                                    ci_nodeID);
-#ifdef EXTRA_STAR_LOOPS
-            scheduler_activate_send(s, cj->mpi.send, task_subtype_part_prep1,
-                                    ci_nodeID);
-#endif
-
-            /* 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);
-          }
-
-        } else if (cj_nodeID != nodeID) {
-
-          /* If the local cell is active, receive data from the foreign cell. */
-          if (ci_active_stars) {
-            scheduler_activate_recv(s, cj->mpi.recv, task_subtype_xv);
-            scheduler_activate_recv(s, cj->mpi.recv, task_subtype_rho);
-#ifdef EXTRA_STAR_LOOPS
-            scheduler_activate_recv(s, cj->mpi.recv, task_subtype_part_prep1);
-#endif
-
-            /* If the local cell is active, more stuff will be needed. */
-            scheduler_activate_send(s, ci->mpi.send, task_subtype_spart_density,
-                                    cj_nodeID);
-#ifdef EXTRA_STAR_LOOPS
-            scheduler_activate_send(s, ci->mpi.send, task_subtype_spart_prep2,
-                                    cj_nodeID);
-#endif
-            cell_activate_drift_spart(ci, s);
-          }
-
-          if (cj_active_stars) {
-            scheduler_activate_recv(s, cj->mpi.recv,
-                                    task_subtype_spart_density);
-#ifdef EXTRA_STAR_LOOPS
-            scheduler_activate_recv(s, cj->mpi.recv, task_subtype_spart_prep2);
-#endif
-
-            /* Is the foreign cell active and will need stuff from us? */
-            scheduler_activate_send(s, ci->mpi.send, task_subtype_xv,
-                                    cj_nodeID);
-            scheduler_activate_send(s, ci->mpi.send, task_subtype_rho,
-                                    cj_nodeID);
-#ifdef EXTRA_STAR_LOOPS
-            scheduler_activate_send(s, ci->mpi.send, task_subtype_part_prep1,
-                                    cj_nodeID);
-#endif
-
-            /* 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);
-          }
-        }
-#endif
-      }
-
-      /* Only interested in sink_swallow tasks as of here. */
-      else if (t->subtype == task_subtype_sink_swallow) {
-
-        /* Too much particle movement? */
-        if (cell_need_rebuild_for_sinks_pair(ci, cj)) *rebuild_space = 1;
-        if (cell_need_rebuild_for_sinks_pair(cj, ci)) *rebuild_space = 1;
-
-#ifdef WITH_MPI
-        error("TODO");
-#endif
-      }
-
-      /* Only interested in black hole density tasks as of here. */
-      else if (t->subtype == task_subtype_bh_density) {
-
-        /* Too much particle movement? */
-        if (cell_need_rebuild_for_black_holes_pair(ci, cj)) *rebuild_space = 1;
-        if (cell_need_rebuild_for_black_holes_pair(cj, ci)) *rebuild_space = 1;
-
-        scheduler_activate(s, ci->hydro.super->black_holes.swallow_ghost_0);
-        scheduler_activate(s, cj->hydro.super->black_holes.swallow_ghost_0);
-
-#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);
-          scheduler_activate_recv(s, ci->mpi.recv, task_subtype_bpart_merger);
-
-          /* 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 */
-          cell_activate_drift_bpart(cj, s);
-
-          /* 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);
-
-          /* 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. */
-          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);
-          scheduler_activate_recv(s, cj->mpi.recv, task_subtype_bpart_merger);
-
-          /* 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 */
-          cell_activate_drift_bpart(ci, s);
-
-          /* 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);
-
-          /* 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. */
-          cell_activate_drift_part(ci, s);
-        }
-#endif
-      }
-
-      /* Only interested in gravity tasks as of here. */
-      else if (t_subtype == task_subtype_grav) {
-
-#ifdef WITH_MPI
-        /* Activate the send/recv tasks. */
-        if (ci_nodeID != nodeID) {
-
-          /* If the local cell is active, receive data from the foreign cell. */
-          if (cj_active_gravity)
-            scheduler_activate_recv(s, ci->mpi.recv, task_subtype_gpart);
-
-          /* Is the foreign cell active and will need stuff from us? */
-          if (ci_active_gravity) {
-
-            struct link *l = scheduler_activate_send(
-                s, cj->mpi.send, task_subtype_gpart, ci_nodeID);
-
-            /* Drift the cell which will be sent at the level at which it is
-               sent, i.e. drift the cell specified in the send task (l->t)
-               itself. */
-            cell_activate_drift_gpart(l->t->ci, s);
-          }
-
-        } else if (cj_nodeID != nodeID) {
-
-          /* If the local cell is active, receive data from the foreign cell. */
-          if (ci_active_gravity)
-            scheduler_activate_recv(s, cj->mpi.recv, task_subtype_gpart);
-
-          /* Is the foreign cell active and will need stuff from us? */
-          if (cj_active_gravity) {
-
-            struct link *l = scheduler_activate_send(
-                s, ci->mpi.send, task_subtype_gpart, cj_nodeID);
-
-            /* Drift the cell which will be sent at the level at which it is
-               sent, i.e. drift the cell specified in the send task (l->t)
-               itself. */
-            cell_activate_drift_gpart(l->t->ci, s);
-          }
-        }
-#endif
-      } /* Only interested in RT tasks as of here. */
-
-      else if (t->subtype == task_subtype_rt_gradient) {
-
-#ifdef WITH_MPI
-        /* Activate the send/recv tasks. */
-        if (ci_nodeID != nodeID) {
-
-          /* If the local cell is active, receive data from the foreign cell. */
-          if (cj_active_rt) {
-
-            scheduler_activate_recv(s, ci->mpi.recv, task_subtype_rt_gradient);
-
-            if (ci_active_rt) {
-              /* We only need updates later on if the other cell is active too
-               */
-              scheduler_activate_recv(s, ci->mpi.recv,
-                                      task_subtype_rt_transport);
-            }
-          } else if (ci_active_rt) {
-#ifdef MPI_SYMMETRIC_FORCE_INTERACTION
-            /* 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. */
-            scheduler_activate_recv(s, ci->mpi.recv, task_subtype_rt_gradient);
-            scheduler_activate_recv(s, ci->mpi.recv, task_subtype_rt_transport);
-#endif
-          }
-
-          /* Is the foreign cell active and will need stuff from us? */
-          if (ci_active_rt) {
-
-            scheduler_activate_send(s, cj->mpi.send, task_subtype_rt_gradient,
-                                    ci_nodeID);
-
-            if (cj_active_rt) {
-              scheduler_activate_send(s, cj->mpi.send,
-                                      task_subtype_rt_transport, ci_nodeID);
-            }
-          } else if (cj_active_rt) {
-#ifdef MPI_SYMMETRIC_FORCE_INTERACTION
-            /* 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 */
-            scheduler_activate_send(s, cj->mpi.send, task_subtype_rt_gradient,
-                                    ci_nodeID);
-            scheduler_activate_send(s, cj->mpi.send, task_subtype_rt_transport,
-                                    ci_nodeID);
-#endif
-          }
-
-        } else if (cj_nodeID != nodeID) {
-
-          /* If the local cell is active, receive data from the foreign cell. */
-          if (ci_active_rt) {
-
-            scheduler_activate_recv(s, cj->mpi.recv, task_subtype_rt_gradient);
-
-            if (cj_active_rt) {
-              /* We only need updates later on if the other cell is active too
-               */
-              scheduler_activate_recv(s, cj->mpi.recv,
-                                      task_subtype_rt_transport);
-            }
-          } else if (cj_active_rt) {
-#ifdef MPI_SYMMETRIC_FORCE_INTERACTION
-            /* 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. */
-            scheduler_activate_recv(s, cj->mpi.recv, task_subtype_rt_gradient);
-            scheduler_activate_recv(s, cj->mpi.recv, task_subtype_rt_transport);
-#endif
-          }
-
-          /* Is the foreign cell active and will need stuff from us? */
-          if (cj_active_rt) {
-
-            scheduler_activate_send(s, ci->mpi.send, task_subtype_rt_gradient,
-                                    cj_nodeID);
-
-            if (ci_active_rt) {
-              /* We only need updates later on if the other cell is active too
-               */
-              scheduler_activate_send(s, ci->mpi.send,
-                                      task_subtype_rt_transport, cj_nodeID);
-            }
-          } else if (ci_active_rt) {
-#ifdef MPI_SYMMETRIC_FORCE_INTERACTION
-            /* 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 */
-            scheduler_activate_send(s, ci->mpi.send, task_subtype_rt_gradient,
-                                    cj_nodeID);
-            scheduler_activate_send(s, ci->mpi.send, task_subtype_rt_transport,
-                                    cj_nodeID);
-#endif
-          }
-        }
-#endif
-      }
-    }
-
-    /* End force for hydro ? */
-    else if (t_type == task_type_end_hydro_force) {
-
-      if (cell_is_active_hydro(t->ci, e)) scheduler_activate(s, t);
-    }
-
-    /* End force for gravity ? */
-    else if (t_type == task_type_end_grav_force) {
-
-      if (cell_is_active_gravity(t->ci, e)) scheduler_activate(s, t);
-    }
-
-    /* Activate the weighting task for neutrinos */
-    else if (t_type == task_type_neutrino_weight) {
-      if (cell_is_active_gravity(t->ci, e)) {
-        scheduler_activate(s, t);
-      }
-    }
-
-    /* Kick ? */
-    else if (t_type == task_type_kick1 || t_type == task_type_kick2) {
-
-      if (cell_is_active_hydro(t->ci, e) || cell_is_active_gravity(t->ci, e) ||
-          cell_is_active_stars(t->ci, e) || cell_is_active_sinks(t->ci, e) ||
-          cell_is_active_black_holes(t->ci, e))
-        scheduler_activate(s, t);
-    }
-
-    /* Hydro ghost tasks ? */
-    else if (t_type == task_type_ghost || t_type == task_type_extra_ghost ||
-             t_type == task_type_ghost_in || t_type == task_type_ghost_out) {
-      if (cell_is_active_hydro(t->ci, e)) scheduler_activate(s, t);
-    }
-
-    /* csds tasks ? */
-    else if (t->type == task_type_csds) {
-      if (cell_is_active_hydro(t->ci, e) || cell_is_active_gravity(t->ci, e) ||
-          cell_is_active_stars(t->ci, e))
-        scheduler_activate(s, t);
-    }
-
-    /* Gravity stuff ? */
-    else if (t_type == task_type_grav_down ||
-             t_type == task_type_grav_long_range ||
-             t_type == task_type_init_grav ||
-             t_type == task_type_init_grav_out ||
-             t_type == task_type_drift_gpart_out ||
-             t_type == task_type_grav_down_in) {
-      if (cell_is_active_gravity(t->ci, e)) scheduler_activate(s, t);
-    }
-
-    /* Multipole - Multipole interaction task */
-    else if (t_type == task_type_grav_mm) {
-
-      /* Local pointers. */
-      const struct cell *ci = t->ci;
-      const 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_gravity = cell_is_active_gravity_mm(ci, e);
-      const int cj_active_gravity = cell_is_active_gravity_mm(cj, e);
-
-      if ((ci_active_gravity && ci_nodeID == nodeID) ||
-          (cj_active_gravity && cj_nodeID == nodeID))
-        scheduler_activate(s, t);
-    }
-
-    /* Star ghost tasks ? */
-    else if (t_type == task_type_stars_ghost ||
-             t_type == task_type_stars_prep_ghost1 ||
-             t_type == task_type_hydro_prep_ghost1 ||
-             t_type == task_type_stars_prep_ghost2) {
-      if (cell_need_activating_stars(t->ci, e, with_star_formation,
-                                     with_star_formation_sink))
-        scheduler_activate(s, t);
-    }
-
-    /* Feedback implicit tasks? */
-    else if (t_type == task_type_stars_in || t_type == task_type_stars_out) {
-      if (cell_need_activating_stars(t->ci, e, with_star_formation,
-                                     with_star_formation_sink))
-        scheduler_activate(s, t);
-    }
-
-    /* Sink implicit tasks? */
-    else if (t_type == task_type_sink_in || t_type == task_type_sink_out ||
-             t_type == task_type_sink_ghost1) {
-      if (cell_is_active_sinks(t->ci, e) || cell_is_active_hydro(t->ci, e))
-        scheduler_activate(s, t);
-    }
-
-    /* Black hole ghost tasks ? */
-    else if (t_type == task_type_bh_density_ghost ||
-             t_type == task_type_bh_swallow_ghost1 ||
-             t_type == task_type_bh_swallow_ghost2 ||
-             t_type == task_type_bh_swallow_ghost3) {
-      if (cell_is_active_black_holes(t->ci, e)) scheduler_activate(s, t);
-    }
-
-    /* Black holes implicit tasks? */
-    else if (t_type == task_type_bh_in || t_type == task_type_bh_out) {
-      if (cell_is_active_black_holes(t->ci, e)) scheduler_activate(s, t);
-    }
-
-    /* Time-step collection? */
-    else if (t_type == task_type_timestep) {
-      t->ci->hydro.updated = 0;
-      t->ci->grav.updated = 0;
-      t->ci->stars.updated = 0;
-      t->ci->sinks.updated = 0;
-      t->ci->black_holes.updated = 0;
-
-      if (!cell_is_empty(t->ci)) {
-        if (cell_is_active_hydro(t->ci, e) ||
-            cell_is_active_gravity(t->ci, e) ||
-            cell_is_active_stars(t->ci, e) || cell_is_active_sinks(t->ci, e) ||
-            cell_is_active_black_holes(t->ci, e))
-          scheduler_activate(s, t);
-      }
-    }
-
-    /* Time-step collection? */
-    else if (t_type == task_type_collect) {
-      t->ci->hydro.updated = 0;
-      t->ci->grav.updated = 0;
-      t->ci->stars.updated = 0;
-      t->ci->sinks.updated = 0;
-      t->ci->black_holes.updated = 0;
-      t->ci->rt.updated = 0; /* this is different from time-step */
-
-      if (!cell_is_empty(t->ci)) {
-        if (cell_is_active_hydro(t->ci, e) ||
-            cell_is_active_gravity(t->ci, e) ||
-            cell_is_active_stars(t->ci, e) || cell_is_active_sinks(t->ci, e) ||
-            cell_is_active_black_holes(t->ci, e) || cell_is_rt_active(t->ci, e))
-          /* this is different from time-step ----^*/
-          scheduler_activate(s, t);
-      }
-    }
-
-    else if ((t_type == task_type_send && t_subtype == task_subtype_tend) ||
-             (t_type == task_type_recv && t_subtype == task_subtype_tend)) {
-      if (!cell_is_empty(t->ci)) {
-        scheduler_activate(s, t);
-      }
-    }
-
-    /* Subgrid tasks: cooling */
-    else if (t_type == task_type_cooling || t_type == task_type_cooling_in ||
-             t_type == task_type_cooling_out) {
-      if (cell_is_active_hydro(t->ci, e)) scheduler_activate(s, t);
-    }
-
-    /* Subgrid tasks: star formation */
-    else if (t_type == task_type_star_formation) {
-      if (cell_is_active_hydro(t->ci, e)) {
-        cell_activate_star_formation_tasks(t->ci, s, with_feedback);
-        cell_activate_super_spart_drifts(t->ci, s);
-      }
-    }
-
-    /* Subgrid tasks: star formation from sinks */
-    else if (t_type == task_type_star_formation_sink) {
-      if (cell_is_active_hydro(t->ci, e) || cell_is_active_sinks(t->ci, e)) {
-        cell_activate_star_formation_sink_tasks(t->ci, s, with_feedback);
-        cell_activate_super_spart_drifts(t->ci, s);
-      }
-    }
-
-    /* Radiative transfer implicit tasks */
-    else if (t->type == task_type_rt_in) {
-      if (cell_is_rt_active(t->ci, e) ||
-          cell_need_activating_stars(t->ci, e, with_star_formation,
-                                     with_star_formation_sink))
-        scheduler_activate(s, t);
-    }
-
-    else if (t->type == task_type_rt_ghost1 || t->type == task_type_rt_ghost2 ||
-             t->type == task_type_rt_transport_out ||
-             t->type == task_type_rt_tchem ||
-             t->type == task_type_rt_advance_cell_time ||
-             t->type == task_type_rt_out) {
-      if (cell_is_rt_active(t->ci, e)) scheduler_activate(s, t);
-      /* Note that rt_collect_times never needs to be active on main steps,
-       * which is always what follows engine_marktasks().*/
-    }
-
-    /* Subgrid tasks: sink formation */
-    else if (t_type == task_type_sink_formation) {
-      if (with_star_formation_sink && t->ci->hydro.count > 0 &&
-          cell_is_active_hydro(t->ci, e)) {
-        cell_activate_sink_formation_tasks(t->ci, s);
-        cell_activate_super_sink_drifts(t->ci, s);
-      }
-    }
-  }
-}
-
-/**
- * @brief Mark tasks to be un-skipped and set the sort flags accordingly.
- *
- * @return 1 if the space has to be rebuilt, 0 otherwise.
- */
-int engine_marktasks(struct engine *e) {
-
-  struct scheduler *s = &e->sched;
-  const ticks tic = getticks();
-  int rebuild_space = 0;
-
-  /* Run through the tasks and mark as skip or not. */
-  size_t extra_data[3] = {(size_t)e, (size_t)rebuild_space, (size_t)&e->sched};
-  threadpool_map(&e->threadpool, engine_marktasks_mapper, s->tasks, s->nr_tasks,
-                 sizeof(struct task), threadpool_auto_chunk_size, extra_data);
-  rebuild_space = extra_data[1];
-
-  if (e->verbose)
-    message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
-            clocks_getunit());
-
-  /* All is well... */
-  return rebuild_space;
-}

src/engine_redistribute.c

@@ -264,8 +264,8 @@ struct redist_mapper_data {
     int *dest =                                                            \
         mydata->dest + (ptrdiff_t)(parts - (struct TYPE *)mydata->base);   \
     int *lcounts = NULL;                                                   \
-    if ((lcounts = (int *)calloc(                                          \
-             sizeof(int), mydata->nr_nodes * mydata->nr_nodes)) == NULL)   \
+    if ((lcounts = (int *)calloc(mydata->nr_nodes * mydata->nr_nodes,      \
+                                 sizeof(int))) == NULL)                    \
       error("Failed to allocate counts thread-specific buffer");           \
     for (int k = 0; k < num_elements; k++) {                               \
       for (int j = 0; j < 3; j++) {                                        \
@@ -319,6 +319,14 @@ void ENGINE_REDISTRIBUTE_DEST_MAPPER(gpart);
  */
 void ENGINE_REDISTRIBUTE_DEST_MAPPER(bpart);
 
+/**
+ * @brief Accumulate the counts of sink particles per cell.
+ * Threadpool helper for accumulating the counts of particles per cell.
+ *
+ * sink version.
+ */
+void ENGINE_REDISTRIBUTE_DEST_MAPPER(sink);
+
 #endif /* redist_mapper_data */
 
 #ifdef WITH_MPI /* savelink_mapper_data */
@@ -399,12 +407,22 @@ void ENGINE_REDISTRIBUTE_SAVELINK_MAPPER(bpart, 1);
 void ENGINE_REDISTRIBUTE_SAVELINK_MAPPER(bpart, 0);
 #endif
 
+/**
+ * @brief Save position of sink-gpart links.
+ * Threadpool helper for accumulating the counts of particles per cell.
+ */
+#ifdef SWIFT_DEBUG_CHECKS
+void ENGINE_REDISTRIBUTE_SAVELINK_MAPPER(sink, 1);
+#else
+void ENGINE_REDISTRIBUTE_SAVELINK_MAPPER(sink, 0);
+#endif
+
 #endif /* savelink_mapper_data */
 
 #ifdef WITH_MPI /* relink_mapper_data */
 
-/* Support for relinking parts, gparts, sparts and bparts after moving between
- * nodes. */
+/* Support for relinking parts, gparts, sparts, bparts and sinks after moving
+ * between nodes. */
 struct relink_mapper_data {
   int nodeID;
   int nr_nodes;
@@ -412,6 +430,7 @@ struct relink_mapper_data {
   int *s_counts;
   int *g_counts;
   int *b_counts;
+  int *sink_counts;
   struct space *s;
 };
 
@@ -435,6 +454,7 @@ void engine_redistribute_relink_mapper(void *map_data, int num_elements,
   int *g_counts = mydata->g_counts;
   int *s_counts = mydata->s_counts;
   int *b_counts = mydata->b_counts;
+  int *sink_counts = mydata->sink_counts;
   struct space *s = mydata->s;
 
   for (int i = 0; i < num_elements; i++) {
@@ -446,12 +466,14 @@ void engine_redistribute_relink_mapper(void *map_data, int num_elements,
     size_t offset_gparts = 0;
     size_t offset_sparts = 0;
     size_t offset_bparts = 0;
+    size_t offset_sinks = 0;
     for (int n = 0; n < node; n++) {
       int ind_recv = n * nr_nodes + nodeID;
       offset_parts += counts[ind_recv];
       offset_gparts += g_counts[ind_recv];
       offset_sparts += s_counts[ind_recv];
       offset_bparts += b_counts[ind_recv];
+      offset_sinks += sink_counts[ind_recv];
     }
 
     /* Number of gparts sent from this node. */
@@ -493,6 +515,17 @@ void engine_redistribute_relink_mapper(void *map_data, int num_elements,
         s->gparts[k].id_or_neg_offset = -partner_index;
         s->bparts[partner_index].gpart = &s->gparts[k];
       }
+
+      /* Does this gpart have a sink partner ? */
+      else if (s->gparts[k].type == swift_type_sink) {
+
+        const ptrdiff_t partner_index =
+            offset_sinks - s->gparts[k].id_or_neg_offset;
+
+        /* Re-link */
+        s->gparts[k].id_or_neg_offset = -partner_index;
+        s->sinks[partner_index].gpart = &s->gparts[k];
+      }
     }
   }
 }
@@ -519,13 +552,6 @@ void engine_redistribute_relink_mapper(void *map_data, int num_elements,
 void engine_redistribute(struct engine *e) {
 
 #ifdef WITH_MPI
-#ifdef SWIFT_DEBUG_CHECKS
-  const int nr_sinks_new = 0;
-#endif
-  if (e->policy & engine_policy_sinks) {
-    error("Not implemented yet");
-  }
-
   const int nr_nodes = e->nr_nodes;
   const int nodeID = e->nodeID;
   struct space *s = e->s;
@@ -536,12 +562,14 @@ void engine_redistribute(struct engine *e) {
   struct gpart *gparts = s->gparts;
   struct spart *sparts = s->sparts;
   struct bpart *bparts = s->bparts;
+  struct sink *sinks = s->sinks;
   ticks tic = getticks();
 
   size_t nr_parts = s->nr_parts;
   size_t nr_gparts = s->nr_gparts;
   size_t nr_sparts = s->nr_sparts;
   size_t nr_bparts = s->nr_bparts;
+  size_t nr_sinks = s->nr_sinks;
 
   /* Start by moving inhibited particles to the end of the arrays */
   for (size_t k = 0; k < nr_parts; /* void */) {
@@ -609,6 +637,27 @@ void engine_redistribute(struct engine *e) {
     }
   }
 
+  /* Now move inhibited sink particles to the end of the arrays */
+  for (size_t k = 0; k < nr_sinks; /* void */) {
+    if (sinks[k].time_bin == time_bin_inhibited ||
+        sinks[k].time_bin == time_bin_not_created) {
+      nr_sinks -= 1;
+
+      /* Swap the particle */
+      memswap(&s->sinks[k], &s->sinks[nr_sinks], sizeof(struct sink));
+
+      /* Swap the link with the gpart */
+      if (s->sinks[k].gpart != NULL) {
+        s->sinks[k].gpart->id_or_neg_offset = -k;
+      }
+      if (s->sinks[nr_sinks].gpart != NULL) {
+        s->sinks[nr_sinks].gpart->id_or_neg_offset = -nr_sinks;
+      }
+    } else {
+      k++;
+    }
+  }
+
   /* Finally do the same with the gravity particles */
   for (size_t k = 0; k < nr_gparts; /* void */) {
     if (gparts[k].time_bin == time_bin_inhibited ||
@@ -626,6 +675,8 @@ void engine_redistribute(struct engine *e) {
         s->sparts[-s->gparts[k].id_or_neg_offset].gpart = &s->gparts[k];
       } else if (s->gparts[k].type == swift_type_black_hole) {
         s->bparts[-s->gparts[k].id_or_neg_offset].gpart = &s->gparts[k];
+      } else if (s->gparts[k].type == swift_type_sink) {
+        s->sinks[-s->gparts[k].id_or_neg_offset].gpart = &s->gparts[k];
       }
 
       if (s->gparts[nr_gparts].type == swift_type_gas) {
@@ -637,6 +688,9 @@ void engine_redistribute(struct engine *e) {
       } else if (s->gparts[nr_gparts].type == swift_type_black_hole) {
         s->bparts[-s->gparts[nr_gparts].id_or_neg_offset].gpart =
             &s->gparts[nr_gparts];
+      } else if (s->gparts[nr_gparts].type == swift_type_sink) {
+        s->sinks[-s->gparts[nr_sinks].id_or_neg_offset].gpart =
+            &s->gparts[nr_gparts];
       }
     } else {
       k++;
@@ -648,7 +702,7 @@ void engine_redistribute(struct engine *e) {
   /* Allocate temporary arrays to store the counts of particles to be sent
    * and the destination of each particle */
   int *counts;
-  if ((counts = (int *)calloc(sizeof(int), nr_nodes * nr_nodes)) == NULL)
+  if ((counts = (int *)calloc(nr_nodes * nr_nodes, sizeof(int))) == NULL)
     error("Failed to allocate counts temporary buffer.");
 
   int *dest;
@@ -727,7 +781,7 @@ void engine_redistribute(struct engine *e) {
 
   /* Get destination of each s-particle */
   int *s_counts;
-  if ((s_counts = (int *)calloc(sizeof(int), nr_nodes * nr_nodes)) == NULL)
+  if ((s_counts = (int *)calloc(nr_nodes * nr_nodes, sizeof(int))) == NULL)
     error("Failed to allocate s_counts temporary buffer.");
 
   int *s_dest;
@@ -793,7 +847,7 @@ void engine_redistribute(struct engine *e) {
 
   /* Get destination of each b-particle */
   int *b_counts;
-  if ((b_counts = (int *)calloc(sizeof(int), nr_nodes * nr_nodes)) == NULL)
+  if ((b_counts = (int *)calloc(nr_nodes * nr_nodes, sizeof(int))) == NULL)
     error("Failed to allocate b_counts temporary buffer.");
 
   int *b_dest;
@@ -857,9 +911,76 @@ void engine_redistribute(struct engine *e) {
   }
   swift_free("b_dest", b_dest);
 
+  /* Get destination of each sink-particle */
+  int *sink_counts;
+  if ((sink_counts = (int *)calloc(nr_nodes * nr_nodes, sizeof(int))) == NULL)
+    error("Failed to allocate sink_counts temporary buffer.");
+
+  int *sink_dest;
+  if ((sink_dest = (int *)swift_malloc("sink_dest", sizeof(int) * nr_sinks)) ==
+      NULL)
+    error("Failed to allocate sink_dest temporary buffer.");
+
+  redist_data.counts = sink_counts;
+  redist_data.dest = sink_dest;
+  redist_data.base = (void *)sinks;
+
+  threadpool_map(&e->threadpool, engine_redistribute_dest_mapper_sink, sinks,
+                 nr_sinks, sizeof(struct sink), threadpool_auto_chunk_size,
+                 &redist_data);
+
+  /* Sort the particles according to their cell index. */
+  if (nr_sinks > 0)
+    space_sinks_sort(s->sinks, sink_dest, &sink_counts[nodeID * nr_nodes],
+                     nr_nodes, 0);
+
+#ifdef SWIFT_DEBUG_CHECKS
+  /* Verify that the sink have been sorted correctly. */
+  for (size_t k = 0; k < nr_sinks; k++) {
+    const struct sink *sink = &s->sinks[k];
+
+    if (sink->time_bin == time_bin_inhibited)
+      error("Inhibited particle found after sorting!");
+
+    if (sink->time_bin == time_bin_not_created)
+      error("Inhibited particle found after sorting!");
+
+    /* New cell index */
+    const int new_cid =
+        cell_getid(s->cdim, sink->x[0] * s->iwidth[0],
+                   sink->x[1] * s->iwidth[1], sink->x[2] * s->iwidth[2]);
+
+    /* New cell of this sink */
+    const struct cell *c = &s->cells_top[new_cid];
+    const int new_node = c->nodeID;
+
+    if (sink_dest[k] != new_node)
+      error("sink's new node index not matching sorted index.");
+
+    if (sink->x[0] < c->loc[0] || sink->x[0] > c->loc[0] + c->width[0] ||
+        sink->x[1] < c->loc[1] || sink->x[1] > c->loc[1] + c->width[1] ||
+        sink->x[2] < c->loc[2] || sink->x[2] > c->loc[2] + c->width[2])
+      error("sink not sorted into the right top-level cell!");
+  }
+#endif
+
+  /* We need to re-link the gpart partners of sinks. */
+  if (nr_sinks > 0) {
+
+    struct savelink_mapper_data savelink_data;
+    savelink_data.nr_nodes = nr_nodes;
+    savelink_data.counts = sink_counts;
+    savelink_data.parts = (void *)sinks;
+    savelink_data.nodeID = nodeID;
+    threadpool_map(&e->threadpool, engine_redistribute_savelink_mapper_sink,
+                   nodes, nr_nodes, sizeof(int), threadpool_auto_chunk_size,
+                   &savelink_data);
+  }
+  swift_free("sink_dest", sink_dest);
+
   /* Get destination of each g-particle */
   int *g_counts;
-  if ((g_counts = (int *)calloc(sizeof(int), nr_nodes * nr_nodes)) == NULL)
+  if ((g_counts = (int *)calloc(nr_nodes * nr_nodes, sizeof(int))) == NULL)
     error("Failed to allocate g_gcount temporary buffer.");
 
   int *g_dest;
@@ -932,22 +1053,30 @@ void engine_redistribute(struct engine *e) {
                     MPI_SUM, MPI_COMM_WORLD) != MPI_SUCCESS)
     error("Failed to allreduce bparticle transfer counts.");
 
+  /* Get all the sink_counts from all the nodes. */
+  if (MPI_Allreduce(MPI_IN_PLACE, sink_counts, nr_nodes * nr_nodes, MPI_INT,
+                    MPI_SUM, MPI_COMM_WORLD) != MPI_SUCCESS)
+    error("Failed to allreduce sink particle transfer counts.");
+
   /* Report how many particles will be moved. */
   if (e->verbose) {
     if (e->nodeID == 0) {
-      size_t total = 0, g_total = 0, s_total = 0, b_total = 0;
-      size_t unmoved = 0, g_unmoved = 0, s_unmoved = 0, b_unmoved = 0;
+      size_t total = 0, g_total = 0, s_total = 0, b_total = 0, sink_total = 0;
+      size_t unmoved = 0, g_unmoved = 0, s_unmoved = 0, b_unmoved = 0,
+             sink_unmoved = 0;
       for (int p = 0, r = 0; p < nr_nodes; p++) {
         for (int n = 0; n < nr_nodes; n++) {
           total += counts[r];
           g_total += g_counts[r];
           s_total += s_counts[r];
           b_total += b_counts[r];
+          sink_total += sink_counts[r];
           if (p == n) {
             unmoved += counts[r];
             g_unmoved += g_counts[r];
             s_unmoved += s_counts[r];
             b_unmoved += b_counts[r];
+            sink_unmoved += sink_counts[r];
           }
           r++;
         }
@@ -967,14 +1096,19 @@ void engine_redistribute(struct engine *e) {
         message("%ld of %ld (%.2f%%) of b-particles moved", b_total - b_unmoved,
                 b_total,
                 100.0 * (double)(b_total - b_unmoved) / (double)b_total);
+      if (sink_total > 0)
+        message(
+            "%ld of %ld (%.2f%%) of sink-particles moved",
+            sink_total - sink_unmoved, sink_total,
+            100.0 * (double)(sink_total - sink_unmoved) / (double)sink_total);
     }
   }
 
-  /* Now each node knows how many parts, sparts, bparts, and gparts will be
-   * transferred to every other node. Get the new numbers of particles for this
-   * node. */
+  /* Now each node knows how many parts, sparts, bparts, sinks and gparts will
+   * be transferred to every other node. Get the new numbers of particles for
+   * this node. */
   size_t nr_parts_new = 0, nr_gparts_new = 0, nr_sparts_new = 0,
-         nr_bparts_new = 0;
+         nr_bparts_new = 0, nr_sinks_new = 0;
   for (int k = 0; k < nr_nodes; k++)
     nr_parts_new += counts[k * nr_nodes + nodeID];
   for (int k = 0; k < nr_nodes; k++)
@@ -983,6 +1117,8 @@ void engine_redistribute(struct engine *e) {
     nr_sparts_new += s_counts[k * nr_nodes + nodeID];
   for (int k = 0; k < nr_nodes; k++)
     nr_bparts_new += b_counts[k * nr_nodes + nodeID];
+  for (int k = 0; k < nr_nodes; k++)
+    nr_sinks_new += sink_counts[k * nr_nodes + nodeID];
 
 #ifdef WITH_CSDS
   const int initial_redistribute = e->ti_current == 0;
@@ -992,6 +1128,7 @@ void engine_redistribute(struct engine *e) {
     size_t spart_offset = 0;
     size_t gpart_offset = 0;
     size_t bpart_offset = 0;
+    size_t sink_offset = 0;
 
     for (int i = 0; i < nr_nodes; i++) {
       const size_t c_ind = engine_rank * nr_nodes + i;
@@ -1002,6 +1139,7 @@ void engine_redistribute(struct engine *e) {
         spart_offset += s_counts[c_ind];
         gpart_offset += g_counts[c_ind];
         bpart_offset += b_counts[c_ind];
+        sink_offset += sink_counts[c_ind];
         continue;
       }
 
@@ -1023,11 +1161,17 @@ void engine_redistribute(struct engine *e) {
         error("TODO");
       }
 
+      /* Log the sinks */
+      if (sink_counts[c_ind] > 0) {
+        error("TODO");
+      }
+
       /* Update the counters */
       part_offset += counts[c_ind];
       spart_offset += s_counts[c_ind];
       gpart_offset += g_counts[c_ind];
       bpart_offset += b_counts[c_ind];
+      sink_offset += sink_counts[c_ind];
     }
   }
 #endif
@@ -1081,6 +1225,15 @@ void engine_redistribute(struct engine *e) {
   s->nr_bparts = nr_bparts_new;
   s->size_bparts = engine_redistribute_alloc_margin * nr_bparts_new;
 
+  /* Sink particles. */
+  new_parts = engine_do_redistribute(
+      "sinks", sink_counts, (char *)s->sinks, nr_sinks_new, sizeof(struct sink),
+      sink_align, sink_mpi_type, nr_nodes, nodeID, e->syncredist);
+  swift_free("sinks", s->sinks);
+  s->sinks = (struct sink *)new_parts;
+  s->nr_sinks = nr_sinks_new;
+  s->size_sinks = engine_redistribute_alloc_margin * nr_sinks_new;
+
   /* All particles have now arrived. Time for some final operations on the
      stuff we just received */
 
@@ -1090,6 +1243,7 @@ void engine_redistribute(struct engine *e) {
     size_t spart_offset = 0;
     size_t gpart_offset = 0;
     size_t bpart_offset = 0;
+    size_t sink_offset = 0;
 
     for (int i = 0; i < nr_nodes; i++) {
       const size_t c_ind = i * nr_nodes + engine_rank;
@@ -1100,6 +1254,7 @@ void engine_redistribute(struct engine *e) {
         spart_offset += s_counts[c_ind];
         gpart_offset += g_counts[c_ind];
         bpart_offset += b_counts[c_ind];
+        sink_offset += sink_counts[c_ind];
         continue;
       }
 
@@ -1121,11 +1276,17 @@ void engine_redistribute(struct engine *e) {
         error("TODO");
       }
 
+      /* Log the sinks */
+      if (sink_counts[c_ind] > 0) {
+        error("TODO");
+      }
+
       /* Update the counters */
       part_offset += counts[c_ind];
       spart_offset += s_counts[c_ind];
       gpart_offset += g_counts[c_ind];
       bpart_offset += b_counts[c_ind];
+      sink_offset += sink_counts[c_ind];
     }
   }
 #endif
@@ -1139,6 +1300,7 @@ void engine_redistribute(struct engine *e) {
   relink_data.g_counts = g_counts;
   relink_data.s_counts = s_counts;
   relink_data.b_counts = b_counts;
+  relink_data.sink_counts = sink_counts;
   relink_data.nodeID = nodeID;
   relink_data.nr_nodes = nr_nodes;
 
@@ -1151,6 +1313,7 @@ void engine_redistribute(struct engine *e) {
   free(g_counts);
   free(s_counts);
   free(b_counts);
+  free(sink_counts);
 
 #ifdef SWIFT_DEBUG_CHECKS
   /* Verify that all parts are in the right place. */
@@ -1186,6 +1349,15 @@ void engine_redistribute(struct engine *e) {
       error("Received b-particle (%zu) that does not belong here (nodeID=%i).",
             k, cells[cid].nodeID);
   }
+  for (size_t k = 0; k < nr_sinks_new; k++) {
+    const int cid = cell_getid(s->cdim, s->sinks[k].x[0] * s->iwidth[0],
+                               s->sinks[k].x[1] * s->iwidth[1],
+                               s->sinks[k].x[2] * s->iwidth[2]);
+    if (cells[cid].nodeID != nodeID)
+      error(
+          "Received sink-particle (%zu) that does not belong here (nodeID=%i).",
+          k, cells[cid].nodeID);
+  }
 
   /* Verify that the links are correct */
   part_verify_links(s->parts, s->gparts, s->sinks, s->sparts, s->bparts,
@@ -1200,10 +1372,11 @@ void engine_redistribute(struct engine *e) {
     for (int k = 0; k < nr_cells; k++)
       if (cells[k].nodeID == nodeID) my_cells += 1;
     message(
-        "node %i now has %zu parts, %zu sparts, %zu bparts and %zu gparts in "
+        "node %i now has %zu parts, %zu sparts, %zu bparts, %zu sinks and %zu "
+        "gparts in "
         "%i cells.",
-        nodeID, nr_parts_new, nr_sparts_new, nr_bparts_new, nr_gparts_new,
-        my_cells);
+        nodeID, nr_parts_new, nr_sparts_new, nr_bparts_new, nr_sinks_new,
+        nr_gparts_new, my_cells);
   }
 
   /* Flag that we do not have any extra particles any more */
@@ -1211,6 +1384,7 @@ void engine_redistribute(struct engine *e) {
   s->nr_extra_gparts = 0;
   s->nr_extra_sparts = 0;
   s->nr_extra_bparts = 0;
+  s->nr_extra_sinks = 0;
 
   /* Flag that a redistribute has taken place */
   e->step_props |= engine_step_prop_redistribute;

src/engine_split_particles.c

@@ -33,6 +33,7 @@
 #include "random.h"
 #include "rt.h"
 #include "star_formation.h"
+#include "tools.h"
 #include "tracers.h"
 
 const int particle_split_factor = 2;
@@ -60,6 +61,7 @@ struct data_split {
   long long offset_id;
   long long *count_id;
   swift_lock_type lock;
+  FILE *extra_split_logger;
 };
 
 /**
@@ -172,19 +174,20 @@ void engine_split_gas_particle_split_mapper(void *restrict map_data, int count,
         memcpy(&global_gparts[k_gparts], gp, sizeof(struct gpart));
       }
 
-      /* Update splitting tree */
-      particle_splitting_update_binary_tree(&global_xparts[k_parts].split_data,
-                                            &xp->split_data);
-
       /* Update the IDs. */
       if (generate_random_ids) {
         /* The gas IDs are always odd, so we multiply by two here to
-         * repsect the parity. */
+         * respect the parity. */
         global_parts[k_parts].id += 2 * (long long)rand_r(&seedp);
       } else {
         global_parts[k_parts].id = offset_id + 2 * atomic_inc(count_id);
       }
 
+      /* Update splitting tree */
+      particle_splitting_update_binary_tree(
+          &xp->split_data, &global_xparts[k_parts].split_data, p->id,
+          global_parts[k_parts].id, data->extra_split_logger, &data->lock);
+
       /* Re-link everything */
       if (with_gravity) {
         global_parts[k_parts].gpart = &global_gparts[k_gparts];
@@ -312,7 +315,10 @@ void engine_split_gas_particles(struct engine *e) {
 
   /* Verify that nothing wrong happened with the IDs */
   if (data_count.max_id > e->max_parts_id) {
-    error("Found a gas particle with an ID larger than the current max!");
+    error(
+        "Found a gas particle with an ID (%lld) larger than the current max "
+        "(%lld)!",
+        data_count.max_id, e->max_parts_id);
   }
 
   /* Be verbose about this. This is an important event */
@@ -432,12 +438,22 @@ void engine_split_gas_particles(struct engine *e) {
   size_t k_parts = s->nr_parts;
   size_t k_gparts = s->nr_gparts;
 
+  FILE *extra_split_logger = NULL;
+  if (e->hydro_properties->log_extra_splits_in_file) {
+    char extra_split_logger_filename[256];
+    sprintf(extra_split_logger_filename, "splits/splits_%04d.txt", engine_rank);
+    extra_split_logger = fopen(extra_split_logger_filename, "a");
+    if (extra_split_logger == NULL) error("Error opening split logger file!");
+  }
+
   /* Loop over the particles again to split them */
   long long local_count_id = 0;
   struct data_split data_split = {
-      e,         mass_threshold, generate_random_ids, &k_parts,
-      &k_gparts, offset_id,      &local_count_id,     0};
+      e,          mass_threshold,    generate_random_ids, &k_parts,
+      &k_gparts,  offset_id,         &local_count_id,
+      /*lock=*/0, extra_split_logger};
   lock_init(&data_split.lock);
+
   threadpool_map(&e->threadpool, engine_split_gas_particle_split_mapper,
                  s->parts, nr_parts_old, sizeof(struct part), 0, &data_split);
   if (lock_destroy(&data_split.lock) != 0) error("Error destroying lock");
@@ -459,7 +475,33 @@ void engine_split_gas_particles(struct engine *e) {
   }
 #endif
 
+  /* Close the logger file */
+  if (e->hydro_properties->log_extra_splits_in_file) fclose(extra_split_logger);
+
   if (e->verbose)
     message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
             clocks_getunit());
 }
+
+void engine_init_split_gas_particles(struct engine *e) {
+
+  if (e->hydro_properties->log_extra_splits_in_file) {
+
+    /* Create the directory to host the logs */
+    if (engine_rank == 0) safe_checkdir("splits", /*create=*/1);
+
+#ifdef WITH_MPI
+    MPI_Barrier(MPI_COMM_WORLD);
+#endif
+
+    /* Create the logger files and add a header */
+    char extra_split_logger_filename[256];
+    sprintf(extra_split_logger_filename, "splits/splits_%04d.txt", engine_rank);
+    FILE *extra_split_logger = fopen(extra_split_logger_filename, "w");
+    fprintf(extra_split_logger, "# %12s %20s %20s %20s %20s\n", "Step", "ID",
+            "Progenitor", "Count", "Tree");
+
+    /* Close everything for now */
+    fclose(extra_split_logger);
+  }
+}