diff --git a/out.txt b/out.txt
deleted file mode 100644
index 9b33fac26edf8c8436b76ba570e5e4f89cc72161..0000000000000000000000000000000000000000
--- a/out.txt
+++ /dev/null
@@ -1,4500 +0,0 @@
-diff --git a/examples/main.c b/examples/main.c
-index 7babf14..2d3ae6a 100644
---- a/examples/main.c
-+++ b/examples/main.c
-@@ -153,7 +153,6 @@ int main(int argc, char *argv[]) {
-   int with_stars = 0;
-   int with_star_formation = 0;
-   int with_feedback = 0;
--  int with_limiter = 0;
-   int with_fp_exceptions = 0;
-   int with_drift_all = 0;
-   int with_mpole_reconstruction = 0;
-@@ -203,8 +202,6 @@ int main(int argc, char *argv[]) {
-       OPT_BOOLEAN('S', "stars", &with_stars, "Run with stars.", NULL, 0, 0),
-       OPT_BOOLEAN('x', "velociraptor", &with_structure_finding,
-                   "Run with structure finding.", NULL, 0, 0),
--      OPT_BOOLEAN(0, "limiter", &with_limiter, "Run with time-step limiter.",
--                  NULL, 0, 0),
- 
-       OPT_GROUP("  Control options:\n"),
-       OPT_BOOLEAN('a', "pin", &with_aff,
-@@ -456,10 +453,8 @@ int main(int argc, char *argv[]) {
- #ifdef WITH_MPI
-   if (with_mpole_reconstruction && nr_nodes > 1)
-     error("Cannot reconstruct m-poles every step over MPI (yet).");
--  if (with_feedback) error("Can't run with feedback over MPI (yet).");
-   if (with_star_formation)
-     error("Can't run with star formation over MPI (yet)");
--  if (with_limiter) error("Can't run with time-step limiter over MPI (yet)");
- #endif
- 
- #if defined(WITH_MPI) && defined(HAVE_VELOCIRAPTOR)
-@@ -696,7 +691,7 @@ int main(int argc, char *argv[]) {
-     /* Initialise the stars properties */
-     if (with_stars)
-       stars_props_init(&stars_properties, &prog_const, &us, params,
--                       &hydro_properties, &cosmo);
-+                       &hydro_properties);
-     else
-       bzero(&stars_properties, sizeof(struct stars_props));
- 
-@@ -846,18 +841,6 @@ int main(int argc, char *argv[]) {
-       fflush(stdout);
-     }
- 
--    /* Verify that we are not using basic modes incorrectly */
--    if (with_hydro && N_total[0] == 0) {
--      error(
--          "ERROR: Running with hydrodynamics but no gas particles found in the "
--          "ICs!");
--    }
--    if ((with_self_gravity || with_external_gravity) && N_total[1] == 0) {
--      error(
--          "ERROR: Running with gravity but no gravity particles found in "
--          "the ICs!");
--    }
--
-     /* Verify that each particle is in it's proper cell. */
-     if (talking && !dry_run) {
-       int icount = 0;
-@@ -900,7 +883,6 @@ int main(int argc, char *argv[]) {
-       engine_policies |= engine_policy_external_gravity;
-     if (with_cosmology) engine_policies |= engine_policy_cosmology;
-     if (with_temperature) engine_policies |= engine_policy_temperature;
--    if (with_limiter) engine_policies |= engine_policy_limiter;
-     if (with_cooling) engine_policies |= engine_policy_cooling;
-     if (with_stars) engine_policies |= engine_policy_stars;
-     if (with_star_formation) engine_policies |= engine_policy_star_formation;
-diff --git a/src/cell.c b/src/cell.c
-index de63e58..4165325 100644
---- a/src/cell.c
-+++ b/src/cell.c
-@@ -181,6 +181,7 @@ int cell_pack(struct cell *restrict c, struct pcell *restrict pc,
- 
-   /* Start by packing the data of the current cell. */
-   pc->hydro.h_max = c->hydro.h_max;
-+  pc->stars.h_max = c->stars.h_max;
-   pc->hydro.ti_end_min = c->hydro.ti_end_min;
-   pc->hydro.ti_end_max = c->hydro.ti_end_max;
-   pc->grav.ti_end_min = c->grav.ti_end_min;
-@@ -285,6 +286,7 @@ int cell_unpack(struct pcell *restrict pc, struct cell *restrict c,
- 
-   /* Unpack the current pcell. */
-   c->hydro.h_max = pc->hydro.h_max;
-+  c->stars.h_max = pc->stars.h_max;
-   c->hydro.ti_end_min = pc->hydro.ti_end_min;
-   c->hydro.ti_end_max = pc->hydro.ti_end_max;
-   c->grav.ti_end_min = pc->grav.ti_end_min;
-@@ -1232,11 +1234,10 @@ void cell_clean_links(struct cell *c, void *data) {
-   c->hydro.density = NULL;
-   c->hydro.gradient = NULL;
-   c->hydro.force = NULL;
--  c->hydro.limiter = NULL;
-   c->grav.grav = NULL;
-   c->grav.mm = NULL;
--  c->stars.density = NULL;
-   c->stars.feedback = NULL;
-+  c->stars.density = NULL;
- }
- 
- /**
-@@ -1603,14 +1604,6 @@ void cell_clear_drift_flags(struct cell *c, void *data) {
- }
- 
- /**
-- * @brief Clear the limiter flags on the given cell.
-- */
--void cell_clear_limiter_flags(struct cell *c, void *data) {
--  c->hydro.do_limiter = 0;
--  c->hydro.do_sub_limiter = 0;
--}
--
--/**
-  * @brief Activate the #part drifts on the given cell.
-  */
- void cell_activate_drift_part(struct cell *c, struct scheduler *s) {
-@@ -1633,10 +1626,7 @@ void cell_activate_drift_part(struct cell *c, struct scheduler *s) {
-     for (struct cell *parent = c->parent;
-          parent != NULL && !parent->hydro.do_sub_drift;
-          parent = parent->parent) {
--
--      /* Mark this cell for drifting */
-       parent->hydro.do_sub_drift = 1;
--
-       if (parent == c->hydro.super) {
- #ifdef SWIFT_DEBUG_CHECKS
-         if (parent->hydro.drift == NULL)
-@@ -1701,45 +1691,6 @@ void cell_activate_drift_spart(struct cell *c, struct scheduler *s) {
- }
- 
- /**
-- * @brief Activate the drifts on the given cell.
-- */
--void cell_activate_limiter(struct cell *c, struct scheduler *s) {
--
--  /* If this cell is already marked for drift, quit early. */
--  if (c->hydro.do_limiter) return;
--
--  /* Mark this cell for limiting. */
--  c->hydro.do_limiter = 1;
--
--  /* Set the do_sub_limiter all the way up and activate the super limiter
--     if this has not yet been done. */
--  if (c == c->super) {
--#ifdef SWIFT_DEBUG_CHECKS
--    if (c->timestep_limiter == NULL)
--      error("Trying to activate un-existing c->timestep_limiter");
--#endif
--    scheduler_activate(s, c->timestep_limiter);
--  } else {
--    for (struct cell *parent = c->parent;
--         parent != NULL && !parent->hydro.do_sub_limiter;
--         parent = parent->parent) {
--
--      /* Mark this cell for limiting */
--      parent->hydro.do_sub_limiter = 1;
--
--      if (parent == c->super) {
--#ifdef SWIFT_DEBUG_CHECKS
--        if (parent->timestep_limiter == NULL)
--          error("Trying to activate un-existing parent->timestep_limiter");
--#endif
--        scheduler_activate(s, parent->timestep_limiter);
--        break;
--      }
--    }
--  }
--}
--
--/**
-  * @brief Activate the sorts up a cell hierarchy.
-  */
- void cell_activate_hydro_sorts_up(struct cell *c, struct scheduler *s) {
-@@ -1802,15 +1753,19 @@ void cell_activate_stars_sorts_up(struct cell *c, struct scheduler *s) {
- 
-   if (c == c->super) {
- #ifdef SWIFT_DEBUG_CHECKS
--    if (c->stars.sorts == NULL)
-+    if ((c->nodeID == engine_rank && c->stars.sorts_local == NULL) ||
-+        (c->nodeID != engine_rank && c->stars.sorts_foreign == NULL))
-       error("Trying to activate un-existing c->stars.sorts");
- #endif
--    scheduler_activate(s, c->stars.sorts);
-     if (c->nodeID == engine_rank) {
-+      scheduler_activate(s, c->stars.sorts_local);
-       // MATTHIEU: to do: do we actually need both drifts here?
-       cell_activate_drift_part(c, s);
-       cell_activate_drift_spart(c, s);
-     }
-+    if (c->nodeID != engine_rank) {
-+      scheduler_activate(s, c->stars.sorts_foreign);
-+    }
-   } else {
- 
-     for (struct cell *parent = c->parent;
-@@ -1819,14 +1774,18 @@ void cell_activate_stars_sorts_up(struct cell *c, struct scheduler *s) {
-       parent->stars.do_sub_sort = 1;
-       if (parent == c->super) {
- #ifdef SWIFT_DEBUG_CHECKS
--        if (parent->stars.sorts == NULL)
-+        if ((c->nodeID == engine_rank && parent->stars.sorts_local == NULL) ||
-+            (c->nodeID != engine_rank && parent->stars.sorts_foreign == NULL))
-           error("Trying to activate un-existing parents->stars.sorts");
- #endif
--        scheduler_activate(s, parent->stars.sorts);
--        if (parent->nodeID == engine_rank) {
-+        if (c->nodeID == engine_rank) {
-+          scheduler_activate(s, parent->stars.sorts_local);
-           cell_activate_drift_part(parent, s);
-           cell_activate_drift_spart(parent, s);
-         }
-+        if (c->nodeID != engine_rank) {
-+          scheduler_activate(s, parent->stars.sorts_foreign);
-+        }
-         break;
-       }
-     }
-@@ -1869,7 +1828,6 @@ void cell_activate_stars_sorts(struct cell *c, int sid, struct scheduler *s) {
- void cell_activate_subcell_hydro_tasks(struct cell *ci, struct cell *cj,
-                                        struct scheduler *s) {
-   const struct engine *e = s->space->e;
--  const int with_limiter = (e->policy & engine_policy_limiter);
- 
-   /* Store the current dx_max and h_max values. */
-   ci->hydro.dx_max_part_old = ci->hydro.dx_max_part;
-@@ -1903,7 +1861,6 @@ void cell_activate_subcell_hydro_tasks(struct cell *ci, struct cell *cj,
- 
-       /* We have reached the bottom of the tree: activate drift */
-       cell_activate_drift_part(ci, s);
--      if (with_limiter) cell_activate_limiter(ci, s);
-     }
-   }
- 
-@@ -2209,12 +2166,6 @@ void cell_activate_subcell_hydro_tasks(struct cell *ci, struct cell *cj,
-       if (ci->nodeID == engine_rank) cell_activate_drift_part(ci, s);
-       if (cj->nodeID == engine_rank) cell_activate_drift_part(cj, s);
- 
--      /* Also activate the time-step limiter */
--      if (ci->nodeID == engine_rank && with_limiter)
--        cell_activate_limiter(ci, s);
--      if (cj->nodeID == engine_rank && with_limiter)
--        cell_activate_limiter(cj, s);
--
-       /* Do we need to sort the cells? */
-       cell_activate_hydro_sorts(ci, sid, s);
-       cell_activate_hydro_sorts(cj, sid, s);
-@@ -2779,7 +2730,6 @@ int cell_unskip_hydro_tasks(struct cell *c, struct scheduler *s) {
- 
-   struct engine *e = s->space->e;
-   const int nodeID = e->nodeID;
--  const int with_limiter = (e->policy & engine_policy_limiter);
-   int rebuild = 0;
- 
-   /* Un-skip the density tasks involved with this cell. */
-@@ -2805,7 +2755,6 @@ int cell_unskip_hydro_tasks(struct cell *c, struct scheduler *s) {
-       /* Activate hydro drift */
-       if (t->type == task_type_self) {
-         if (ci_nodeID == nodeID) cell_activate_drift_part(ci, s);
--        if (ci_nodeID == nodeID && with_limiter) cell_activate_limiter(ci, s);
-       }
- 
-       /* Set the correct sorting flags and activate hydro drifts */
-@@ -2820,10 +2769,6 @@ int cell_unskip_hydro_tasks(struct cell *c, struct scheduler *s) {
-         if (ci_nodeID == nodeID) cell_activate_drift_part(ci, s);
-         if (cj_nodeID == nodeID) cell_activate_drift_part(cj, s);
- 
--        /* Activate the limiter tasks. */
--        if (ci_nodeID == nodeID && with_limiter) cell_activate_limiter(ci, s);
--        if (cj_nodeID == nodeID && with_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);
-@@ -2858,11 +2803,7 @@ int cell_unskip_hydro_tasks(struct cell *c, struct scheduler *s) {
-         }
- 
-         /* If the foreign cell is active, we want its ti_end values. */
--        if (ci_active || with_limiter) scheduler_activate(s, ci->mpi.recv_ti);
--
--        if (with_limiter) scheduler_activate(s, ci->mpi.limiter.recv);
--        if (with_limiter)
--          scheduler_activate_send(s, cj->mpi.limiter.send, ci->nodeID);
-+        if (ci_active) scheduler_activate(s, ci->mpi.recv_ti);
- 
-         /* Is the foreign cell active and will need stuff from us? */
-         if (ci_active) {
-@@ -2872,7 +2813,6 @@ int cell_unskip_hydro_tasks(struct cell *c, struct scheduler *s) {
-           /* 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);
--          if (with_limiter) cell_activate_limiter(cj, s);
- 
-           /* If the local cell is also active, more stuff will be needed. */
-           if (cj_active) {
-@@ -2885,8 +2825,7 @@ int cell_unskip_hydro_tasks(struct cell *c, struct scheduler *s) {
-         }
- 
-         /* If the local cell is active, send its ti_end values. */
--        if (cj_active || with_limiter)
--          scheduler_activate_send(s, cj->mpi.send_ti, ci_nodeID);
-+        if (cj_active) scheduler_activate_send(s, cj->mpi.send_ti, ci_nodeID);
- 
-       } else if (cj_nodeID != nodeID) {
- 
-@@ -2903,11 +2842,7 @@ int cell_unskip_hydro_tasks(struct cell *c, struct scheduler *s) {
-         }
- 
-         /* If the foreign cell is active, we want its ti_end values. */
--        if (cj_active || with_limiter) scheduler_activate(s, cj->mpi.recv_ti);
--
--        if (with_limiter) scheduler_activate(s, cj->mpi.limiter.recv);
--        if (with_limiter)
--          scheduler_activate_send(s, ci->mpi.limiter.send, cj->nodeID);
-+        if (cj_active) scheduler_activate(s, cj->mpi.recv_ti);
- 
-         /* Is the foreign cell active and will need stuff from us? */
-         if (cj_active) {
-@@ -2917,7 +2852,6 @@ int cell_unskip_hydro_tasks(struct cell *c, struct scheduler *s) {
-           /* 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);
--          if (with_limiter) cell_activate_limiter(ci, s);
- 
-           /* If the local cell is also active, more stuff will be needed. */
-           if (ci_active) {
-@@ -2931,8 +2865,7 @@ int cell_unskip_hydro_tasks(struct cell *c, struct scheduler *s) {
-         }
- 
-         /* If the local cell is active, send its ti_end values. */
--        if (ci_active || with_limiter)
--          scheduler_activate_send(s, ci->mpi.send_ti, cj_nodeID);
-+        if (ci_active) scheduler_activate_send(s, ci->mpi.send_ti, cj_nodeID);
-       }
- #endif
-     }
-@@ -2945,8 +2878,6 @@ int cell_unskip_hydro_tasks(struct cell *c, struct scheduler *s) {
-       scheduler_activate(s, l->t);
-     for (struct link *l = c->hydro.force; l != NULL; l = l->next)
-       scheduler_activate(s, l->t);
--    for (struct link *l = c->hydro.limiter; l != NULL; l = l->next)
--      scheduler_activate(s, l->t);
- 
-     if (c->hydro.extra_ghost != NULL)
-       scheduler_activate(s, c->hydro.extra_ghost);
-@@ -2960,6 +2891,7 @@ int cell_unskip_hydro_tasks(struct cell *c, struct scheduler *s) {
-     if (c->hydro.cooling != NULL) scheduler_activate(s, c->hydro.cooling);
-     if (c->hydro.star_formation != NULL)
-       scheduler_activate(s, c->hydro.star_formation);
-+    if (c->sourceterms != NULL) scheduler_activate(s, c->sourceterms);
-     if (c->logger != NULL) scheduler_activate(s, c->logger);
-   }
- 
-@@ -3144,6 +3076,10 @@ int cell_unskip_stars_tasks(struct cell *c, struct scheduler *s) {
-     struct cell *cj = t->cj;
-     const int ci_active = cell_is_active_stars(ci, e);
-     const int cj_active = (cj != NULL) ? cell_is_active_stars(cj, e) : 0;
-+#ifdef WITH_MPI
-+    const int ci_nodeID = ci->nodeID;
-+    const int cj_nodeID = (cj != NULL) ? cj->nodeID : -1;
-+#endif
- 
-     /* Only activate tasks that involve a local active cell. */
-     if ((ci_active && ci->nodeID == nodeID) ||
-@@ -3161,38 +3097,42 @@ int cell_unskip_stars_tasks(struct cell *c, struct scheduler *s) {
-       /* Set the correct sorting flags and activate hydro drifts */
-       else if (t->type == task_type_pair) {
-         /* Do ci */
--        /* 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);
--
--        /* Check the sorts and activate them if needed. */
--        cell_activate_stars_sorts(ci, t->flags, s);
--        cell_activate_hydro_sorts(cj, t->flags, s);
-+        if (ci_active && cj->hydro.count != 0 && ci->stars.count != 0) {
-+          /* 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);
-+
-+          /* Check the sorts and activate them if needed. */
-+          cell_activate_stars_sorts(ci, t->flags, s);
-+          cell_activate_hydro_sorts(cj, t->flags, s);
-+        }
- 
-         /* Do cj */
--        /* hydro for ci */
--        atomic_or(&ci->hydro.requires_sorts, 1 << t->flags);
--        ci->hydro.dx_max_sort_old = ci->hydro.dx_max_sort;
--
--        /* stars for cj */
--        atomic_or(&cj->stars.requires_sorts, 1 << t->flags);
--        cj->stars.dx_max_sort_old = cj->stars.dx_max_sort;
--
--        /* Activate the drift tasks. */
--        if (cj->nodeID == nodeID) cell_activate_drift_spart(cj, s);
--        if (ci->nodeID == nodeID) cell_activate_drift_part(ci, s);
--
--        /* Check the sorts and activate them if needed. */
--        cell_activate_hydro_sorts(ci, t->flags, s);
--        cell_activate_stars_sorts(cj, t->flags, s);
-+        if (cj_active && ci->hydro.count != 0 && cj->stars.count != 0) {
-+          /* hydro for ci */
-+          atomic_or(&ci->hydro.requires_sorts, 1 << t->flags);
-+          ci->hydro.dx_max_sort_old = ci->hydro.dx_max_sort;
-+
-+          /* stars for cj */
-+          atomic_or(&cj->stars.requires_sorts, 1 << t->flags);
-+          cj->stars.dx_max_sort_old = cj->stars.dx_max_sort;
-+
-+          /* Activate the drift tasks. */
-+          if (cj->nodeID == nodeID) cell_activate_drift_spart(cj, s);
-+          if (ci->nodeID == nodeID) cell_activate_drift_part(ci, s);
-+
-+          /* 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. */
-@@ -3207,110 +3147,77 @@ int cell_unskip_stars_tasks(struct cell *c, struct scheduler *s) {
-       /* Check whether there was too much particle motion, i.e. the
-          cell neighbour conditions were violated. */
-       if (cell_need_rebuild_for_stars_pair(ci, cj)) rebuild = 1;
-+      if (cell_need_rebuild_for_hydro_pair(ci, cj)) rebuild = 1;
- 
- #ifdef WITH_MPI
--      error("MPI with stars not implemented");
--      /* /\* Activate the send/recv tasks. *\/ */
--      /* if (ci->nodeID != nodeID) { */
--
--      /*   /\* If the local cell is active, receive data from the foreign cell.
--       * *\/ */
--      /*   if (cj_active) { */
--      /*     scheduler_activate(s, ci->hydro.recv_xv); */
--      /*     if (ci_active) { */
--      /*       scheduler_activate(s, ci->hydro.recv_rho); */
--
--      /*     } */
--      /*   } */
--
--      /*   /\* If the foreign cell is active, we want its ti_end values. *\/ */
--      /*   if (ci_active) scheduler_activate(s, ci->mpi.recv_ti); */
--
--      /*   /\* Is the foreign cell active and will need stuff from us? *\/ */
--      /*   if (ci_active) { */
--
--      /*     scheduler_activate_send(s, cj->hydro.send_xv, ci->nodeID); */
-+      /* Activate the send/recv tasks. */
-+      if (ci_nodeID != 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); */
-+        if (cj_active) {
-+          scheduler_activate(s, ci->mpi.hydro.recv_xv);
-+          if (ci_active) {
-+            scheduler_activate(s, ci->mpi.stars.recv);
-+          }
-+        }
- 
--      /*     /\* If the local cell is also active, more stuff will be needed.
--       * *\/ */
--      /*     if (cj_active) { */
--      /*       scheduler_activate_send(s, cj->hydro.send_rho, ci->nodeID); */
-+        /* If the foreign cell is active, we want its ti_end values. */
-+        if (ci_active) scheduler_activate(s, ci->mpi.recv_ti);
- 
--      /*     } */
--      /*   } */
-+        /* Is the foreign cell active and will need stuff from us? */
-+        if (ci_active) {
- 
--      /*   /\* If the local cell is active, send its ti_end values. *\/ */
--      /*   if (cj_active) scheduler_activate_send(s, cj->mpi.send_ti,
--       * ci->nodeID);
--       */
-+          scheduler_activate_send(s, cj->mpi.hydro.send_xv, ci_nodeID);
- 
--      /* } else if (cj->nodeID != 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);
- 
--      /*   /\* If the local cell is active, receive data from the foreign cell.
--       * *\/ */
--      /*   if (ci_active) { */
--      /*     scheduler_activate(s, cj->hydro.recv_xv); */
--      /*     if (cj_active) { */
--      /*       scheduler_activate(s, cj->hydro.recv_rho); */
-+          /* If the local cell is also active, more stuff will be needed.
-+           */
-+          if (cj_active) {
-+            scheduler_activate_send(s, cj->mpi.stars.send, ci_nodeID);
-+          }
-+        }
- 
--      /*     } */
--      /*   } */
-+        /* If the local cell is active, send its ti_end values. */
-+        if (cj_active) scheduler_activate_send(s, cj->mpi.send_ti, ci_nodeID);
- 
--      /*   /\* If the foreign cell is active, we want its ti_end values. *\/ */
--      /*   if (cj_active) scheduler_activate(s, cj->mpi.recv_ti); */
-+      } else if (cj_nodeID != nodeID) {
- 
--      /*   /\* Is the foreign cell active and will need stuff from us? *\/ */
--      /*   if (cj_active) { */
-+        /* If the local cell is active, receive data from the foreign cell. */
-+        if (ci_active) {
-+          scheduler_activate(s, cj->mpi.hydro.recv_xv);
-+          if (cj_active) {
-+            scheduler_activate(s, cj->mpi.stars.recv);
-+          }
-+        }
- 
--      /*     scheduler_activate_send(s, ci->hydro.send_xv, cj->nodeID); */
-+        /* If the foreign cell is active, we want its ti_end values. */
-+        if (cj_active) scheduler_activate(s, cj->mpi.recv_ti);
- 
--      /*     /\* 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); */
-+        /* Is the foreign cell active and will need stuff from us? */
-+        if (cj_active) {
- 
--      /*     /\* If the local cell is also active, more stuff will be needed.
--       * *\/ */
--      /*     if (ci_active) { */
-+          scheduler_activate_send(s, ci->mpi.hydro.send_xv, cj_nodeID);
- 
--      /*       scheduler_activate_send(s, ci->hydro.send_rho, 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);
- 
--      /*     } */
--      /*   } */
-+          /* If the local cell is also active, more stuff will be needed.
-+           */
-+          if (ci_active) {
-+            scheduler_activate_send(s, ci->mpi.stars.send, cj_nodeID);
-+          }
-+        }
- 
--      /*   /\* If the local cell is active, send its ti_end values. *\/ */
--      /*   if (ci_active) scheduler_activate_send(s, ci->mpi.send_ti,
--       * cj->nodeID);
--       */
--      /* } */
-+        /* If the local cell is active, send its ti_end values. */
-+        if (ci_active) scheduler_activate_send(s, ci->mpi.send_ti, cj_nodeID);
-+      }
- #endif
-     }
-   }
- 
--  /* Un-skip the feedback tasks involved with this cell. */
--  for (struct link *l = c->stars.feedback; l != NULL; l = l->next) {
--    struct task *t = l->t;
--    struct cell *ci = t->ci;
--    struct cell *cj = t->cj;
--    const int ci_active = cell_is_active_stars(ci, e);
--    const int cj_active = (cj != NULL) ? cell_is_active_stars(cj, e) : 0;
--
--    /* Only activate tasks that involve a local active cell. */
--    if ((ci_active && ci->nodeID == nodeID) ||
--        (cj_active && cj->nodeID == nodeID)) {
--      scheduler_activate(s, t);
--    }
--    if (t->type == task_type_pair || t->type == task_type_sub_pair) {
--
--      /* Check whether there was too much particle motion, i.e. the
--         cell neighbour conditions were violated. */
--      if (cell_need_rebuild_for_stars_pair(ci, cj)) rebuild = 1;
--    }
--  }
--
-   /* Unskip all the other task types. */
-   if (c->nodeID == nodeID && cell_is_active_stars(c, e)) {
- 
-diff --git a/src/common_io.c b/src/common_io.c
-index e6c0d01..24e7401 100644
---- a/src/common_io.c
-+++ b/src/common_io.c
-@@ -142,7 +142,7 @@ void io_read_attribute(hid_t grp, const char* name, enum IO_DATA_TYPE type,
-  * Calls #error() if an error occurs.
-  */
- void io_write_attribute(hid_t grp, const char* name, enum IO_DATA_TYPE type,
--                        const void* data, int num) {
-+                        void* data, int num) {
- 
-   const hid_t h_space = H5Screate(H5S_SIMPLE);
-   if (h_space < 0)
-@@ -387,332 +387,6 @@ void io_write_engine_policy(hid_t h_file, const struct engine* e) {
-   H5Gclose(h_grp);
- }
- 
--void io_write_cell_offsets(hid_t h_grp, const int cdim[3],
--                           const struct cell* cells_top, const int nr_cells,
--                           const double width[3], const int nodeID,
--                           const long long global_counts[swift_type_count],
--                           const long long global_offsets[swift_type_count],
--                           const struct unit_system* internal_units,
--                           const struct unit_system* snapshot_units) {
--
--  double cell_width[3] = {width[0], width[1], width[2]};
--
--  /* Temporary memory for the cell-by-cell information */
--  double* centres = NULL;
--  centres = (double*)malloc(3 * nr_cells * sizeof(double));
--
--  /* Count of particles in each cell */
--  long long *count_part = NULL, *count_gpart = NULL, *count_spart = NULL;
--  count_part = (long long*)malloc(nr_cells * sizeof(long long));
--  count_gpart = (long long*)malloc(nr_cells * sizeof(long long));
--  count_spart = (long long*)malloc(nr_cells * sizeof(long long));
--
--  /* Global offsets of particles in each cell */
--  long long *offset_part = NULL, *offset_gpart = NULL, *offset_spart = NULL;
--  offset_part = (long long*)malloc(nr_cells * sizeof(long long));
--  offset_gpart = (long long*)malloc(nr_cells * sizeof(long long));
--  offset_spart = (long long*)malloc(nr_cells * sizeof(long long));
--
--  /* Offsets of the 0^th element */
--  offset_part[0] = 0;
--  offset_gpart[0] = 0;
--  offset_spart[0] = 0;
--
--  /* Collect the cell information of *local* cells */
--  long long local_offset_part = 0;
--  long long local_offset_gpart = 0;
--  long long local_offset_spart = 0;
--  for (int i = 0; i < nr_cells; ++i) {
--
--    if (cells_top[i].nodeID == nodeID) {
--
--      /* Centre of each cell */
--      centres[i * 3 + 0] = cells_top[i].loc[0] + cell_width[0] * 0.5;
--      centres[i * 3 + 1] = cells_top[i].loc[1] + cell_width[1] * 0.5;
--      centres[i * 3 + 2] = cells_top[i].loc[2] + cell_width[2] * 0.5;
--
--      /* Count real particles that will be written */
--      count_part[i] = cells_top[i].hydro.count - cells_top[i].hydro.inhibited;
--      count_gpart[i] = cells_top[i].grav.count - cells_top[i].grav.inhibited;
--      count_spart[i] = cells_top[i].stars.count - cells_top[i].stars.inhibited;
--
--      /* Only count DM gpart (gpart without friends) */
--      count_gpart[i] -= count_part[i];
--      count_gpart[i] -= count_spart[i];
--
--      /* Offsets including the global offset of all particles on this MPI rank
--       */
--      offset_part[i] = local_offset_part + global_offsets[swift_type_gas];
--      offset_gpart[i] =
--          local_offset_gpart + global_offsets[swift_type_dark_matter];
--      offset_spart[i] = local_offset_spart + global_offsets[swift_type_stars];
--
--      local_offset_part += count_part[i];
--      local_offset_gpart += count_gpart[i];
--      local_offset_spart += count_spart[i];
--
--    } else {
--
--      /* Just zero everything for the foregin cells */
--
--      centres[i * 3 + 0] = 0.;
--      centres[i * 3 + 1] = 0.;
--      centres[i * 3 + 2] = 0.;
--
--      count_part[i] = 0;
--      count_gpart[i] = 0;
--      count_spart[i] = 0;
--
--      offset_part[i] = 0;
--      offset_gpart[i] = 0;
--      offset_spart[i] = 0;
--    }
--  }
--
--#ifdef WITH_MPI
--  /* Now, reduce all the arrays. Note that we use a bit-wise OR here. This
--     is safe as we made sure only local cells have non-zero values. */
--  if (nodeID == 0) {
--    MPI_Reduce(MPI_IN_PLACE, count_part, nr_cells, MPI_LONG_LONG_INT, MPI_BOR,
--               0, MPI_COMM_WORLD);
--  } else {
--    MPI_Reduce(count_part, NULL, nr_cells, MPI_LONG_LONG_INT, MPI_BOR, 0,
--               MPI_COMM_WORLD);
--  }
--  if (nodeID == 0) {
--    MPI_Reduce(MPI_IN_PLACE, count_gpart, nr_cells, MPI_LONG_LONG_INT, MPI_BOR,
--               0, MPI_COMM_WORLD);
--  } else {
--    MPI_Reduce(count_gpart, NULL, nr_cells, MPI_LONG_LONG_INT, MPI_BOR, 0,
--               MPI_COMM_WORLD);
--  }
--  if (nodeID == 0) {
--    MPI_Reduce(MPI_IN_PLACE, count_spart, nr_cells, MPI_LONG_LONG_INT, MPI_BOR,
--               0, MPI_COMM_WORLD);
--  } else {
--    MPI_Reduce(count_spart, NULL, nr_cells, MPI_LONG_LONG_INT, MPI_BOR, 0,
--               MPI_COMM_WORLD);
--  }
--  if (nodeID == 0) {
--    MPI_Reduce(MPI_IN_PLACE, offset_part, nr_cells, MPI_LONG_LONG_INT, MPI_BOR,
--               0, MPI_COMM_WORLD);
--  } else {
--    MPI_Reduce(offset_part, NULL, nr_cells, MPI_LONG_LONG_INT, MPI_BOR, 0,
--               MPI_COMM_WORLD);
--  }
--  if (nodeID == 0) {
--    MPI_Reduce(MPI_IN_PLACE, offset_gpart, nr_cells, MPI_LONG_LONG_INT, MPI_BOR,
--               0, MPI_COMM_WORLD);
--  } else {
--    MPI_Reduce(offset_gpart, NULL, nr_cells, MPI_LONG_LONG_INT, MPI_BOR, 0,
--               MPI_COMM_WORLD);
--  }
--  if (nodeID == 0) {
--    MPI_Reduce(MPI_IN_PLACE, offset_spart, nr_cells, MPI_LONG_LONG_INT, MPI_BOR,
--               0, MPI_COMM_WORLD);
--  } else {
--    MPI_Reduce(offset_spart, NULL, nr_cells, MPI_LONG_LONG_INT, MPI_BOR, 0,
--               MPI_COMM_WORLD);
--  }
--
--  /* For the centres we use a sum as MPI does not like bit-wise operations
--     on floating point numbers */
--  if (nodeID == 0) {
--    MPI_Reduce(MPI_IN_PLACE, centres, 3 * nr_cells, MPI_DOUBLE, MPI_SUM, 0,
--               MPI_COMM_WORLD);
--  } else {
--    MPI_Reduce(centres, NULL, 3 * nr_cells, MPI_DOUBLE, MPI_SUM, 0,
--               MPI_COMM_WORLD);
--  }
--#endif
--
--  /* Only rank 0 actually writes */
--  if (nodeID == 0) {
--
--    /* Unit conversion if necessary */
--    const double factor = units_conversion_factor(
--        internal_units, snapshot_units, UNIT_CONV_LENGTH);
--    if (factor != 1.) {
--
--      /* Convert the cell centres */
--      for (int i = 0; i < nr_cells; ++i) {
--        centres[i * 3 + 0] *= factor;
--        centres[i * 3 + 1] *= factor;
--        centres[i * 3 + 2] *= factor;
--      }
--
--      /* Convert the cell widths */
--      cell_width[0] *= factor;
--      cell_width[1] *= factor;
--      cell_width[2] *= factor;
--    }
--
--    /* Write some meta-information first */
--    hid_t h_subgrp =
--        H5Gcreate(h_grp, "Meta-data", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
--    if (h_subgrp < 0) error("Error while creating meta-data sub-group");
--    io_write_attribute(h_subgrp, "nr_cells", INT, &nr_cells, 1);
--    io_write_attribute(h_subgrp, "size", DOUBLE, cell_width, 3);
--    io_write_attribute(h_subgrp, "dimension", INT, cdim, 3);
--    H5Gclose(h_subgrp);
--
--    /* Write the centres to the group */
--    hsize_t shape[2] = {nr_cells, 3};
--    hid_t h_space = H5Screate(H5S_SIMPLE);
--    if (h_space < 0) error("Error while creating data space for cell centres");
--    hid_t h_err = H5Sset_extent_simple(h_space, 2, shape, shape);
--    if (h_err < 0)
--      error("Error while changing shape of gas offsets data space.");
--    hid_t h_data = H5Dcreate(h_grp, "Centres", io_hdf5_type(DOUBLE), h_space,
--                             H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
--    if (h_data < 0) error("Error while creating dataspace for gas offsets.");
--    h_err = H5Dwrite(h_data, io_hdf5_type(DOUBLE), h_space, H5S_ALL,
--                     H5P_DEFAULT, centres);
--    if (h_err < 0) error("Error while writing centres.");
--    H5Dclose(h_data);
--    H5Sclose(h_space);
--
--    /* Group containing the offsets for each particle type */
--    h_subgrp =
--        H5Gcreate(h_grp, "Offsets", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
--    if (h_subgrp < 0) error("Error while creating offsets sub-group");
--
--    if (global_counts[swift_type_gas] > 0) {
--
--      shape[0] = nr_cells;
--      shape[1] = 1;
--      h_space = H5Screate(H5S_SIMPLE);
--      if (h_space < 0) error("Error while creating data space for gas offsets");
--      h_err = H5Sset_extent_simple(h_space, 1, shape, shape);
--      if (h_err < 0)
--        error("Error while changing shape of gas offsets data space.");
--      h_data = H5Dcreate(h_subgrp, "PartType0", io_hdf5_type(LONGLONG), h_space,
--                         H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
--      if (h_data < 0) error("Error while creating dataspace for gas offsets.");
--      h_err = H5Dwrite(h_data, io_hdf5_type(LONGLONG), h_space, H5S_ALL,
--                       H5P_DEFAULT, offset_part);
--      if (h_err < 0) error("Error while writing gas offsets.");
--      H5Dclose(h_data);
--      H5Sclose(h_space);
--    }
--
--    if (global_counts[swift_type_dark_matter] > 0) {
--
--      shape[0] = nr_cells;
--      shape[1] = 1;
--      h_space = H5Screate(H5S_SIMPLE);
--      if (h_space < 0) error("Error while creating data space for DM offsets");
--      h_err = H5Sset_extent_simple(h_space, 1, shape, shape);
--      if (h_err < 0)
--        error("Error while changing shape of DM offsets data space.");
--      h_data = H5Dcreate(h_subgrp, "PartType1", io_hdf5_type(LONGLONG), h_space,
--                         H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
--      if (h_data < 0) error("Error while creating dataspace for DM offsets.");
--      h_err = H5Dwrite(h_data, io_hdf5_type(LONGLONG), h_space, H5S_ALL,
--                       H5P_DEFAULT, offset_gpart);
--      if (h_err < 0) error("Error while writing DM offsets.");
--      H5Dclose(h_data);
--      H5Sclose(h_space);
--    }
--
--    if (global_counts[swift_type_stars] > 0) {
--
--      shape[0] = nr_cells;
--      shape[1] = 1;
--      h_space = H5Screate(H5S_SIMPLE);
--      if (h_space < 0)
--        error("Error while creating data space for stars offsets");
--      h_err = H5Sset_extent_simple(h_space, 1, shape, shape);
--      if (h_err < 0)
--        error("Error while changing shape of stars offsets data space.");
--      h_data = H5Dcreate(h_subgrp, "PartType4", io_hdf5_type(LONGLONG), h_space,
--                         H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
--      if (h_data < 0) error("Error while creating dataspace for star offsets.");
--      h_err = H5Dwrite(h_data, io_hdf5_type(LONGLONG), h_space, H5S_ALL,
--                       H5P_DEFAULT, offset_spart);
--      if (h_err < 0) error("Error while writing star offsets.");
--      H5Dclose(h_data);
--      H5Sclose(h_space);
--    }
--
--    H5Gclose(h_subgrp);
--
--    /* Group containing the counts for each particle type */
--    h_subgrp =
--        H5Gcreate(h_grp, "Counts", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
--    if (h_subgrp < 0) error("Error while creating counts sub-group");
--
--    if (global_counts[swift_type_gas] > 0) {
--
--      shape[0] = nr_cells;
--      shape[1] = 1;
--      h_space = H5Screate(H5S_SIMPLE);
--      if (h_space < 0) error("Error while creating data space for gas counts");
--      h_err = H5Sset_extent_simple(h_space, 1, shape, shape);
--      if (h_err < 0)
--        error("Error while changing shape of gas counts data space.");
--      h_data = H5Dcreate(h_subgrp, "PartType0", io_hdf5_type(LONGLONG), h_space,
--                         H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
--      if (h_data < 0) error("Error while creating dataspace for gas counts.");
--      h_err = H5Dwrite(h_data, io_hdf5_type(LONGLONG), h_space, H5S_ALL,
--                       H5P_DEFAULT, count_part);
--      if (h_err < 0) error("Error while writing gas counts.");
--      H5Dclose(h_data);
--      H5Sclose(h_space);
--    }
--
--    if (global_counts[swift_type_dark_matter] > 0) {
--
--      shape[0] = nr_cells;
--      shape[1] = 1;
--      h_space = H5Screate(H5S_SIMPLE);
--      if (h_space < 0) error("Error while creating data space for DM counts");
--      h_err = H5Sset_extent_simple(h_space, 1, shape, shape);
--      if (h_err < 0)
--        error("Error while changing shape of DM counts data space.");
--      h_data = H5Dcreate(h_subgrp, "PartType1", io_hdf5_type(LONGLONG), h_space,
--                         H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
--      if (h_data < 0) error("Error while creating dataspace for DM counts.");
--      h_err = H5Dwrite(h_data, io_hdf5_type(LONGLONG), h_space, H5S_ALL,
--                       H5P_DEFAULT, count_gpart);
--      if (h_err < 0) error("Error while writing DM counts.");
--      H5Dclose(h_data);
--      H5Sclose(h_space);
--    }
--
--    if (global_counts[swift_type_stars] > 0) {
--
--      shape[0] = nr_cells;
--      shape[1] = 1;
--      h_space = H5Screate(H5S_SIMPLE);
--      if (h_space < 0)
--        error("Error while creating data space for stars counts");
--      h_err = H5Sset_extent_simple(h_space, 1, shape, shape);
--      if (h_err < 0)
--        error("Error while changing shape of stars counts data space.");
--      h_data = H5Dcreate(h_subgrp, "PartType4", io_hdf5_type(LONGLONG), h_space,
--                         H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
--      if (h_data < 0) error("Error while creating dataspace for star counts.");
--      h_err = H5Dwrite(h_data, io_hdf5_type(LONGLONG), h_space, H5S_ALL,
--                       H5P_DEFAULT, count_spart);
--      if (h_err < 0) error("Error while writing star counts.");
--      H5Dclose(h_data);
--      H5Sclose(h_space);
--    }
--
--    H5Gclose(h_subgrp);
--  }
--
--  /* Free everything we allocated */
--  free(centres);
--  free(count_part);
--  free(count_gpart);
--  free(count_spart);
--  free(offset_part);
--  free(offset_gpart);
--  free(offset_spart);
--}
--
- #endif /* HAVE_HDF5 */
- 
- /**
-diff --git a/src/cooling/EAGLE/cooling.c b/src/cooling/EAGLE/cooling.c
-index a08a427..b94d60f 100644
---- a/src/cooling/EAGLE/cooling.c
-+++ b/src/cooling/EAGLE/cooling.c
-@@ -79,38 +79,32 @@ __attribute__((always_inline)) INLINE void get_redshift_index(
-     float z, int *z_index, float *dz,
-     struct cooling_function_data *restrict cooling) {
- 
--  /* Before the earliest redshift or before hydrogen reionization, flag for
-+  /* before the earliest redshift or before hydrogen reionization, flag for
-    * collisional cooling */
-   if (z > cooling->H_reion_z) {
-     *z_index = eagle_cooling_N_redshifts;
-     *dz = 0.0;
-   }
--
--  /* From reionization use the cooling tables */
-+  /* from reionization use the cooling tables */
-   else if (z > cooling->Redshifts[eagle_cooling_N_redshifts - 1] &&
-            z <= cooling->H_reion_z) {
-     *z_index = eagle_cooling_N_redshifts + 1;
-     *dz = 0.0;
-   }
--
--  /* At the end, just use the last value */
-+  /* at the end, just use the last value */
-   else if (z <= cooling->Redshifts[0]) {
-     *z_index = 0;
-     *dz = 0.0;
--  }
--
--  /* Normal case: search... */
--  else {
-+  } else {
- 
-     /* start at the previous index and search */
--    for (int i = cooling->previous_z_index; i >= 0; i--) {
--      if (z > cooling->Redshifts[i]) {
--
--        *z_index = i;
--        cooling->previous_z_index = i;
-+    for (int iz = cooling->previous_z_index; iz >= 0; iz--) {
-+      if (z > cooling->Redshifts[iz]) {
- 
--        *dz = (z - cooling->Redshifts[i]) /
--              (cooling->Redshifts[i + 1] - cooling->Redshifts[i]);
-+        *z_index = iz;
-+        cooling->previous_z_index = iz;
-+        *dz = (z - cooling->Redshifts[iz]) /
-+              (cooling->Redshifts[iz + 1] - cooling->Redshifts[iz]);
-         break;
-       }
-     }
-@@ -133,40 +127,20 @@ void cooling_update(const struct cosmology *cosmo,
-   /* Current redshift */
-   const float redshift = cosmo->z;
- 
--  /* What is the current table index along the redshift axis? */
-+  /* Get index along the redshift index of the tables */
-   int z_index = -1;
-   float dz = 0.f;
--  get_redshift_index(redshift, &z_index, &dz, cooling);
--  cooling->dz = dz;
--
--  /* Do we already have the correct tables loaded? */
--  if (cooling->z_index == z_index) return;
--
--  /* Which table should we load ? */
--  if (z_index >= eagle_cooling_N_redshifts) {
--
--    if (z_index == eagle_cooling_N_redshifts + 1) {
--
--      /* Bewtween re-ionization and first table */
--      get_redshift_invariant_table(cooling, /* photodis=*/0);
--
--    } else {
--
--      /* Above re-ionization */
--      get_redshift_invariant_table(cooling, /* photodis=*/1);
--    }
--
-+  if (redshift > cooling->H_reion_z) {
-+    z_index = -2;
-+  } else if (redshift > cooling->Redshifts[eagle_cooling_N_redshifts - 1]) {
-+    z_index = -1;
-   } else {
--
--    /* Normal case: two tables bracketing the current z */
--    const int low_z_index = z_index;
--    const int high_z_index = z_index + 1;
--
--    get_cooling_table(cooling, low_z_index, high_z_index);
-+    get_redshift_index(redshift, &z_index, &dz, cooling);
-   }
--
--  /* Store the currently loaded index */
-   cooling->z_index = z_index;
-+  cooling->dz = dz;
-+
-+  eagle_check_cooling_tables(cooling, restart_flag);
- }
- 
- /**
-@@ -492,10 +466,7 @@ void cooling_cool_part(const struct phys_const *restrict phys_const,
-      (See cosmology theory document for the derivation) */
-   const double delta_redshift = -dt * cosmo->H * cosmo->a_inv;
- 
--  /* Get this particle's abundance ratios compared to solar
--   * Note that we need to add S and Ca that are in the tables but not tracked
--   * by the particles themselves.
--   * The order is [H, He, C, N, O, Ne, Mg, Si, S, Ca, Fe] */
-+  /* Get this particle's abundance ratios */
-   float abundance_ratio[chemistry_element_count + 2];
-   abundance_ratio_to_solar(p, cooling, abundance_ratio);
- 
-@@ -764,32 +735,28 @@ void cooling_init_backend(struct swift_params *parameter_file,
-                           struct cooling_function_data *cooling) {
- 
-   /* read some parameters */
--  parser_get_param_string(parameter_file, "EAGLECooling:dir_name",
-+  parser_get_param_string(parameter_file, "EagleCooling:filename",
-                           cooling->cooling_table_path);
--  cooling->H_reion_z =
--      parser_get_param_float(parameter_file, "EAGLECooling:H_reion_z");
-+  cooling->H_reion_z = parser_get_param_float(
-+      parameter_file, "EagleCooling:reionisation_redshift");
-+  cooling->calcium_over_silicon_ratio = parser_get_param_float(
-+      parameter_file, "EAGLEChemistry:CalciumOverSilicon");
-+  cooling->sulphur_over_silicon_ratio = parser_get_param_float(
-+      parameter_file, "EAGLEChemistry:SulphurOverSilicon");
-   cooling->He_reion_z_centre =
--      parser_get_param_float(parameter_file, "EAGLECooling:He_reion_z_centre");
-+      parser_get_param_float(parameter_file, "EagleCooling:He_reion_z_centre");
-   cooling->He_reion_z_sigma =
--      parser_get_param_float(parameter_file, "EAGLECooling:He_reion_z_sigma");
-+      parser_get_param_float(parameter_file, "EagleCooling:He_reion_z_sigma");
-   cooling->He_reion_heat_cgs =
--      parser_get_param_float(parameter_file, "EAGLECooling:He_reion_ev_p_H");
--
--  /* Optional parameters to correct the abundances */
--  cooling->Ca_over_Si_ratio_in_solar = parser_get_opt_param_float(
--      parameter_file, "EAGLECooling:Ca_over_Si_in_solar", 1.f);
--  cooling->S_over_Si_ratio_in_solar = parser_get_opt_param_float(
--      parameter_file, "EAGLECooling:S_over_Si_in_solar", 1.f);
-+      parser_get_param_float(parameter_file, "EagleCooling:He_reion_ev_pH");
- 
--  /* Convert to cgs (units used internally by the cooling routines) */
-+  /* convert to cgs */
-   cooling->He_reion_heat_cgs *=
-       phys_const->const_electron_volt *
-       units_cgs_conversion_factor(us, UNIT_CONV_ENERGY);
- 
--  /* Read in the list of redshifts */
-+  /* read in cooling table header */
-   get_cooling_redshifts(cooling);
--
--  /* Read in cooling table header */
-   char fname[eagle_table_path_name_length + 12];
-   sprintf(fname, "%sz_0.000.hdf5", cooling->cooling_table_path);
-   read_cooling_header(fname, cooling);
-@@ -797,7 +764,7 @@ void cooling_init_backend(struct swift_params *parameter_file,
-   /* Allocate space for cooling tables */
-   allocate_cooling_tables(cooling);
- 
--  /* Compute conversion factors */
-+  /* compute conversion factors */
-   cooling->internal_energy_to_cgs =
-       units_cgs_conversion_factor(us, UNIT_CONV_ENERGY_PER_UNIT_MASS);
-   cooling->internal_energy_from_cgs = 1. / cooling->internal_energy_to_cgs;
-@@ -839,16 +806,14 @@ void cooling_init_backend(struct swift_params *parameter_file,
-                               cooling->T_CMB_0 * cooling->T_CMB_0 *
-                               cooling->T_CMB_0;
- 
--  /* Set the redshift indices to invalid values */
--  cooling->z_index = -10;
--  cooling->previous_z_index = eagle_cooling_N_redshifts + 2;
--
-+  /* set low_z_index to -10 to indicate we haven't read any tables yet */
-+  cooling->low_z_index = -10;
-   /* set previous_z_index and to last value of redshift table*/
-   cooling->previous_z_index = eagle_cooling_N_redshifts - 2;
- 
-   /* Check if we are running with the newton scheme */
-   cooling->newton_flag = parser_get_opt_param_int(
--      parameter_file, "EAGLECooling:newton_integration", 0);
-+      parameter_file, "EagleCooling:newton_integration", 0);
- }
- 
- /**
-@@ -869,13 +834,10 @@ void cooling_restore_tables(struct cooling_function_data *cooling,
-   sprintf(fname, "%sz_0.000.hdf5", cooling->cooling_table_path);
-   read_cooling_header(fname, cooling);
- 
--  /* Allocate memory for the tables */
-+  /* Read relevant cooling tables.
-+   * Third variable in cooling_update flag to mark restart*/
-   allocate_cooling_tables(cooling);
--
--  /* Force a re-read of the cooling tables */
--  cooling->z_index = -10;
--  cooling->previous_z_index = eagle_cooling_N_redshifts + 2;
--  cooling_update(cosmo, cooling, /*restart_flag=*/1);
-+  cooling_update(cosmo, cooling, /*restart=*/1);
- }
- 
- /**
-@@ -904,7 +866,6 @@ void cooling_clean(struct cooling_function_data *cooling) {
-   free(cooling->HeFrac);
-   free(cooling->Therm);
-   free(cooling->SolarAbundances);
--  free(cooling->SolarAbundances_inv);
- 
-   /* Free the tables */
-   free(cooling->table.metal_heating);
-diff --git a/src/cooling/EAGLE/cooling_tables.c b/src/cooling/EAGLE/cooling_tables.c
-index c66b7eb..4f3aed0 100644
---- a/src/cooling/EAGLE/cooling_tables.c
-+++ b/src/cooling/EAGLE/cooling_tables.c
-@@ -31,7 +31,6 @@
- #include <string.h>
- 
- /* Local includes. */
--#include "chemistry_struct.h"
- #include "cooling_struct.h"
- #include "cooling_tables.h"
- #include "error.h"
-@@ -40,7 +39,7 @@
- /**
-  * @brief Names of the elements in the order they are stored in the files
-  */
--static const char *eagle_tables_element_names[eagle_cooling_N_metal] = {
-+static const char *eagle_tables_element_names[9] = {
-     "Carbon",  "Nitrogen", "Oxygen",  "Neon", "Magnesium",
-     "Silicon", "Sulphur",  "Calcium", "Iron"};
- 
-@@ -212,10 +211,6 @@ void read_cooling_header(const char *fname,
-   if (N_SolarAbundances != eagle_cooling_N_abundances)
-     error("Invalid solar abundances array length.");
- 
--  /* Check value */
--  if (N_SolarAbundances != chemistry_element_count + 2)
--    error("Number of abundances not compatible with the chemistry model.");
--
-   dataset = H5Dopen(tempfile_id, "/Header/Number_of_metals", H5P_DEFAULT);
-   status = H5Dread(dataset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
-                    &N_Elements);
-@@ -242,10 +237,6 @@ void read_cooling_header(const char *fname,
-   if (posix_memalign((void **)&cooling->SolarAbundances, SWIFT_STRUCT_ALIGNMENT,
-                      N_SolarAbundances * sizeof(float)) != 0)
-     error("Failed to allocate Solar abundances table");
--  if (posix_memalign((void **)&cooling->SolarAbundances_inv,
--                     SWIFT_STRUCT_ALIGNMENT,
--                     N_SolarAbundances * sizeof(float)) != 0)
--    error("Failed to allocate Solar abundances inverses table");
- 
-   /* read in values for each of the arrays */
-   dataset = H5Dopen(tempfile_id, "/Solar/Temperature_bins", H5P_DEFAULT);
-@@ -294,10 +285,6 @@ void read_cooling_header(const char *fname,
-   for (int i = 0; i < N_nH; i++) {
-     cooling->nH[i] = log10(cooling->nH[i]);
-   }
--  /* Compute inverse of solar mass fractions */
--  for (int i = 0; i < N_SolarAbundances; ++i) {
--    cooling->SolarAbundances_inv[i] = 1.f / cooling->SolarAbundances[i];
--  }
- 
- #else
-   error("Need HDF5 to read cooling tables");
-@@ -359,10 +346,9 @@ void allocate_cooling_tables(struct cooling_function_data *restrict cooling) {
-  * used to obtain temperature of particle)
-  *
-  * @param cooling #cooling_function_data structure
-- * @param photodis Are we loading the photo-dissociation table?
-  */
--void get_redshift_invariant_table(
--    struct cooling_function_data *restrict cooling, const int photodis) {
-+static void get_redshift_invariant_table(
-+    struct cooling_function_data *restrict cooling) {
- #ifdef HAVE_HDF5
- 
-   /* Temporary tables */
-@@ -391,12 +377,10 @@ void get_redshift_invariant_table(
- 
-   /* Decide which high redshift table to read. Indices set in cooling_update */
-   char filename[eagle_table_path_name_length + 21];
--  if (photodis) {
--    sprintf(filename, "%sz_photodis.hdf5", cooling->cooling_table_path);
--    message("Reading cooling table 'z_photodis.hdf5'");
--  } else {
-+  if (cooling->low_z_index == -1) {
-     sprintf(filename, "%sz_8.989nocompton.hdf5", cooling->cooling_table_path);
--    message("Reading cooling table 'z_8.989nocompton.hdf5'");
-+  } else if (cooling->low_z_index == -2) {
-+    sprintf(filename, "%sz_photodis.hdf5", cooling->cooling_table_path);
-   }
- 
-   hid_t file_id = H5Fopen(filename, H5F_ACC_RDONLY, H5P_DEFAULT);
-@@ -557,11 +541,8 @@ void get_redshift_invariant_table(
-  * used to obtain temperature of particle)
-  *
-  * @param cooling #cooling_function_data structure
-- * @param low_z_index Index of the lowest redshift table to load.
-- * @param high_z_index Index of the highest redshift table to load.
-  */
--void get_cooling_table(struct cooling_function_data *restrict cooling,
--                       const int low_z_index, const int high_z_index) {
-+static void get_cooling_table(struct cooling_function_data *restrict cooling) {
- 
- #ifdef HAVE_HDF5
- 
-@@ -591,10 +572,11 @@ void get_cooling_table(struct cooling_function_data *restrict cooling,
- 
-   /* Read in tables, transpose so that values for indices which vary most are
-    * adjacent. Repeat for redshift above and redshift below current value.  */
--  for (int z_index = low_z_index; z_index <= high_z_index; z_index++) {
-+  for (int z_index = cooling->low_z_index; z_index <= cooling->high_z_index;
-+       z_index++) {
- 
-     /* Index along redhsift dimension for the subset of tables we read */
--    const int local_z_index = z_index - low_z_index;
-+    const int local_z_index = z_index - cooling->low_z_index;
- 
- #ifdef SWIFT_DEBUG_CHECKS
-     if (local_z_index >= eagle_cooling_N_loaded_redshifts)
-@@ -602,12 +584,9 @@ void get_cooling_table(struct cooling_function_data *restrict cooling,
- #endif
- 
-     /* Open table for this redshift index */
--    char fname[eagle_table_path_name_length + 12];
-+    char fname[eagle_table_path_name_length + 32];
-     sprintf(fname, "%sz_%1.3f.hdf5", cooling->cooling_table_path,
-             cooling->Redshifts[z_index]);
--    message("Reading cooling table 'z_%1.3f.hdf5'",
--            cooling->Redshifts[z_index]);
--
-     hid_t file_id = H5Fopen(fname, H5F_ACC_RDONLY, H5P_DEFAULT);
-     if (file_id < 0) error("unable to open file %s", fname);
- 
-@@ -748,10 +727,48 @@ void get_cooling_table(struct cooling_function_data *restrict cooling,
-   free(he_electron_abundance);
- 
- #ifdef SWIFT_DEBUG_CHECKS
--  message("Done reading in general cooling table");
-+  message("done reading in general cooling table");
- #endif
- 
- #else
-   error("Need HDF5 to read cooling tables");
- #endif
- }
-+
-+/**
-+ * @brief Constructs the data structure containting the relevant cooling tables
-+ * for the redshift index (set in cooling_update)
-+ *
-+ * @param cooling #cooling_function_data structure
-+ */
-+static void eagle_readtable(struct cooling_function_data *restrict cooling) {
-+
-+  if (cooling->z_index < 0) {
-+    /* z_index is set to < 0 in cooling_update if need
-+     * to read any of the high redshift tables */
-+    get_redshift_invariant_table(cooling);
-+  } else {
-+    get_cooling_table(cooling);
-+  }
-+}
-+
-+/**
-+ * @brief Checks the tables that are currently loaded in memory and read
-+ * new ones if necessary.
-+ *
-+ * @param cooling The #cooling_function_data we play with.
-+ * @param restart_flag Flag indicating if we are restarting a run
-+ */
-+void eagle_check_cooling_tables(struct cooling_function_data *restrict cooling,
-+                                const int restart_flag) {
-+
-+  /* Do we already have the right table in memory? */
-+  if (cooling->low_z_index == cooling->z_index && restart_flag != 0) return;
-+
-+  /* Record the table indices */
-+  cooling->low_z_index = cooling->z_index;
-+  cooling->high_z_index = cooling->z_index + 1;
-+
-+  /* Load the damn thing */
-+  eagle_readtable(cooling);
-+}
-diff --git a/src/cosmology.c b/src/cosmology.c
-index be23343..4718ed5 100644
---- a/src/cosmology.c
-+++ b/src/cosmology.c
-@@ -576,8 +576,6 @@ void cosmology_init_no_cosmo(struct cosmology *c) {
-   c->a_dot = 0.;
-   c->time = 0.;
-   c->universe_age_at_present_day = 0.;
--  c->Hubble_time = 0.;
--  c->lookback_time = 0.;
- 
-   /* Initialise the interpolation tables */
-   c->drift_fac_interp_table = NULL;
-diff --git a/src/debug.c b/src/debug.c
-index 809d704..8b2012e 100644
---- a/src/debug.c
-+++ b/src/debug.c
-@@ -179,6 +179,14 @@ int checkSpacehmax(struct space *s) {
-     }
-   }
- 
-+  float cell_stars_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].stars.h_max > cell_stars_h_max) {
-+      cell_stars_h_max = s->cells_top[k].stars.h_max;
-+    }
-+  }
-+
-   /* Now all particles. */
-   float part_h_max = 0.0f;
-   for (size_t k = 0; k < s->nr_parts; k++) {
-@@ -187,10 +195,21 @@ int checkSpacehmax(struct space *s) {
-     }
-   }
- 
-+  /* Now all the sparticles. */
-+  float spart_h_max = 0.0f;
-+  for (size_t k = 0; k < s->nr_sparts; k++) {
-+    if (s->sparts[k].h > spart_h_max) {
-+      spart_h_max = s->sparts[k].h;
-+    }
-+  }
-+
-   /*  If within some epsilon we are OK. */
--  if (fabsf(cell_h_max - part_h_max) <= FLT_EPSILON) return 1;
-+  if (fabsf(cell_h_max - part_h_max) <= FLT_EPSILON &&
-+      fabsf(cell_stars_h_max - spart_h_max) <= FLT_EPSILON)
-+    return 1;
- 
-   /* There is a problem. Hunt it down. */
-+  /* part */
-   for (int k = 0; k < s->nr_cells; k++) {
-     if (s->cells_top[k].nodeID == s->e->nodeID) {
-       if (s->cells_top[k].hydro.h_max > part_h_max) {
-@@ -207,6 +226,23 @@ int checkSpacehmax(struct space *s) {
-     }
-   }
- 
-+  /* spart */
-+  for (int k = 0; k < s->nr_cells; k++) {
-+    if (s->cells_top[k].nodeID == s->e->nodeID) {
-+      if (s->cells_top[k].stars.h_max > spart_h_max) {
-+        message("cell %d is inconsistent (%f > %f)", k,
-+                s->cells_top[k].stars.h_max, spart_h_max);
-+      }
-+    }
-+  }
-+
-+  for (size_t k = 0; k < s->nr_sparts; k++) {
-+    if (s->sparts[k].h > cell_stars_h_max) {
-+      message("spart %lld is inconsistent (%f > %f)", s->sparts[k].id,
-+              s->sparts[k].h, cell_stars_h_max);
-+    }
-+  }
-+
-   return 0;
- }
- 
-@@ -225,6 +261,8 @@ int checkCellhdxmax(const struct cell *c, int *depth) {
- 
-   float h_max = 0.0f;
-   float dx_max = 0.0f;
-+  float stars_h_max = 0.0f;
-+  float stars_dx_max = 0.0f;
-   int result = 1;
- 
-   const double loc_min[3] = {c->loc[0], c->loc[1], c->loc[2]};
-@@ -260,6 +298,33 @@ int checkCellhdxmax(const struct cell *c, int *depth) {
-     dx_max = max(dx_max, sqrt(dx2));
-   }
- 
-+  const size_t nr_sparts = c->stars.count;
-+  struct spart *sparts = c->stars.parts;
-+  for (size_t k = 0; k < nr_sparts; k++) {
-+
-+    struct spart *const sp = &sparts[k];
-+
-+    if (sp->x[0] < loc_min[0] || sp->x[0] >= loc_max[0] ||
-+        sp->x[1] < loc_min[1] || sp->x[1] >= loc_max[1] ||
-+        sp->x[2] < loc_min[2] || sp->x[2] >= loc_max[2]) {
-+
-+      message(
-+          "Inconsistent part position p->x=[%e %e %e], c->loc=[%e %e %e] "
-+          "c->width=[%e %e %e]",
-+          sp->x[0], sp->x[1], sp->x[2], c->loc[0], c->loc[1], c->loc[2],
-+          c->width[0], c->width[1], c->width[2]);
-+
-+      result = 0;
-+    }
-+
-+    const float dx2 = sp->x_diff[0] * sp->x_diff[0] +
-+                      sp->x_diff[1] * sp->x_diff[1] +
-+                      sp->x_diff[2] * sp->x_diff[2];
-+
-+    stars_h_max = max(stars_h_max, sp->h);
-+    stars_dx_max = max(stars_dx_max, sqrt(dx2));
-+  }
-+
-   if (c->split) {
-     for (int k = 0; k < 8; k++) {
-       if (c->progeny[k] != NULL) {
-@@ -283,6 +348,19 @@ int checkCellhdxmax(const struct cell *c, int *depth) {
-     result = 0;
-   }
- 
-+  if (c->stars.h_max != stars_h_max) {
-+    message("%d Inconsistent stars_h_max: cell %f != parts %f", *depth,
-+            c->stars.h_max, stars_h_max);
-+    message("location: %f %f %f", c->loc[0], c->loc[1], c->loc[2]);
-+    result = 0;
-+  }
-+  if (c->stars.dx_max_part != stars_dx_max) {
-+    message("%d Inconsistent stars_dx_max: %f != %f", *depth,
-+            c->stars.dx_max_part, stars_dx_max);
-+    message("location: %f %f %f", c->loc[0], c->loc[1], c->loc[2]);
-+    result = 0;
-+  }
-+
-   return result;
- }
- 
-diff --git a/src/engine.c b/src/engine.c
-index 2c5582f..964c43c 100644
---- a/src/engine.c
-+++ b/src/engine.c
-@@ -113,8 +113,7 @@ const char *engine_policy_names[] = {"none",
-                                      "stars",
-                                      "structure finding",
-                                      "star formation",
--                                     "feedback",
--                                     "time-step limiter"};
-+                                     "feedback"};
- 
- /** The rank of the engine as a global variable (for messages). */
- int engine_rank;
-@@ -1908,10 +1907,6 @@ int engine_estimate_nr_tasks(struct engine *e) {
- #endif
- #endif
-   }
--  if (e->policy & engine_policy_limiter) {
--    n1 += 18;
--    n2 += 1;
--  }
-   if (e->policy & engine_policy_self_gravity) {
-     n1 += 125;
-     n2 += 8;
-@@ -2590,20 +2585,15 @@ void engine_skip_force_and_kick(struct engine *e) {
-     /* Skip everything that updates the particles */
-     if (t->type == task_type_drift_part || t->type == task_type_drift_gpart ||
-         t->type == task_type_kick1 || t->type == task_type_kick2 ||
--        t->type == task_type_timestep ||
--        t->type == task_type_timestep_limiter ||
--        t->subtype == task_subtype_force ||
--        t->subtype == task_subtype_limiter || t->subtype == task_subtype_grav ||
--        t->type == task_type_end_force ||
-+        t->type == task_type_timestep || t->subtype == task_subtype_force ||
-+        t->subtype == task_subtype_grav || t->type == task_type_end_force ||
-         t->type == task_type_grav_long_range || t->type == task_type_grav_mm ||
--        t->type == task_type_grav_down || t->type == task_type_cooling ||
--	t->subtype == task_subtype_stars_feedback)
-+        t->type == task_type_grav_down || t->type == task_type_cooling)
-       t->skip = 1;
-   }
- 
-   /* Run through the cells and clear some flags. */
-   space_map_cells_pre(e->s, 1, cell_clear_drift_flags, NULL);
--  space_map_cells_pre(e->s, 1, cell_clear_limiter_flags, NULL);
- }
- 
- /**
-@@ -2750,7 +2740,7 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs,
-   if (!flag_entropy_ICs) {
- 
-     if (e->nodeID == 0) message("Converting internal energy variable.");
--    
-+
-     space_convert_quantities(e->s, e->verbose);
- 
-     /* Correct what we did (e.g. in PE-SPH, need to recompute rho_bar) */
-@@ -2781,7 +2771,7 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs,
- 
-   /* Prepare all the tasks again for a new round */
-   engine_marktasks(e);
--    
-+
-   /* No drift this time */
-   engine_skip_drift(e);
- 
-@@ -2813,11 +2803,6 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs,
-     gravity_exact_force_check(e->s, e, 1e-1);
- #endif
- 
--#ifdef SWIFT_DEBUG_CHECKS
--  /* Make sure all woken-up particles have been processed */
--  space_check_limiter(e->s);
--#endif
--
-   /* Recover the (integer) end of the next time-step */
-   engine_collect_end_of_step(e, 1);
- 
-@@ -3075,11 +3060,6 @@ void engine_step(struct engine *e) {
-     gravity_exact_force_check(e->s, e, 1e-1);
- #endif
- 
--#ifdef SWIFT_DEBUG_CHECKS
--  /* Make sure all woken-up particles have been processed */
--  space_check_limiter(e->s);
--#endif
--
-   /* Collect information about the next time-step */
-   engine_collect_end_of_step(e, 1);
-   e->forcerebuild = e->collect_group1.forcerebuild;
-diff --git a/src/engine_maketasks.c b/src/engine_maketasks.c
-index c39659b..4d03d80 100644
---- a/src/engine_maketasks.c
-+++ b/src/engine_maketasks.c
-@@ -204,19 +204,94 @@ void engine_addtasks_send_hydro(struct engine *e, struct cell *ci,
- }
- 
- /**
-+ * @brief Add send tasks for the stars pairs to a hierarchy of cells.
-+ *
-+ * @param e The #engine.
-+ * @param ci The sending #cell.
-+ * @param cj Dummy cell containing the nodeID of the receiving node.
-+ * @param t_xv The send_xv #task, if it has already been created.
-+ * @param t_feed The send_feed #task, if it has already been created.
-+ */
-+void engine_addtasks_send_stars(struct engine *e, struct cell *ci,
-+                                struct cell *cj, struct task *t_xv,
-+                                struct task *t_feed) {
-+
-+#ifdef WITH_MPI
-+  struct link *l = NULL;
-+  struct scheduler *s = &e->sched;
-+  const int nodeID = cj->nodeID;
-+
-+  /* Check if any of the density tasks are for the target node. */
-+  for (l = ci->stars.density; l != NULL; l = l->next)
-+    if (l->t->ci->nodeID == nodeID ||
-+        (l->t->cj != NULL && l->t->cj->nodeID == nodeID))
-+      break;
-+
-+  /* If so, attach send tasks. */
-+  if (l != NULL) {
-+    /* Get the task if created in hydro part */
-+    struct link *hydro = NULL;
-+    for (hydro = ci->mpi.hydro.send_xv; hydro != NULL; hydro = hydro->next) {
-+      if (hydro->t->ci->nodeID == nodeID ||
-+          (hydro->t->cj != NULL && hydro->t->cj->nodeID == nodeID)) {
-+        break;
-+      }
-+    }
-+
-+    if (t_xv == NULL) {
-+      /* Make sure this cell is tagged. */
-+      cell_ensure_tagged(ci);
-+
-+      /* Already exists, just need to get it */
-+      if (hydro != NULL) {
-+        t_xv = hydro->t;
-+
-+        /* This task does not exists, need to create it */
-+      } else {
-+
-+        /* Create the tasks and their dependencies? */
-+        t_xv = scheduler_addtask(s, task_type_send, task_subtype_xv,
-+                                 ci->mpi.tag, 0, ci, cj);
-+
-+        /* Drift before you send */
-+        scheduler_addunlock(s, ci->hydro.super->hydro.drift, t_xv);
-+      }
-+
-+      /* Create the tasks and their dependencies? */
-+      t_feed = scheduler_addtask(s, task_type_send, task_subtype_spart,
-+                                 ci->mpi.tag, 0, ci, cj);
-+
-+      /* Ghost before you send */
-+      scheduler_addunlock(s, ci->super->stars.ghost_out, t_feed);
-+    }
-+
-+    if (hydro == NULL) {
-+      engine_addlink(e, &ci->mpi.hydro.send_xv, t_xv);
-+    }
-+    engine_addlink(e, &ci->mpi.stars.send, t_feed);
-+  }
-+
-+  /* Recurse? */
-+  if (ci->split)
-+    for (int k = 0; k < 8; k++)
-+      if (ci->progeny[k] != NULL)
-+        engine_addtasks_send_stars(e, ci->progeny[k], cj, t_xv, t_feed);
-+
-+#else
-+  error("SWIFT was not compiled with MPI support.");
-+#endif
-+}
-+
-+/**
-  * @brief Add send tasks for the time-step to a hierarchy of cells.
-  *
-  * @param e The #engine.
-  * @param ci The sending #cell.
-  * @param cj Dummy cell containing the nodeID of the receiving node.
-  * @param t_ti The send_ti #task, if it has already been created.
-- * @param t_limiter The send_limiter #task, if already created.
-- * @param with_limiter Are we running with the time-step limiter?
-  */
- void engine_addtasks_send_timestep(struct engine *e, struct cell *ci,
--                                   struct cell *cj, struct task *t_ti,
--                                   struct task *t_limiter,
--                                   const int with_limiter) {
-+                                   struct cell *cj, struct task *t_ti) {
- 
- #ifdef WITH_MPI
-   struct link *l = NULL;
-@@ -236,6 +311,12 @@ void engine_addtasks_send_timestep(struct engine *e, struct cell *ci,
-           (l->t->cj != NULL && l->t->cj->nodeID == nodeID))
-         break;
- 
-+  if (l == NULL)
-+    for (l = ci->stars.density; l != NULL; l = l->next)
-+      if (l->t->ci->nodeID == nodeID ||
-+          (l->t->cj != NULL && l->t->cj->nodeID == nodeID))
-+        break;
-+
-   /* If found anything, attach send tasks. */
-   if (l != NULL) {
- 
-@@ -248,31 +329,19 @@ void engine_addtasks_send_timestep(struct engine *e, struct cell *ci,
-       t_ti = scheduler_addtask(s, task_type_send, task_subtype_tend,
-                                ci->mpi.tag, 0, ci, cj);
- 
--      if (with_limiter)
--        t_limiter = scheduler_addtask(s, task_type_send, task_subtype_limiter,
--                                      ci->mpi.tag, 0, ci, cj);
--
-       /* The super-cell's timestep task should unlock the send_ti task. */
-       scheduler_addunlock(s, ci->super->timestep, t_ti);
--      if (with_limiter) scheduler_addunlock(s, t_limiter, ci->super->timestep);
--      if (with_limiter)
--        scheduler_addunlock(s, t_limiter, ci->super->timestep_limiter);
--      if (with_limiter) scheduler_addunlock(s, ci->super->kick2, t_limiter);
--      if (with_limiter)
--        scheduler_addunlock(s, ci->super->timestep_limiter, t_ti);
-     }
- 
-     /* Add them to the local cell. */
-     engine_addlink(e, &ci->mpi.send_ti, t_ti);
--    if (with_limiter) engine_addlink(e, &ci->mpi.limiter.send, t_limiter);
-   }
- 
-   /* Recurse? */
-   if (ci->split)
-     for (int k = 0; k < 8; k++)
-       if (ci->progeny[k] != NULL)
--        engine_addtasks_send_timestep(e, ci->progeny[k], cj, t_ti, t_limiter,
--                                      with_limiter);
-+        engine_addtasks_send_timestep(e, ci->progeny[k], cj, t_ti);
- 
- #else
-   error("SWIFT was not compiled with MPI support.");
-@@ -349,6 +418,72 @@ void engine_addtasks_recv_hydro(struct engine *e, struct cell *c,
- }
- 
- /**
-+ * @brief Add recv tasks for stars pairs to a hierarchy of cells.
-+ *
-+ * @param e The #engine.
-+ * @param c The foreign #cell.
-+ * @param t_xv The recv_xv #task, if it has already been created.
-+ * @param t_feed The recv_feed #task, if it has already been created.
-+ */
-+void engine_addtasks_recv_stars(struct engine *e, struct cell *c,
-+                                struct task *t_xv, struct task *t_feed) {
-+
-+#ifdef WITH_MPI
-+  struct scheduler *s = &e->sched;
-+  int new_task = 0;
-+
-+  /* Have we reached a level where there are any stars (or hydro) tasks ? */
-+  if (t_xv == NULL && (c->stars.density != NULL || c->hydro.density != NULL)) {
-+
-+#ifdef SWIFT_DEBUG_CHECKS
-+    /* Make sure this cell has a valid tag. */
-+    if (c->mpi.tag < 0) error("Trying to receive from untagged cell.");
-+#endif  // SWIFT_DEBUG_CHECKS
-+
-+    /* Create the tasks. */
-+    if (c->mpi.hydro.recv_xv == NULL) {
-+      new_task = 1;
-+      t_xv = scheduler_addtask(s, task_type_recv, task_subtype_xv, c->mpi.tag,
-+                               0, c, NULL);
-+    } else {
-+      t_xv = c->mpi.hydro.recv_xv;
-+    }
-+    t_feed = scheduler_addtask(s, task_type_recv, task_subtype_spart,
-+                               c->mpi.tag, 0, c, NULL);
-+
-+    /* Need to sort task before feedback loop */
-+    scheduler_addunlock(s, t_feed, c->super->stars.sorts_foreign);
-+  }
-+
-+  c->mpi.hydro.recv_xv = t_xv;
-+  c->mpi.stars.recv = t_feed;
-+
-+  /* Add dependencies. */
-+  if (c->hydro.sorts != NULL && new_task) {
-+    scheduler_addunlock(s, t_xv, c->hydro.sorts);
-+  }
-+
-+  for (struct link *l = c->stars.density; l != NULL; l = l->next) {
-+    scheduler_addunlock(s, t_xv, l->t);
-+    scheduler_addunlock(s, l->t, t_feed);
-+  }
-+
-+  for (struct link *l = c->stars.feedback; l != NULL; l = l->next) {
-+    scheduler_addunlock(s, t_feed, l->t);
-+  }
-+
-+  /* Recurse? */
-+  if (c->split)
-+    for (int k = 0; k < 8; k++)
-+      if (c->progeny[k] != NULL)
-+        engine_addtasks_recv_stars(e, c->progeny[k], t_xv, t_feed);
-+
-+#else
-+  error("SWIFT was not compiled with MPI support.");
-+#endif
-+}
-+
-+/**
-  * @brief Add recv tasks for gravity pairs to a hierarchy of cells.
-  *
-  * @param e The #engine.
-@@ -396,12 +531,9 @@ void engine_addtasks_recv_gravity(struct engine *e, struct cell *c,
-  * @param e The #engine.
-  * @param c The foreign #cell.
-  * @param t_ti The recv_ti #task, if already been created.
-- * @param t_limiter The recv_limiter #task, if already created.
-- * @param with_limiter Are we running with the time-step limiter?
-  */
- void engine_addtasks_recv_timestep(struct engine *e, struct cell *c,
--                                   struct task *t_ti, struct task *t_limiter,
--                                   const int with_limiter) {
-+                                   struct task *t_ti) {
- 
- #ifdef WITH_MPI
-   struct scheduler *s = &e->sched;
-@@ -416,42 +548,24 @@ void engine_addtasks_recv_timestep(struct engine *e, struct cell *c,
- 
-     t_ti = scheduler_addtask(s, task_type_recv, task_subtype_tend, c->mpi.tag,
-                              0, c, NULL);
--
--    if (with_limiter)
--      t_limiter = scheduler_addtask(s, task_type_recv, task_subtype_limiter,
--                                    c->mpi.tag, 0, c, NULL);
-   }
- 
-   c->mpi.recv_ti = t_ti;
- 
--  for (struct link *l = c->grav.grav; l != NULL; l = l->next) {
-+  for (struct link *l = c->grav.grav; l != NULL; l = l->next)
-     scheduler_addunlock(s, l->t, t_ti);
--  }
--
--  if (with_limiter) {
--
--    for (struct link *l = c->hydro.force; l != NULL; l = l->next) {
--      scheduler_addunlock(s, l->t, t_limiter);
--    }
--
--    for (struct link *l = c->hydro.limiter; l != NULL; l = l->next) {
--      scheduler_addunlock(s, t_limiter, l->t);
--      scheduler_addunlock(s, l->t, t_ti);
--    }
- 
--  } else {
-+  for (struct link *l = c->hydro.force; l != NULL; l = l->next)
-+    scheduler_addunlock(s, l->t, t_ti);
- 
--    for (struct link *l = c->hydro.force; l != NULL; l = l->next) {
--      scheduler_addunlock(s, l->t, t_ti);
--    }
--  }
-+  for (struct link *l = c->stars.feedback; l != NULL; l = l->next)
-+    scheduler_addunlock(s, l->t, t_ti);
- 
-   /* Recurse? */
-   if (c->split)
-     for (int k = 0; k < 8; k++)
-       if (c->progeny[k] != NULL)
--        engine_addtasks_recv_timestep(e, c->progeny[k], t_ti, t_limiter,
--                                      with_limiter);
-+        engine_addtasks_recv_timestep(e, c->progeny[k], t_ti);
- 
- #else
-   error("SWIFT was not compiled with MPI support.");
-@@ -475,7 +589,6 @@ void engine_make_hierarchical_tasks_common(struct engine *e, struct cell *c) {
-   struct scheduler *s = &e->sched;
-   const int is_with_cooling = (e->policy & engine_policy_cooling);
-   const int is_with_star_formation = (e->policy & engine_policy_star_formation);
--  const int with_limiter = (e->policy & engine_policy_limiter);
- 
-   /* Are we in a super-cell ? */
-   if (c->super == c) {
-@@ -504,7 +617,7 @@ void engine_make_hierarchical_tasks_common(struct engine *e, struct cell *c) {
-                                        task_subtype_none, 0, 0, c, NULL);
- 
-       /* Subgrid tasks */
--      if (is_with_cooling && c->hydro.count_total > 0) {
-+      if (is_with_cooling) {
- 
-         c->hydro.cooling = scheduler_addtask(s, task_type_cooling,
-                                              task_subtype_none, 0, 0, c, NULL);
-@@ -516,7 +629,7 @@ void engine_make_hierarchical_tasks_common(struct engine *e, struct cell *c) {
-         scheduler_addunlock(s, c->end_force, c->kick2);
-       }
- 
--      if (is_with_star_formation && c->hydro.count_total > 0) {
-+      if (is_with_star_formation) {
- 
-         c->hydro.star_formation = scheduler_addtask(
-             s, task_type_star_formation, task_subtype_none, 0, 0, c, NULL);
-@@ -530,16 +643,6 @@ void engine_make_hierarchical_tasks_common(struct engine *e, struct cell *c) {
- 
-       scheduler_addunlock(s, c->timestep, c->kick1);
- 
--      /* Time-step limiting */
--      if (with_limiter) {
--        c->timestep_limiter = scheduler_addtask(
--            s, task_type_timestep_limiter, task_subtype_none, 0, 0, c, NULL);
--
--        /* Make sure it is not run before kick2 */
--        scheduler_addunlock(s, c->timestep, c->timestep_limiter);
--        scheduler_addunlock(s, c->timestep_limiter, c->kick1);
--      }
--
- #if defined(WITH_LOGGER)
-       scheduler_addunlock(s, c->kick1, c->logger);
- #endif
-@@ -661,11 +764,8 @@ void engine_add_stars_ghosts(struct engine *e, struct cell *c,
-                              struct task *stars_ghost_in,
-                              struct task *stars_ghost_out) {
- 
--  /* Abort as there are no star particles here? */
--  if (c->stars.count_total == 0) return;
--
-   /* If we have reached the leaf OR have to few particles to play with*/
--  if (!c->split || c->stars.count_total < engine_max_sparts_per_ghost) {
-+  if (!c->split || c->stars.count < engine_max_sparts_per_ghost) {
- 
-     /* Add the ghost task and its dependencies */
-     struct scheduler *s = &e->sched;
-@@ -688,11 +788,8 @@ void engine_add_stars_ghosts(struct engine *e, struct cell *c,
- void engine_add_ghosts(struct engine *e, struct cell *c, struct task *ghost_in,
-                        struct task *ghost_out) {
- 
--  /* Abort as there are no hydro particles here? */
--  if (c->hydro.count_total == 0) return;
--
-   /* If we have reached the leaf OR have to few particles to play with*/
--  if (!c->split || c->hydro.count_total < engine_max_parts_per_ghost) {
-+  if (!c->split || c->hydro.count < engine_max_parts_per_ghost) {
- 
-     /* Add the ghost task and its dependencies */
-     struct scheduler *s = &e->sched;
-@@ -700,7 +797,6 @@ void engine_add_ghosts(struct engine *e, struct cell *c, struct task *ghost_in,
-         scheduler_addtask(s, task_type_ghost, task_subtype_none, 0, 0, c, NULL);
-     scheduler_addunlock(s, ghost_in, c->hydro.ghost);
-     scheduler_addunlock(s, c->hydro.ghost, ghost_out);
--
-   } else {
-     /* Keep recursing */
-     for (int k = 0; k < 8; k++)
-@@ -783,13 +879,16 @@ void engine_make_hierarchical_tasks_stars(struct engine *e, struct cell *c) {
- 
-   /* Are we in a super-cell ? */
-   if (c->super == c) {
--
--    /* Add the sort task. */
--    c->stars.sorts = scheduler_addtask(s, task_type_stars_sort,
--                                       task_subtype_none, 0, 0, c, NULL);
-+    /* Foreign tasks only */
-+    if (c->nodeID != e->nodeID) {
-+      c->stars.sorts_foreign = scheduler_addtask(
-+          s, task_type_stars_sort_foreign, task_subtype_none, 0, 0, c, NULL);
-+    }
- 
-     /* Local tasks only... */
-     if (c->nodeID == e->nodeID) {
-+      c->stars.sorts_local = scheduler_addtask(
-+          s, task_type_stars_sort_local, task_subtype_none, 0, 0, c, NULL);
- 
-       /* Generate the ghost tasks. */
-       c->stars.ghost_in =
-@@ -1054,11 +1153,19 @@ void engine_count_and_link_tasks_mapper(void *map_data, int num_elements,
-     }
- 
-     /* Link stars sort tasks to all the higher sort task. */
--    if (t_type == task_type_stars_sort) {
-+    if (t_type == task_type_stars_sort_local) {
-       for (struct cell *finger = t->ci->parent; finger != NULL;
--           finger = finger->parent)
--        if (finger->stars.sorts != NULL)
--          scheduler_addunlock(sched, t, finger->stars.sorts);
-+           finger = finger->parent) {
-+        if (finger->stars.sorts_local != NULL)
-+          scheduler_addunlock(sched, t, finger->stars.sorts_local);
-+      }
-+    }
-+    if (t_type == task_type_stars_sort_foreign) {
-+      for (struct cell *finger = t->ci->parent; finger != NULL;
-+           finger = finger->parent) {
-+        if (finger->stars.sorts_foreign != NULL)
-+          scheduler_addunlock(sched, t, finger->stars.sorts_foreign);
-+      }
-     }
- 
-     /* Link self tasks to cells. */
-@@ -1332,8 +1439,7 @@ void engine_link_gravity_tasks(struct engine *e) {
-  */
- static inline void engine_make_hydro_loops_dependencies(
-     struct scheduler *sched, struct task *density, struct task *gradient,
--    struct task *force, struct task *limiter, struct cell *c, int with_cooling,
--    int with_limiter) {
-+    struct task *force, struct cell *c, int with_cooling) {
- 
-   /* density loop --> ghost --> gradient loop --> extra_ghost */
-   /* extra_ghost --> force loop  */
-@@ -1351,15 +1457,14 @@ static inline void engine_make_hydro_loops_dependencies(
-  * @param sched The #scheduler.
-  * @param density The density task to link.
-  * @param force The force task to link.
-- * @param limiter The limiter task to link.
-  * @param c The cell.
-- * @param with_cooling Are we running with cooling switched on?
-- * @param with_limiter Are we running with limiter switched on?
-+ * @param with_cooling Are we running with cooling switched on ?
-  */
--static inline void engine_make_hydro_loops_dependencies(
--    struct scheduler *sched, struct task *density, struct task *force,
--    struct task *limiter, struct cell *c, int with_cooling, int with_limiter) {
--
-+static inline void engine_make_hydro_loops_dependencies(struct scheduler *sched,
-+                                                        struct task *density,
-+                                                        struct task *force,
-+                                                        struct cell *c,
-+                                                        int with_cooling) {
-   /* density loop --> ghost --> force loop */
-   scheduler_addunlock(sched, density, c->hydro.super->hydro.ghost_in);
-   scheduler_addunlock(sched, c->hydro.super->hydro.ghost_out, force);
-@@ -1381,7 +1486,6 @@ static inline void engine_make_stars_loops_dependencies(struct scheduler *sched,
-   /* density loop --> ghost --> feedback loop*/
-   scheduler_addunlock(sched, density, c->super->stars.ghost_in);
-   scheduler_addunlock(sched, c->super->stars.ghost_out, feedback);
--  scheduler_addunlock(sched, c->super->hydro.ghost_out, density);
- }
- 
- /**
-@@ -1401,12 +1505,6 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
-   struct scheduler *sched = &e->sched;
-   const int nodeID = e->nodeID;
-   const int with_cooling = (e->policy & engine_policy_cooling);
--  const int with_limiter = (e->policy & engine_policy_limiter);
--#ifdef EXTRA_HYDRO_LOOP
--  struct task *t_gradient = NULL;
--#endif
--  struct task *t_force = NULL;
--  struct task *t_limiter = NULL;
- 
-   for (int ind = 0; ind < num_elements; ind++) {
-     struct task *t = &((struct task *)map_data)[ind];
-@@ -1424,53 +1522,31 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
- 
- #ifdef EXTRA_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, 0, 0, t->ci, NULL);
--      t_force = scheduler_addtask(sched, task_type_self, task_subtype_force, 0,
--                                  0, t->ci, NULL);
--
--      /* and the task for the time-step limiter */
--      if (with_limiter)
--        t_limiter = scheduler_addtask(sched, task_type_self,
--                                      task_subtype_limiter, 0, 0, t->ci, NULL);
-+      struct task *t2 = scheduler_addtask(
-+          sched, task_type_self, task_subtype_gradient, 0, 0, t->ci, NULL);
-+      struct task *t3 = scheduler_addtask(
-+          sched, task_type_self, task_subtype_force, 0, 0, t->ci, NULL);
- 
-       /* Add the link between the new loops and the cell */
--      engine_addlink(e, &t->ci->hydro.gradient, t_gradient);
--      engine_addlink(e, &t->ci->hydro.force, t_force);
--      if (with_limiter) engine_addlink(e, &t->ci->hydro.limiter, t_limiter);
-+      engine_addlink(e, &t->ci->hydro.gradient, t2);
-+      engine_addlink(e, &t->ci->hydro.force, t3);
- 
-       /* Now, build all the dependencies for the hydro */
--      engine_make_hydro_loops_dependencies(sched, t, t_gradient, t_force,
--                                           t_limiter, t->ci, with_cooling,
--                                           with_limiter);
--      scheduler_addunlock(sched, t_force, t->ci->super->end_force);
--      if (with_limiter)
--        scheduler_addunlock(sched, t->ci->super->kick2, t_limiter);
--      if (with_limiter)
--        scheduler_addunlock(sched, t_limiter, t->ci->super->timestep);
-+      engine_make_hydro_loops_dependencies(sched, t, t2, t3, t->ci,
-+                                           with_cooling);
-+      scheduler_addunlock(sched, t3, t->ci->super->end_force);
- #else
- 
-       /* Start by constructing the task for the second hydro loop */
--      t_force = scheduler_addtask(sched, task_type_self, task_subtype_force, 0,
--                                  0, t->ci, NULL);
--
--      /* and the task for the time-step limiter */
--      if (with_limiter)
--        t_limiter = scheduler_addtask(sched, task_type_self,
--                                      task_subtype_limiter, 0, 0, t->ci, NULL);
-+      struct task *t2 = scheduler_addtask(
-+          sched, task_type_self, task_subtype_force, 0, 0, t->ci, NULL);
- 
-       /* Add the link between the new loop and the cell */
--      engine_addlink(e, &t->ci->hydro.force, t_force);
--      if (with_limiter) engine_addlink(e, &t->ci->hydro.limiter, t_limiter);
-+      engine_addlink(e, &t->ci->hydro.force, t2);
- 
-       /* Now, build all the dependencies for the hydro */
--      engine_make_hydro_loops_dependencies(sched, t, t_force, t_limiter, t->ci,
--                                           with_cooling, with_limiter);
--      scheduler_addunlock(sched, t_force, t->ci->super->end_force);
--      if (with_limiter)
--        scheduler_addunlock(sched, t->ci->super->kick2, t_limiter);
--      if (with_limiter)
--        scheduler_addunlock(sched, t_limiter, t->ci->super->timestep);
-+      engine_make_hydro_loops_dependencies(sched, t, t2, t->ci, with_cooling);
-+      scheduler_addunlock(sched, t2, t->ci->super->end_force);
- #endif
-     }
- 
-@@ -1489,103 +1565,54 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
- 
- #ifdef EXTRA_HYDRO_LOOP
-       /* Start by constructing the task for the second and third hydro loop */
--      t_gradient = scheduler_addtask(sched, task_type_pair,
--                                     task_subtype_gradient, 0, 0, t->ci, t->cj);
--      t_force = scheduler_addtask(sched, task_type_pair, task_subtype_force, 0,
--                                  0, t->ci, t->cj);
--
--      /* and the task for the time-step limiter */
--      if (with_limiter)
--        t_limiter = scheduler_addtask(sched, task_type_pair,
--                                      task_subtype_limiter, 0, 0, t->ci, t->cj);
-+      struct task *t2 = scheduler_addtask(
-+          sched, task_type_pair, task_subtype_gradient, 0, 0, t->ci, t->cj);
-+      struct task *t3 = scheduler_addtask(
-+          sched, task_type_pair, task_subtype_force, 0, 0, t->ci, t->cj);
- 
-       /* Add the link between the new loop and both cells */
--      engine_addlink(e, &t->ci->hydro.gradient, t_gradient);
--      engine_addlink(e, &t->cj->hydro.gradient, t_gradient);
--      engine_addlink(e, &t->ci->hydro.force, t_force);
--      engine_addlink(e, &t->cj->hydro.force, t_force);
--      if (with_limiter) engine_addlink(e, &t->ci->hydro.limiter, t_limiter);
--      if (with_limiter) engine_addlink(e, &t->cj->hydro.limiter, t_limiter);
-+      engine_addlink(e, &t->ci->hydro.gradient, t2);
-+      engine_addlink(e, &t->cj->hydro.gradient, t2);
-+      engine_addlink(e, &t->ci->hydro.force, t3);
-+      engine_addlink(e, &t->cj->hydro.force, t3);
- 
-       /* 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 (t->ci->nodeID == nodeID) {
--        engine_make_hydro_loops_dependencies(sched, t, t_gradient, t_force,
--                                             t_limiter, t->ci, with_cooling,
--                                             with_limiter);
--        scheduler_addunlock(sched, t_force, t->ci->super->end_force);
--        if (with_limiter)
--          scheduler_addunlock(sched, t->ci->super->kick2, t_limiter);
--        if (with_limiter)
--          scheduler_addunlock(sched, t_limiter, t->ci->super->timestep);
--        if (with_limiter)
--          scheduler_addunlock(sched, t_limiter, t->ci->super->timestep_limiter);
-+        engine_make_hydro_loops_dependencies(sched, t, t2, t3, t->ci,
-+                                             with_cooling);
-+        scheduler_addunlock(sched, t3, t->ci->super->end_force);
-       }
-       if (t->cj->nodeID == nodeID) {
--        if (t->ci->hydro.super != t->cj->hydro.super) {
--          engine_make_hydro_loops_dependencies(sched, t, t_gradient, t_force,
--                                               t_limiter, t->cj, with_cooling,
--                                               with_limiter);
--        }
--
--        if (t->ci->super != t->cj->super) {
--          scheduler_addunlock(sched, t_force, t->cj->super->end_force);
--          if (with_limiter)
--            scheduler_addunlock(sched, t->cj->super->kick2, t_limiter);
--          if (with_limiter)
--            scheduler_addunlock(sched, t_limiter, t->cj->super->timestep);
--          if (with_limiter)
--            scheduler_addunlock(sched, t_limiter,
--                                t->cj->super->timestep_limiter);
--        }
-+        if (t->ci->hydro.super != t->cj->hydro.super)
-+          engine_make_hydro_loops_dependencies(sched, t, t2, t3, t->cj,
-+                                               with_cooling);
-+        if (t->ci->super != t->cj->super)
-+          scheduler_addunlock(sched, t3, t->cj->super->end_force);
-       }
- 
- #else
- 
-       /* Start by constructing the task for the second hydro loop */
--      t_force = scheduler_addtask(sched, task_type_pair, task_subtype_force, 0,
--                                  0, t->ci, t->cj);
--
--      /* and the task for the time-step limiter */
--      if (with_limiter)
--        t_limiter = scheduler_addtask(sched, task_type_pair,
--                                      task_subtype_limiter, 0, 0, t->ci, t->cj);
-+      struct task *t2 = scheduler_addtask(
-+          sched, task_type_pair, task_subtype_force, 0, 0, t->ci, t->cj);
- 
-       /* Add the link between the new loop and both cells */
--      engine_addlink(e, &t->ci->hydro.force, t_force);
--      engine_addlink(e, &t->cj->hydro.force, t_force);
--      if (with_limiter) engine_addlink(e, &t->ci->hydro.limiter, t_limiter);
--      if (with_limiter) engine_addlink(e, &t->cj->hydro.limiter, t_limiter);
-+      engine_addlink(e, &t->ci->hydro.force, t2);
-+      engine_addlink(e, &t->cj->hydro.force, t2);
- 
-       /* 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 (t->ci->nodeID == nodeID) {
--        engine_make_hydro_loops_dependencies(sched, t, t_force, t_limiter,
--                                             t->ci, with_cooling, with_limiter);
--        scheduler_addunlock(sched, t_force, t->ci->super->end_force);
--        if (with_limiter)
--          scheduler_addunlock(sched, t->ci->super->kick2, t_limiter);
--        if (with_limiter)
--          scheduler_addunlock(sched, t_limiter, t->ci->super->timestep);
--        if (with_limiter)
--          scheduler_addunlock(sched, t_limiter, t->ci->super->timestep_limiter);
-+        engine_make_hydro_loops_dependencies(sched, t, t2, t->ci, with_cooling);
-+        scheduler_addunlock(sched, t2, t->ci->super->end_force);
-       }
-       if (t->cj->nodeID == nodeID) {
--        if (t->ci->hydro.super != t->cj->hydro.super) {
--          engine_make_hydro_loops_dependencies(
--              sched, t, t_force, t_limiter, t->cj, with_cooling, with_limiter);
--        }
--
--        if (t->ci->super != t->cj->super) {
--          scheduler_addunlock(sched, t_force, t->cj->super->end_force);
--          if (with_limiter)
--            scheduler_addunlock(sched, t->cj->super->kick2, t_limiter);
--          if (with_limiter)
--            scheduler_addunlock(sched, t_limiter, t->cj->super->timestep);
--          if (with_limiter)
--            scheduler_addunlock(sched, t_limiter,
--                                t->cj->super->timestep_limiter);
--        }
-+        if (t->ci->hydro.super != t->cj->hydro.super)
-+          engine_make_hydro_loops_dependencies(sched, t, t2, t->cj,
-+                                               with_cooling);
-+        if (t->ci->super != t->cj->super)
-+          scheduler_addunlock(sched, t2, t->cj->super->end_force);
-       }
- 
- #endif
-@@ -1603,65 +1630,39 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
- #ifdef EXTRA_HYDRO_LOOP
- 
-       /* Start by constructing the task for the second and third hydro loop */
--      t_gradient =
-+      struct task *t2 =
-           scheduler_addtask(sched, task_type_sub_self, task_subtype_gradient,
--                            t->flags, 0, t->ci, NULL);
--      t_force = scheduler_addtask(sched, task_type_sub_self, task_subtype_force,
--                                  t->flags, 0, t->ci, NULL);
--
--      /* and the task for the time-step limiter */
--      if (with_limiter)
--        t_limiter =
--            scheduler_addtask(sched, task_type_sub_self, task_subtype_limiter,
--                              t->flags, 0, t->ci, NULL);
-+                            t->flags, 0, t->ci, t->cj);
-+      struct task *t3 =
-+          scheduler_addtask(sched, task_type_sub_self, task_subtype_force,
-+                            t->flags, 0, t->ci, t->cj);
- 
-       /* Add the link between the new loop and the cell */
--      engine_addlink(e, &t->ci->hydro.gradient, t_gradient);
--      engine_addlink(e, &t->ci->hydro.force, t_force);
--      if (with_limiter) engine_addlink(e, &t->ci->hydro.limiter, t_limiter);
-+      engine_addlink(e, &t->ci->hydro.gradient, t2);
-+      engine_addlink(e, &t->ci->hydro.force, t3);
- 
-       /* 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 (t->ci->nodeID == nodeID) {
--        engine_make_hydro_loops_dependencies(sched, t, t_gradient, t_force,
--                                             t_limiter, t->ci, with_cooling,
--                                             with_limiter);
--        scheduler_addunlock(sched, t_force, t->ci->super->end_force);
--        if (with_limiter)
--          scheduler_addunlock(sched, t->ci->super->kick2, t_limiter);
--        if (with_limiter)
--          scheduler_addunlock(sched, t_limiter, t->ci->super->timestep);
--        if (with_limiter)
--          scheduler_addunlock(sched, t_limiter, t->ci->super->timestep_limiter);
-+        engine_make_hydro_loops_dependencies(sched, t, t2, t3, t->ci,
-+                                             with_cooling);
-+        scheduler_addunlock(sched, t3, t->ci->super->end_force);
-       }
- 
- #else
-       /* Start by constructing the task for the second hydro loop */
--      t_force = scheduler_addtask(sched, task_type_sub_self, task_subtype_force,
--                                  t->flags, 0, t->ci, NULL);
--
--      /* and the task for the time-step limiter */
--      if (with_limiter)
--        t_limiter =
--            scheduler_addtask(sched, task_type_sub_self, task_subtype_limiter,
--                              t->flags, 0, t->ci, NULL);
-+      struct task *t2 =
-+          scheduler_addtask(sched, task_type_sub_self, task_subtype_force,
-+                            t->flags, 0, t->ci, t->cj);
- 
-       /* Add the link between the new loop and the cell */
--      engine_addlink(e, &t->ci->hydro.force, t_force);
--      if (with_limiter) engine_addlink(e, &t->ci->hydro.limiter, t_limiter);
-+      engine_addlink(e, &t->ci->hydro.force, t2);
- 
-       /* 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 (t->ci->nodeID == nodeID) {
--        engine_make_hydro_loops_dependencies(sched, t, t_force, t_limiter,
--                                             t->ci, with_cooling, with_limiter);
--        scheduler_addunlock(sched, t_force, t->ci->super->end_force);
--        if (with_limiter)
--          scheduler_addunlock(sched, t->ci->super->kick2, t_limiter);
--        if (with_limiter)
--          scheduler_addunlock(sched, t_limiter, t->ci->super->timestep);
--        if (with_limiter)
--          scheduler_addunlock(sched, t_limiter, t->ci->super->timestep_limiter);
-+        engine_make_hydro_loops_dependencies(sched, t, t2, t->ci, with_cooling);
-+        scheduler_addunlock(sched, t2, t->ci->super->end_force);
-       }
- #endif
-     }
-@@ -1683,106 +1684,56 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
- #ifdef EXTRA_HYDRO_LOOP
- 
-       /* Start by constructing the task for the second and third hydro loop */
--      t_gradient =
-+      struct task *t2 =
-           scheduler_addtask(sched, task_type_sub_pair, task_subtype_gradient,
-                             t->flags, 0, t->ci, t->cj);
--      t_force = scheduler_addtask(sched, task_type_sub_pair, task_subtype_force,
--                                  t->flags, 0, t->ci, t->cj);
--
--      /* and the task for the time-step limiter */
--      if (with_limiter)
--        t_limiter =
--            scheduler_addtask(sched, task_type_sub_pair, task_subtype_limiter,
--                              t->flags, 0, t->ci, t->cj);
-+      struct task *t3 =
-+          scheduler_addtask(sched, task_type_sub_pair, task_subtype_force,
-+                            t->flags, 0, t->ci, t->cj);
- 
-       /* Add the link between the new loop and both cells */
--      engine_addlink(e, &t->ci->hydro.gradient, t_gradient);
--      engine_addlink(e, &t->cj->hydro.gradient, t_gradient);
--      engine_addlink(e, &t->ci->hydro.force, t_force);
--      engine_addlink(e, &t->cj->hydro.force, t_force);
--      if (with_limiter) engine_addlink(e, &t->ci->hydro.limiter, t_limiter);
--      if (with_limiter) engine_addlink(e, &t->cj->hydro.limiter, t_limiter);
-+      engine_addlink(e, &t->ci->hydro.gradient, t2);
-+      engine_addlink(e, &t->cj->hydro.gradient, t2);
-+      engine_addlink(e, &t->ci->hydro.force, t3);
-+      engine_addlink(e, &t->cj->hydro.force, t3);
- 
-       /* 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 (t->ci->nodeID == nodeID) {
--        engine_make_hydro_loops_dependencies(sched, t, t_gradient, t_force,
--                                             t_limiter, t->ci, with_cooling,
--                                             with_limiter);
--        scheduler_addunlock(sched, t_force, t->ci->super->end_force);
--        if (with_limiter)
--          scheduler_addunlock(sched, t->ci->super->kick2, t_limiter);
--        if (with_limiter)
--          scheduler_addunlock(sched, t_limiter, t->ci->super->timestep);
--        if (with_limiter)
--          scheduler_addunlock(sched, t_limiter, t->ci->super->timestep_limiter);
-+        engine_make_hydro_loops_dependencies(sched, t, t2, t3, t->ci,
-+                                             with_cooling);
-+        scheduler_addunlock(sched, t3, t->ci->super->end_force);
-       }
-       if (t->cj->nodeID == nodeID) {
--        if (t->ci->hydro.super != t->cj->hydro.super) {
--          engine_make_hydro_loops_dependencies(sched, t, t_gradient, t_force,
--                                               t_limiter, t->cj, with_cooling,
--                                               with_limiter);
--        }
--
--        if (t->ci->super != t->cj->super) {
--          scheduler_addunlock(sched, t_force, t->cj->super->end_force);
--          if (with_limiter)
--            scheduler_addunlock(sched, t->cj->super->kick2, t_limiter);
--          if (with_limiter)
--            scheduler_addunlock(sched, t_limiter, t->cj->super->timestep);
--          if (with_limiter)
--            scheduler_addunlock(sched, t_limiter,
--                                t->cj->super->timestep_limiter);
--        }
-+        if (t->ci->hydro.super != t->cj->hydro.super)
-+          engine_make_hydro_loops_dependencies(sched, t, t2, t3, t->cj,
-+                                               with_cooling);
-+        if (t->ci->super != t->cj->super)
-+          scheduler_addunlock(sched, t3, t->cj->super->end_force);
-       }
- 
- #else
-       /* Start by constructing the task for the second hydro loop */
--      t_force = scheduler_addtask(sched, task_type_sub_pair, task_subtype_force,
--                                  t->flags, 0, t->ci, t->cj);
--
--      /* and the task for the time-step limiter */
--      if (with_limiter)
--        t_limiter =
--            scheduler_addtask(sched, task_type_sub_pair, task_subtype_limiter,
--                              t->flags, 0, t->ci, t->cj);
-+      struct task *t2 =
-+          scheduler_addtask(sched, task_type_sub_pair, task_subtype_force,
-+                            t->flags, 0, t->ci, t->cj);
- 
-       /* Add the link between the new loop and both cells */
--      engine_addlink(e, &t->ci->hydro.force, t_force);
--      engine_addlink(e, &t->cj->hydro.force, t_force);
--      if (with_limiter) engine_addlink(e, &t->ci->hydro.limiter, t_limiter);
--      if (with_limiter) engine_addlink(e, &t->cj->hydro.limiter, t_limiter);
-+      engine_addlink(e, &t->ci->hydro.force, t2);
-+      engine_addlink(e, &t->cj->hydro.force, t2);
- 
-       /* 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 (t->ci->nodeID == nodeID) {
--        engine_make_hydro_loops_dependencies(sched, t, t_force, t_limiter,
--                                             t->ci, with_cooling, with_limiter);
--
--        scheduler_addunlock(sched, t_force, t->ci->super->end_force);
--        if (with_limiter)
--          scheduler_addunlock(sched, t->ci->super->kick2, t_limiter);
--        if (with_limiter)
--          scheduler_addunlock(sched, t_limiter, t->ci->super->timestep);
--        if (with_limiter)
--          scheduler_addunlock(sched, t_limiter, t->ci->super->timestep_limiter);
-+        engine_make_hydro_loops_dependencies(sched, t, t2, t->ci, with_cooling);
-+        scheduler_addunlock(sched, t2, t->ci->super->end_force);
-       }
-       if (t->cj->nodeID == nodeID) {
--        if (t->ci->hydro.super != t->cj->hydro.super) {
--          engine_make_hydro_loops_dependencies(
--              sched, t, t_force, t_limiter, t->cj, with_cooling, with_limiter);
--        }
--
--        if (t->ci->super != t->cj->super) {
--          scheduler_addunlock(sched, t_force, t->cj->super->end_force);
--          if (with_limiter)
--            scheduler_addunlock(sched, t->cj->super->kick2, t_limiter);
--          if (with_limiter)
--            scheduler_addunlock(sched, t_limiter, t->cj->super->timestep);
--          if (with_limiter)
--            scheduler_addunlock(sched, t_limiter,
--                                t->cj->super->timestep_limiter);
--        }
-+        if (t->ci->hydro.super != t->cj->hydro.super)
-+          engine_make_hydro_loops_dependencies(sched, t, t2, t->cj,
-+                                               with_cooling);
-+        if (t->ci->super != t->cj->super)
-+          scheduler_addunlock(sched, t2, t->cj->super->end_force);
-       }
- #endif
-     }
-@@ -1810,7 +1761,7 @@ void engine_make_extra_starsloop_tasks_mapper(void *map_data, int num_elements,
-     struct task *t = &((struct task *)map_data)[ind];
- 
-     /* Sort tasks depend on the drift and gravity drift of the cell. */
--    if (t->type == task_type_stars_sort && t->ci->nodeID == engine_rank) {
-+    if (t->type == task_type_stars_sort_local) {
-       scheduler_addunlock(sched, t->ci->hydro.super->hydro.drift, t);
-       scheduler_addunlock(sched, t->ci->super->grav.drift, t);
-     }
-@@ -1845,19 +1796,25 @@ void engine_make_extra_starsloop_tasks_mapper(void *map_data, int num_elements,
-       /* Make all stars density tasks depend on the hydro drift and sorts,
-        * gravity drift and star sorts. */
-       if (t->ci->nodeID == engine_rank)
--        scheduler_addunlock(sched, t->ci->super->hydro.drift, t);
--      scheduler_addunlock(sched, t->ci->super->hydro.sorts, t);
--      if (t->cj->nodeID == engine_rank)
--        scheduler_addunlock(sched, t->cj->super->grav.drift, t);
--      scheduler_addunlock(sched, t->ci->super->stars.sorts, t);
-+        scheduler_addunlock(sched, t->ci->hydro.super->hydro.drift, t);
-+      scheduler_addunlock(sched, t->ci->hydro.super->hydro.sorts, t);
- 
--      if (t->ci->super != t->cj->super) {
-+      if (t->ci->nodeID == engine_rank) {
-+        scheduler_addunlock(sched, t->ci->super->grav.drift, t);
-+        scheduler_addunlock(sched, t->ci->super->stars.sorts_local, t);
-+      }
-+
-+      if (t->ci->hydro.super != t->cj->hydro.super) {
-         if (t->cj->nodeID == engine_rank)
--          scheduler_addunlock(sched, t->cj->super->hydro.drift, t);
--        scheduler_addunlock(sched, t->cj->super->hydro.sorts, t);
--        if (t->ci->nodeID == engine_rank)
--          scheduler_addunlock(sched, t->ci->super->grav.drift, t);
--        scheduler_addunlock(sched, t->cj->super->stars.sorts, t);
-+          scheduler_addunlock(sched, t->cj->hydro.super->hydro.drift, t);
-+        scheduler_addunlock(sched, t->cj->hydro.super->hydro.sorts, t);
-+      }
-+
-+      if (t->ci->super != t->cj->super) {
-+        if (t->cj->nodeID == engine_rank) {
-+          scheduler_addunlock(sched, t->cj->super->grav.drift, t);
-+          scheduler_addunlock(sched, t->cj->super->stars.sorts_local, t);
-+        }
-       }
- 
-       /* Start by constructing the task for the second stars loop */
-@@ -1865,6 +1822,16 @@ void engine_make_extra_starsloop_tasks_mapper(void *map_data, int num_elements,
-           scheduler_addtask(sched, task_type_pair, task_subtype_stars_feedback,
-                             0, 0, t->ci, t->cj);
- 
-+      /* Add sort before feedback loop */
-+      if (t->ci->nodeID != engine_rank) {
-+        scheduler_addunlock(sched, t->ci->super->stars.sorts_foreign, t2);
-+      }
-+      if (t->ci->super != t->cj->super) {
-+        if (t->cj->nodeID != engine_rank) {
-+          scheduler_addunlock(sched, t->cj->super->stars.sorts_foreign, t2);
-+        }
-+      }
-+
-       /* Add the link between the new loop and both cells */
-       engine_addlink(e, &t->ci->stars.feedback, t2);
-       engine_addlink(e, &t->cj->stars.feedback, t2);
-@@ -1875,10 +1842,10 @@ void engine_make_extra_starsloop_tasks_mapper(void *map_data, int num_elements,
-         scheduler_addunlock(sched, t2, t->ci->super->end_force);
-       }
-       if (t->cj->nodeID == nodeID) {
--        if (t->ci->super != t->cj->super)
-+        if (t->ci->super != t->cj->super) {
-           engine_make_stars_loops_dependencies(sched, t, t2, t->cj);
--        if (t->ci->super != t->cj->super)
-           scheduler_addunlock(sched, t2, t->cj->super->end_force);
-+        }
-       }
-     }
- 
-@@ -1888,10 +1855,10 @@ void engine_make_extra_starsloop_tasks_mapper(void *map_data, int num_elements,
- 
-       /* Make all stars density tasks depend on the hydro drift and sorts,
-        * gravity drift and star sorts. */
--      scheduler_addunlock(sched, t->ci->super->hydro.drift, t);
--      scheduler_addunlock(sched, t->ci->super->hydro.sorts, t);
-+      scheduler_addunlock(sched, t->ci->hydro.super->hydro.drift, t);
-+      scheduler_addunlock(sched, t->ci->hydro.super->hydro.sorts, t);
-       scheduler_addunlock(sched, t->ci->super->grav.drift, t);
--      scheduler_addunlock(sched, t->ci->super->stars.sorts, t);
-+      scheduler_addunlock(sched, t->ci->super->stars.sorts_local, t);
- 
-       /* Start by constructing the task for the second stars loop */
-       struct task *t2 = scheduler_addtask(sched, task_type_sub_self,
-@@ -1915,19 +1882,25 @@ void engine_make_extra_starsloop_tasks_mapper(void *map_data, int num_elements,
-       /* Make all stars density tasks depend on the hydro drift and sorts,
-        * gravity drift and star sorts. */
-       if (t->cj->nodeID == engine_rank)
--        scheduler_addunlock(sched, t->cj->super->hydro.drift, t);
--      scheduler_addunlock(sched, t->cj->super->hydro.sorts, t);
--      if (t->cj->nodeID == engine_rank)
-+        scheduler_addunlock(sched, t->cj->hydro.super->hydro.drift, t);
-+      scheduler_addunlock(sched, t->cj->hydro.super->hydro.sorts, t);
-+
-+      if (t->cj->nodeID == engine_rank) {
-         scheduler_addunlock(sched, t->cj->super->grav.drift, t);
--      scheduler_addunlock(sched, t->ci->super->stars.sorts, t);
-+        scheduler_addunlock(sched, t->cj->super->stars.sorts_local, t);
-+      }
- 
--      if (t->ci->super != t->cj->super) {
--        if (t->ci->nodeID == engine_rank)
--          scheduler_addunlock(sched, t->ci->super->hydro.drift, t);
--        scheduler_addunlock(sched, t->ci->super->hydro.sorts, t);
-+      if (t->ci->hydro.super != t->cj->hydro.super) {
-         if (t->ci->nodeID == engine_rank)
-+          scheduler_addunlock(sched, t->ci->hydro.super->hydro.drift, t);
-+        scheduler_addunlock(sched, t->ci->hydro.super->hydro.sorts, t);
-+      }
-+
-+      if (t->ci->super != t->cj->super) {
-+        if (t->ci->nodeID == engine_rank) {
-           scheduler_addunlock(sched, t->ci->super->grav.drift, t);
--        scheduler_addunlock(sched, t->cj->super->stars.sorts, t);
-+          scheduler_addunlock(sched, t->ci->super->stars.sorts_local, t);
-+        }
-       }
- 
-       /* Start by constructing the task for the second stars loop */
-@@ -1935,6 +1908,16 @@ void engine_make_extra_starsloop_tasks_mapper(void *map_data, int num_elements,
-                                           task_subtype_stars_feedback, t->flags,
-                                           0, t->ci, t->cj);
- 
-+      /* Add the sort before feedback */
-+      if (t->cj->nodeID != engine_rank) {
-+        scheduler_addunlock(sched, t->cj->super->stars.sorts_foreign, t2);
-+      }
-+      if (t->ci->super != t->cj->super) {
-+        if (t->ci->nodeID != engine_rank) {
-+          scheduler_addunlock(sched, t->ci->super->stars.sorts_foreign, t2);
-+        }
-+      }
-+
-       /* Add the link between the new loop and both cells */
-       engine_addlink(e, &t->ci->stars.feedback, t2);
-       engine_addlink(e, &t->cj->stars.feedback, t2);
-@@ -2020,7 +2003,9 @@ void engine_make_starsloop_tasks_mapper(void *map_data, int num_elements,
-           struct cell *cj = &cells[cjd];
- 
-           /* Is that neighbour local and does it have particles ? */
--          if (cid >= cjd || (cj->stars.count == 0 && cj->hydro.count == 0) ||
-+          if (cid >= cjd ||
-+              ((cj->stars.count == 0 || ci->hydro.count == 0) &&
-+               (cj->hydro.count == 0 || ci->stars.count == 0)) ||
-               (ci->nodeID != nodeID && cj->nodeID != nodeID))
-             continue;
- 
-@@ -2028,6 +2013,50 @@ void engine_make_starsloop_tasks_mapper(void *map_data, int num_elements,
-           const int sid = sortlistID[(kk + 1) + 3 * ((jj + 1) + 3 * (ii + 1))];
-           scheduler_addtask(sched, task_type_pair, task_subtype_stars_density,
-                             sid, 0, ci, cj);
-+
-+#ifdef SWIFT_DEBUG_CHECKS
-+#ifdef WITH_MPI
-+
-+          /* Let's cross-check that we had a proxy for that cell */
-+          if (ci->nodeID == nodeID && cj->nodeID != engine_rank) {
-+
-+            /* Find the proxy for this node */
-+            const int proxy_id = e->proxy_ind[cj->nodeID];
-+            if (proxy_id < 0)
-+              error("No proxy exists for that foreign node %d!", cj->nodeID);
-+
-+            const struct proxy *p = &e->proxies[proxy_id];
-+
-+            /* Check whether the cell exists in the proxy */
-+            int n = 0;
-+            for (n = 0; n < p->nr_cells_in; n++)
-+              if (p->cells_in[n] == cj) break;
-+            if (n == p->nr_cells_in)
-+              error(
-+                  "Cell %d not found in the proxy but trying to construct "
-+                  "stars task!",
-+                  cjd);
-+          } else if (cj->nodeID == nodeID && ci->nodeID != engine_rank) {
-+
-+            /* Find the proxy for this node */
-+            const int proxy_id = e->proxy_ind[ci->nodeID];
-+            if (proxy_id < 0)
-+              error("No proxy exists for that foreign node %d!", ci->nodeID);
-+
-+            const struct proxy *p = &e->proxies[proxy_id];
-+
-+            /* Check whether the cell exists in the proxy */
-+            int n = 0;
-+            for (n = 0; n < p->nr_cells_in; n++)
-+              if (p->cells_in[n] == ci) break;
-+            if (n == p->nr_cells_in)
-+              error(
-+                  "Cell %d not found in the proxy but trying to construct "
-+                  "stars task!",
-+                  cid);
-+          }
-+#endif /* WITH_MPI */
-+#endif /* SWIFT_DEBUG_CHECKS */
-         }
-       }
-     }
-@@ -2168,7 +2197,6 @@ struct cell_type_pair {
- void engine_addtasks_send_mapper(void *map_data, int num_elements,
-                                  void *extra_data) {
-   struct engine *e = (struct engine *)extra_data;
--  const int with_limiter = (e->policy & engine_policy_limiter);
-   struct cell_type_pair *cell_type_pairs = (struct cell_type_pair *)map_data;
- 
-   for (int k = 0; k < num_elements; k++) {
-@@ -2177,7 +2205,7 @@ void engine_addtasks_send_mapper(void *map_data, int num_elements,
-     const int type = cell_type_pairs[k].type;
- 
-     /* Add the send task for the particle timesteps. */
--    engine_addtasks_send_timestep(e, ci, cj, NULL, NULL, with_limiter);
-+    engine_addtasks_send_timestep(e, ci, cj, NULL);
- 
-     /* Add the send tasks for the cells in the proxy that have a hydro
-      * connection. */
-@@ -2185,6 +2213,12 @@ void engine_addtasks_send_mapper(void *map_data, int num_elements,
-       engine_addtasks_send_hydro(e, ci, cj, /*t_xv=*/NULL,
-                                  /*t_rho=*/NULL, /*t_gradient=*/NULL);
- 
-+    /* Add the send tasks for the cells in the proxy that have a stars
-+     * connection. */
-+    if ((e->policy & engine_policy_feedback) && (type & proxy_cell_type_hydro))
-+      engine_addtasks_send_stars(e, ci, cj, /*t_xv=*/NULL,
-+                                 /*t_rho=*/NULL);
-+
-     /* Add the send tasks for the cells in the proxy that have a gravity
-      * connection. */
-     if ((e->policy & engine_policy_self_gravity) &&
-@@ -2196,7 +2230,6 @@ void engine_addtasks_send_mapper(void *map_data, int num_elements,
- void engine_addtasks_recv_mapper(void *map_data, int num_elements,
-                                  void *extra_data) {
-   struct engine *e = (struct engine *)extra_data;
--  const int with_limiter = (e->policy & engine_policy_limiter);
-   struct cell_type_pair *cell_type_pairs = (struct cell_type_pair *)map_data;
- 
-   for (int k = 0; k < num_elements; k++) {
-@@ -2204,13 +2237,18 @@ void engine_addtasks_recv_mapper(void *map_data, int num_elements,
-     const int type = cell_type_pairs[k].type;
- 
-     /* Add the recv task for the particle timesteps. */
--    engine_addtasks_recv_timestep(e, ci, NULL, NULL, with_limiter);
-+    engine_addtasks_recv_timestep(e, ci, NULL);
- 
-     /* Add the recv tasks for the cells in the proxy that have a hydro
-      * connection. */
-     if ((e->policy & engine_policy_hydro) && (type & proxy_cell_type_hydro))
-       engine_addtasks_recv_hydro(e, ci, NULL, NULL, NULL);
- 
-+    /* Add the recv tasks for the cells in the proxy that have a stars
-+     * connection. */
-+    if ((e->policy & engine_policy_feedback) && (type & proxy_cell_type_hydro))
-+      engine_addtasks_recv_stars(e, ci, NULL, NULL);
-+
-     /* Add the recv tasks for the cells in the proxy that have a gravity
-      * connection. */
-     if ((e->policy & engine_policy_self_gravity) &&
-@@ -2294,7 +2332,6 @@ void engine_maketasks(struct engine *e) {
-   const size_t self_grav_tasks_per_cell = 125;
-   const size_t ext_grav_tasks_per_cell = 1;
-   const size_t stars_tasks_per_cell = 27;
--  const size_t limiter_tasks_per_cell = 27;
- 
-   if (e->policy & engine_policy_hydro)
-     e->size_links += s->tot_cells * hydro_tasks_per_cell;
-@@ -2304,8 +2341,6 @@ void engine_maketasks(struct engine *e) {
-     e->size_links += s->tot_cells * self_grav_tasks_per_cell;
-   if (e->policy & engine_policy_stars)
-     e->size_links += s->tot_cells * stars_tasks_per_cell;
--  if (e->policy & engine_policy_limiter)
--    e->size_links += s->tot_cells * limiter_tasks_per_cell;
- 
-   /* Allocate the new link list */
-   if ((e->links = (struct link *)malloc(sizeof(struct link) * e->size_links)) ==
-@@ -2402,9 +2437,6 @@ void engine_maketasks(struct engine *e) {
-   tic2 = getticks();
- 
- #ifdef WITH_MPI
--  if (e->policy & engine_policy_feedback)
--    error("Cannot run stellar feedback with MPI (yet).");
--
-   /* Add the communication tasks if MPI is being used. */
-   if (e->policy & engine_policy_mpi) {
- 
-diff --git a/src/engine_marktasks.c b/src/engine_marktasks.c
-index 3a26dbb..55dddba 100644
---- a/src/engine_marktasks.c
-+++ b/src/engine_marktasks.c
-@@ -53,6 +53,94 @@
- #include "proxy.h"
- #include "timers.h"
- 
-+#ifdef WITH_MPI
-+/**
-+ * @brief Activate the MPI for the stars
-+ */
-+void engine_activate_stars_mpi(struct engine *e, struct scheduler *s,
-+                               struct cell *ci, struct cell *cj) {
-+
-+  const int nodeID = e->nodeID;
-+  const int ci_nodeID = ci->nodeID;
-+  const int cj_nodeID = cj->nodeID;
-+  const int ci_active_stars = cell_is_active_stars(ci, e) &&
-+                              ci->stars.count != 0 && cj->hydro.count != 0;
-+  const int cj_active_stars = cell_is_active_stars(cj, e) &&
-+                              cj->stars.count != 0 && ci->hydro.count != 0;
-+
-+  /* Activate the send/recv tasks. */
-+  if (ci_nodeID != nodeID) {
-+
-+    if (cj_active_stars) {
-+      scheduler_activate(s, ci->mpi.hydro.recv_xv);
-+      if (ci_active_stars) {
-+        scheduler_activate(s, ci->mpi.stars.recv);
-+      }
-+    }
-+
-+    /* If the foreign cell is active, we want its ti_end values. */
-+    if (ci_active_stars) scheduler_activate(s, ci->mpi.recv_ti);
-+
-+    /* Is the foreign cell active and will need stuff from us? */
-+    if (ci_active_stars) {
-+
-+      struct link *l =
-+          scheduler_activate_send(s, cj->mpi.hydro.send_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_stars) {
-+        scheduler_activate_send(s, cj->mpi.stars.send, ci_nodeID);
-+      }
-+    }
-+
-+    /* If the local cell is active, send its ti_end values. */
-+    if (cj_active_stars) {
-+      scheduler_activate_send(s, cj->mpi.send_ti, ci_nodeID);
-+    }
-+
-+  } else if (cj_nodeID != nodeID) {
-+    /* If the local cell is active, receive data from the foreign cell. */
-+    if (ci_active_stars) {
-+
-+      scheduler_activate(s, cj->mpi.hydro.recv_xv);
-+      if (cj_active_stars) {
-+        scheduler_activate(s, cj->mpi.stars.recv);
-+      }
-+    }
-+
-+    /* If the foreign cell is active, we want its ti_end values. */
-+    if (cj_active_stars) scheduler_activate(s, cj->mpi.recv_ti);
-+
-+    /* Is the foreign cell active and will need stuff from us? */
-+    if (cj_active_stars) {
-+
-+      struct link *l =
-+          scheduler_activate_send(s, ci->mpi.hydro.send_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_stars) {
-+        scheduler_activate_send(s, ci->mpi.stars.send, cj_nodeID);
-+      }
-+    }
-+
-+    /* If the local cell is active, send its ti_end values. */
-+    if (ci_active_stars) {
-+      scheduler_activate_send(s, ci->mpi.send_ti, cj_nodeID);
-+    }
-+  }
-+}
-+#endif
-+
- /**
-  * @brief Mark tasks to be un-skipped and set the sort flags accordingly.
-  *        Threadpool mapper function.
-@@ -69,7 +157,6 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
-   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_limiter = e->policy & engine_policy_limiter;
- 
-   for (int ind = 0; ind < num_elements; ind++) {
- 
-@@ -91,7 +178,6 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
-         if (cell_is_active_hydro(ci, e)) {
-           scheduler_activate(s, t);
-           cell_activate_drift_part(ci, s);
--          if (with_limiter) cell_activate_limiter(ci, s);
-         }
-       }
- 
-@@ -101,7 +187,6 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
-         if (cell_is_active_hydro(ci, e)) {
-           scheduler_activate(s, t);
-           cell_activate_subcell_hydro_tasks(ci, NULL, s);
--          if (with_limiter) cell_activate_limiter(ci, s);
-         }
-       }
- 
-@@ -114,16 +199,6 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
-         if (cell_is_active_hydro(ci, e)) scheduler_activate(s, t);
-       }
- 
--      else if (t->type == task_type_self &&
--               t->subtype == task_subtype_limiter) {
--        if (cell_is_active_hydro(ci, e)) scheduler_activate(s, t);
--      }
--
--      else if (t->type == task_type_sub_self &&
--               t->subtype == task_subtype_limiter) {
--        if (cell_is_active_hydro(ci, e)) scheduler_activate(s, t);
--      }
--
- #ifdef EXTRA_HYDRO_LOOP
-       else if (t_type == task_type_self && t_subtype == task_subtype_gradient) {
-         if (cell_is_active_hydro(ci, e)) scheduler_activate(s, t);
-@@ -220,7 +295,6 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
-       /* 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))) {
-@@ -240,10 +314,6 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
-           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_limiter) cell_activate_limiter(ci, s);
--          if (cj_nodeID == nodeID && with_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);
-@@ -273,39 +343,43 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
-         if (t_type == task_type_pair) {
- 
-           /* Do ci */
--          /* Store some values. */
--          atomic_or(&cj->hydro.requires_sorts, 1 << t->flags);
--          atomic_or(&ci->stars.requires_sorts, 1 << t->flags);
-+          if (ci_active_stars && ci->nodeID == engine_rank) {
-+            /* Store some values. */
-+            atomic_or(&cj->hydro.requires_sorts, 1 << t->flags);
-+            atomic_or(&ci->stars.requires_sorts, 1 << t->flags);
- 
--          cj->hydro.dx_max_sort_old = cj->hydro.dx_max_sort;
--          ci->stars.dx_max_sort_old = ci->stars.dx_max_sort;
-+            cj->hydro.dx_max_sort_old = cj->hydro.dx_max_sort;
-+            ci->stars.dx_max_sort_old = ci->stars.dx_max_sort;
- 
--          /* Activate the hydro drift tasks. */
--          if (ci_nodeID == nodeID) cell_activate_drift_spart(ci, s);
-+            /* Activate the hydro 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) cell_activate_drift_part(cj, s);
- 
--          /* Check the sorts and activate them if needed. */
--          cell_activate_hydro_sorts(cj, t->flags, 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);
-+            cell_activate_stars_sorts(ci, t->flags, s);
-+          }
- 
-           /* Do cj */
--          /* Store some values. */
--          atomic_or(&ci->hydro.requires_sorts, 1 << t->flags);
--          atomic_or(&cj->stars.requires_sorts, 1 << t->flags);
-+          if (ci_active_stars && ci->nodeID == engine_rank) {
-+            /* Store some values. */
-+            atomic_or(&ci->hydro.requires_sorts, 1 << t->flags);
-+            atomic_or(&cj->stars.requires_sorts, 1 << t->flags);
- 
--          ci->hydro.dx_max_sort_old = ci->hydro.dx_max_sort;
--          cj->stars.dx_max_sort_old = cj->stars.dx_max_sort;
-+            ci->hydro.dx_max_sort_old = ci->hydro.dx_max_sort;
-+            cj->stars.dx_max_sort_old = cj->stars.dx_max_sort;
- 
--          /* Activate the hydro drift tasks. */
--          if (ci_nodeID == nodeID) cell_activate_drift_part(ci, s);
-+            /* Activate the hydro drift tasks. */
-+            if (ci_nodeID == nodeID) cell_activate_drift_part(ci, s);
- 
--          if (cj_nodeID == nodeID) cell_activate_drift_spart(cj, s);
-+            if (cj_nodeID == nodeID) cell_activate_drift_spart(cj, 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);
-+            /* 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. */
-@@ -368,7 +442,6 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
- 
-           /* Is the foreign cell active and will need stuff from us? */
-           if (ci_active_hydro) {
--
-             struct link *l =
-                 scheduler_activate_send(s, cj->mpi.hydro.send_xv, ci_nodeID);
- 
-@@ -396,6 +469,7 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
- 
-           /* If the local cell is active, receive data from the foreign cell. */
-           if (ci_active_hydro) {
-+
-             scheduler_activate(s, cj->mpi.hydro.recv_xv);
-             if (cj_active_hydro) {
-               scheduler_activate(s, cj->mpi.hydro.recv_rho);
-@@ -443,8 +517,11 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
- 
-         /* Too much particle movement? */
-         if (cell_need_rebuild_for_stars_pair(ci, cj)) *rebuild_space = 1;
-+        if (cell_need_rebuild_for_hydro_pair(ci, cj)) *rebuild_space = 1;
- 
--        // LOIC: Need implementing MPI case
-+#ifdef WITH_MPI
-+        engine_activate_stars_mpi(e, s, ci, cj);
-+#endif
-       }
- 
-       /* Only interested in gravity tasks as of here. */
-diff --git a/src/gravity_properties.c b/src/gravity_properties.c
-index e548e30..fffbf22 100644
---- a/src/gravity_properties.c
-+++ b/src/gravity_properties.c
-@@ -170,22 +170,20 @@ void gravity_props_print_snapshot(hid_t h_grpgrav,
-   io_write_attribute_s(h_grpgrav, "Softening style",
-                        kernel_gravity_softening_name);
-   io_write_attribute_f(
--      h_grpgrav, "Comoving softening length [internal units]",
-+      h_grpgrav, "Comoving softening length",
-       p->epsilon_comoving * kernel_gravity_softening_plummer_equivalent);
-+  io_write_attribute_f(h_grpgrav,
-+                       "Comoving Softening length (Plummer equivalent)",
-+                       p->epsilon_comoving);
-   io_write_attribute_f(
--      h_grpgrav,
--      "Comoving Softening length (Plummer equivalent)  [internal units]",
--      p->epsilon_comoving);
--  io_write_attribute_f(
--      h_grpgrav, "Maximal physical softening length  [internal units]",
-+      h_grpgrav, "Maximal physical softening length",
-       p->epsilon_max_physical * kernel_gravity_softening_plummer_equivalent);
-   io_write_attribute_f(h_grpgrav,
--                       "Maximal physical softening length (Plummer equivalent) "
--                       " [internal units]",
-+                       "Maximal physical softening length (Plummer equivalent)",
-                        p->epsilon_max_physical);
-   io_write_attribute_f(h_grpgrav, "Opening angle", p->theta_crit);
-   io_write_attribute_s(h_grpgrav, "Scheme", GRAVITY_IMPLEMENTATION);
--  io_write_attribute_i(h_grpgrav, "MM order", SELF_GRAVITY_MULTIPOLE_ORDER);
-+  io_write_attribute_d(h_grpgrav, "MM order", SELF_GRAVITY_MULTIPOLE_ORDER);
-   io_write_attribute_f(h_grpgrav, "Mesh a_smooth", p->a_smooth);
-   io_write_attribute_f(h_grpgrav, "Mesh r_cut_max ratio", p->r_cut_max_ratio);
-   io_write_attribute_f(h_grpgrav, "Mesh r_cut_min ratio", p->r_cut_min_ratio);
-diff --git a/src/hydro_properties.c b/src/hydro_properties.c
-index ac226c9..85f88d4 100644
---- a/src/hydro_properties.c
-+++ b/src/hydro_properties.c
-@@ -239,8 +239,7 @@ void hydro_props_print_snapshot(hid_t h_grpsph, const struct hydro_props *p) {
-   io_write_attribute_f(h_grpsph, "Kernel delta N_ngb", p->delta_neighbours);
-   io_write_attribute_f(h_grpsph, "Kernel eta", p->eta_neighbours);
-   io_write_attribute_f(h_grpsph, "Smoothing length tolerance", p->h_tolerance);
--  io_write_attribute_f(h_grpsph, "Maximal smoothing length [internal units]",
--                       p->h_max);
-+  io_write_attribute_f(h_grpsph, "Maximal smoothing length", p->h_max);
-   io_write_attribute_f(h_grpsph, "CFL parameter", p->CFL_condition);
-   io_write_attribute_f(h_grpsph, "Volume log(max(delta h))",
-                        p->log_max_h_change);
-@@ -249,12 +248,8 @@ void hydro_props_print_snapshot(hid_t h_grpsph, const struct hydro_props *p) {
-   io_write_attribute_i(h_grpsph, "Max ghost iterations",
-                        p->max_smoothing_iterations);
-   io_write_attribute_f(h_grpsph, "Minimal temperature", p->minimal_temperature);
--  io_write_attribute_f(h_grpsph,
--                       "Minimal energy per unit mass [internal units]",
--                       p->minimal_internal_energy);
-   io_write_attribute_f(h_grpsph, "Initial temperature", p->initial_temperature);
--  io_write_attribute_f(h_grpsph,
--                       "Initial energy per unit mass [internal units]",
-+  io_write_attribute_f(h_grpsph, "Initial energy per unit mass",
-                        p->initial_internal_energy);
-   io_write_attribute_f(h_grpsph, "Hydrogen mass fraction",
-                        p->hydrogen_mass_fraction);
-@@ -265,11 +260,8 @@ void hydro_props_print_snapshot(hid_t h_grpsph, const struct hydro_props *p) {
-                        p->viscosity.alpha_max);
-   io_write_attribute_f(h_grpsph, "Alpha viscosity (min)",
-                        p->viscosity.alpha_min);
--  io_write_attribute_f(h_grpsph, "Viscosity decay length [internal units]",
--                       p->viscosity.length);
-+  io_write_attribute_f(h_grpsph, "Viscosity decay length", p->viscosity.length);
-   io_write_attribute_f(h_grpsph, "Beta viscosity", const_viscosity_beta);
--  io_write_attribute_f(h_grpsph, "Max v_sig ratio (limiter)",
--                       const_limiter_max_v_sig_ratio);
- }
- #endif
- 
-diff --git a/src/parallel_io.c b/src/parallel_io.c
-index 0c6514b..611e4af 100644
---- a/src/parallel_io.c
-+++ b/src/parallel_io.c
-@@ -957,6 +957,7 @@ void prepare_file(struct engine* e, const char* baseName, long long N_total[6],
-   const int with_cosmology = e->policy & engine_policy_cosmology;
- 
-   FILE* xmfFile = 0;
-+  int periodic = e->s->periodic;
-   int numFiles = 1;
- 
-   /* First time, we need to create the XMF file */
-@@ -983,25 +984,27 @@ void prepare_file(struct engine* e, const char* baseName, long long N_total[6],
-   xmf_write_outputheader(xmfFile, fileName, e->time);
- 
-   /* Open header to write simulation properties */
--  /* message("Writing file header..."); */
-+  /* message("Writing runtime parameters..."); */
-   hid_t h_grp =
--      H5Gcreate(h_file, "/Header", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
--  if (h_grp < 0) error("Error while creating file header\n");
-+      H5Gcreate(h_file, "/RuntimePars", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
-+  if (h_grp < 0) error("Error while creating runtime parameters group\n");
- 
--  /* Convert basic output information to snapshot units */
--  const double factor_time =
--      units_conversion_factor(internal_units, snapshot_units, UNIT_CONV_TIME);
--  const double factor_length =
--      units_conversion_factor(internal_units, snapshot_units, UNIT_CONV_LENGTH);
--  const double dblTime = e->time * factor_time;
--  const double dim[3] = {e->s->dim[0] * factor_length,
--                         e->s->dim[1] * factor_length,
--                         e->s->dim[2] * factor_length};
-+  /* Write the relevant information */
-+  io_write_attribute(h_grp, "PeriodicBoundariesOn", INT, &periodic, 1);
-+
-+  /* Close runtime parameters */
-+  H5Gclose(h_grp);
-+
-+  /* Open header to write simulation properties */
-+  /* message("Writing file header..."); */
-+  h_grp = H5Gcreate(h_file, "/Header", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
-+  if (h_grp < 0) error("Error while creating file header\n");
- 
-   /* Print the relevant information and print status */
--  io_write_attribute(h_grp, "BoxSize", DOUBLE, dim, 3);
-+  io_write_attribute(h_grp, "BoxSize", DOUBLE, e->s->dim, 3);
-+  double dblTime = e->time;
-   io_write_attribute(h_grp, "Time", DOUBLE, &dblTime, 1);
--  const int dimension = (int)hydro_dimension;
-+  int dimension = (int)hydro_dimension;
-   io_write_attribute(h_grp, "Dimension", INT, &dimension, 1);
-   io_write_attribute(h_grp, "Redshift", DOUBLE, &e->cosmology->z, 1);
-   io_write_attribute(h_grp, "Scale-factor", DOUBLE, &e->cosmology->a, 1);
-@@ -1279,32 +1282,6 @@ void write_output_parallel(struct engine* e, const char* baseName,
-     snprintf(fileName, FILENAME_BUFFER_SIZE, "%s_%04i.hdf5", baseName,
-              e->snapshot_output_count);
- 
--  /* Now write the top-level cell structure */
--  hid_t h_file_cells = 0, h_grp_cells = 0;
--  if (mpi_rank == 0) {
--
--    /* Open the snapshot on rank 0 */
--    h_file_cells = H5Fopen(fileName, H5F_ACC_RDWR, H5P_DEFAULT);
--    if (h_file_cells < 0)
--      error("Error while opening file '%s' on rank %d.", fileName, mpi_rank);
--
--    /* Create the group we want in the file */
--    h_grp_cells = H5Gcreate(h_file_cells, "/Cells", H5P_DEFAULT, H5P_DEFAULT,
--                            H5P_DEFAULT);
--    if (h_grp_cells < 0) error("Error while creating cells group");
--  }
--
--  /* Write the location of the particles in the arrays */
--  io_write_cell_offsets(h_grp_cells, e->s->cdim, e->s->cells_top,
--                        e->s->nr_cells, e->s->width, mpi_rank, N_total, offset,
--                        internal_units, snapshot_units);
--
--  /* Close everything */
--  if (mpi_rank == 0) {
--    H5Gclose(h_grp_cells);
--    H5Fclose(h_file_cells);
--  }
--
-   /* Prepare some file-access properties */
-   hid_t plist_id = H5Pcreate(H5P_FILE_ACCESS);
- 
-diff --git a/src/runner.c b/src/runner.c
-index 5552782..5b2e745 100644
---- a/src/runner.c
-+++ b/src/runner.c
-@@ -64,7 +64,6 @@
- #include "task.h"
- #include "timers.h"
- #include "timestep.h"
--#include "timestep_limiter.h"
- #include "tracers.h"
- 
- #define TASK_LOOP_DENSITY 0
-@@ -95,19 +94,14 @@
- #undef FUNCTION
- #undef FUNCTION_TASK_LOOP
- 
--/* Import the limiter loop functions. */
--#define FUNCTION limiter
--#define FUNCTION_TASK_LOOP TASK_LOOP_LIMITER
--#include "runner_doiact.h"
--#undef FUNCTION
--#undef FUNCTION_TASK_LOOP
--
- /* Import the gravity loop functions. */
- #include "runner_doiact_grav.h"
- 
- /* Import the stars density loop functions. */
- #define FUNCTION density
-+#define ONLY_LOCAL 1
- #include "runner_doiact_stars.h"
-+#undef ONLY_LOCAL
- #undef FUNCTION
- 
- /* Import the stars feedback loop functions. */
-@@ -128,14 +122,11 @@ void runner_do_stars_ghost(struct runner *r, struct cell *c, int timer) {
-   struct spart *restrict sparts = c->stars.parts;
-   const struct engine *e = r->e;
-   const struct cosmology *cosmo = e->cosmology;
--  const int with_cosmology = (e->policy & engine_policy_cosmology);
-   const struct stars_props *stars_properties = e->stars_properties;
-   const float stars_h_max = stars_properties->h_max;
-   const float eps = stars_properties->h_tolerance;
-   const float stars_eta_dim = pow_dimension(stars_properties->eta_neighbours);
-   const int max_smoothing_iter = stars_properties->max_smoothing_iterations;
--  const integertime_t ti_current = e->ti_current;
--  const double time_base = e->time_base;
-   int redo = 0, scount = 0;
- 
-   TIMER_TIC;
-@@ -162,7 +153,7 @@ void runner_do_stars_ghost(struct runner *r, struct cell *c, int timer) {
-     /* While there are particles that need to be updated... */
-     for (int num_reruns = 0; scount > 0 && num_reruns < max_smoothing_iter;
-          num_reruns++) {
--        
-+
-       /* Reset the redo-count. */
-       redo = 0;
- 
-@@ -172,18 +163,6 @@ void runner_do_stars_ghost(struct runner *r, struct cell *c, int timer) {
-         /* Get a direct pointer on the part. */
-         struct spart *sp = &sparts[sid[i]];
- 
--	/* get particle timestep */
--        double dt;
--        if (with_cosmology) {
--          const integertime_t ti_step = get_integer_timestep(sp->time_bin);
--          const integertime_t ti_begin =
--              get_integer_time_begin(ti_current - 1, sp->time_bin);
--          dt = cosmology_get_therm_kick_factor(e->cosmology, ti_begin,
--                                                     ti_begin + ti_step);
--        } else {
--          dt = get_timestep(sp->time_bin, time_base);
--        }
--
- #ifdef SWIFT_DEBUG_CHECKS
-         /* Is this part within the timestep? */
-         if (!spart_is_active(sp, e))
-@@ -217,6 +196,15 @@ void runner_do_stars_ghost(struct runner *r, struct cell *c, int timer) {
-               sp->density.wcount_dh * h_old_dim +
-               hydro_dimension * sp->density.wcount * h_old_dim_minus_one;
- 
-+          /* Skip if h is already h_max and we don't have enough neighbours */
-+          if ((sp->h >= stars_h_max) && (f < 0.f)) {
-+
-+            stars_reset_acceleration(sp);
-+            /* Ok, we are done with this particle */
-+            continue;
-+          }
-+
-+          /* Normal case: Use Newton-Raphson to get a better value of h */
-           /* Avoid floating point exception from f_prime = 0 */
-           h_new = h_old - f / (f_prime + FLT_MIN);
- #ifdef SWIFT_DEBUG_CHECKS
-@@ -262,9 +250,10 @@ void runner_do_stars_ghost(struct runner *r, struct cell *c, int timer) {
-         }
- 
-         /* We now have a particle whose smoothing length has converged */
-+        stars_reset_acceleration(sp);
- 
-         /* Compute the stellar evolution  */
--        stars_evolve_spart(sp, stars_properties, cosmo, dt);
-+        stars_evolve_spart(sp, stars_properties, cosmo);
-       }
- 
-       /* We now need to treat the particles whose smoothing length had not
-@@ -893,8 +882,9 @@ void runner_do_stars_sort(struct runner *r, struct cell *c, int flags,
-   if (flags == 0 && !c->stars.do_sub_sort) return;
- 
-   /* Check that the particles have been moved to the current time */
--  if (flags && !cell_are_spart_drifted(c, r->e))
-+  if (flags && !cell_are_spart_drifted(c, r->e)) {
-     error("Sorting un-drifted cell c->nodeID=%d", c->nodeID);
-+  }
- 
- #ifdef SWIFT_DEBUG_CHECKS
-   /* Make sure the sort flags are consistent (downward). */
-@@ -984,6 +974,9 @@ void runner_do_stars_sort(struct runner *r, struct cell *c, int flags,
- 
-         /* Copy the minimum into the new sort array. */
-         finger[ind].d = buff[inds[0]];
-+        if (c->progeny[inds[0]] == NULL ||
-+            c->progeny[inds[0]]->stars.count == 0)
-+          error("should not do it");
-         finger[ind].i = fingers[inds[0]]->i + off[inds[0]];
- 
-         /* Update the buffer. */
-@@ -1696,26 +1689,19 @@ void runner_do_kick1(struct runner *r, struct cell *c, int timer) {
-       /* If particle needs to be kicked */
-       if (part_is_starting(p, e)) {
- 
--#ifdef SWIFT_DEBUG_CHECKS
--        if (p->wakeup == time_bin_awake)
--          error("Woken-up particle that has not been processed in kick1");
--#endif
--
--        /* Skip particles that have been woken up and treated by the limiter. */
--        if (p->wakeup != time_bin_not_awake) continue;
--
-         const integertime_t ti_step = get_integer_timestep(p->time_bin);
-         const integertime_t ti_begin =
-             get_integer_time_begin(ti_current + 1, p->time_bin);
- 
- #ifdef SWIFT_DEBUG_CHECKS
--        const integertime_t ti_end = ti_begin + ti_step;
-+        const integertime_t ti_end =
-+            get_integer_time_end(ti_current + 1, p->time_bin);
- 
-         if (ti_begin != ti_current)
-           error(
-               "Particle in wrong time-bin, ti_end=%lld, ti_begin=%lld, "
--              "ti_step=%lld time_bin=%d wakeup=%d ti_current=%lld",
--              ti_end, ti_begin, ti_step, p->time_bin, p->wakeup, ti_current);
-+              "ti_step=%lld time_bin=%d ti_current=%lld",
-+              ti_end, ti_begin, ti_step, p->time_bin, ti_current);
- #endif
- 
-         /* Time interval for this half-kick */
-@@ -1877,60 +1863,39 @@ void runner_do_kick2(struct runner *r, struct cell *c, int timer) {
-       /* If particle needs to be kicked */
-       if (part_is_active(p, e)) {
- 
--        integertime_t ti_begin, ti_end, ti_step;
--
--#ifdef SWIFT_DEBUG_CHECKS
--        if (p->wakeup == time_bin_awake)
--          error("Woken-up particle that has not been processed in kick1");
--#endif
--
--        if (p->wakeup == time_bin_not_awake) {
--
--          /* Time-step from a regular kick */
--          ti_step = get_integer_timestep(p->time_bin);
--          ti_begin = get_integer_time_begin(ti_current, p->time_bin);
--          ti_end = ti_begin + ti_step;
--
--        } else {
--
--          /* Time-step that follows a wake-up call */
--          ti_begin = get_integer_time_begin(ti_current, p->wakeup);
--          ti_end = get_integer_time_end(ti_current, p->time_bin);
--          ti_step = ti_end - ti_begin;
--
--          /* Reset the flag. Everything is back to normal from now on. */
--          p->wakeup = time_bin_awake;
--        }
-+        const integertime_t ti_step = get_integer_timestep(p->time_bin);
-+        const integertime_t ti_begin =
-+            get_integer_time_begin(ti_current, p->time_bin);
- 
- #ifdef SWIFT_DEBUG_CHECKS
-         if (ti_begin + ti_step != ti_current)
-           error(
-               "Particle in wrong time-bin, ti_begin=%lld, ti_step=%lld "
--              "time_bin=%d wakeup=%d ti_current=%lld",
--              ti_begin, ti_step, p->time_bin, p->wakeup, ti_current);
-+              "time_bin=%d ti_current=%lld",
-+              ti_begin, ti_step, p->time_bin, ti_current);
- #endif
-         /* Time interval for this half-kick */
-         double dt_kick_grav, dt_kick_hydro, dt_kick_therm, dt_kick_corr;
-         if (with_cosmology) {
-           dt_kick_hydro = cosmology_get_hydro_kick_factor(
--              cosmo, ti_begin + ti_step / 2, ti_end);
-+              cosmo, ti_begin + ti_step / 2, ti_begin + ti_step);
-           dt_kick_grav = cosmology_get_grav_kick_factor(
--              cosmo, ti_begin + ti_step / 2, ti_end);
-+              cosmo, ti_begin + ti_step / 2, ti_begin + ti_step);
-           dt_kick_therm = cosmology_get_therm_kick_factor(
--              cosmo, ti_begin + ti_step / 2, ti_end);
-+              cosmo, ti_begin + ti_step / 2, ti_begin + ti_step);
-           dt_kick_corr = cosmology_get_corr_kick_factor(
--              cosmo, ti_begin + ti_step / 2, ti_end);
-+              cosmo, ti_begin + ti_step / 2, ti_begin + ti_step);
-         } else {
--          dt_kick_hydro = (ti_end - (ti_begin + ti_step / 2)) * time_base;
--          dt_kick_grav = (ti_end - (ti_begin + ti_step / 2)) * time_base;
--          dt_kick_therm = (ti_end - (ti_begin + ti_step / 2)) * time_base;
--          dt_kick_corr = (ti_end - (ti_begin + ti_step / 2)) * time_base;
-+          dt_kick_hydro = (ti_step / 2) * time_base;
-+          dt_kick_grav = (ti_step / 2) * time_base;
-+          dt_kick_therm = (ti_step / 2) * time_base;
-+          dt_kick_corr = (ti_step / 2) * time_base;
-         }
- 
-         /* Finish the time-step with a second half-kick */
-         kick_part(p, xp, dt_kick_hydro, dt_kick_grav, dt_kick_therm,
-                   dt_kick_corr, cosmo, hydro_props, ti_begin + ti_step / 2,
--                  ti_end);
-+                  ti_begin + ti_step);
- 
- #ifdef SWIFT_DEBUG_CHECKS
-         /* Check that kick and the drift are synchronized */
-@@ -2333,144 +2298,6 @@ void runner_do_timestep(struct runner *r, struct cell *c, int timer) {
- }
- 
- /**
-- * @brief Apply the time-step limiter to all awaken particles in a cell
-- * hierarchy.
-- *
-- * @param r The task #runner.
-- * @param c The #cell.
-- * @param force Limit the particles irrespective of the #cell flags.
-- * @param timer Are we timing this ?
-- */
--void runner_do_limiter(struct runner *r, struct cell *c, int force, int timer) {
--
--  const struct engine *e = r->e;
--  const integertime_t ti_current = e->ti_current;
--  const int count = c->hydro.count;
--  struct part *restrict parts = c->hydro.parts;
--  struct xpart *restrict xparts = c->hydro.xparts;
--
--  TIMER_TIC;
--
--#ifdef SWIFT_DEBUG_CHECKS
--  /* Check that we only limit local cells. */
--  if (c->nodeID != engine_rank) error("Limiting dt of a foreign cell is nope.");
--#endif
--
--  integertime_t ti_hydro_end_min = max_nr_timesteps, ti_hydro_end_max = 0,
--                ti_hydro_beg_max = 0;
--  integertime_t ti_gravity_end_min = max_nr_timesteps, ti_gravity_end_max = 0,
--                ti_gravity_beg_max = 0;
--
--  /* Limit irrespective of cell flags? */
--  force |= c->hydro.do_limiter;
--
--  /* Early abort? */
--  if (c->hydro.count == 0) {
--
--    /* Clear the limiter flags. */
--    c->hydro.do_limiter = 0;
--    c->hydro.do_sub_limiter = 0;
--    return;
--  }
--
--  /* Loop over the progeny ? */
--  if (c->split && (force || c->hydro.do_sub_limiter)) {
--    for (int k = 0; k < 8; k++) {
--      if (c->progeny[k] != NULL) {
--        struct cell *restrict cp = c->progeny[k];
--
--        /* Recurse */
--        runner_do_limiter(r, cp, force, 0);
--
--        /* And aggregate */
--        ti_hydro_end_min = min(cp->hydro.ti_end_min, ti_hydro_end_min);
--        ti_hydro_end_max = max(cp->hydro.ti_end_max, ti_hydro_end_max);
--        ti_hydro_beg_max = max(cp->hydro.ti_beg_max, ti_hydro_beg_max);
--        ti_gravity_end_min = min(cp->grav.ti_end_min, ti_gravity_end_min);
--        ti_gravity_end_max = max(cp->grav.ti_end_max, ti_gravity_end_max);
--        ti_gravity_beg_max = max(cp->grav.ti_beg_max, ti_gravity_beg_max);
--      }
--    }
--
--    /* Store the updated values */
--    c->hydro.ti_end_min = min(c->hydro.ti_end_min, ti_hydro_end_min);
--    c->hydro.ti_end_max = max(c->hydro.ti_end_max, ti_hydro_end_max);
--    c->hydro.ti_beg_max = max(c->hydro.ti_beg_max, ti_hydro_beg_max);
--    c->grav.ti_end_min = min(c->grav.ti_end_min, ti_gravity_end_min);
--    c->grav.ti_end_max = max(c->grav.ti_end_max, ti_gravity_end_max);
--    c->grav.ti_beg_max = max(c->grav.ti_beg_max, ti_gravity_beg_max);
--
--  } else if (!c->split && force) {
--
--    ti_hydro_end_min = c->hydro.ti_end_min;
--    ti_hydro_end_max = c->hydro.ti_end_max;
--    ti_hydro_beg_max = c->hydro.ti_beg_max;
--    ti_gravity_end_min = c->grav.ti_end_min;
--    ti_gravity_end_max = c->grav.ti_end_max;
--    ti_gravity_beg_max = c->grav.ti_beg_max;
--
--    /* Loop over the gas particles in this cell. */
--    for (int k = 0; k < count; k++) {
--
--      /* Get a handle on the part. */
--      struct part *restrict p = &parts[k];
--      struct xpart *restrict xp = &xparts[k];
--
--      /* Avoid inhibited particles */
--      if (part_is_inhibited(p, e)) continue;
--
--      /* If the particle will be active no need to wake it up */
--      if (part_is_active(p, e) && p->wakeup != time_bin_not_awake)
--        p->wakeup = time_bin_not_awake;
--
--      /* Bip, bip, bip... wake-up time */
--      if (p->wakeup == time_bin_awake) {
--
--        /* Apply the limiter and get the new time-step size */
--        const integertime_t ti_new_step = timestep_limit_part(p, xp, e);
--
--        /* What is the next sync-point ? */
--        ti_hydro_end_min = min(ti_current + ti_new_step, ti_hydro_end_min);
--        ti_hydro_end_max = max(ti_current + ti_new_step, ti_hydro_end_max);
--
--        /* What is the next starting point for this cell ? */
--        ti_hydro_beg_max = max(ti_current, ti_hydro_beg_max);
--
--        /* Also limit the gpart counter-part */
--        if (p->gpart != NULL) {
--
--          /* Register the time-bin */
--          p->gpart->time_bin = p->time_bin;
--
--          /* What is the next sync-point ? */
--          ti_gravity_end_min =
--              min(ti_current + ti_new_step, ti_gravity_end_min);
--          ti_gravity_end_max =
--              max(ti_current + ti_new_step, ti_gravity_end_max);
--
--          /* What is the next starting point for this cell ? */
--          ti_gravity_beg_max = max(ti_current, ti_gravity_beg_max);
--        }
--      }
--    }
--
--    /* Store the updated values */
--    c->hydro.ti_end_min = min(c->hydro.ti_end_min, ti_hydro_end_min);
--    c->hydro.ti_end_max = max(c->hydro.ti_end_max, ti_hydro_end_max);
--    c->hydro.ti_beg_max = max(c->hydro.ti_beg_max, ti_hydro_beg_max);
--    c->grav.ti_end_min = min(c->grav.ti_end_min, ti_gravity_end_min);
--    c->grav.ti_end_max = max(c->grav.ti_end_max, ti_gravity_end_max);
--    c->grav.ti_beg_max = max(c->grav.ti_beg_max, ti_gravity_beg_max);
--  }
--
--  /* Clear the limiter flags. */
--  c->hydro.do_limiter = 0;
--  c->hydro.do_sub_limiter = 0;
--
--  if (timer) TIMER_TOC(timer_do_limiter);
--}
--
--/**
-  * @brief End the force calculation of all active particles in a cell
-  * by multiplying the acccelerations by the relevant constants
-  *
-@@ -2779,9 +2606,11 @@ void runner_do_recv_gpart(struct runner *r, struct cell *c, int timer) {
-  *
-  * @param r The runner thread.
-  * @param c The cell.
-+ * @param clear_sorts Should we clear the sort flag and hence trigger a sort ?
-  * @param timer Are we timing this ?
-  */
--void runner_do_recv_spart(struct runner *r, struct cell *c, int timer) {
-+void runner_do_recv_spart(struct runner *r, struct cell *c, int clear_sorts,
-+                          int timer) {
- 
- #ifdef WITH_MPI
- 
-@@ -2789,19 +2618,22 @@ void runner_do_recv_spart(struct runner *r, struct cell *c, int timer) {
-   const size_t nr_sparts = c->stars.count;
-   const integertime_t ti_current = r->e->ti_current;
- 
--  error("Need to add h_max computation");
--
-   TIMER_TIC;
- 
-   integertime_t ti_gravity_end_min = max_nr_timesteps;
-   integertime_t ti_gravity_end_max = 0;
-+  integertime_t ti_stars_end_min = max_nr_timesteps;
-   timebin_t time_bin_min = num_time_bins;
-   timebin_t time_bin_max = 0;
-+  float h_max = 0.f;
- 
- #ifdef SWIFT_DEBUG_CHECKS
-   if (c->nodeID == engine_rank) error("Updating a local cell!");
- #endif
- 
-+  /* Clear this cell's sorted mask. */
-+  if (clear_sorts) c->stars.sorted = 0;
-+
-   /* If this cell is a leaf, collect the particle data. */
-   if (!c->split) {
- 
-@@ -2810,22 +2642,27 @@ void runner_do_recv_spart(struct runner *r, struct cell *c, int timer) {
-       if (sparts[k].time_bin == time_bin_inhibited) continue;
-       time_bin_min = min(time_bin_min, sparts[k].time_bin);
-       time_bin_max = max(time_bin_max, sparts[k].time_bin);
-+      h_max = max(h_max, sparts[k].h);
-     }
- 
-     /* Convert into a time */
-     ti_gravity_end_min = get_integer_time_end(ti_current, time_bin_min);
-     ti_gravity_end_max = get_integer_time_end(ti_current, time_bin_max);
-+    ti_stars_end_min = get_integer_time_end(ti_current, time_bin_min);
-   }
- 
-   /* Otherwise, recurse and collect. */
-   else {
-     for (int k = 0; k < 8; k++) {
-       if (c->progeny[k] != NULL && c->progeny[k]->stars.count > 0) {
--        runner_do_recv_spart(r, c->progeny[k], 0);
-+        runner_do_recv_spart(r, c->progeny[k], clear_sorts, 0);
-         ti_gravity_end_min =
-             min(ti_gravity_end_min, c->progeny[k]->grav.ti_end_min);
-         ti_gravity_end_max =
-             max(ti_gravity_end_max, c->progeny[k]->grav.ti_end_max);
-+        ti_stars_end_min =
-+            min(ti_stars_end_min, c->progeny[k]->stars.ti_end_min);
-+        h_max = max(h_max, c->progeny[k]->stars.h_max);
-       }
-     }
-   }
-@@ -2836,12 +2673,19 @@ void runner_do_recv_spart(struct runner *r, struct cell *c, int timer) {
-         "Received a cell at an incorrect time c->ti_end_min=%lld, "
-         "e->ti_current=%lld.",
-         ti_gravity_end_min, ti_current);
-+
-+  if (ti_stars_end_min < ti_current)
-+    error(
-+        "Received a cell at an incorrect time c->ti_end_min=%lld, "
-+        "e->ti_current=%lld.",
-+        ti_stars_end_min, ti_current);
- #endif
- 
-   /* ... and store. */
-   // c->grav.ti_end_min = ti_gravity_end_min;
-   // c->grav.ti_end_max = ti_gravity_end_max;
-   c->grav.ti_old_part = ti_current;
-+  c->stars.h_max = h_max;
- 
-   if (timer) TIMER_TOC(timer_dorecv_spart);
- 
-@@ -2922,13 +2766,11 @@ void *runner_main(void *data) {
- #endif
-           else if (t->subtype == task_subtype_force)
-             runner_doself2_branch_force(r, ci);
--          else if (t->subtype == task_subtype_limiter)
--            runner_doself2_branch_limiter(r, ci);
-           else if (t->subtype == task_subtype_grav)
-             runner_doself_recursive_grav(r, ci, 1);
-           else if (t->subtype == task_subtype_external_grav)
-             runner_do_grav_external(r, ci, 1);
--          else if (t->subtype == task_subtype_stars_density) 
-+          else if (t->subtype == task_subtype_stars_density)
-             runner_doself_stars_density(r, ci, 1);
-           else if (t->subtype == task_subtype_stars_feedback)
-             runner_doself_stars_feedback(r, ci, 1);
-@@ -2945,8 +2787,6 @@ void *runner_main(void *data) {
- #endif
-           else if (t->subtype == task_subtype_force)
-             runner_dopair2_branch_force(r, ci, cj);
--          else if (t->subtype == task_subtype_limiter)
--            runner_dopair2_branch_limiter(r, ci, cj);
-           else if (t->subtype == task_subtype_grav)
-             runner_dopair_recursive_grav(r, ci, cj, 1);
-           else if (t->subtype == task_subtype_stars_density)
-@@ -2966,8 +2806,6 @@ void *runner_main(void *data) {
- #endif
-           else if (t->subtype == task_subtype_force)
-             runner_dosub_self2_force(r, ci, 1);
--          else if (t->subtype == task_subtype_limiter)
--            runner_dosub_self2_limiter(r, ci, 1);
-           else if (t->subtype == task_subtype_stars_density)
-             runner_dosub_self_stars_density(r, ci, 1);
-           else if (t->subtype == task_subtype_stars_feedback)
-@@ -2985,8 +2823,6 @@ void *runner_main(void *data) {
- #endif
-           else if (t->subtype == task_subtype_force)
-             runner_dosub_pair2_force(r, ci, cj, t->flags, 1);
--          else if (t->subtype == task_subtype_limiter)
--            runner_dosub_pair2_limiter(r, ci, cj, t->flags, 1);
-           else if (t->subtype == task_subtype_stars_density)
-             runner_dosub_pair_stars_density(r, ci, cj, t->flags, 1);
-           else if (t->subtype == task_subtype_stars_feedback)
-@@ -3003,7 +2839,8 @@ void *runner_main(void *data) {
-           /* Reset the sort flags as our work here is done. */
-           t->flags = 0;
-           break;
--        case task_type_stars_sort:
-+        case task_type_stars_sort_local:
-+        case task_type_stars_sort_foreign:
-           /* Cleanup only if any of the indices went stale. */
-           runner_do_stars_sort(
-               r, ci, t->flags,
-@@ -3046,9 +2883,6 @@ void *runner_main(void *data) {
-         case task_type_timestep:
-           runner_do_timestep(r, ci, 1);
-           break;
--        case task_type_timestep_limiter:
--          runner_do_limiter(r, ci, 0, 1);
--          break;
- #ifdef WITH_MPI
-         case task_type_send:
-           if (t->subtype == task_subtype_tend) {
-@@ -3065,12 +2899,10 @@ void *runner_main(void *data) {
-             runner_do_recv_part(r, ci, 0, 1);
-           } else if (t->subtype == task_subtype_gradient) {
-             runner_do_recv_part(r, ci, 0, 1);
--          } else if (t->subtype == task_subtype_limiter) {
--            runner_do_recv_part(r, ci, 0, 1);
-           } else if (t->subtype == task_subtype_gpart) {
-             runner_do_recv_gpart(r, ci, 1);
-           } else if (t->subtype == task_subtype_spart) {
--            runner_do_recv_spart(r, ci, 1);
-+            runner_do_recv_spart(r, ci, 1, 1);
-           } else if (t->subtype == task_subtype_multipole) {
-             cell_unpack_multipoles(ci, (struct gravity_tensors *)t->buff);
-             free(t->buff);
-diff --git a/src/scheduler.c b/src/scheduler.c
-index ad6af73..d5fe7b6 100644
---- a/src/scheduler.c
-+++ b/src/scheduler.c
-@@ -60,42 +60,11 @@
- void scheduler_clear_active(struct scheduler *s) { s->active_count = 0; }
- 
- /**
-- * @brief Increase the space available for unlocks. Only call when
-- *        current index == s->size_unlock;
-- */
--static void scheduler_extend_unlocks(struct scheduler *s) {
--
--  /* Allocate the new buffer. */
--  const int size_unlocks_new = s->size_unlocks * 2;
--  struct task **unlocks_new =
--      (struct task **)malloc(sizeof(struct task *) * size_unlocks_new);
--  int *unlock_ind_new = (int *)malloc(sizeof(int) * size_unlocks_new);
--  if (unlocks_new == NULL || unlock_ind_new == NULL)
--    error("Failed to re-allocate unlocks.");
--
--  /* Wait for all writes to the old buffer to complete. */
--  while (s->completed_unlock_writes < s->size_unlocks)
--    ;
--
--  /* Copy the buffers. */
--  memcpy(unlocks_new, s->unlocks, sizeof(struct task *) * s->size_unlocks);
--  memcpy(unlock_ind_new, s->unlock_ind, sizeof(int) * s->size_unlocks);
--  free(s->unlocks);
--  free(s->unlock_ind);
--  s->unlocks = unlocks_new;
--  s->unlock_ind = unlock_ind_new;
--
--  /* Publish the new buffer size. */
--  s->size_unlocks = size_unlocks_new;
--}
--
--/**
-  * @brief Add an unlock_task to the given task.
-  *
-  * @param s The #scheduler.
-  * @param ta The unlocking #task.
-  * @param tb The #task that will be unlocked.
--
-  */
- void scheduler_addunlock(struct scheduler *s, struct task *ta,
-                          struct task *tb) {
-@@ -108,21 +77,37 @@ void scheduler_addunlock(struct scheduler *s, struct task *ta,
-   const int ind = atomic_inc(&s->nr_unlocks);
- 
-   /* Does the buffer need to be grown? */
--  if (ind == s->size_unlocks) scheduler_extend_unlocks(s);
--
--#ifdef SWIFT_DEBUG_CHECKS
--  if (ind > s->size_unlocks * 2)
--    message("unlocks guard enabled: %d / %d", ind, s->size_unlocks);
--#endif
-+  if (ind == s->size_unlocks) {
-+    /* Allocate the new buffer. */
-+    struct task **unlocks_new;
-+    int *unlock_ind_new;
-+    const int size_unlocks_new = s->size_unlocks * 2;
-+    if ((unlocks_new = (struct task **)malloc(sizeof(struct task *) *
-+                                              size_unlocks_new)) == NULL ||
-+        (unlock_ind_new = (int *)malloc(sizeof(int) * size_unlocks_new)) ==
-+            NULL)
-+      error("Failed to re-allocate unlocks.");
-+
-+    /* Wait for all writes to the old buffer to complete. */
-+    while (s->completed_unlock_writes < ind)
-+      ;
-+
-+    /* Copy the buffers. */
-+    memcpy(unlocks_new, s->unlocks, sizeof(struct task *) * ind);
-+    memcpy(unlock_ind_new, s->unlock_ind, sizeof(int) * ind);
-+    free(s->unlocks);
-+    free(s->unlock_ind);
-+    s->unlocks = unlocks_new;
-+    s->unlock_ind = unlock_ind_new;
-+
-+    /* Publish the new buffer size. */
-+    s->size_unlocks = size_unlocks_new;
-+  }
- 
-   /* Wait for there to actually be space at my index. */
-   while (ind > s->size_unlocks)
-     ;
- 
--  /* Guard against case when more than (old) s->size_unlocks unlocks
--   * are now pending. */
--  if (ind == s->size_unlocks) scheduler_extend_unlocks(s);
--
-   /* Write the unlock to the scheduler. */
-   s->unlocks[ind] = tb;
-   s->unlock_ind[ind] = ta - s->tasks;
-@@ -130,7 +115,7 @@ void scheduler_addunlock(struct scheduler *s, struct task *ta,
- }
- 
- /**
-- * @brief compute the number of similar dependencies
-+ * @brief compute the number of same dependencies
-  *
-  * @param s The #scheduler
-  * @param ta The #task
-@@ -198,11 +183,10 @@ struct task_dependency {
- };
- 
- #ifdef WITH_MPI
--
- /**
-  * @brief Define the #task_dependency for MPI
-  *
-- * @param tstype The MPI_Datatype to initialize
-+ * @param tstype The #MPI_Datatype to initialize
-  */
- void task_dependency_define(MPI_Datatype *tstype) {
- 
-@@ -528,7 +512,7 @@ void scheduler_write_dependencies(struct scheduler *s, int verbose) {
-   /* Be clean */
-   free(task_dep);
- 
--  if (verbose)
-+  if (verbose && s->nodeID == 0)
-     message("Printing task graph took %.3f %s.",
-             clocks_from_ticks(getticks() - tic), clocks_getunit());
- }
-@@ -1030,8 +1014,17 @@ static void scheduler_splittask_stars(struct task *t, struct scheduler *s) {
-     redo = 0;
- 
-     /* Empty task? */
-+    /* Need defines in order to evaluate after check for t->ci == NULL */
-+#define pair_no_part                                          \
-+  t->type == task_type_pair &&                                \
-+      (t->ci->stars.count == 0 || t->cj->hydro.count == 0) && \
-+      (t->cj->stars.count == 0 || t->ci->hydro.count == 0)
-+#define self_no_part           \
-+  t->type == task_type_self && \
-+      (t->ci->stars.count == 0 || t->ci->hydro.count == 0)
-+
-     if ((t->ci == NULL) || (t->type == task_type_pair && t->cj == NULL) ||
--        t->ci->stars.count == 0 || (t->cj != NULL && t->cj->stars.count == 0)) {
-+        (self_no_part) || (pair_no_part)) {
-       t->type = task_type_none;
-       t->subtype = task_subtype_none;
-       t->cj = NULL;
-@@ -1079,9 +1072,11 @@ static void scheduler_splittask_stars(struct task *t, struct scheduler *s) {
- 
-           /* Make a task for each pair of progeny */
-           for (int j = 0; j < 8; j++)
--            if (ci->progeny[j] != NULL && ci->progeny[j]->stars.count)
-+            if (ci->progeny[j] != NULL &&
-+                (ci->progeny[j]->stars.count || ci->progeny[j]->hydro.count))
-               for (int k = j + 1; k < 8; k++)
--                if (ci->progeny[k] != NULL && ci->progeny[k]->stars.count)
-+                if (ci->progeny[k] != NULL && (ci->progeny[k]->stars.count ||
-+                                               ci->progeny[k]->hydro.count))
-                   scheduler_splittask_stars(
-                       scheduler_addtask(s, task_type_pair, t->subtype,
-                                         sub_sid_flag[j][k], 0, ci->progeny[j],
-@@ -1110,14 +1105,25 @@ static void scheduler_splittask_stars(struct task *t, struct scheduler *s) {
-       double shift[3];
-       const int sid = space_getsid(s->space, &ci, &cj, shift);
- 
-+#ifdef SWIFT_DEBUG_CHECKS
-+      if (sid != t->flags)
-+        error("Got pair task with incorrect flags: sid=%d flags=%lld", sid,
-+              t->flags);
-+#endif
-+
-+      /* Compute number of interactions */
-+      const int ci_interaction = (ci->stars.count * cj->hydro.count);
-+      const int cj_interaction = (cj->stars.count * ci->hydro.count);
-+
-+      const int number_interactions = (ci_interaction + cj_interaction);
-+
-       /* Should this task be split-up? */
-       if (cell_can_split_pair_stars_task(ci) &&
-           cell_can_split_pair_stars_task(cj)) {
- 
-         /* Replace by a single sub-task? */
--        if (scheduler_dosub && /* Use division to avoid integer overflow. */
--            ci->stars.count * sid_scale[sid] <
--                space_subsize_pair_stars / cj->stars.count &&
-+        if (scheduler_dosub &&
-+            number_interactions * sid_scale[sid] < space_subsize_pair_stars &&
-             !sort_is_corner(sid)) {
- 
-           /* Make this task a sub task. */
-@@ -1466,15 +1472,17 @@ static void scheduler_splittask_stars(struct task *t, struct scheduler *s) {
- 
-         /* Otherwise, break it up if it is too large? */
-       } else if (scheduler_doforcesplit && ci->split && cj->split &&
--                 (ci->stars.count > space_maxsize / cj->stars.count)) {
-+                 (number_interactions > space_maxsize)) {
- 
-         /* Replace the current task. */
-         t->type = task_type_none;
- 
-         for (int j = 0; j < 8; j++)
--          if (ci->progeny[j] != NULL && ci->progeny[j]->stars.count)
-+          if (ci->progeny[j] != NULL &&
-+              (ci->progeny[j]->stars.count || ci->progeny[j]->hydro.count))
-             for (int k = 0; k < 8; k++)
--              if (cj->progeny[k] != NULL && cj->progeny[k]->stars.count) {
-+              if (cj->progeny[k] != NULL && (cj->progeny[k]->stars.count ||
-+                                             cj->progeny[k]->hydro.count)) {
-                 struct task *tl =
-                     scheduler_addtask(s, task_type_pair, t->subtype, 0, 0,
-                                       ci->progeny[j], cj->progeny[k]);
-@@ -1981,7 +1989,8 @@ void scheduler_reweight(struct scheduler *s, int verbose) {
-                (sizeof(int) * 8 - intrinsics_clz(t->ci->hydro.count));
-         break;
- 
--      case task_type_stars_sort:
-+      case task_type_stars_sort_local:
-+      case task_type_stars_sort_foreign:
-         cost = wscale * intrinsics_popcount(t->flags) * scount_i *
-                (sizeof(int) * 8 - intrinsics_clz(t->ci->stars.count));
-         break;
-@@ -2275,7 +2284,8 @@ void scheduler_enqueue(struct scheduler *s, struct task *t) {
-       case task_type_kick2:
-       case task_type_stars_ghost:
-       case task_type_logger:
--      case task_type_stars_sort:
-+      case task_type_stars_sort_local:
-+      case task_type_stars_sort_foreign:
-       case task_type_timestep:
-         qid = t->ci->super->owner;
-         break;
-diff --git a/src/serial_io.c b/src/serial_io.c
-index c05c43a..a273d03 100644
---- a/src/serial_io.c
-+++ b/src/serial_io.c
-@@ -776,6 +776,7 @@ void write_output_serial(struct engine* e, const char* baseName,
-                          int mpi_size, MPI_Comm comm, MPI_Info info) {
- 
-   hid_t h_file = 0, h_grp = 0;
-+  int periodic = e->s->periodic;
-   int numFiles = 1;
-   const struct part* parts = e->s->parts;
-   const struct xpart* xparts = e->s->xparts;
-@@ -845,24 +846,27 @@ void write_output_serial(struct engine* e, const char* baseName,
-     if (h_file < 0) error("Error while opening file '%s'.", fileName);
- 
-     /* Open header to write simulation properties */
-+    /* message("Writing runtime parameters..."); */
-+    h_grp = H5Gcreate(h_file, "/RuntimePars", H5P_DEFAULT, H5P_DEFAULT,
-+                      H5P_DEFAULT);
-+    if (h_grp < 0) error("Error while creating runtime parameters group\n");
-+
-+    /* Write the relevant information */
-+    io_write_attribute(h_grp, "PeriodicBoundariesOn", INT, &periodic, 1);
-+
-+    /* Close runtime parameters */
-+    H5Gclose(h_grp);
-+
-+    /* Open header to write simulation properties */
-     /* message("Writing file header..."); */
-     h_grp = H5Gcreate(h_file, "/Header", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
-     if (h_grp < 0) error("Error while creating file header\n");
- 
--    /* Convert basic output information to snapshot units */
--    const double factor_time =
--        units_conversion_factor(internal_units, snapshot_units, UNIT_CONV_TIME);
--    const double factor_length = units_conversion_factor(
--        internal_units, snapshot_units, UNIT_CONV_LENGTH);
--    const double dblTime = e->time * factor_time;
--    const double dim[3] = {e->s->dim[0] * factor_length,
--                           e->s->dim[1] * factor_length,
--                           e->s->dim[2] * factor_length};
--
-     /* Print the relevant information and print status */
--    io_write_attribute(h_grp, "BoxSize", DOUBLE, dim, 3);
-+    io_write_attribute(h_grp, "BoxSize", DOUBLE, e->s->dim, 3);
-+    double dblTime = e->time;
-     io_write_attribute(h_grp, "Time", DOUBLE, &dblTime, 1);
--    const int dimension = (int)hydro_dimension;
-+    int dimension = (int)hydro_dimension;
-     io_write_attribute(h_grp, "Dimension", INT, &dimension, 1);
-     io_write_attribute(h_grp, "Redshift", DOUBLE, &e->cosmology->z, 1);
-     io_write_attribute(h_grp, "Scale-factor", DOUBLE, &e->cosmology->a, 1);
-@@ -1024,32 +1028,6 @@ void write_output_serial(struct engine* e, const char* baseName,
-     H5Fclose(h_file);
-   }
- 
--  /* Now write the top-level cell structure */
--  hid_t h_file_cells = 0, h_grp_cells = 0;
--  if (mpi_rank == 0) {
--
--    /* Open the snapshot on rank 0 */
--    h_file_cells = H5Fopen(fileName, H5F_ACC_RDWR, H5P_DEFAULT);
--    if (h_file_cells < 0)
--      error("Error while opening file '%s' on rank %d.", fileName, mpi_rank);
--
--    /* Create the group we want in the file */
--    h_grp_cells = H5Gcreate(h_file_cells, "/Cells", H5P_DEFAULT, H5P_DEFAULT,
--                            H5P_DEFAULT);
--    if (h_grp_cells < 0) error("Error while creating cells group");
--  }
--
--  /* Write the location of the particles in the arrays */
--  io_write_cell_offsets(h_grp_cells, e->s->cdim, e->s->cells_top,
--                        e->s->nr_cells, e->s->width, mpi_rank, N_total, offset,
--                        internal_units, snapshot_units);
--
--  /* Close everything */
--  if (mpi_rank == 0) {
--    H5Gclose(h_grp_cells);
--    H5Fclose(h_file_cells);
--  }
--
-   /* Now loop over ranks and write the data */
-   for (int rank = 0; rank < mpi_size; ++rank) {
- 
-diff --git a/src/single_io.c b/src/single_io.c
-index 178f59a..93d5f7f 100644
---- a/src/single_io.c
-+++ b/src/single_io.c
-@@ -639,6 +639,7 @@ void write_output_single(struct engine* e, const char* baseName,
-                          const struct unit_system* snapshot_units) {
- 
-   hid_t h_file = 0, h_grp = 0;
-+  int periodic = e->s->periodic;
-   int numFiles = 1;
-   const struct part* parts = e->s->parts;
-   const struct xpart* xparts = e->s->xparts;
-@@ -698,24 +699,27 @@ void write_output_single(struct engine* e, const char* baseName,
-   if (h_file < 0) error("Error while opening file '%s'.", fileName);
- 
-   /* Open header to write simulation properties */
-+  /* message("Writing runtime parameters..."); */
-+  h_grp =
-+      H5Gcreate(h_file, "/RuntimePars", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
-+  if (h_grp < 0) error("Error while creating runtime parameters group\n");
-+
-+  /* Write the relevant information */
-+  io_write_attribute(h_grp, "PeriodicBoundariesOn", INT, &periodic, 1);
-+
-+  /* Close runtime parameters */
-+  H5Gclose(h_grp);
-+
-+  /* Open header to write simulation properties */
-   /* message("Writing file header..."); */
-   h_grp = H5Gcreate(h_file, "/Header", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
-   if (h_grp < 0) error("Error while creating file header\n");
- 
--  /* Convert basic output information to snapshot units */
--  const double factor_time =
--      units_conversion_factor(internal_units, snapshot_units, UNIT_CONV_TIME);
--  const double factor_length =
--      units_conversion_factor(internal_units, snapshot_units, UNIT_CONV_LENGTH);
--  const double dblTime = e->time * factor_time;
--  const double dim[3] = {e->s->dim[0] * factor_length,
--                         e->s->dim[1] * factor_length,
--                         e->s->dim[2] * factor_length};
--
-   /* Print the relevant information and print status */
--  io_write_attribute(h_grp, "BoxSize", DOUBLE, dim, 3);
-+  io_write_attribute(h_grp, "BoxSize", DOUBLE, e->s->dim, 3);
-+  double dblTime = e->time;
-   io_write_attribute(h_grp, "Time", DOUBLE, &dblTime, 1);
--  const int dimension = (int)hydro_dimension;
-+  int dimension = (int)hydro_dimension;
-   io_write_attribute(h_grp, "Dimension", INT, &dimension, 1);
-   io_write_attribute(h_grp, "Redshift", DOUBLE, &e->cosmology->z, 1);
-   io_write_attribute(h_grp, "Scale-factor", DOUBLE, &e->cosmology->a, 1);
-@@ -822,17 +826,6 @@ void write_output_single(struct engine* e, const char* baseName,
-   /* Print the system of Units used internally */
-   io_write_unit_system(h_file, internal_units, "InternalCodeUnits");
- 
--  /* Now write the top-level cell structure */
--  long long global_offsets[swift_type_count] = {0};
--  h_grp = H5Gcreate(h_file, "/Cells", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
--  if (h_grp < 0) error("Error while creating cells group");
--
--  /* Write the location of the particles in the arrays */
--  io_write_cell_offsets(h_grp, e->s->cdim, e->s->cells_top, e->s->nr_cells,
--                        e->s->width, e->nodeID, N_total, global_offsets,
--                        internal_units, snapshot_units);
--  H5Gclose(h_grp);
--
-   /* Tell the user if a conversion will be needed */
-   if (e->verbose) {
-     if (units_are_equal(snapshot_units, internal_units)) {
-diff --git a/src/space.c b/src/space.c
-index 1b5796f..a420eab 100644
---- a/src/space.c
-+++ b/src/space.c
-@@ -183,13 +183,13 @@ void space_rebuild_recycle_mapper(void *map_data, int num_elements,
-       space_recycle_list(s, cell_rec_begin, cell_rec_end, multipole_rec_begin,
-                          multipole_rec_end);
-     c->hydro.sorts = NULL;
--    c->stars.sorts = NULL;
-+    c->stars.sorts_local = NULL;
-+    c->stars.sorts_foreign = NULL;
-     c->nr_tasks = 0;
-     c->grav.nr_mm_tasks = 0;
-     c->hydro.density = NULL;
-     c->hydro.gradient = NULL;
-     c->hydro.force = NULL;
--    c->hydro.limiter = NULL;
-     c->grav.grav = NULL;
-     c->grav.mm = NULL;
-     c->hydro.dx_max_part = 0.0f;
-@@ -224,12 +224,12 @@ void space_rebuild_recycle_mapper(void *map_data, int num_elements,
-     c->kick1 = NULL;
-     c->kick2 = NULL;
-     c->timestep = NULL;
--    c->timestep_limiter = NULL;
-     c->end_force = NULL;
-     c->hydro.drift = NULL;
-     c->grav.drift = NULL;
-     c->grav.drift_out = NULL;
-     c->hydro.cooling = NULL;
-+    c->sourceterms = NULL;
-     c->grav.long_range = NULL;
-     c->grav.down_in = NULL;
-     c->grav.down = NULL;
-@@ -245,8 +245,6 @@ void space_rebuild_recycle_mapper(void *map_data, int num_elements,
-     c->stars.do_sub_sort = 0;
-     c->grav.do_sub_drift = 0;
-     c->hydro.do_sub_drift = 0;
--    c->hydro.do_sub_limiter = 0;
--    c->hydro.do_limiter = 0;
-     c->hydro.ti_end_min = -1;
-     c->hydro.ti_end_max = -1;
-     c->grav.ti_end_min = -1;
-@@ -274,15 +272,15 @@ void space_rebuild_recycle_mapper(void *map_data, int num_elements,
-     c->mpi.hydro.recv_rho = NULL;
-     c->mpi.hydro.recv_gradient = NULL;
-     c->mpi.grav.recv = NULL;
-+    c->mpi.stars.recv = NULL;
-     c->mpi.recv_ti = NULL;
--    c->mpi.limiter.recv = NULL;
- 
-     c->mpi.hydro.send_xv = NULL;
-     c->mpi.hydro.send_rho = NULL;
-     c->mpi.hydro.send_gradient = NULL;
-     c->mpi.grav.send = NULL;
-+    c->mpi.stars.send = NULL;
-     c->mpi.send_ti = NULL;
--    c->mpi.limiter.send = NULL;
- #endif
-   }
- }
-@@ -548,6 +546,8 @@ void space_regrid(struct space *s, int verbose) {
-           c->mpi.hydro.send_xv = NULL;
-           c->mpi.hydro.send_rho = NULL;
-           c->mpi.hydro.send_gradient = NULL;
-+          c->mpi.stars.send = NULL;
-+          c->mpi.stars.recv = NULL;
-           c->mpi.grav.recv = NULL;
-           c->mpi.grav.send = NULL;
- #endif  // WITH_MPI
-@@ -2712,8 +2712,6 @@ void space_split_recursive(struct space *s, struct cell *c,
-       cp->stars.do_sub_sort = 0;
-       cp->grav.do_sub_drift = 0;
-       cp->hydro.do_sub_drift = 0;
--      cp->hydro.do_sub_limiter = 0;
--      cp->hydro.do_limiter = 0;
- #ifdef WITH_MPI
-       cp->mpi.tag = -1;
- #endif  // WITH_MPI
-@@ -3511,9 +3509,6 @@ void space_first_init_sparts_mapper(void *restrict map_data, int count,
-     sp[k].v[1] *= a_factor_vel;
-     sp[k].v[2] *= a_factor_vel;
- 
--    // MATTHIEU: TO DO: Read smoothing lengths from ICs!!!!
--    //sp[k].h = s->dim[0] / 10. / kernel_gamma;
--
- #ifdef HYDRO_DIMENSION_2D
-     sp[k].x[2] = 0.f;
-     sp[k].v[2] = 0.f;
-@@ -4312,49 +4307,6 @@ void space_check_timesteps(struct space *s) {
- }
- 
- /**
-- * @brief #threadpool mapper function for the limiter debugging check
-- */
--void space_check_limiter_mapper(void *map_data, int nr_parts,
--                                void *extra_data) {
--#ifdef SWIFT_DEBUG_CHECKS
--  /* Unpack the data */
--  struct part *restrict parts = (struct part *)map_data;
--
--  /* Verify that all limited particles have been treated */
--  for (int k = 0; k < nr_parts; k++) {
--
--    if (parts[k].time_bin == time_bin_inhibited) continue;
--
--    if (parts[k].wakeup == time_bin_awake)
--      error("Particle still woken up! id=%lld", parts[k].id);
--
--    if (parts[k].gpart != NULL)
--      if (parts[k].time_bin != parts[k].gpart->time_bin)
--        error("Gpart not on the same time-bin as part");
--  }
--#else
--  error("Calling debugging code without debugging flag activated.");
--#endif
--}
--
--/**
-- * @brief Checks that all particles have their wakeup flag in a correct state.
-- *
-- * Should only be used for debugging purposes.
-- *
-- * @param s The #space to check.
-- */
--void space_check_limiter(struct space *s) {
--#ifdef SWIFT_DEBUG_CHECKS
--
--  threadpool_map(&s->e->threadpool, space_check_limiter_mapper, s->parts,
--                 s->nr_parts, sizeof(struct part), 1000, NULL);
--#else
--  error("Calling debugging code without debugging flag activated.");
--#endif
--}
--
--/**
-  * @brief Resets all the individual cell task counters to 0.
-  *
-  * Should only be used for debugging purposes.
-diff --git a/src/task.c b/src/task.c
-index 4d5695f..097a6f1 100644
---- a/src/task.c
-+++ b/src/task.c
-@@ -66,7 +66,6 @@ const char *taskID_names[task_type_count] = {"none",
-                                              "kick1",
-                                              "kick2",
-                                              "timestep",
--                                             "timestep_limiter",
-                                              "send",
-                                              "recv",
-                                              "grav_long_range",
-@@ -80,14 +79,15 @@ const char *taskID_names[task_type_count] = {"none",
-                                              "stars_ghost_in",
-                                              "stars_ghost",
-                                              "stars_ghost_out",
--                                             "stars_sort"};
-+                                             "stars_sort_local",
-+                                             "stars_sort_foreign"};
- 
- /* Sub-task type names. */
- const char *subtaskID_names[task_subtype_count] = {
--    "none",    "density",       "gradient",      "force",
--    "limiter", "grav",          "external_grav", "tend",
--    "xv",      "rho",           "gpart",         "multipole",
--    "spart",   "stars_density", "stars_feedback"};
-+    "none",          "density",       "gradient",  "force",
-+    "grav",          "external_grav", "tend",      "xv",
-+    "rho",           "gpart",         "multipole", "spart",
-+    "stars_density", "stars_feedback"};
- 
- #ifdef WITH_MPI
- /* MPI communicators for the subtypes. */
-@@ -141,7 +141,6 @@ __attribute__((always_inline)) INLINE static enum task_actions task_acts_on(
-     case task_type_sort:
-     case task_type_ghost:
-     case task_type_extra_ghost:
--    case task_type_timestep_limiter:
-     case task_type_cooling:
-       return task_action_part;
-       break;
-@@ -150,7 +149,8 @@ __attribute__((always_inline)) INLINE static enum task_actions task_acts_on(
-       return task_action_all;
- 
-     case task_type_stars_ghost:
--    case task_type_stars_sort:
-+    case task_type_stars_sort_local:
-+    case task_type_stars_sort_foreign:
-       return task_action_spart;
-       break;
- 
-@@ -163,7 +163,6 @@ __attribute__((always_inline)) INLINE static enum task_actions task_acts_on(
-         case task_subtype_density:
-         case task_subtype_gradient:
-         case task_subtype_force:
--        case task_subtype_limiter:
-           return task_action_part;
-           break;
- 
-@@ -340,8 +339,6 @@ void task_unlock(struct task *t) {
- 
-     case task_type_drift_part:
-     case task_type_sort:
--    case task_type_ghost:
--    case task_type_timestep_limiter:
-       cell_unlocktree(ci);
-       break;
- 
-@@ -350,7 +347,8 @@ void task_unlock(struct task *t) {
-       cell_gunlocktree(ci);
-       break;
- 
--    case task_type_stars_sort:
-+    case task_type_stars_sort_local:
-+    case task_type_stars_sort_foreign:
-       cell_sunlocktree(ci);
-       break;
- 
-@@ -467,13 +465,12 @@ int task_lock(struct task *t) {
- 
-     case task_type_drift_part:
-     case task_type_sort:
--    case task_type_ghost:
--    case task_type_timestep_limiter:
-       if (ci->hydro.hold) return 0;
-       if (cell_locktree(ci) != 0) return 0;
-       break;
- 
--    case task_type_stars_sort:
-+    case task_type_stars_sort_local:
-+    case task_type_stars_sort_foreign:
-       if (ci->stars.hold) return 0;
-       if (cell_slocktree(ci) != 0) return 0;
-       break;
-@@ -637,7 +634,7 @@ void task_print(const struct task *t) {
-  * graph.
-  *
-  * @param type The #task type.
-- * @param subtype The #task subtype.
-+ * @param subtype The #subtask type.
-  * @param cluster (return) The group name (should be allocated)
-  */
- void task_get_group_name(int type, int subtype, char *cluster) {
-@@ -662,11 +659,11 @@ void task_get_group_name(int type, int subtype, char *cluster) {
-     case task_subtype_grav:
-       strcpy(cluster, "Gravity");
-       break;
--    case task_subtype_limiter:
--      strcpy(cluster, "Timestep_limiter");
--      break;
-     case task_subtype_stars_density:
--      strcpy(cluster, "Stars");
-+      strcpy(cluster, "StarsDensity");
-+      break;
-+    case task_subtype_stars_feedback:
-+      strcpy(cluster, "StarsFeedback");
-       break;
-     default:
-       strcpy(cluster, "None");
-diff --git a/src/tools.c b/src/tools.c
-index b73fdff..c0400aa 100644
---- a/src/tools.c
-+++ b/src/tools.c
-@@ -343,7 +343,6 @@ void pairs_all_stars_density(struct runner *r, struct cell *ci,
-   const double dim[3] = {r->e->s->dim[0], r->e->s->dim[1], r->e->s->dim[2]};
-   const struct engine *e = r->e;
-   const struct cosmology *cosmo = e->cosmology;
--  const struct stars_props *stars_properties = e->stars_properties;
-   const float a = cosmo->a;
-   const float H = cosmo->H;
- 
-@@ -360,7 +359,6 @@ void pairs_all_stars_density(struct runner *r, struct cell *ci,
-     for (int j = 0; j < cj->hydro.count; ++j) {
- 
-       struct part *pj = &cj->hydro.parts[j];
--      struct xpart *xpj = &cj->hydro.xparts[j];
- 
-       /* Pairwise distance */
-       r2 = 0.0f;
-@@ -373,7 +371,7 @@ void pairs_all_stars_density(struct runner *r, struct cell *ci,
-       /* Hit or miss? */
-       if (r2 < hig2) {
-         /* Interact */
--        runner_iact_nonsym_stars_density(r2, dx, hi, pj->h, spi, pj, a, H, cosmo, stars_properties, xpj);
-+        runner_iact_nonsym_stars_density(r2, dx, hi, pj->h, spi, pj, a, H);
-       }
-     }
-   }
-@@ -391,7 +389,6 @@ void pairs_all_stars_density(struct runner *r, struct cell *ci,
-     for (int i = 0; i < ci->hydro.count; ++i) {
- 
-       struct part *pi = &ci->hydro.parts[i];
--      struct xpart *xpi = &ci->hydro.xparts[i];
- 
-       /* Pairwise distance */
-       r2 = 0.0f;
-@@ -404,7 +401,7 @@ void pairs_all_stars_density(struct runner *r, struct cell *ci,
-       /* Hit or miss? */
-       if (r2 < hjg2) {
-         /* Interact */
--        runner_iact_nonsym_stars_density(r2, dx, hj, pi->h, spj, pi, a, H, cosmo, stars_properties, xpi);
-+        runner_iact_nonsym_stars_density(r2, dx, hj, pi->h, spj, pi, a, H);
-       }
-     }
-   }
-@@ -507,10 +504,8 @@ void self_all_stars_density(struct runner *r, struct cell *ci) {
-   float r2, hi, hj, hig2, dxi[3];
-   struct spart *spi;
-   struct part *pj;
--  struct xpart *xpj;
-   const struct engine *e = r->e;
-   const struct cosmology *cosmo = e->cosmology;
--  const struct stars_props *stars_properties = e->stars_properties;
-   const float a = cosmo->a;
-   const float H = cosmo->H;
- 
-@@ -526,7 +521,6 @@ void self_all_stars_density(struct runner *r, struct cell *ci) {
-     for (int j = 0; j < ci->hydro.count; ++j) {
- 
-       pj = &ci->hydro.parts[j];
--      xpj = &ci->hydro.xparts[j];
-       hj = pj->h;
- 
-       /* Pairwise distance */
-@@ -539,7 +533,7 @@ void self_all_stars_density(struct runner *r, struct cell *ci) {
-       /* Hit or miss? */
-       if (r2 > 0.f && r2 < hig2) {
-         /* Interact */
--        runner_iact_nonsym_stars_density(r2, dxi, hi, hj, spi, pj, a, H, cosmo, stars_properties, xpj);
-+        runner_iact_nonsym_stars_density(r2, dxi, hi, hj, spi, pj, a, H);
-       }
-     }
-   }
-diff --git a/tests/testStellarEvolution.c b/tests/testStellarEvolution.c
-deleted file mode 100644
-index 1f5088a..0000000
---- a/tests/testStellarEvolution.c
-+++ /dev/null
-@@ -1,92 +0,0 @@
--/*******************************************************************************
-- * This file is part of SWIFT.
-- * Copyright (C) 2015 Matthieu Schaller (matthieu.schaller@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/>.
-- *
-- ******************************************************************************/
--#include "../config.h"
--
--/* Local headers. */
--#include "swift.h"
--
--/*
-- * @brief 
-- * 
-- */
--int main(int argc, char **argv) {
--  // Declare relevant structs
--  struct swift_params *params = malloc(sizeof(struct swift_params));
--  struct unit_system us;
--  struct chemistry_global_data chem_data;
--  struct part p;
--  struct xpart xp;
--  struct spart sp;
--  struct phys_const phys_const;
--  struct cosmology cosmo;
--  struct hydro_props hydro_properties;
--  struct stars_props stars_properties;
--  char *parametersFileName = "./testStellarEvolution.yml";
--
--  /* Read the parameter file */
--  if (params == NULL) error("Error allocating memory for the parameter file.");
--  message("Reading runtime parameters from file '%s'", parametersFileName);
--  parser_read_file(parametersFileName, params);
--
--  /* Init units */
--  units_init_from_params(&us, params, "InternalUnitSystem");
--  phys_const_init(&us, params, &phys_const);
--
--  /* Init chemistry */
--  chemistry_init(params, &us, &phys_const, &chem_data);
--  chemistry_first_init_part(&phys_const, &us, &cosmo, &chem_data, &p, &xp);
--  chemistry_print(&chem_data);
--
--  /* Init cosmology */
--  cosmology_init(params, &us, &phys_const, &cosmo);
--  cosmology_print(&cosmo);
--
--  /* Init hydro properties */
--  hydro_props_init(&hydro_properties, &phys_const, &us, params);
--
--  /* Init star properties */
--  stars_props_init(&stars_properties, &phys_const, &us, params, &hydro_properties);
--
--  /* Init spart */
--  stars_first_init_spart(&sp);
--
--  /* Evolve spart */
--  float dt = 1.0e-6;
--  stars_evolve_spart(&sp, &stars_properties, &cosmo, dt);
--
--  for (int i = 0; i < 9; i++) {
--    message("element %d to distribute fraction %.5e", i, sp.to_distribute.chemistry_data.metal_mass_fraction[i]);
--  }
--  message("to distribute mass %.5e",sp.to_distribute.mass);
--  message("to distribute num_SNIa %.5e",sp.to_distribute.num_SNIa);
--  message("to distribute ejecta_specific_thermal_energy %.5e",sp.to_distribute.ejecta_specific_thermal_energy);
--  message("to distribute metal_mass_fraction_total %.5e",sp.to_distribute.chemistry_data.metal_mass_fraction_total);
--  message("to distribute mass_from_AGB %.5e",sp.to_distribute.chemistry_data.mass_from_AGB);
--  message("to distribute metal_mass_fraction_from_AGB %.5e",sp.to_distribute.chemistry_data.metal_mass_fraction_from_AGB);
--  message("to distribute mass_from_SNII %.5e",sp.to_distribute.chemistry_data.mass_from_SNII);
--  message("to distribute metal_mass_fraction_from_SNII %.5e",sp.to_distribute.chemistry_data.metal_mass_fraction_from_SNII);
--  message("to distribute mass_from_SNIa %.5e",sp.to_distribute.chemistry_data.mass_from_SNIa);
--  message("to distribute metal_mass_fraction_from_SNIa %.5e",sp.to_distribute.chemistry_data.metal_mass_fraction_from_SNIa);
--  message("to distribute iron_mass_fraction_from_SNIa %.5e",sp.to_distribute.chemistry_data.iron_mass_fraction_from_SNIa);
--
--  message("done test");
--
--  free(params);
--  return 0;
--}
diff --git a/out_stars.txt b/out_stars.txt
deleted file mode 100644
index 28ef2fc90b520b1ee823e0c2096ba862310d3d58..0000000000000000000000000000000000000000
--- a/out_stars.txt
+++ /dev/null
@@ -1,987 +0,0 @@
-diff --git a/src/stars/Default/stars.h b/src/stars/Default/stars.h
-index 8632057..a8faf43 100644
---- a/src/stars/Default/stars.h
-+++ b/src/stars/Default/stars.h
-@@ -56,8 +56,6 @@ __attribute__((always_inline)) INLINE static void stars_init_spart(
-     struct spart* sp) {
- 
- #ifdef DEBUG_INTERACTIONS_STARS
--  for (int i = 0; i < MAX_NUM_OF_NEIGHBOURS_STARS; ++i)
--    sp->ids_ngbs_density[i] = -1;
-   sp->num_ngb_density = 0;
- #endif
- 
-@@ -145,6 +143,20 @@ __attribute__((always_inline)) INLINE static void stars_spart_has_no_neighbours(
-  */
- __attribute__((always_inline)) INLINE static void stars_evolve_spart(
-     struct spart* restrict sp, const struct stars_props* stars_properties,
--    const struct cosmology* cosmo, float dt) {}
-+    const struct cosmology* cosmo) {}
- 
-+/**
-+ * @brief Reset acceleration fields of a particle
-+ *
-+ * This is the equivalent of hydro_reset_acceleration.
-+ * We do not compute the acceleration on star, therefore no need to use it.
-+ *
-+ * @param p The particle to act upon
-+ */
-+__attribute__((always_inline)) INLINE static void stars_reset_acceleration(
-+    struct spart* restrict p) {
-+#ifdef DEBUG_INTERACTIONS_STARS
-+  p->num_ngb_force = 0;
-+#endif
-+}
- #endif /* SWIFT_DEFAULT_STARS_H */
-diff --git a/src/stars/Default/stars_iact.h b/src/stars/Default/stars_iact.h
-index a0da318..7e36d4c 100644
---- a/src/stars/Default/stars_iact.h
-+++ b/src/stars/Default/stars_iact.h
-@@ -14,9 +14,7 @@ __attribute__((always_inline)) INLINE static void
- runner_iact_nonsym_stars_density(float r2, const float *dx, float hi, float hj,
-                                  struct spart *restrict si,
-                                  const struct part *restrict pj, float a,
--                                 float H, const struct cosmology *restrict cosmo,
--                                 const struct stars_props *restrict stars_properties,
--                                 struct xpart *restrict xp) {
-+                                 float H) {
- 
-   float wi, wi_dx;
- 
-@@ -35,8 +33,6 @@ runner_iact_nonsym_stars_density(float r2, const float *dx, float hi, float hj,
- 
- #ifdef DEBUG_INTERACTIONS_STARS
-   /* Update ngb counters */
--  if (si->num_ngb_density < MAX_NUM_OF_NEIGHBOURS_STARS)
--    si->ids_ngbs_density[si->num_ngb_density] = pj->id;
-   ++si->num_ngb_density;
- #endif
- }
-@@ -56,9 +52,10 @@ runner_iact_nonsym_stars_density(float r2, const float *dx, float hi, float hj,
- __attribute__((always_inline)) INLINE static void
- runner_iact_nonsym_stars_feedback(float r2, const float *dx, float hi, float hj,
-                                   struct spart *restrict si,
--                                  struct part *restrict pj, float a, float H,
--                                  const struct cosmology *restrict cosmo,
--                                  const struct stars_props *restrict stars_properties,
--                                  struct xpart *restrict xp) {
--  si->to_distribute.mass = 2;
-+                                  struct part *restrict pj, float a, float H) {
-+
-+#ifdef DEBUG_INTERACTIONS_STARS
-+  /* Update ngb counters */
-+  ++si->num_ngb_force;
-+#endif
- }
-diff --git a/src/stars/Default/stars_io.h b/src/stars/Default/stars_io.h
-index 8b0c75f..fbed04a 100644
---- a/src/stars/Default/stars_io.h
-+++ b/src/stars/Default/stars_io.h
-@@ -60,10 +60,10 @@ INLINE static void stars_write_particles(const struct spart *sparts,
-                                          struct io_props *list,
-                                          int *num_fields) {
- 
--  /* Say how much we want to read */
-+  /* Say how much we want to write */
-   *num_fields = 5;
- 
--  /* List what we want to read */
-+  /* List what we want to write */
-   list[0] = io_make_output_field("Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH,
-                                  sparts, x);
-   list[1] =
-@@ -74,6 +74,17 @@ INLINE static void stars_write_particles(const struct spart *sparts,
-                                  sparts, id);
-   list[4] = io_make_output_field("SmoothingLength", FLOAT, 1, UNIT_CONV_LENGTH,
-                                  sparts, h);
-+
-+#ifdef DEBUG_INTERACTIONS_STARS
-+
-+  list += *num_fields;
-+  *num_fields += 2;
-+
-+  list[0] = io_make_output_field("Num_ngb_density", INT, 1, UNIT_CONV_NO_UNITS,
-+                                 sparts, num_ngb_density);
-+  list[1] = io_make_output_field("Num_ngb_force", INT, 1, UNIT_CONV_NO_UNITS,
-+                                 sparts, num_ngb_force);
-+#endif
- }
- 
- /**
-@@ -91,8 +102,7 @@ INLINE static void stars_props_init(struct stars_props *sp,
-                                     const struct phys_const *phys_const,
-                                     const struct unit_system *us,
-                                     struct swift_params *params,
--                                    const struct hydro_props *p,
--				    const struct cosmology *cosmo) {
-+                                    const struct hydro_props *p) {
- 
-   /* Kernel properties */
-   sp->eta_neighbours = parser_get_opt_param_float(
-diff --git a/src/stars/Default/stars_part.h b/src/stars/Default/stars_part.h
-index 089f9c2..c487528 100644
---- a/src/stars/Default/stars_part.h
-+++ b/src/stars/Default/stars_part.h
-@@ -65,12 +65,6 @@ struct spart {
- 
-   } density;
- 
--  struct {
--    /* Mass of ejecta */
--    float mass;
--
--  } to_distribute;
--
- #ifdef SWIFT_DEBUG_CHECKS
- 
-   /* Time of the last drift */
-@@ -82,11 +76,11 @@ struct spart {
- #endif
- 
- #ifdef DEBUG_INTERACTIONS_STARS
--  /*! List of interacting particles in the density SELF and PAIR */
--  long long ids_ngbs_density[MAX_NUM_OF_NEIGHBOURS_STARS];
--
-   /*! Number of interactions in the density SELF and PAIR */
-   int num_ngb_density;
-+
-+  /*! Number of interactions in the force SELF and PAIR */
-+  int num_ngb_force;
- #endif
- 
- } SWIFT_STRUCT_ALIGN;
-diff --git a/src/stars/const/stars.h b/src/stars/const/stars.h
-deleted file mode 100644
-index d043b72..0000000
---- a/src/stars/const/stars.h
-+++ /dev/null
-@@ -1,177 +0,0 @@
--/*******************************************************************************
-- * This file is part of SWIFT.
-- * Coypright (c) 2016 Matthieu Schaller (matthieu.schaller@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/>.
-- *
-- ******************************************************************************/
--#ifndef SWIFT_CONST_STARS_H
--#define SWIFT_CONST_STARS_H
--
--#include <float.h>
--#include "minmax.h"
--
--/**
-- * @brief Computes the gravity time-step of a given star particle.
-- *
-- * @param sp Pointer to the s-particle data.
-- */
--__attribute__((always_inline)) INLINE static float stars_compute_timestep(
--    const struct spart* const sp) {
--
--  return FLT_MAX;
--}
--
--/**
-- * @brief Initialises the s-particles for the first time
-- *
-- * This function is called only once just after the ICs have been
-- * read in to do some conversions.
-- *
-- * @param sp The particle to act upon
-- */
--__attribute__((always_inline)) INLINE static void stars_first_init_spart(
--    struct spart* sp) {
--
--  sp->time_bin = 0;
--}
--
--/**
-- * @brief Prepares a s-particle for its interactions
-- *
-- * @param sp The particle to act upon
-- */
--__attribute__((always_inline)) INLINE static void stars_init_spart(
--    struct spart* sp) {
--
--#ifdef DEBUG_INTERACTIONS_STARS
--  for (int i = 0; i < MAX_NUM_OF_NEIGHBOURS_STARS; ++i)
--    sp->ids_ngbs_density[i] = -1;
--  sp->num_ngb_density = 0;
--#endif
--
--  sp->density.wcount = 0.f;
--  sp->density.wcount_dh = 0.f;
--
--  sp->omega_normalisation_inv = 0.f;
--  sp->ngb_mass = 0.f;
--}
--
--/**
-- * @brief Sets the values to be predicted in the drifts to their values at a
-- * kick time
-- *
-- * @param sp The particle.
-- */
--__attribute__((always_inline)) INLINE static void stars_reset_predicted_values(
--    struct spart* restrict sp) {}
--
--/**
-- * @brief Finishes the calculation of (non-gravity) forces acting on stars
-- *
-- * Multiplies the forces and accelerations by the appropiate constants
-- *
-- * @param sp The particle to act upon
-- */
--__attribute__((always_inline)) INLINE static void stars_end_force(
--    struct spart* sp) {}
--
--/**
-- * @brief Kick the additional variables
-- *
-- * @param sp The particle to act upon
-- * @param dt The time-step for this kick
-- */
--__attribute__((always_inline)) INLINE static void stars_kick_extra(
--    struct spart* sp, float dt) {}
--
--/**
-- * @brief Finishes the calculation of density on stars
-- *
-- * @param sp The particle to act upon
-- * @param cosmo The current cosmological model.
-- */
--__attribute__((always_inline)) INLINE static void stars_end_density(
--    struct spart* sp, const struct cosmology* cosmo) {
--
--  /* Some smoothing length multiples. */
--  const float h = sp->h;
--  const float h_inv = 1.0f / h;                       /* 1/h */
--  const float h_inv_dim = pow_dimension(h_inv);       /* 1/h^d */
--  const float h_inv_dim_plus_one = h_inv_dim * h_inv; /* 1/h^(d+1) */
--
--  /* Finish the calculation by inserting the missing h-factors */
--  sp->density.wcount *= h_inv_dim;
--  sp->density.wcount_dh *= h_inv_dim_plus_one;
--}
--
--/**
-- * @brief Sets all particle fields to sensible values when the #spart has 0
-- * ngbs.
-- *
-- * @param sp The particle to act upon
-- * @param cosmo The current cosmological model.
-- */
--__attribute__((always_inline)) INLINE static void stars_spart_has_no_neighbours(
--    struct spart* restrict sp, const struct cosmology* cosmo) {
--
--  /* Some smoothing length multiples. */
--  const float h = sp->h;
--  const float h_inv = 1.0f / h;                 /* 1/h */
--  const float h_inv_dim = pow_dimension(h_inv); /* 1/h^d */
--
--  /* Re-set problematic values */
--  sp->density.wcount = kernel_root * h_inv_dim;
--  sp->density.wcount_dh = 0.f;
--}
--
--/**
-- * @brief Evolve the stellar properties of a #spart.
-- *
-- * This function allows for example to compute the SN rate before sending
-- * this information to a different MPI rank.
-- *
-- * @param sp The particle to act upon
-- * @param cosmo The current cosmological model.
-- * @param stars_properties The #stars_props
-- * @param dt Timestep over which the particle evolves.
-- */
--__attribute__((always_inline)) INLINE static void stars_evolve_spart(
--    struct spart* restrict sp, const struct stars_props* stars_properties,
--    const struct cosmology* cosmo, double dt) {
--  
--  /* Proportion of quantities to be released each timestep */
--  float feedback_factor = dt/stars_properties->feedback_timescale;
--
--  /* Amount of mass to distribute in one step */
--  sp->to_distribute.mass = sp->mass_init * feedback_factor;
--
--  /* Set all enrichment quantities to constant values */
--  for (int i = 0; i < chemistry_element_count; i++) sp->to_distribute.chemistry_data.metal_mass_fraction[i] = 1.f/chemistry_element_count;
--  sp->to_distribute.chemistry_data.metal_mass_fraction_total = 1.f - 2.f/chemistry_element_count;
--  sp->to_distribute.chemistry_data.mass_from_AGB = 1.0e-2 * sp->to_distribute.mass;
--  sp->to_distribute.chemistry_data.metal_mass_fraction_from_AGB = 1.0e-2;
--  sp->to_distribute.chemistry_data.mass_from_SNII = 1.0e-2 * sp->to_distribute.mass;
--  sp->to_distribute.chemistry_data.metal_mass_fraction_from_SNII = 1.0e-2;
--  sp->to_distribute.chemistry_data.mass_from_SNIa = 1.0e-2 * sp->to_distribute.mass;
--  sp->to_distribute.chemistry_data.metal_mass_fraction_from_SNIa = 1.0e-2;
--  sp->to_distribute.chemistry_data.iron_mass_fraction_from_SNIa = 1.0e-2;
--
--  /* Set feedback to constant values */
--  sp->to_distribute.num_SNIa = 1.0e11 * dt;
--  sp->to_distribute.ejecta_specific_thermal_energy = 1.0e-3;
--  
--}
--
--#endif /* SWIFT_CONST_STARS_H */
-diff --git a/src/stars/const/stars_debug.h b/src/stars/const/stars_debug.h
-deleted file mode 100644
-index 53a54b0..0000000
---- a/src/stars/const/stars_debug.h
-+++ /dev/null
-@@ -1,31 +0,0 @@
--/*******************************************************************************
-- * This file is part of SWIFT.
-- * Coypright (c) 2016 Matthieu Schaller (matthieu.schaller@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/>.
-- *
-- ******************************************************************************/
--#ifndef SWIFT_CONST_STARS_DEBUG_H
--#define SWIFT_CONST_STARS_DEBUG_H
--
--__attribute__((always_inline)) INLINE static void stars_debug_particle(
--    const struct spart* p) {
--  printf(
--      "x=[%.3e,%.3e,%.3e], "
--      "v_full=[%.3e,%.3e,%.3e] p->mass=%.3e \n t_begin=%d, t_end=%d\n",
--      p->x[0], p->x[1], p->x[2], p->v_full[0], p->v_full[1], p->v_full[2],
--      p->mass, p->ti_begin, p->ti_end);
--}
--
--#endif /* SWIFT_CONST_STARS_DEBUG_H */
-diff --git a/src/stars/const/stars_iact.h b/src/stars/const/stars_iact.h
-deleted file mode 100644
-index 3bce69f..0000000
---- a/src/stars/const/stars_iact.h
-+++ /dev/null
-@@ -1,212 +0,0 @@
--/**
-- * @brief Density interaction between two particles (non-symmetric).
-- *
-- * @param r2 Comoving square distance between the two particles.
-- * @param dx Comoving vector separating both particles (pi - pj).
-- * @param hi Comoving smoothing-length of particle i.
-- * @param hj Comoving smoothing-length of particle j.
-- * @param si First sparticle.
-- * @param pj Second particle (not updated).
-- * @param a Current scale factor.
-- * @param H Current Hubble parameter.
-- */
--__attribute__((always_inline)) INLINE static void
--runner_iact_nonsym_stars_density(float r2, const float *dx, float hi, float hj,
--                                 struct spart *restrict si,
--                                 const struct part *restrict pj, float a,
--                                 float H, const struct cosmology *restrict cosmo,
--				 const struct stars_props *restrict stars_properties,
--				 struct xpart *restrict xp) {
--
--  float wi, wi_dx;
--
--  /* Get r and 1/r. */
--  const float r_inv = 1.0f / sqrtf(r2);
--  const float r = r2 * r_inv;
--
--  /* Compute the kernel function */
--  const float hi_inv = 1.0f / hi;
--  const float ui = r * hi_inv;
--  kernel_deval(ui, &wi, &wi_dx);
--  
--  float wj, wj_dx;
--  const float hj_inv = 1.0f / hj;
--  const float uj = r * hj_inv;
--  kernel_deval(uj, &wj, &wj_dx);
--
--  /* Compute contribution to the number of neighbours */
--  si->density.wcount += wi;
--  si->density.wcount_dh -= (hydro_dimension * wi + ui * wi_dx);
--
--  /* Add mass of pj to neighbour mass of si  */
--  si->ngb_mass += hydro_get_mass(pj);
--
--  /* Add contribution of pj to normalisation of kernel (IMPROVE COMMENT?) */
--  si->omega_normalisation_inv += wj / hydro_get_physical_density(pj,cosmo);
--
--#ifdef DEBUG_INTERACTIONS_STARS
--  /* Update ngb counters */
--  if (si->num_ngb_density < MAX_NUM_OF_NEIGHBOURS_STARS)
--    si->ids_ngbs_density[si->num_ngb_density] = pj->id;
--  ++si->num_ngb_density;
--#endif
--}
--
--/**
-- * @brief Increases thermal energy of particle due 
-- * to feedback by specified amount
-- *
-- * @param du change in internal energy
-- * @param p Particle we're acting on
-- * @param cosmo Cosmology struct
-- */
--static inline void thermal_feedback(float du, struct part * restrict p, 
--				    struct xpart * restrict xp,
--				    const struct cosmology * restrict cosmo) {
--
--  float u = hydro_get_physical_internal_energy(p, xp, cosmo);
--  hydro_set_physical_internal_energy(p, cosmo, u + du);
--  // Just setting p->entropy is not enough because xp->entropy_full gets updated with p->entropy_dt
--  // TODO: ADD HYDRO FUNCTIONS FOR UPDATING DRIFTED AND NON DRIFTED INTERNAL ENERGY AND GET RID OF 
--  // THE ENTROPY UPDATE HERE.
--  xp->entropy_full = p->entropy;
--}
--
--/**
-- * @brief Feedback interaction between two particles (non-symmetric).
-- * Used for updating properties of gas particles neighbouring a star particle
-- *
-- * @param r2 Comoving square distance between the two particles.
-- * @param dx Comoving vector separating both particles (si - pj).
-- * @param hi Comoving smoothing-length of particle i.
-- * @param hj Comoving smoothing-length of particle j.
-- * @param si First (star) particle.
-- * @param pj Second (gas) particle.
-- * @param a Current scale factor.
-- * @param H Current Hubble parameter.
-- */
--__attribute__((always_inline)) INLINE static void
--runner_iact_nonsym_stars_feedback(float r2, const float *dx, float hi, float hj,
--                                  struct spart *restrict si,
--                                  struct part *restrict pj, float a, float H,
--				  const struct cosmology *restrict cosmo,
--				  const struct stars_props *restrict stars_properties,
--				  struct xpart *restrict xp) {
--  // ALEXEI: GET RID OF a AND H IN SIGNATURE SINCE THESE CAN BE DERIVED FROM COSMO?
--
--  float wj;
--
--  /* Get r and 1/r. */
--  const float r_inv = 1.0f / sqrtf(r2);
--  const float r = r2 * r_inv;
--
--  /* Compute the kernel function */
--  const float hj_inv = 1.0f / hj;
--  const float uj = r * hj_inv;
--  kernel_eval(uj, &wj);
--
--  /* Compute weighting for distributing various properties (TODO: better comment?) */
--  // ALEXEI: come up with better name for omega_frac?
--  float omega_frac = wj/hydro_get_physical_density(pj,cosmo)/si->omega_normalisation_inv;
--
--  /* Update particle mass */
--  float current_mass = hydro_get_mass(pj);
--  //float new_mass = current_mass + si->to_distribute.mass*omega_frac;
--  float new_mass = current_mass;
--  hydro_set_mass(pj,new_mass);
--
--  // ALEXEI: do we want to use the smoothed mass fraction?
--  /* Update total metallicity */
--  float current_metal_mass_total = pj->chemistry_data.metal_mass_fraction_total * current_mass;
--  float new_metal_mass_total = current_metal_mass_total + si->to_distribute.chemistry_data.metal_mass_fraction_total * 
--    			       si->to_distribute.mass * omega_frac;
--  pj->chemistry_data.metal_mass_fraction_total = new_metal_mass_total/new_mass;
--  
--  /* Update mass fraction of each tracked element  */
--  for (int elem = 0; elem < chemistry_element_count; elem++) {
--    float current_metal_mass = pj->chemistry_data.metal_mass_fraction[elem] * current_mass;
--    float new_metal_mass = current_metal_mass + si->to_distribute.chemistry_data.metal_mass_fraction[elem] * 
--    			   si->to_distribute.mass * omega_frac;
--    pj->chemistry_data.metal_mass_fraction[elem] = new_metal_mass/new_mass;
--  }
--
--  /* Update iron mass fraction from SNIa  */
--  float current_iron_from_SNIa_mass = pj->chemistry_data.iron_mass_fraction_from_SNIa * current_mass;
--  float new_iron_from_SNIa_mass = current_iron_from_SNIa_mass + si->to_distribute.chemistry_data.iron_mass_fraction_from_SNIa * 
--    			       si->to_distribute.mass * omega_frac;
--  pj->chemistry_data.iron_mass_fraction_from_SNIa = new_iron_from_SNIa_mass/new_mass;
--
--  /* Update mass from SNIa */
--  pj->chemistry_data.mass_from_SNIa += si->to_distribute.chemistry_data.mass_from_SNIa * omega_frac;
--
--  /* Update metal mass fraction from SNIa */
--  float current_metal_mass_fraction_from_SNIa = pj->chemistry_data.metal_mass_fraction_from_SNIa * current_mass;
--  float new_metal_mass_fraction_from_SNIa = current_metal_mass_fraction_from_SNIa + si->to_distribute.chemistry_data.metal_mass_fraction_from_SNIa * 
--    			       si->to_distribute.mass * omega_frac;
--  pj->chemistry_data.metal_mass_fraction_from_SNIa = new_metal_mass_fraction_from_SNIa/new_mass;
--
--  /* Update mass from SNII */
--  pj->chemistry_data.mass_from_SNII += si->to_distribute.chemistry_data.mass_from_SNII * omega_frac;
--
--  /* Update metal mass fraction from SNII */
--  float current_metal_mass_fraction_from_SNII = pj->chemistry_data.metal_mass_fraction_from_SNII * current_mass;
--  float new_metal_mass_fraction_from_SNII = current_metal_mass_fraction_from_SNII + si->to_distribute.chemistry_data.metal_mass_fraction_from_SNII * 
--    			       si->to_distribute.mass * omega_frac;
--  pj->chemistry_data.metal_mass_fraction_from_SNII = new_metal_mass_fraction_from_SNII/new_mass;
--
--  /* Update mass from AGB */
--  pj->chemistry_data.mass_from_AGB += si->to_distribute.chemistry_data.mass_from_AGB * omega_frac;
--
--  /* Update metal mass fraction from AGB */
--  float current_metal_mass_fraction_from_AGB = pj->chemistry_data.metal_mass_fraction_from_AGB * current_mass;
--  float new_metal_mass_fraction_from_AGB = current_metal_mass_fraction_from_AGB + si->to_distribute.chemistry_data.metal_mass_fraction_from_AGB * 
--    			       si->to_distribute.mass * omega_frac;
--  pj->chemistry_data.metal_mass_fraction_from_AGB = new_metal_mass_fraction_from_AGB/new_mass;
--
--  /* Update momentum */
--  for (int i = 0; i < 3; i++) {
--    // Do we need to calculate relative velocities here?
--    pj->v[i] += si->to_distribute.mass * omega_frac * si->v[i];
--  }
--
--  /* Energy feedback */
--  float heating_probability = -1.f, du = 0.f, d_energy = 0.f;
--  d_energy = si->to_distribute.mass * (si->to_distribute.ejecta_specific_thermal_energy 
--     + 0.5*(si->v[0]*si->v[0] + si->v[1]*si->v[1] + si->v[2]*si->v[2]) * cosmo->a2_inv);
--
--  if (stars_properties->continuous_heating) {
--    // We're doing ONLY continuous heating 
--    d_energy += si->to_distribute.num_SNIa * stars_properties->total_energy_SNe * omega_frac * si->mass_init;
--    du = d_energy/hydro_get_mass(pj);
--    thermal_feedback(du,pj,xp,cosmo);
--  } else {
--    // We're doing stochastic heating
--    heating_probability = stars_properties->units_factor1 * si->to_distribute.num_SNIa *
--                          stars_properties->SNIa_energy_fraction /
--                          (stars_properties->deltaT_desired * si->ngb_mass);
--    //message("probability %.5e factor1 %.5e num_snia %.5e energy_fraction %.5e deltaT_desired %.5e ngb_mass %.5e temp_to_u_factor %.5e",heating_probability, stars_properties->units_factor1, si->to_distribute.num_SNIa,
--    //                      stars_properties->SNIa_energy_fraction, stars_properties->deltaT_desired, si->ngb_mass, stars_properties->temp_to_u_factor);
--    // ALEXEI: CHECK UNITS HERE. Eagle does this update in cgs, we should probably keep it in internal units.
--    du = stars_properties->deltaT_desired * stars_properties->temp_to_u_factor;
--    if (heating_probability >= 1) {
--      du = stars_properties->units_factor2 * si->to_distribute.num_SNIa / si->ngb_mass;
--      heating_probability = 1; 
--    }
--  }
--
--  /* pick random number to see if we do stochastic heating */
--  unsigned int seed = pj->id;
--  double random_num = rand_r(&seed) * stars_properties->inv_rand_max;
--  if (random_num < heating_probability) {
--    // ALEXEI: As above, check units
--    thermal_feedback(du, pj, xp, cosmo);
--    // Debugging...
--    message("we did some heating! id %llu probability %.5e random_num %.5e old energy %.5e energy change %.5e", pj->id, heating_probability, random_num, hydro_get_physical_internal_energy(pj,xp,cosmo)*hydro_get_mass(pj), du*hydro_get_mass(pj));
--  } else {
--    // Debugging...
--    //message("we missed heating... id %llu probability %.5e random_num %.5e du %.5e", pj->id, heating_probability, random_num, du);
--  }
--
--  /* Decrease the mass of star particle (TO CHECK: WHAT ABOUT INTERNAL ENERGY?); */
--  si->mass -= si->to_distribute.mass;
--}
-diff --git a/src/stars/const/stars_io.h b/src/stars/const/stars_io.h
-deleted file mode 100644
-index cdd317f..0000000
---- a/src/stars/const/stars_io.h
-+++ /dev/null
-@@ -1,225 +0,0 @@
--/*******************************************************************************
-- * This file is part of SWIFT.
-- * Coypright (c) 2016 Matthieu Schaller (matthieu.schaller@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/>.
-- *
-- ******************************************************************************/
--#ifndef SWIFT_CONST_STARS_IO_H
--#define SWIFT_CONST_STARS_IO_H
--
--#include "io_properties.h"
--#include "stars_part.h"
--
--/**
-- * @brief Specifies which s-particle fields to read from a dataset
-- *
-- * @param sparts The s-particle array.
-- * @param list The list of i/o properties to read.
-- * @param num_fields The number of i/o fields to read.
-- */
--INLINE static void stars_read_particles(struct spart *sparts,
--                                        struct io_props *list,
--                                        int *num_fields) {
--
--  /* Say how much we want to read */
--  *num_fields = 6;
--
--  /* List what we want to read */
--  list[0] = io_make_input_field("Coordinates", DOUBLE, 3, COMPULSORY,
--                                UNIT_CONV_LENGTH, sparts, x);
--  list[1] = io_make_input_field("Velocities", FLOAT, 3, COMPULSORY,
--                                UNIT_CONV_SPEED, sparts, v);
--  list[2] = io_make_input_field("Masses", FLOAT, 1, COMPULSORY, UNIT_CONV_MASS,
--                                sparts, mass);
--  /* Temporary addition, discuss with Folkert what to do about initial mass
--   * when sparts are read from ICs. */
--  list[3] = io_make_input_field("Masses", FLOAT, 1, COMPULSORY, UNIT_CONV_MASS,
--                                sparts, mass_init);
--  list[4] = io_make_input_field("ParticleIDs", LONGLONG, 1, COMPULSORY,
--                                UNIT_CONV_NO_UNITS, sparts, id);
--  list[5] = io_make_input_field("SmoothingLength", FLOAT, 1, OPTIONAL,
--                                UNIT_CONV_LENGTH, sparts, h);
--}
--
--/**
-- * @brief Specifies which s-particle fields to write to a dataset
-- *
-- * @param sparts The s-particle array.
-- * @param list The list of i/o properties to write.
-- * @param num_fields The number of i/o fields to write.
-- */
--INLINE static void stars_write_particles(const struct spart *sparts,
--                                         struct io_props *list,
--                                         int *num_fields) {
--
--  /* Say how much we want to read */
--  *num_fields = 5;
--
--  /* List what we want to read */
--  list[0] = io_make_output_field("Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH,
--                                 sparts, x);
--  list[1] =
--      io_make_output_field("Velocities", FLOAT, 3, UNIT_CONV_SPEED, sparts, v);
--  list[2] =
--      io_make_output_field("Masses", FLOAT, 1, UNIT_CONV_MASS, sparts, mass);
--  list[3] = io_make_output_field("ParticleIDs", LONGLONG, 1, UNIT_CONV_NO_UNITS,
--                                 sparts, id);
--  list[4] = io_make_output_field("SmoothingLength", FLOAT, 1, UNIT_CONV_LENGTH,
--                                 sparts, h);
--}
--
--/**
-- * @brief Initialize the global properties of the stars scheme.
-- *
-- * By default, takes the values provided by the hydro.
-- *
-- * @param sp The #stars_props.
-- * @param phys_const The physical constants in the internal unit system.
-- * @param us The internal unit system.
-- * @param params The parsed parameters.
-- * @param p The already read-in properties of the hydro scheme.
-- */
--INLINE static void stars_props_init(struct stars_props *sp,
--                                    const struct phys_const *phys_const,
--                                    const struct unit_system *us,
--                                    struct swift_params *params,
--                                    const struct hydro_props *p,
--				    const struct cosmology *cosmo) {
--
--  /* Kernel properties */
--  sp->eta_neighbours = parser_get_opt_param_float(
--      params, "Stars:resolution_eta", p->eta_neighbours);
--
--  /* Tolerance for the smoothing length Newton-Raphson scheme */
--  sp->h_tolerance =
--      parser_get_opt_param_float(params, "Stars:h_tolerance", p->h_tolerance);
--
--  /* Get derived properties */
--  sp->target_neighbours = pow_dimension(sp->eta_neighbours) * kernel_norm;
--  const float delta_eta = sp->eta_neighbours * (1.f + sp->h_tolerance);
--  sp->delta_neighbours =
--      (pow_dimension(delta_eta) - pow_dimension(sp->eta_neighbours)) *
--      kernel_norm;
--
--  /* Maximal smoothing length */
--  sp->h_max = parser_get_opt_param_float(params, "Stars:h_max", p->h_max);
--
--  /* Number of iterations to converge h */
--  sp->max_smoothing_iterations = parser_get_opt_param_int(
--      params, "Stars:max_ghost_iterations", p->max_smoothing_iterations);
--
--  /* Time integration properties */
--  const float max_volume_change =
--      parser_get_opt_param_float(params, "Stars:max_volume_change", -1);
--  if (max_volume_change == -1)
--    sp->log_max_h_change = p->log_max_h_change;
--  else
--    sp->log_max_h_change = logf(powf(max_volume_change, hydro_dimension_inv));
--
--  /* Set stellar evolution parameters */
--  sp->deltaT_desired = 3.16228e7;
--  // BE CAREFUL OF UNITS HERE...
--  sp->temp_to_u_factor = phys_const->const_boltzmann_k / (p->mu_ionised * (hydro_gamma_minus_one) * phys_const->const_proton_mass);
--  // Why is this one here? copied from EAGLE...
--  sp->SNIa_energy_fraction = 1.0e0;
--  sp->total_energy_SNe = 1.0e51/units_cgs_conversion_factor(us,UNIT_CONV_ENERGY);
--  
--  /* Set to 1 if using continuous heating, 0 for stochastic  */
--  sp->continuous_heating = parser_get_opt_param_int(params, "Stars:continuous_heating", 0);
--
--  // CHANGE NAME. CURRENTLY SAME AS IN EAGLE
--  sp->units_factor2 = sp->total_energy_SNe * cosmo->h;
--  sp->units_factor1 = sp->units_factor2 / sp->temp_to_u_factor;
--
--  sp->feedback_timescale = parser_get_param_float(params, "Stars:feedback_timescale");
--
--  // CHANGE THIS TO BE CONSISTENT WITH RAND MAX USED IN STAR FORMATION
--  sp->inv_rand_max = 1.0/RAND_MAX;
--
--}
--
--/**
-- * @brief Print the global properties of the stars scheme.
-- *
-- * @param sp The #stars_props.
-- */
--INLINE static void stars_props_print(const struct stars_props *sp) {
--
--  /* Now stars */
--  message("Stars kernel: %s with eta=%f (%.2f neighbours).", kernel_name,
--          sp->eta_neighbours, sp->target_neighbours);
--
--  message("Stars relative tolerance in h: %.5f (+/- %.4f neighbours).",
--          sp->h_tolerance, sp->delta_neighbours);
--
--  message(
--      "Stars integration: Max change of volume: %.2f "
--      "(max|dlog(h)/dt|=%f).",
--      pow_dimension(expf(sp->log_max_h_change)), sp->log_max_h_change);
--
--  if (sp->h_max != FLT_MAX)
--    message("Maximal smoothing length allowed: %.4f", sp->h_max);
--
--  message("Maximal iterations in ghost task set to %d",
--          sp->max_smoothing_iterations);
--}
--
--#if defined(HAVE_HDF5)
--INLINE static void stars_props_print_snapshot(hid_t h_grpstars,
--                                              const struct stars_props *sp) {
--
--  io_write_attribute_s(h_grpstars, "Kernel function", kernel_name);
--  io_write_attribute_f(h_grpstars, "Kernel target N_ngb",
--                       sp->target_neighbours);
--  io_write_attribute_f(h_grpstars, "Kernel delta N_ngb", sp->delta_neighbours);
--  io_write_attribute_f(h_grpstars, "Kernel eta", sp->eta_neighbours);
--  io_write_attribute_f(h_grpstars, "Smoothing length tolerance",
--                       sp->h_tolerance);
--  io_write_attribute_f(h_grpstars, "Maximal smoothing length", sp->h_max);
--  io_write_attribute_f(h_grpstars, "Volume log(max(delta h))",
--                       sp->log_max_h_change);
--  io_write_attribute_f(h_grpstars, "Volume max change time-step",
--                       pow_dimension(expf(sp->log_max_h_change)));
--  io_write_attribute_i(h_grpstars, "Max ghost iterations",
--                       sp->max_smoothing_iterations);
--}
--#endif
--
--/**
-- * @brief Write a #stars_props struct to the given FILE as a stream of bytes.
-- *
-- * @param p the struct
-- * @param stream the file stream
-- */
--INLINE static void stars_props_struct_dump(const struct stars_props *p,
--                                           FILE *stream) {
--  restart_write_blocks((void *)p, sizeof(struct stars_props), 1, stream,
--                       "starsprops", "stars props");
--}
--
--/**
-- * @brief Restore a stars_props struct from the given FILE as a stream of
-- * bytes.
-- *
-- * @param p the struct
-- * @param stream the file stream
-- */
--INLINE static void stars_props_struct_restore(const struct stars_props *p,
--                                              FILE *stream) {
--  restart_read_blocks((void *)p, sizeof(struct stars_props), 1, stream, NULL,
--                      "stars props");
--}
--
--#endif /* SWIFT_DEFAULT_STAR_IO_H */
-diff --git a/src/stars/const/stars_part.h b/src/stars/const/stars_part.h
-deleted file mode 100644
-index 942a73c..0000000
---- a/src/stars/const/stars_part.h
-+++ /dev/null
-@@ -1,160 +0,0 @@
--/*******************************************************************************
-- * This file is part of SWIFT.
-- * Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@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/>.
-- *
-- ******************************************************************************/
--#ifndef SWIFT_DEFAULT_STAR_PART_H
--#define SWIFT_DEFAULT_STAR_PART_H
--
--/* Some standard headers. */
--#include <stdlib.h>
--
--/* Read chemistry */
--#include "chemistry_struct.h" 
--
--/**
-- * @brief Particle fields for the star particles.
-- *
-- * All quantities related to gravity are stored in the associate #gpart.
-- */
--struct spart {
--
--  /*! Particle ID. */
--  long long id;
--
--  /*! Pointer to corresponding gravity part. */
--  struct gpart* gpart;
--
--  /*! Particle position. */
--  double x[3];
--
--  /* Offset between current position and position at last tree rebuild. */
--  float x_diff[3];
--
--  /* Offset between current position and position at last tree rebuild. */
--  float x_diff_sort[3];
--
--  /*! Particle velocity. */
--  float v[3];
--
--  /*! Star mass */
--  float mass;
--
--  /*! Initial star mass */
--  float mass_init;
--
--  /* Particle cutoff radius. */
--  float h;
--
--  /*! Particle time bin */
--  timebin_t time_bin;
--
--  struct {
--    /* Number of neighbours. */
--    float wcount;
--
--    /* Number of neighbours spatial derivative. */
--    float wcount_dh;
--
--  } density;
--
--  struct {
--    /* Mass of ejecta */
--    float mass;
--
--    /* Mass fractions of ejecta */
--    struct chemistry_part_data chemistry_data;
--  
--    float ejecta_specific_thermal_energy;
--
--    /* Number of type 1a SNe per unit mass */
--    float num_SNIa;
--
--    /* total mass of neighbouring gas particles (TODO: CHECK IF THIS IS DIFFERENT THAN ngb_mass DECLARED BELOW) */
--    float ngb_mass;
--
--  } to_distribute;
--
--  /* kernel normalisation factor (equivalent to metalweight_norm in eagle_enrich.c:811, IMPROVE COMMENT) */
--  float omega_normalisation_inv;
--
--  /* total mass of neighbouring gas particles */
--  float ngb_mass;
--
--#ifdef SWIFT_DEBUG_CHECKS
--
--  /* Time of the last drift */
--  integertime_t ti_drift;
--
--  /* Time of the last kick */
--  integertime_t ti_kick;
--
--#endif
--
--#ifdef DEBUG_INTERACTIONS_STARS
--  /*! List of interacting particles in the density SELF and PAIR */
--  long long ids_ngbs_density[MAX_NUM_OF_NEIGHBOURS_STARS];
--
--  /*! Number of interactions in the density SELF and PAIR */
--  int num_ngb_density;
--#endif
--
--} SWIFT_STRUCT_ALIGN;
--
--/**
-- * @brief Contains all the constants and parameters of the stars scheme
-- */
--struct stars_props {
--
--  /*! Resolution parameter */
--  float eta_neighbours;
--
--  /*! Target weightd number of neighbours (for info only)*/
--  float target_neighbours;
--
--  /*! Smoothing length tolerance */
--  float h_tolerance;
--
--  /*! Tolerance on neighbour number  (for info only)*/
--  float delta_neighbours;
--
--  /*! Maximal smoothing length */
--  float h_max;
--
--  /*! Maximal number of iterations to converge h */
--  int max_smoothing_iterations;
--
--  /*! Maximal change of h over one time-step */
--  float log_max_h_change;
--
--  int continuous_heating;
--
--  float SNIa_energy_fraction;
--  float deltaT_desired;
--  float temp_to_u_factor;
--  float total_energy_SNe;
--
--  // Conversion factors copied from EAGLE. CHANGE NAME TO BE MORE DESCRIPTIVE
--  float units_factor1, units_factor2;
--
--  float feedback_timescale;
--
--  // CHANGE THIS TO BE CONSISTENT WITH RAND MAX USED IN STAR FORMATION
--  double inv_rand_max;
--
--};
--
--#endif /* SWIFT_DEFAULT_STAR_PART_H */
diff --git a/out_stars_default.txt b/out_stars_default.txt
deleted file mode 100644
index 5d7f579016988caa3262e2dd5c7d72c390472426..0000000000000000000000000000000000000000
--- a/out_stars_default.txt
+++ /dev/null
@@ -1,152 +0,0 @@
-diff --git a/src/stars/Default/stars.h b/src/stars/Default/stars.h
-index 8632057..a8faf43 100644
---- a/src/stars/Default/stars.h
-+++ b/src/stars/Default/stars.h
-@@ -56,8 +56,6 @@ __attribute__((always_inline)) INLINE static void stars_init_spart(
-     struct spart* sp) {
- 
- #ifdef DEBUG_INTERACTIONS_STARS
--  for (int i = 0; i < MAX_NUM_OF_NEIGHBOURS_STARS; ++i)
--    sp->ids_ngbs_density[i] = -1;
-   sp->num_ngb_density = 0;
- #endif
- 
-@@ -145,6 +143,20 @@ __attribute__((always_inline)) INLINE static void stars_spart_has_no_neighbours(
-  */
- __attribute__((always_inline)) INLINE static void stars_evolve_spart(
-     struct spart* restrict sp, const struct stars_props* stars_properties,
--    const struct cosmology* cosmo, float dt) {}
-+    const struct cosmology* cosmo) {}
- 
-+/**
-+ * @brief Reset acceleration fields of a particle
-+ *
-+ * This is the equivalent of hydro_reset_acceleration.
-+ * We do not compute the acceleration on star, therefore no need to use it.
-+ *
-+ * @param p The particle to act upon
-+ */
-+__attribute__((always_inline)) INLINE static void stars_reset_acceleration(
-+    struct spart* restrict p) {
-+#ifdef DEBUG_INTERACTIONS_STARS
-+  p->num_ngb_force = 0;
-+#endif
-+}
- #endif /* SWIFT_DEFAULT_STARS_H */
-diff --git a/src/stars/Default/stars_iact.h b/src/stars/Default/stars_iact.h
-index a0da318..7e36d4c 100644
---- a/src/stars/Default/stars_iact.h
-+++ b/src/stars/Default/stars_iact.h
-@@ -14,9 +14,7 @@ __attribute__((always_inline)) INLINE static void
- runner_iact_nonsym_stars_density(float r2, const float *dx, float hi, float hj,
-                                  struct spart *restrict si,
-                                  const struct part *restrict pj, float a,
--                                 float H, const struct cosmology *restrict cosmo,
--                                 const struct stars_props *restrict stars_properties,
--                                 struct xpart *restrict xp) {
-+                                 float H) {
- 
-   float wi, wi_dx;
- 
-@@ -35,8 +33,6 @@ runner_iact_nonsym_stars_density(float r2, const float *dx, float hi, float hj,
- 
- #ifdef DEBUG_INTERACTIONS_STARS
-   /* Update ngb counters */
--  if (si->num_ngb_density < MAX_NUM_OF_NEIGHBOURS_STARS)
--    si->ids_ngbs_density[si->num_ngb_density] = pj->id;
-   ++si->num_ngb_density;
- #endif
- }
-@@ -56,9 +52,10 @@ runner_iact_nonsym_stars_density(float r2, const float *dx, float hi, float hj,
- __attribute__((always_inline)) INLINE static void
- runner_iact_nonsym_stars_feedback(float r2, const float *dx, float hi, float hj,
-                                   struct spart *restrict si,
--                                  struct part *restrict pj, float a, float H,
--                                  const struct cosmology *restrict cosmo,
--                                  const struct stars_props *restrict stars_properties,
--                                  struct xpart *restrict xp) {
--  si->to_distribute.mass = 2;
-+                                  struct part *restrict pj, float a, float H) {
-+
-+#ifdef DEBUG_INTERACTIONS_STARS
-+  /* Update ngb counters */
-+  ++si->num_ngb_force;
-+#endif
- }
-diff --git a/src/stars/Default/stars_io.h b/src/stars/Default/stars_io.h
-index 8b0c75f..fbed04a 100644
---- a/src/stars/Default/stars_io.h
-+++ b/src/stars/Default/stars_io.h
-@@ -60,10 +60,10 @@ INLINE static void stars_write_particles(const struct spart *sparts,
-                                          struct io_props *list,
-                                          int *num_fields) {
- 
--  /* Say how much we want to read */
-+  /* Say how much we want to write */
-   *num_fields = 5;
- 
--  /* List what we want to read */
-+  /* List what we want to write */
-   list[0] = io_make_output_field("Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH,
-                                  sparts, x);
-   list[1] =
-@@ -74,6 +74,17 @@ INLINE static void stars_write_particles(const struct spart *sparts,
-                                  sparts, id);
-   list[4] = io_make_output_field("SmoothingLength", FLOAT, 1, UNIT_CONV_LENGTH,
-                                  sparts, h);
-+
-+#ifdef DEBUG_INTERACTIONS_STARS
-+
-+  list += *num_fields;
-+  *num_fields += 2;
-+
-+  list[0] = io_make_output_field("Num_ngb_density", INT, 1, UNIT_CONV_NO_UNITS,
-+                                 sparts, num_ngb_density);
-+  list[1] = io_make_output_field("Num_ngb_force", INT, 1, UNIT_CONV_NO_UNITS,
-+                                 sparts, num_ngb_force);
-+#endif
- }
- 
- /**
-@@ -91,8 +102,7 @@ INLINE static void stars_props_init(struct stars_props *sp,
-                                     const struct phys_const *phys_const,
-                                     const struct unit_system *us,
-                                     struct swift_params *params,
--                                    const struct hydro_props *p,
--				    const struct cosmology *cosmo) {
-+                                    const struct hydro_props *p) {
- 
-   /* Kernel properties */
-   sp->eta_neighbours = parser_get_opt_param_float(
-diff --git a/src/stars/Default/stars_part.h b/src/stars/Default/stars_part.h
-index 089f9c2..c487528 100644
---- a/src/stars/Default/stars_part.h
-+++ b/src/stars/Default/stars_part.h
-@@ -65,12 +65,6 @@ struct spart {
- 
-   } density;
- 
--  struct {
--    /* Mass of ejecta */
--    float mass;
--
--  } to_distribute;
--
- #ifdef SWIFT_DEBUG_CHECKS
- 
-   /* Time of the last drift */
-@@ -82,11 +76,11 @@ struct spart {
- #endif
- 
- #ifdef DEBUG_INTERACTIONS_STARS
--  /*! List of interacting particles in the density SELF and PAIR */
--  long long ids_ngbs_density[MAX_NUM_OF_NEIGHBOURS_STARS];
--
-   /*! Number of interactions in the density SELF and PAIR */
-   int num_ngb_density;
-+
-+  /*! Number of interactions in the force SELF and PAIR */
-+  int num_ngb_force;
- #endif
- 
- } SWIFT_STRUCT_ALIGN;
diff --git a/out_stars_runner.txt b/out_stars_runner.txt
deleted file mode 100644
index 452857dca3bd5e5bb15b2890177be1c7d8a7696d..0000000000000000000000000000000000000000
--- a/out_stars_runner.txt
+++ /dev/null
@@ -1,238 +0,0 @@
-diff --git a/src/runner_doiact_stars.h b/src/runner_doiact_stars.h
-index 80996e1..9cf9ea1 100644
---- a/src/runner_doiact_stars.h
-+++ b/src/runner_doiact_stars.h
-@@ -74,9 +74,13 @@
-  * @param timer 1 if the time is to be recorded.
-  */
- void DOSELF1_STARS(struct runner *r, struct cell *c, int timer) {
-+
-+#ifdef SWIFT_DEBUG_CHECKS
-+  if (c->nodeID != engine_rank) error("Should be run on a different node");
-+#endif
-+
-   const struct engine *e = r->e;
-   const struct cosmology *cosmo = e->cosmology;
--  const struct stars_props *stars_properties = e->stars_properties;
- 
-   /* Anything to do here? */
-   if (!cell_is_active_stars(c, e)) return;
-@@ -90,7 +94,6 @@ void DOSELF1_STARS(struct runner *r, struct cell *c, int timer) {
-   const int count = c->hydro.count;
-   struct spart *restrict sparts = c->stars.parts;
-   struct part *restrict parts = c->hydro.parts;
--  struct xpart *restrict xparts = c->hydro.xparts;
- 
-   /* Loop over the sparts in ci. */
-   for (int sid = 0; sid < scount; sid++) {
-@@ -108,7 +111,6 @@ void DOSELF1_STARS(struct runner *r, struct cell *c, int timer) {
- 
-       /* Get a pointer to the jth particle. */
-       struct part *restrict pj = &parts[pjd];
--      struct xpart *restrict xpj = &xparts[pjd];
-       const float hj = pj->h;
- 
-       /* Compute the pairwise distance. */
-@@ -125,7 +127,7 @@ void DOSELF1_STARS(struct runner *r, struct cell *c, int timer) {
- #endif
- 
-       if (r2 > 0.f && r2 < hig2) {
--        IACT_STARS(r2, dx, hi, hj, si, pj, a, H, cosmo, stars_properties, xpj);
-+        IACT_STARS(r2, dx, hi, hj, si, pj, a, H);
-       }
-     } /* loop over the parts in ci. */
-   }   /* loop over the sparts in ci. */
-@@ -141,9 +143,14 @@ void DOSELF1_STARS(struct runner *r, struct cell *c, int timer) {
- void DO_NONSYM_PAIR1_STARS(struct runner *r, struct cell *restrict ci,
-                            struct cell *restrict cj) {
- 
-+#ifdef SWIFT_DEBUG_CHECKS
-+#ifdef ONLY_LOCAL
-+  if (ci->nodeID != engine_rank) error("Should be run on a different node");
-+#endif
-+#endif
-+
-   const struct engine *e = r->e;
-   const struct cosmology *cosmo = e->cosmology;
--  const struct stars_props *stars_properties = e->stars_properties;
- 
-   /* Anything to do here? */
-   if (!cell_is_active_stars(ci, e)) return;
-@@ -152,7 +159,6 @@ void DO_NONSYM_PAIR1_STARS(struct runner *r, struct cell *restrict ci,
-   const int count_j = cj->hydro.count;
-   struct spart *restrict sparts_i = ci->stars.parts;
-   struct part *restrict parts_j = cj->hydro.parts;
--  struct xpart *restrict xparts_j = cj->hydro.xparts;
- 
-   /* Cosmological terms */
-   const float a = cosmo->a;
-@@ -183,7 +189,6 @@ void DO_NONSYM_PAIR1_STARS(struct runner *r, struct cell *restrict ci,
- 
-       /* Get a pointer to the jth particle. */
-       struct part *restrict pj = &parts_j[pjd];
--      struct xpart *restrict xpj = &xparts_j[pjd];
-       const float hj = pj->h;
- 
-       /* Compute the pairwise distance. */
-@@ -199,7 +204,7 @@ void DO_NONSYM_PAIR1_STARS(struct runner *r, struct cell *restrict ci,
-         error("Particle pj not drifted to current time");
- #endif
- 
--      if (r2 < hig2) IACT_STARS(r2, dx, hi, hj, si, pj, a, H, cosmo, stars_properties, xpj);
-+      if (r2 < hig2) IACT_STARS(r2, dx, hi, hj, si, pj, a, H);
- 
-     } /* loop over the parts in cj. */
-   }   /* loop over the parts in ci. */
-@@ -208,9 +213,16 @@ void DO_NONSYM_PAIR1_STARS(struct runner *r, struct cell *restrict ci,
- void DOPAIR1_STARS(struct runner *r, struct cell *restrict ci,
-                    struct cell *restrict cj, int timer) {
- 
--  if (ci->stars.count != 0 && cj->hydro.count != 0)
-+#ifdef ONLY_LOCAL
-+  const int ci_local = ci->nodeID == engine_rank;
-+  const int cj_local = cj->nodeID == engine_rank;
-+#else
-+  const int ci_local = 1;
-+  const int cj_local = 1;
-+#endif
-+  if (ci_local && ci->stars.count != 0 && cj->hydro.count != 0)
-     DO_NONSYM_PAIR1_STARS(r, ci, cj);
--  if (cj->stars.count != 0 && ci->hydro.count != 0)
-+  if (cj_local && cj->stars.count != 0 && ci->hydro.count != 0)
-     DO_NONSYM_PAIR1_STARS(r, cj, ci);
- }
- 
-@@ -233,13 +245,14 @@ void DOPAIR1_SUBSET_STARS(struct runner *r, struct cell *restrict ci,
-                           int scount, struct cell *restrict cj,
-                           const double *shift) {
- 
-+#ifdef SWIFT_DEBUG_CHECKS
-+  if (ci->nodeID != engine_rank) error("Should be run on a different node");
-+#endif
-   const struct engine *e = r->e;
-   const struct cosmology *cosmo = e->cosmology;
--  const struct stars_props *stars_properties = e->stars_properties;
- 
-   const int count_j = cj->hydro.count;
-   struct part *restrict parts_j = cj->hydro.parts;
--  struct xpart *restrict xparts_j = cj->hydro.xparts;
- 
-   /* Cosmological terms */
-   const float a = cosmo->a;
-@@ -265,7 +278,6 @@ void DOPAIR1_SUBSET_STARS(struct runner *r, struct cell *restrict ci,
- 
-       /* Get a pointer to the jth particle. */
-       struct part *restrict pj = &parts_j[pjd];
--      struct xpart *restrict xpj = &xparts_j[pjd];
- 
-       /* Compute the pairwise distance. */
-       float r2 = 0.0f;
-@@ -282,7 +294,7 @@ void DOPAIR1_SUBSET_STARS(struct runner *r, struct cell *restrict ci,
- #endif
-       /* Hit or miss? */
-       if (r2 < hig2) {
--        IACT_STARS(r2, dx, hi, pj->h, spi, pj, a, H, cosmo, stars_properties, xpj);
-+        IACT_STARS(r2, dx, hi, pj->h, spi, pj, a, H);
-       }
-     } /* loop over the parts in cj. */
-   }   /* loop over the parts in ci. */
-@@ -302,9 +314,12 @@ void DOSELF1_SUBSET_STARS(struct runner *r, struct cell *restrict ci,
-                           struct spart *restrict sparts, int *restrict ind,
-                           int scount) {
- 
-+#ifdef SWIFT_DEBUG_CHECKS
-+  if (ci->nodeID != engine_rank) error("Should be run on a different node");
-+#endif
-+
-   const struct engine *e = r->e;
-   const struct cosmology *cosmo = e->cosmology;
--  const struct stars_props *stars_properties = e->stars_properties;
- 
-   /* Cosmological terms */
-   const float a = cosmo->a;
-@@ -312,7 +327,6 @@ void DOSELF1_SUBSET_STARS(struct runner *r, struct cell *restrict ci,
- 
-   const int count_i = ci->hydro.count;
-   struct part *restrict parts_j = ci->hydro.parts;
--  struct xpart *restrict xparts_j = ci->hydro.xparts;
- 
-   /* Loop over the parts in ci. */
-   for (int spid = 0; spid < scount; spid++) {
-@@ -335,7 +349,6 @@ void DOSELF1_SUBSET_STARS(struct runner *r, struct cell *restrict ci,
- 
-       /* Get a pointer to the jth particle. */
-       struct part *restrict pj = &parts_j[pjd];
--      struct xpart *restrict xpj = &xparts_j[pjd];
-       const float hj = pj->h;
- 
-       /* Compute the pairwise distance. */
-@@ -353,7 +366,7 @@ void DOSELF1_SUBSET_STARS(struct runner *r, struct cell *restrict ci,
- 
-       /* Hit or miss? */
-       if (r2 > 0.f && r2 < hig2) {
--        IACT_STARS(r2, dx, hi, hj, spi, pj, a, H, cosmo, stars_properties, xpj);
-+        IACT_STARS(r2, dx, hi, hj, spi, pj, a, H);
-       }
-     } /* loop over the parts in cj. */
-   }   /* loop over the parts in ci. */
-@@ -1007,6 +1020,8 @@ void DOSELF1_BRANCH_STARS(struct runner *r, struct cell *c) {
-                                                                             \
-     for (int pjd = 0; pjd < cj->TYPE.count; pjd++) {                        \
-       const struct PART *p = &cj->TYPE.parts[sort_j[pjd].i];                \
-+      if (PART##_is_inhibited(p, e)) continue;                              \
-+                                                                            \
-       const float d = p->x[0] * runner_shift[sid][0] +                      \
-                       p->x[1] * runner_shift[sid][1] +                      \
-                       p->x[2] * runner_shift[sid][2];                       \
-@@ -1019,9 +1034,9 @@ void DOSELF1_BRANCH_STARS(struct runner *r, struct cell *c) {
-             "cj->nodeID=%d "                                                \
-             "ci->nodeID=%d d=%e sort_j[pjd].d=%e cj->" #TYPE                \
-             ".dx_max_sort=%e "                                              \
--            "cj->" #TYPE ".dx_max_sort_old=%e",                             \
-+            "cj->" #TYPE ".dx_max_sort_old=%e, %i",                         \
-             cj->nodeID, ci->nodeID, d, sort_j[pjd].d, cj->TYPE.dx_max_sort, \
--            cj->TYPE.dx_max_sort_old);                                      \
-+            cj->TYPE.dx_max_sort_old, cj->cellID);                          \
-     }                                                                       \
-   })
- 
-@@ -1040,8 +1055,17 @@ void DOPAIR1_BRANCH_STARS(struct runner *r, struct cell *ci, struct cell *cj) {
-   const struct engine *restrict e = r->e;
-   const int ci_active = cell_is_active_stars(ci, e);
-   const int cj_active = cell_is_active_stars(cj, e);
--  const int do_ci = (ci->stars.count != 0 && cj->hydro.count != 0 && ci_active);
--  const int do_cj = (cj->stars.count != 0 && ci->hydro.count != 0 && cj_active);
-+#ifdef ONLY_LOCAL
-+  const int ci_local = ci->nodeID == engine_rank;
-+  const int cj_local = cj->nodeID == engine_rank;
-+#else
-+  const int ci_local = 1;
-+  const int cj_local = 1;
-+#endif
-+  const int do_ci =
-+      (ci->stars.count != 0 && cj->hydro.count != 0 && ci_active && ci_local);
-+  const int do_cj =
-+      (cj->stars.count != 0 && ci->hydro.count != 0 && cj_active && cj_local);
- 
-   /* Anything to do here? */
-   if (!do_ci && !do_cj) return;
-@@ -1326,10 +1350,17 @@ void DOSUB_PAIR1_STARS(struct runner *r, struct cell *ci, struct cell *cj,
-   /* Otherwise, compute the pair directly. */
-   else {
- 
-+#ifdef ONLY_LOCAL
-+    const int ci_local = ci->nodeID == engine_rank;
-+    const int cj_local = cj->nodeID == engine_rank;
-+#else
-+    const int ci_local = 1;
-+    const int cj_local = 1;
-+#endif
-     const int do_ci = ci->stars.count != 0 && cj->hydro.count != 0 &&
--                      cell_is_active_stars(ci, e);
-+                      cell_is_active_stars(ci, e) && ci_local;
-     const int do_cj = cj->stars.count != 0 && ci->hydro.count != 0 &&
--                      cell_is_active_stars(cj, e);
-+                      cell_is_active_stars(cj, e) && cj_local;
- 
-     if (do_ci) {
-