diff --git a/examples/main.c b/examples/main.c
index 8b5200850438aebfbd39bb84d4f9b41602d34aa6..c1f8b99c857b7344e4c290c2e21e4f32d099dd2d 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -947,7 +947,19 @@ int main(int argc, char *argv[]) {
 
     if (with_hydro) {
 #ifdef NONE_SPH
-      error("Can't run with hydro when compiled without hydro model!");
+      error("Can't run with hydro when compiled without a hydro model!");
+#endif
+    }
+    if (with_stars) {
+#ifdef STARS_NONE
+      error("Can't run with stars when compiled without a stellar model!");
+#endif
+    }
+    if (with_black_holes) {
+#ifdef BLACK_HOLES_NONE
+      error(
+          "Can't run with black holes when compiled without a black hole "
+          "model!");
 #endif
     }
 
diff --git a/src/cell.c b/src/cell.c
index f1ee31cb7676de5a68f385d8c5b3eb3688a02a41..018db5f60f8a5b7efb36647930ca505573b9b9e8 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -515,7 +515,9 @@ void cell_sanitize(struct cell *c, int treated) {
   struct part *parts = c->hydro.parts;
   struct spart *sparts = c->stars.parts;
   float h_max = 0.f;
+  float h_max_active = 0.f;
   float stars_h_max = 0.f;
+  float stars_h_max_active = 0.f;
 
   /* Treat cells will <1000 particles */
   if (count < 1000 && !treated) {
@@ -542,19 +544,28 @@ void cell_sanitize(struct cell *c, int treated) {
 
         /* And collect */
         h_max = max(h_max, c->progeny[k]->hydro.h_max);
+        h_max_active = max(h_max_active, c->progeny[k]->hydro.h_max_active);
         stars_h_max = max(stars_h_max, c->progeny[k]->stars.h_max);
+        stars_h_max_active =
+            max(stars_h_max_active, c->progeny[k]->stars.h_max_active);
       }
     }
   } else {
-    /* Get the new value of h_max */
+    /* Get the new value of h_max (note all particles are active) */
     for (int i = 0; i < count; ++i) h_max = max(h_max, parts[i].h);
+    for (int i = 0; i < count; ++i)
+      h_max_active = max(h_max_active, parts[i].h);
     for (int i = 0; i < scount; ++i)
       stars_h_max = max(stars_h_max, sparts[i].h);
+    for (int i = 0; i < scount; ++i)
+      stars_h_max_active = max(stars_h_max_active, sparts[i].h);
   }
 
   /* Record the change */
   c->hydro.h_max = h_max;
+  c->hydro.h_max_active = h_max_active;
   c->stars.h_max = stars_h_max;
+  c->stars.h_max_active = stars_h_max_active;
 }
 
 /**
diff --git a/src/cell.h b/src/cell.h
index aaeb263164e595555026b1c9be1e78c361a12e0f..08fbf13139916b1f57647377fae3ceef5e971975 100644
--- a/src/cell.h
+++ b/src/cell.h
@@ -1101,11 +1101,15 @@ __attribute__((always_inline)) INLINE static void cell_malloc_hydro_sorts(
 __attribute__((always_inline)) INLINE static void cell_free_hydro_sorts(
     struct cell *c) {
 
+#ifdef NONE_SPH
+  /* Nothing to do as we have no particles and hence no sorts */
+#else
   if (c->hydro.sort != NULL) {
     swift_free("hydro.sort", c->hydro.sort);
     c->hydro.sort = NULL;
     c->hydro.sort_allocated = 0;
   }
+#endif
 }
 
 /**
@@ -1214,11 +1218,15 @@ __attribute__((always_inline)) INLINE static void cell_malloc_stars_sorts(
 __attribute__((always_inline)) INLINE static void cell_free_stars_sorts(
     struct cell *c) {
 
+#ifdef STARS_NONE
+  /* Nothing to do as we have no particles and hence no sorts */
+#else
   if (c->stars.sort != NULL) {
     swift_free("stars.sort", c->stars.sort);
     c->stars.sort = NULL;
     c->stars.sort_allocated = 0;
   }
+#endif
 }
 
 /**
diff --git a/src/cell_black_holes.h b/src/cell_black_holes.h
index 1ee1545db8f1df9749080a79fa7e1ab0bb812358..ead886a3752a032f462414345a7189e2799e2cb7 100644
--- a/src/cell_black_holes.h
+++ b/src/cell_black_holes.h
@@ -33,50 +33,75 @@
  */
 struct cell_black_holes {
 
-  /*! Pointer to the #bpart data. */
-  struct bpart *parts;
+  /* If we are not using BHs, compact as much of the unecessary variables
+     into an anonymous union to save memory in the cell structure. */
+#ifdef BLACK_HOLES_NONE
+  union {
+#endif
 
-  /*! The drift task for bparts */
-  struct task *drift;
+    /*! Pointer to the #bpart data. */
+    struct bpart *parts;
 
-  /*! Implicit tasks marking the entry of the BH physics block of tasks */
-  struct task *black_holes_in;
+    /*! The drift task for bparts */
+    struct task *drift;
 
-  /*! Implicit tasks marking the exit of the BH physics block of tasks */
-  struct task *black_holes_out;
+    /*! Implicit tasks marking the entry of the BH physics block of tasks */
+    struct task *black_holes_in;
 
-  /*! The black hole ghost task itself */
-  struct task *density_ghost;
+    /*! Implicit tasks marking the exit of the BH physics block of tasks */
+    struct task *black_holes_out;
 
-  /*! The ghost tasks related to BH swallowing */
-  struct task *swallow_ghost_0;
-  struct task *swallow_ghost_1;
-  struct task *swallow_ghost_2;
+    /*! The black hole ghost task itself */
+    struct task *density_ghost;
 
-  /*! Linked list of the tasks computing this cell's BH density. */
-  struct link *density;
+    /*! The ghost tasks related to BH swallowing */
+    struct task *swallow_ghost_0;
+    struct task *swallow_ghost_1;
+    struct task *swallow_ghost_2;
 
-  /*! Linked list of the tasks computing this cell's BH swallowing and
-   * merging. */
-  struct link *swallow;
+    /*! Linked list of the tasks computing this cell's BH density. */
+    struct link *density;
 
-  /*! Linked list of the tasks processing the particles to swallow */
-  struct link *do_gas_swallow;
+    /*! Linked list of the tasks computing this cell's BH swallowing and
+     * merging. */
+    struct link *swallow;
 
-  /*! Linked list of the tasks processing the particles to swallow */
-  struct link *do_bh_swallow;
+    /*! Linked list of the tasks processing the particles to swallow */
+    struct link *do_gas_swallow;
 
-  /*! Linked list of the tasks computing this cell's BH feedback. */
-  struct link *feedback;
+    /*! Linked list of the tasks processing the particles to swallow */
+    struct link *do_bh_swallow;
 
-  /*! Last (integer) time the cell's bpart were drifted forward in time. */
-  integertime_t ti_old_part;
+    /*! Linked list of the tasks computing this cell's BH feedback. */
+    struct link *feedback;
 
-  /*! Minimum end of (integer) time step in this cell for black tasks. */
-  integertime_t ti_end_min;
+    /*! Last (integer) time the cell's bpart were drifted forward in time. */
+    integertime_t ti_old_part;
+
+    /*! Maximum end of (integer) time step in this cell for black hole tasks. */
+    integertime_t ti_end_max;
+
+    /*! Nr of #bpart this cell can hold after addition of new #bpart. */
+    int count_total;
+
+    /*! Max smoothing length of active particles in this cell. */
+    float h_max_active;
+
+    /*! Values of h_max before the drifts, used for sub-cell tasks. */
+    float h_max_old;
+
+    /*! Maximum part movement in this cell since last construction. */
+    float dx_max_part;
 
-  /*! Maximum end of (integer) time step in this cell for black hole tasks. */
-  integertime_t ti_end_max;
+    /*! Values of dx_max before the drifts, used for sub-cell tasks. */
+    float dx_max_part_old;
+
+#ifdef BLACK_HOLES_NONE
+  };
+#endif
+
+  /*! Maximum end of (integer) time step in this cell for black tasks. */
+  integertime_t ti_end_min;
 
   /*! Maximum beginning of (integer) time step in this cell for black hole
    * tasks. */
@@ -85,29 +110,17 @@ struct cell_black_holes {
   /*! Spin lock for various uses (#bpart case). */
   swift_lock_type lock;
 
-  /*! Number of #bpart updated in this cell. */
-  int updated;
+  /*! Max smoothing length in this cell. */
+  float h_max;
 
   /*! Is the #bpart data of this cell being used in a sub-cell? */
   int hold;
 
+  /*! Number of #bpart updated in this cell. */
+  int updated;
+
   /*! Nr of #bpart in this cell. */
   int count;
-
-  /*! Nr of #bpart this cell can hold after addition of new #bpart. */
-  int count_total;
-
-  /*! Max smoothing length in this cell. */
-  float h_max;
-
-  /*! Values of h_max before the drifts, used for sub-cell tasks. */
-  float h_max_old;
-
-  /*! Maximum part movement in this cell since last construction. */
-  float dx_max_part;
-
-  /*! Values of dx_max before the drifts, used for sub-cell tasks. */
-  float dx_max_part_old;
 };
 
 #endif /* SWIFT_CELL_BLACK_HOLES_H */
diff --git a/src/cell_convert_part.c b/src/cell_convert_part.c
index edd199cf06e2654d5ef69612b2975669615ef09d..b7d27022a694aa6d0af3c68dd4e5a66d37580a98 100644
--- a/src/cell_convert_part.c
+++ b/src/cell_convert_part.c
@@ -860,7 +860,7 @@ struct spart *cell_convert_part_to_spart(struct engine *e, struct cell *c,
 #endif
 
   /* Set a smoothing length */
-  sp->h = max(c->stars.h_max, c->hydro.h_max);
+  sp->h = p->h;
 
   /* Here comes the Sun! */
   return sp;
diff --git a/src/cell_drift.c b/src/cell_drift.c
index 291316b17b67364942ee0c67a9c1f6b0558d43ba..c176c9a6b321ab31e3d06df4c7ccd7f135c330b1 100644
--- a/src/cell_drift.c
+++ b/src/cell_drift.c
@@ -57,6 +57,7 @@ void cell_drift_part(struct cell *c, const struct engine *e, int force) {
   float dx_max = 0.f, dx2_max = 0.f;
   float dx_max_sort = 0.0f, dx2_max_sort = 0.f;
   float cell_h_max = 0.f;
+  float cell_h_max_active = 0.f;
 
   /* Drift irrespective of cell flags? */
   force = (force || cell_get_flag(c, cell_flag_do_hydro_drift));
@@ -97,11 +98,13 @@ void cell_drift_part(struct cell *c, const struct engine *e, int force) {
         dx_max = max(dx_max, cp->hydro.dx_max_part);
         dx_max_sort = max(dx_max_sort, cp->hydro.dx_max_sort);
         cell_h_max = max(cell_h_max, cp->hydro.h_max);
+        cell_h_max_active = max(cell_h_max_active, cp->hydro.h_max_active);
       }
     }
 
     /* Store the values */
     c->hydro.h_max = cell_h_max;
+    c->hydro.h_max_active = cell_h_max_active;
     c->hydro.dx_max_part = dx_max;
     c->hydro.dx_max_sort = dx_max_sort;
 
@@ -233,6 +236,9 @@ void cell_drift_part(struct cell *c, const struct engine *e, int force) {
                            e->cooling_func, e->time);
         rt_init_part(p);
         rt_reset_part(p);
+
+        /* Update the maximal active smoothing length in the cell */
+        cell_h_max_active = max(cell_h_max_active, p->h);
       }
     }
 
@@ -242,6 +248,7 @@ void cell_drift_part(struct cell *c, const struct engine *e, int force) {
 
     /* Store the values */
     c->hydro.h_max = cell_h_max;
+    c->hydro.h_max_active = cell_h_max_active;
     c->hydro.dx_max_part = dx_max;
     c->hydro.dx_max_sort = dx_max_sort;
 
@@ -414,6 +421,7 @@ void cell_drift_spart(struct cell *c, const struct engine *e, int force) {
   float dx_max = 0.f, dx2_max = 0.f;
   float dx_max_sort = 0.0f, dx2_max_sort = 0.f;
   float cell_h_max = 0.f;
+  float cell_h_max_active = 0.f;
 
   /* Drift irrespective of cell flags? */
   force = (force || cell_get_flag(c, cell_flag_do_stars_drift));
@@ -454,11 +462,13 @@ void cell_drift_spart(struct cell *c, const struct engine *e, int force) {
         dx_max = max(dx_max, cp->stars.dx_max_part);
         dx_max_sort = max(dx_max_sort, cp->stars.dx_max_sort);
         cell_h_max = max(cell_h_max, cp->stars.h_max);
+        cell_h_max_active = max(cell_h_max_active, cp->stars.h_max_active);
       }
     }
 
     /* Store the values */
     c->stars.h_max = cell_h_max;
+    c->stars.h_max_active = cell_h_max_active;
     c->stars.dx_max_part = dx_max;
     c->stars.dx_max_sort = dx_max_sort;
 
@@ -553,6 +563,9 @@ void cell_drift_spart(struct cell *c, const struct engine *e, int force) {
         stars_init_spart(sp);
         feedback_init_spart(sp);
         rt_init_spart(sp);
+
+        /* Update the maximal active smoothing length in the cell */
+        cell_h_max_active = max(cell_h_max_active, sp->h);
       }
     }
 
@@ -562,6 +575,7 @@ void cell_drift_spart(struct cell *c, const struct engine *e, int force) {
 
     /* Store the values */
     c->stars.h_max = cell_h_max;
+    c->stars.h_max_active = cell_h_max_active;
     c->stars.dx_max_part = dx_max;
     c->stars.dx_max_sort = dx_max_sort;
 
@@ -593,6 +607,7 @@ void cell_drift_bpart(struct cell *c, const struct engine *e, int force) {
 
   float dx_max = 0.f, dx2_max = 0.f;
   float cell_h_max = 0.f;
+  float cell_h_max_active = 0.f;
 
   /* Drift irrespective of cell flags? */
   force = (force || cell_get_flag(c, cell_flag_do_bh_drift));
@@ -633,11 +648,14 @@ void cell_drift_bpart(struct cell *c, const struct engine *e, int force) {
         /* Update */
         dx_max = max(dx_max, cp->black_holes.dx_max_part);
         cell_h_max = max(cell_h_max, cp->black_holes.h_max);
+        cell_h_max_active =
+            max(cell_h_max_active, cp->black_holes.h_max_active);
       }
     }
 
     /* Store the values */
     c->black_holes.h_max = cell_h_max;
+    c->black_holes.h_max_active = cell_h_max_active;
     c->black_holes.dx_max_part = dx_max;
 
     /* Update the time of the last drift */
@@ -726,6 +744,9 @@ void cell_drift_bpart(struct cell *c, const struct engine *e, int force) {
       /* Get ready for a density calculation */
       if (bpart_is_active(bp, e)) {
         black_holes_init_bpart(bp);
+
+        /* Update the maximal active smoothing length in the cell */
+        cell_h_max_active = max(cell_h_max_active, bp->h);
       }
     }
 
@@ -734,6 +755,7 @@ void cell_drift_bpart(struct cell *c, const struct engine *e, int force) {
 
     /* Store the values */
     c->black_holes.h_max = cell_h_max;
+    c->black_holes.h_max_active = cell_h_max_active;
     c->black_holes.dx_max_part = dx_max;
 
     /* Update the time of the last drift */
@@ -762,6 +784,7 @@ void cell_drift_sink(struct cell *c, const struct engine *e, int force) {
 
   float dx_max = 0.f, dx2_max = 0.f;
   float cell_r_max = 0.f;
+  float cell_r_max_active = 0.f;
 
   /* Drift irrespective of cell flags? */
   force = (force || cell_get_flag(c, cell_flag_do_sink_drift));
@@ -802,11 +825,13 @@ void cell_drift_sink(struct cell *c, const struct engine *e, int force) {
         /* Update */
         dx_max = max(dx_max, cp->sinks.dx_max_part);
         cell_r_max = max(cell_r_max, cp->sinks.r_cut_max);
+        cell_r_max_active = max(cell_r_max_active, cp->sinks.r_cut_max_active);
       }
     }
 
     /* Store the values */
     c->sinks.r_cut_max = cell_r_max;
+    c->sinks.r_cut_max_active = cell_r_max_active;
     c->sinks.dx_max_part = dx_max;
 
     /* Update the time of the last drift */
@@ -891,6 +916,8 @@ void cell_drift_sink(struct cell *c, const struct engine *e, int force) {
       /* Get ready for a density calculation */
       if (sink_is_active(sink, e)) {
         sink_init_sink(sink);
+
+        cell_r_max_active = max(cell_r_max_active, sink->r_cut);
       }
     }
 
@@ -899,6 +926,7 @@ void cell_drift_sink(struct cell *c, const struct engine *e, int force) {
 
     /* Store the values */
     c->sinks.r_cut_max = cell_r_max;
+    c->sinks.r_cut_max_active = cell_r_max_active;
     c->sinks.dx_max_part = dx_max;
 
     /* Update the time of the last drift */
diff --git a/src/cell_hydro.h b/src/cell_hydro.h
index 6a89c5d91d312bbf043887f2085127984dc54128..de158cb25923da549379e5f1176f925895c4cb5d 100644
--- a/src/cell_hydro.h
+++ b/src/cell_hydro.h
@@ -33,161 +33,174 @@
  */
 struct cell_hydro {
 
-  /*! Pointer to the #part data. */
-  struct part *parts;
+  /* If we are not using hydro, compact as much of the unecessary variables
+     into an anonymous union to save memory in the cell structure. */
+#ifdef NONE_SPH
+  union {
+#endif
 
-  /*! Pointer to the #xpart data. */
-  struct xpart *xparts;
+    /*! Pointer to the #part data. */
+    struct part *parts;
 
-  /*! Pointer for the sorted indices. */
-  struct sort_entry *sort;
+    /*! Pointer to the #xpart data. */
+    struct xpart *xparts;
 
-  /*! Super cell, i.e. the highest-level parent cell that has a hydro
-   * pair/self tasks */
-  struct cell *super;
+    /*! Pointer for the sorted indices. */
+    struct sort_entry *sort;
 
-  /*! The task computing this cell's sorts. */
-  struct task *sorts;
+    /*! Super cell, i.e. the highest-level parent cell that has a hydro
+     * pair/self tasks */
+    struct cell *super;
 
-  /*! The drift task for parts */
-  struct task *drift;
+    /*! The task computing this cell's sorts. */
+    struct task *sorts;
 
-  /*! Linked list of the tasks computing this cell's hydro density. */
-  struct link *density;
+    /*! The drift task for parts */
+    struct task *drift;
 
-  /* Linked list of the tasks computing this cell's hydro gradients. */
-  struct link *gradient;
+    /*! Linked list of the tasks computing this cell's hydro density. */
+    struct link *density;
 
-  /*! Linked list of the tasks computing this cell's hydro forces. */
-  struct link *force;
+    /* Linked list of the tasks computing this cell's hydro gradients. */
+    struct link *gradient;
 
-  /*! Linked list of the tasks computing this cell's limiter. */
-  struct link *limiter;
+    /*! Linked list of the tasks computing this cell's hydro forces. */
+    struct link *force;
 
-  /*! Dependency implicit task for the ghost  (in->ghost->out)*/
-  struct task *ghost_in;
+    /*! Linked list of the tasks computing this cell's limiter. */
+    struct link *limiter;
 
-  /*! Dependency implicit task for the ghost  (in->ghost->out)*/
-  struct task *ghost_out;
+    /*! Dependency implicit task for the ghost  (in->ghost->out)*/
+    struct task *ghost_in;
 
-  /*! The ghost task itself */
-  struct task *ghost;
+    /*! Dependency implicit task for the ghost  (in->ghost->out)*/
+    struct task *ghost_out;
 
-  /*! The extra ghost task for complex hydro schemes */
-  struct task *extra_ghost;
+    /*! The ghost task itself */
+    struct task *ghost;
 
-  /*! The task to end the force calculation */
-  struct task *end_force;
+    /*! The extra ghost task for complex hydro schemes */
+    struct task *extra_ghost;
 
-  /*! Dependency implicit task for cooling (in->cooling->out) */
-  struct task *cooling_in;
+    /*! The task to end the force calculation */
+    struct task *end_force;
 
-  /*! Dependency implicit task for cooling (in->cooling->out) */
-  struct task *cooling_out;
+    /*! Dependency implicit task for cooling (in->cooling->out) */
+    struct task *cooling_in;
 
-  /*! Task for cooling */
-  struct task *cooling;
+    /*! Dependency implicit task for cooling (in->cooling->out) */
+    struct task *cooling_out;
 
-  /*! Task for star formation */
-  struct task *star_formation;
+    /*! Task for cooling */
+    struct task *cooling;
 
-  /*! Task for sink formation */
-  struct task *sink_formation;
+    /*! Task for star formation */
+    struct task *star_formation;
 
-  /*! Task for sorting the stars again after a SF event */
-  struct task *stars_resort;
+    /*! Task for sink formation */
+    struct task *sink_formation;
 
-  /*! Radiative transfer ghost in task */
-  struct task *rt_in;
+    /*! Task for sorting the stars again after a SF event */
+    struct task *stars_resort;
 
-  /*! Task for self/pair injection step of radiative transfer */
-  struct link *rt_inject;
+    /*! Radiative transfer ghost in task */
+    struct task *rt_in;
 
-  /*! Radiative transfer ghost1 task (finishes up injection) */
-  struct task *rt_ghost1;
+    /*! Task for self/pair injection step of radiative transfer */
+    struct link *rt_inject;
 
-  /*! Task for self/pair gradient step of radiative transfer */
-  struct link *rt_gradient;
+    /*! Radiative transfer ghost1 task (finishes up injection) */
+    struct task *rt_ghost1;
 
-  /*! Radiative transfer ghost2 task */
-  struct task *rt_ghost2;
+    /*! Task for self/pair gradient step of radiative transfer */
+    struct link *rt_gradient;
 
-  /*! Task for self/pair transport step of radiative transfer */
-  struct link *rt_transport;
+    /*! Radiative transfer ghost2 task */
+    struct task *rt_ghost2;
 
-  /*! Radiative transfer transport out task */
-  struct task *rt_transport_out;
+    /*! Task for self/pair transport step of radiative transfer */
+    struct link *rt_transport;
 
-  /*! Radiative transfer thermochemistry task */
-  struct task *rt_tchem;
+    /*! Radiative transfer transport out task */
+    struct task *rt_transport_out;
 
-  /*! Radiative transfer ghost out task */
-  struct task *rt_out;
+    /*! Radiative transfer thermochemistry task */
+    struct task *rt_tchem;
 
-  /*! Last (integer) time the cell's part were drifted forward in time. */
-  integertime_t ti_old_part;
+    /*! Radiative transfer ghost out task */
+    struct task *rt_out;
 
-  /*! Minimum end of (integer) time step in this cell for hydro tasks. */
-  integertime_t ti_end_min;
+    /*! Last (integer) time the cell's part were drifted forward in time. */
+    integertime_t ti_old_part;
 
-  /*! Maximum end of (integer) time step in this cell for hydro tasks. */
-  integertime_t ti_end_max;
+    /*! Maximum end of (integer) time step in this cell for hydro tasks. */
+    integertime_t ti_end_max;
 
-  /*! Maximum beginning of (integer) time step in this cell for hydro tasks.
-   */
-  integertime_t ti_beg_max;
+    /*! Max smoothing length of active particles in this cell. */
+    float h_max_active;
 
-  /*! Spin lock for various uses (#part case). */
-  swift_lock_type lock;
+    /*! Values of h_max before the drifts, used for sub-cell tasks. */
+    float h_max_old;
 
-  /*! Max smoothing length in this cell. */
-  float h_max;
+    /*! Maximum part movement in this cell since last construction. */
+    float dx_max_part;
 
-  /*! Maximum part movement in this cell since last construction. */
-  float dx_max_part;
+    /*! Maximum particle movement in this cell since the last sort. */
+    float dx_max_sort;
 
-  /*! Maximum particle movement in this cell since the last sort. */
-  float dx_max_sort;
+    /*! Values of dx_max before the drifts, used for sub-cell tasks. */
+    float dx_max_part_old;
 
-  /*! Values of h_max before the drifts, used for sub-cell tasks. */
-  float h_max_old;
+    /*! Values of dx_max_sort before the drifts, used for sub-cell tasks. */
+    float dx_max_sort_old;
 
-  /*! Values of dx_max before the drifts, used for sub-cell tasks. */
-  float dx_max_part_old;
+    /*! Nr of #part this cell can hold after addition of new #part. */
+    int count_total;
 
-  /*! Values of dx_max_sort before the drifts, used for sub-cell tasks. */
-  float dx_max_sort_old;
+    /*! Bit mask of sort directions that will be needed in the next timestep. */
+    uint16_t requires_sorts;
 
-  /*! Nr of #part in this cell. */
-  int count;
+    /*! Bit mask of sorts that need to be computed for this cell. */
+    uint16_t do_sort;
 
-  /*! Nr of #part this cell can hold after addition of new #part. */
-  int count_total;
+    /*! Bit-mask indicating the sorted directions */
+    uint16_t sorted;
 
-  /*! Number of #part updated in this cell. */
-  int updated;
+    /*! Bit-mask indicating the sorted directions */
+    uint16_t sort_allocated;
 
-  /*! Is the #part data of this cell being used in a sub-cell? */
-  int hold;
+#ifdef SWIFT_DEBUG_CHECKS
 
-  /*! Bit mask of sort directions that will be needed in the next timestep. */
-  uint16_t requires_sorts;
+    /*! Last (integer) time the cell's sort arrays were updated. */
+    integertime_t ti_sort;
 
-  /*! Bit mask of sorts that need to be computed for this cell. */
-  uint16_t do_sort;
+#endif
 
-  /*! Bit-mask indicating the sorted directions */
-  uint16_t sorted;
+#ifdef NONE_SPH
+  };
+#endif
 
-  /*! Bit-mask indicating the sorted directions */
-  uint16_t sort_allocated;
+  /*! Minimum end of (integer) time step in this cell for hydro tasks. */
+  integertime_t ti_end_min;
 
-#ifdef SWIFT_DEBUG_CHECKS
+  /*! Maximum beginning of (integer) time step in this cell for hydro tasks.
+   */
+  integertime_t ti_beg_max;
 
-  /*! Last (integer) time the cell's sort arrays were updated. */
-  integertime_t ti_sort;
+  /*! Spin lock for various uses (#part case). */
+  swift_lock_type lock;
 
-#endif
+  /*! Max smoothing length in this cell. */
+  float h_max;
+
+  /*! Number of #part updated in this cell. */
+  int updated;
+
+  /*! Is the #part data of this cell being used in a sub-cell? */
+  int hold;
+
+  /*! Nr of #part in this cell. */
+  int count;
 };
 
 #endif /* SWIFT_CELL_HYDRO_H */
diff --git a/src/cell_sinks.h b/src/cell_sinks.h
index c6ad0478f53e5890cb0c88f4a99209fa31a8a5ab..5500968f3f28b4e56e08cd0cc2089e9ed20e7daf 100644
--- a/src/cell_sinks.h
+++ b/src/cell_sinks.h
@@ -33,30 +33,58 @@
  */
 struct cell_sinks {
 
-  /*! Pointer to the #sink data. */
-  struct sink *parts;
+  /* If we are not using sinks, compact as much of the unecessary variables
+     into an anonymous union to save memory in the cell structure. */
+#ifdef SINK_NONE
+  union {
+#endif
 
-  /*! Linked list of the tasks computing this cell's sink formation checks. */
-  struct link *compute_formation;
+    /*! Pointer to the #sink data. */
+    struct sink *parts;
 
-  /*! The drift task for sinks */
-  struct task *drift;
+    /*! Linked list of the tasks computing this cell's sink formation checks. */
+    struct link *compute_formation;
 
-  /*! Implicit tasks marking the entry of the sink block of tasks */
-  struct task *sink_in;
+    /*! The drift task for sinks */
+    struct task *drift;
 
-  /*! Implicit tasks marking the exit of the sink block of tasks */
-  struct task *sink_out;
+    /*! Implicit tasks marking the entry of the sink block of tasks */
+    struct task *sink_in;
 
-  /*! Last (integer) time the cell's sink were drifted forward in time. */
-  integertime_t ti_old_part;
+    /*! Implicit tasks marking the exit of the sink block of tasks */
+    struct task *sink_out;
+
+    /*! Last (integer) time the cell's sink were drifted forward in time. */
+    integertime_t ti_old_part;
+
+    /*! Maximum end of (integer) time step in this cell for sink tasks. */
+    integertime_t ti_end_max;
+
+    /*! Spin lock for sink formation use. */
+    swift_lock_type sink_formation_lock;
+
+    /*! Nr of #sink this cell can hold after addition of new one. */
+    int count_total;
+
+    /*! Max cut off radius of active particles in this cell. */
+    float r_cut_max_active;
+
+    /*! Values of r_cut_max before the drifts, used for sub-cell tasks. */
+    float r_cut_max_old;
+
+    /*! Maximum part movement in this cell since last construction. */
+    float dx_max_part;
+
+    /*! Values of dx_max before the drifts, used for sub-cell tasks. */
+    float dx_max_part_old;
+
+#ifdef SINK_NONE
+  };
+#endif
 
   /*! Minimum end of (integer) time step in this cell for sink tasks. */
   integertime_t ti_end_min;
 
-  /*! Maximum end of (integer) time step in this cell for sink tasks. */
-  integertime_t ti_end_max;
-
   /*! Maximum beginning of (integer) time step in this cell for sink
    * tasks. */
   integertime_t ti_beg_max;
@@ -64,14 +92,8 @@ struct cell_sinks {
   /*! Spin lock for various uses (#sink case). */
   swift_lock_type lock;
 
-  /*! Spin lock for sink formation use. */
-  swift_lock_type sink_formation_lock;
-
-  /*! Nr of #sink in this cell. */
-  int count;
-
-  /*! Nr of #sink this cell can hold after addition of new one. */
-  int count_total;
+  /*! Max cut off radius in this cell. */
+  float r_cut_max;
 
   /*! Number of #sink updated in this cell. */
   int updated;
@@ -79,17 +101,8 @@ struct cell_sinks {
   /*! Is the #sink data of this cell being used in a sub-cell? */
   int hold;
 
-  /*! Max cut off radius in this cell. */
-  float r_cut_max;
-
-  /*! Values of r_cut_max before the drifts, used for sub-cell tasks. */
-  float r_cut_max_old;
-
-  /*! Maximum part movement in this cell since last construction. */
-  float dx_max_part;
-
-  /*! Values of dx_max before the drifts, used for sub-cell tasks. */
-  float dx_max_part_old;
+  /*! Nr of #sink in this cell. */
+  int count;
 };
 
 #endif /* SWIFT_CELL_SINKS_H */
diff --git a/src/cell_stars.h b/src/cell_stars.h
index bb27b23036be1095f069efac29ceaf7444f6fe02..a4a005e9ebfc5ccfeb0b64dcc20474df1927ab01 100644
--- a/src/cell_stars.h
+++ b/src/cell_stars.h
@@ -34,105 +34,118 @@
  */
 struct cell_stars {
 
-  /*! Pointer to the #spart data. */
-  struct spart *parts;
+  /* If we are not using stars, compact as much of the unecessary variables
+     into an anonymous union to save memory in the cell structure. */
+#ifdef STARS_NONE
+  union {
+#endif
 
-  /*! Pointer to the #spart data at rebuild time. */
-  struct spart *parts_rebuild;
+    /*! Pointer to the #spart data. */
+    struct spart *parts;
 
-  /*! The star ghost task itself */
-  struct task *ghost;
+    /*! Pointer to the #spart data at rebuild time. */
+    struct spart *parts_rebuild;
 
-  /*! Linked list of the tasks computing this cell's star density. */
-  struct link *density;
+    /*! The star ghost task itself */
+    struct task *ghost;
 
-  /*! Linked list of the tasks computing this cell's star feedback. */
-  struct link *feedback;
+    /*! Linked list of the tasks computing this cell's star density. */
+    struct link *density;
 
-  /*! The task computing this cell's sorts before the density. */
-  struct task *sorts;
+    /*! Linked list of the tasks computing this cell's star feedback. */
+    struct link *feedback;
 
-  /*! The drift task for sparts */
-  struct task *drift;
+    /*! The task computing this cell's sorts before the density. */
+    struct task *sorts;
 
-  /*! Implicit tasks marking the entry of the stellar physics block of tasks
-   */
-  struct task *stars_in;
+    /*! The drift task for sparts */
+    struct task *drift;
 
-  /*! Implicit tasks marking the exit of the stellar physics block of tasks */
-  struct task *stars_out;
+    /*! Implicit tasks marking the entry of the stellar physics block of tasks
+     */
+    struct task *stars_in;
 
-  /*! Pointer for the sorted indices. */
-  struct sort_entry *sort;
+    /*! Implicit tasks marking the exit of the stellar physics block of tasks */
+    struct task *stars_out;
 
-  /*! Last (integer) time the cell's spart were drifted forward in time. */
-  integertime_t ti_old_part;
+    /*! Pointer for the sorted indices. */
+    struct sort_entry *sort;
 
-  /*! Minimum end of (integer) time step in this cell for star tasks. */
-  integertime_t ti_end_min;
+    /*! Last (integer) time the cell's spart were drifted forward in time. */
+    integertime_t ti_old_part;
 
-  /*! Maximum end of (integer) time step in this cell for star tasks. */
-  integertime_t ti_end_max;
+    /*! Maximum end of (integer) time step in this cell for star tasks. */
+    integertime_t ti_end_max;
 
-  /*! Maximum beginning of (integer) time step in this cell for star tasks.
-   */
-  integertime_t ti_beg_max;
+    /*! Spin lock for star formation use. */
+    swift_lock_type star_formation_lock;
 
-  /*! Spin lock for various uses (#spart case). */
-  swift_lock_type lock;
+    /*! Nr of #spart this cell can hold after addition of new #spart. */
+    int count_total;
 
-  /*! Spin lock for star formation use. */
-  swift_lock_type star_formation_lock;
+    /*! Max smoothing length of active particles in this cell. */
+    float h_max_active;
 
-  /*! Nr of #spart in this cell. */
-  int count;
+    /*! Values of h_max before the drifts, used for sub-cell tasks. */
+    float h_max_old;
 
-  /*! Nr of #spart this cell can hold after addition of new #spart. */
-  int count_total;
+    /*! Maximum part movement in this cell since last construction. */
+    float dx_max_part;
 
-  /*! Number of #spart updated in this cell. */
-  int updated;
+    /*! Values of dx_max before the drifts, used for sub-cell tasks. */
+    float dx_max_part_old;
 
-  /*! Is the #spart data of this cell being used in a sub-cell? */
-  int hold;
+    /*! Maximum particle movement in this cell since the last sort. */
+    float dx_max_sort;
 
-  /*! Max smoothing length in this cell. */
-  float h_max;
+    /*! Values of dx_max_sort before the drifts, used for sub-cell tasks. */
+    float dx_max_sort_old;
+
+    /*! Bit mask of sort directions that will be needed in the next timestep. */
+    uint16_t requires_sorts;
+
+    /*! Bit-mask indicating the sorted directions */
+    uint16_t sorted;
+
+    /*! Bit-mask indicating the sorted directions */
+    uint16_t sort_allocated;
 
-  /*! Values of h_max before the drifts, used for sub-cell tasks. */
-  float h_max_old;
+    /*! Bit mask of sorts that need to be computed for this cell. */
+    uint16_t do_sort;
 
-  /*! Maximum part movement in this cell since last construction. */
-  float dx_max_part;
+    /*! Star formation history struct */
+    struct star_formation_history sfh;
 
-  /*! Values of dx_max before the drifts, used for sub-cell tasks. */
-  float dx_max_part_old;
+#ifdef SWIFT_DEBUG_CHECKS
+    /*! Last (integer) time the cell's sort arrays were updated. */
+    integertime_t ti_sort;
+#endif
 
-  /*! Maximum particle movement in this cell since the last sort. */
-  float dx_max_sort;
+#ifdef STARS_NONE
+  };
+#endif
 
-  /*! Values of dx_max_sort before the drifts, used for sub-cell tasks. */
-  float dx_max_sort_old;
+  /*! Maximum end of (integer) time step in this cell for star tasks. */
+  integertime_t ti_end_min;
 
-  /*! Bit mask of sort directions that will be needed in the next timestep. */
-  uint16_t requires_sorts;
+  /*! Maximum beginning of (integer) time step in this cell for star tasks.
+   */
+  integertime_t ti_beg_max;
 
-  /*! Bit-mask indicating the sorted directions */
-  uint16_t sorted;
+  /*! Spin lock for various uses (#spart case). */
+  swift_lock_type lock;
 
-  /*! Bit-mask indicating the sorted directions */
-  uint16_t sort_allocated;
+  /*! Max smoothing length in this cell. */
+  float h_max;
 
-  /*! Bit mask of sorts that need to be computed for this cell. */
-  uint16_t do_sort;
+  /*! Number of #spart updated in this cell. */
+  int updated;
 
-  /*! Star formation history struct */
-  struct star_formation_history sfh;
+  /*! Nr of #spart in this cell. */
+  int count;
 
-#ifdef SWIFT_DEBUG_CHECKS
-  /*! Last (integer) time the cell's sort arrays were updated. */
-  integertime_t ti_sort;
-#endif
+  /*! Is the #spart data of this cell being used in a sub-cell? */
+  int hold;
 };
 
 #endif /* SWIFT_CELL_STARS_H */
diff --git a/src/engine.c b/src/engine.c
index 1d00e2091cf823d4f9ecebda4b68605ebca43c5c..337ee530a3c47f2235b2bb15492bb88ffe1eef68 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -1256,10 +1256,12 @@ void engine_rebuild(struct engine *e, const int repartitioned,
   if (clean_smoothing_length_values) space_sanitize(e->s);
 
 /* If in parallel, exchange the cell structure, top-level and neighbouring
- * multipoles. */
+ * multipoles. To achieve this, free the foreign particle buffers first. */
 #ifdef WITH_MPI
   if (e->policy & engine_policy_self_gravity) engine_exchange_top_multipoles(e);
 
+  space_free_foreign_parts(e->s, /*clear_cell_pointers=*/1);
+
   engine_exchange_cells(e);
 #endif
 
diff --git a/src/engine_io.c b/src/engine_io.c
index 9f5896942ad3105d849ed02e1c2f25dcf4dd86f4..e5ed79a04352badefe7890757ec65f1abb542b1b 100644
--- a/src/engine_io.c
+++ b/src/engine_io.c
@@ -378,11 +378,21 @@ void engine_check_for_dumps(struct engine *e) {
         e->force_checks_snapshot_flag = 1;
 #endif
 
+        /* Free the mesh memory to get some breathing space */
+        if ((e->policy & engine_policy_self_gravity) && e->s->periodic)
+          pm_mesh_free(e->mesh);
+
         /* Do we want FoF group IDs in the snapshot? */
         if (with_fof && e->snapshot_invoke_fof) {
           engine_fof(e, /*dump_results=*/0, /*seed_black_holes=*/0);
         }
 
+        /* Free the foreign particles to get more breathing space.
+         * This cannot be done before FOF as comms are used in there. */
+#ifdef WITH_MPI
+        space_free_foreign_parts(e->s, /*clear_cell_pointers=*/1);
+#endif
+
         /* Do we want a corresponding VELOCIraptor output? */
         if (with_stf && e->snapshot_invoke_stf && !e->stf_this_timestep) {
 
@@ -407,6 +417,13 @@ void engine_check_for_dumps(struct engine *e) {
 #endif
         }
 
+        /* Reallocate freed memory */
+        if ((e->policy & engine_policy_self_gravity) && e->s->periodic)
+          pm_mesh_allocate(e->mesh);
+#ifdef WITH_MPI
+        engine_allocate_foreign_particles(e);
+#endif
+
         /* ... and find the next output time */
         engine_compute_next_snapshot_time(e);
         break;
@@ -423,6 +440,13 @@ void engine_check_for_dumps(struct engine *e) {
 
       case output_stf:
 
+        /* Free the mesh memory to get some breathing space */
+        if ((e->policy & engine_policy_self_gravity) && e->s->periodic)
+          pm_mesh_free(e->mesh);
+#ifdef WITH_MPI
+        space_free_foreign_parts(e->s, /*clear_cell_pointers=*/1);
+#endif
+
 #ifdef HAVE_VELOCIRAPTOR
         /* Unleash the raptor! */
         if (!e->stf_this_timestep) {
@@ -430,6 +454,13 @@ void engine_check_for_dumps(struct engine *e) {
           e->step_props |= engine_step_prop_stf;
         }
 
+        /* Reallocate freed memory */
+        if ((e->policy & engine_policy_self_gravity) && e->s->periodic)
+          pm_mesh_allocate(e->mesh);
+#ifdef WITH_MPI
+        engine_allocate_foreign_particles(e);
+#endif
+
         /* ... and find the next output time */
         engine_compute_next_stf_time(e);
 #else
diff --git a/src/part.h b/src/part.h
index cad222c0bc789b2755d55fb96e7748e3cc9c258f..8382fa10157066e59b99549677fe3ddf83185593 100644
--- a/src/part.h
+++ b/src/part.h
@@ -101,9 +101,7 @@ struct threadpool;
 #endif
 
 /* Import the right star particle definition */
-#if defined(FEEDBACK_CONST)
-#include "./stars/const/stars_part.h"
-#elif defined(STARS_NONE)
+#if defined(STARS_NONE)
 #include "./stars/Default/stars_part.h"
 #elif defined(STARS_EAGLE)
 #include "./stars/EAGLE/stars_part.h"
diff --git a/src/runner_ghost.c b/src/runner_ghost.c
index da95ee95e749402a2facdd4e008493e9838b5f99..67de56005712e8988a6ff2fe245fca8d9f10c230 100644
--- a/src/runner_ghost.c
+++ b/src/runner_ghost.c
@@ -89,7 +89,8 @@ void runner_do_stars_ghost(struct runner *r, struct cell *c, int timer) {
   int redo = 0, scount = 0;
 
   /* Running value of the maximal smoothing length */
-  double h_max = c->stars.h_max;
+  float h_max = c->stars.h_max;
+  float h_max_active = c->stars.h_max_active;
 
   TIMER_TIC;
 
@@ -110,6 +111,7 @@ void runner_do_stars_ghost(struct runner *r, struct cell *c, int timer) {
 
         /* Update h_max */
         h_max = max(h_max, c->progeny[k]->stars.h_max);
+        h_max_active = max(h_max_active, c->progeny[k]->stars.h_max_active);
       }
     }
   } else {
@@ -349,6 +351,7 @@ void runner_do_stars_ghost(struct runner *r, struct cell *c, int timer) {
 
         /* Check if h_max has increased */
         h_max = max(h_max, sp->h);
+        h_max_active = max(h_max_active, sp->h);
 
         stars_reset_feedback(sp);
 
@@ -500,12 +503,27 @@ void runner_do_stars_ghost(struct runner *r, struct cell *c, int timer) {
 
   /* Update h_max */
   c->stars.h_max = h_max;
+  c->stars.h_max_active = h_max_active;
+
+#ifdef SWIFT_DEBUG_CHECKS
+  for (int i = 0; i < c->stars.count; ++i) {
+    const struct spart *sp = &c->stars.parts[i];
+    const float h = c->stars.parts[i].h;
+    if (spart_is_inhibited(sp, e)) continue;
+
+    if (h > c->stars.h_max)
+      error("Particle has h larger than h_max (id=%lld)", sp->id);
+    if (spart_is_active(sp, e) && h > c->stars.h_max_active)
+      error("Active particle has h larger than h_max_active (id=%lld)", sp->id);
+  }
+#endif
 
   /* The ghost may not always be at the top level.
    * Therefore we need to update h_max between the super- and top-levels */
   if (c->stars.ghost) {
     for (struct cell *tmp = c->parent; tmp != NULL; tmp = tmp->parent) {
       atomic_max_f(&tmp->stars.h_max, h_max);
+      atomic_max_f(&tmp->stars.h_max_active, h_max_active);
     }
   }
 
@@ -535,7 +553,8 @@ void runner_do_black_holes_density_ghost(struct runner *r, struct cell *c,
   int redo = 0, bcount = 0;
 
   /* Running value of the maximal smoothing length */
-  double h_max = c->black_holes.h_max;
+  float h_max = c->black_holes.h_max;
+  float h_max_active = c->black_holes.h_max_active;
 
   TIMER_TIC;
 
@@ -551,6 +570,8 @@ void runner_do_black_holes_density_ghost(struct runner *r, struct cell *c,
 
         /* Update h_max */
         h_max = max(h_max, c->progeny[k]->black_holes.h_max);
+        h_max_active =
+            max(h_max_active, c->progeny[k]->black_holes.h_max_active);
       }
     }
   } else {
@@ -739,6 +760,7 @@ void runner_do_black_holes_density_ghost(struct runner *r, struct cell *c,
 
         /* Check if h_max has increased */
         h_max = max(h_max, bp->h);
+        h_max_active = max(h_max_active, bp->h);
       }
 
       /* We now need to treat the particles whose smoothing length had not
@@ -811,12 +833,27 @@ void runner_do_black_holes_density_ghost(struct runner *r, struct cell *c,
 
   /* Update h_max */
   c->black_holes.h_max = h_max;
+  c->black_holes.h_max_active = h_max_active;
+
+#ifdef SWIFT_DEBUG_CHECKS
+  for (int i = 0; i < c->black_holes.count; ++i) {
+    const struct bpart *bp = &c->black_holes.parts[i];
+    const float h = c->black_holes.parts[i].h;
+    if (bpart_is_inhibited(bp, e)) continue;
+
+    if (h > c->black_holes.h_max)
+      error("Particle has h larger than h_max (id=%lld)", bp->id);
+    if (bpart_is_active(bp, e) && h > c->black_holes.h_max_active)
+      error("Active particle has h larger than h_max_active (id=%lld)", bp->id);
+  }
+#endif
 
   /* The ghost may not always be at the top level.
    * Therefore we need to update h_max between the super- and top-levels */
   if (c->black_holes.density_ghost) {
     for (struct cell *tmp = c->parent; tmp != NULL; tmp = tmp->parent) {
       atomic_max_f(&tmp->black_holes.h_max, h_max);
+      atomic_max_f(&tmp->black_holes.h_max_active, h_max_active);
     }
   }
 
@@ -1005,7 +1042,8 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
   int redo = 0, count = 0;
 
   /* Running value of the maximal smoothing length */
-  double h_max = c->hydro.h_max;
+  float h_max = c->hydro.h_max;
+  float h_max_active = c->hydro.h_max_active;
 
   TIMER_TIC;
 
@@ -1021,6 +1059,7 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
 
         /* Update h_max */
         h_max = max(h_max, c->progeny[k]->hydro.h_max);
+        h_max_active = max(h_max_active, c->progeny[k]->hydro.h_max_active);
       }
     }
   } else {
@@ -1289,8 +1328,9 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
 
         /* We now have a particle whose smoothing length has converged */
 
-        /* Check if h_max is increased */
+        /* Check if h_max has increased */
         h_max = max(h_max, p->h);
+        h_max_active = max(h_max_active, p->h);
 
 #ifdef EXTRA_HYDRO_LOOP
 
@@ -1409,12 +1449,27 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
 
   /* Update h_max */
   c->hydro.h_max = h_max;
+  c->hydro.h_max_active = h_max_active;
+
+#ifdef SWIFT_DEBUG_CHECKS
+  for (int i = 0; i < c->hydro.count; ++i) {
+    const struct part *p = &c->hydro.parts[i];
+    const float h = c->hydro.parts[i].h;
+    if (part_is_inhibited(p, e)) continue;
+
+    if (h > c->hydro.h_max)
+      error("Particle has h larger than h_max (id=%lld)", p->id);
+    if (part_is_active(p, e) && h > c->hydro.h_max_active)
+      error("Active particle has h larger than h_max_active (id=%lld)", p->id);
+  }
+#endif
 
   /* The ghost may not always be at the top level.
    * Therefore we need to update h_max between the super- and top-levels */
   if (c->hydro.ghost) {
     for (struct cell *tmp = c->parent; tmp != NULL; tmp = tmp->parent) {
       atomic_max_f(&tmp->hydro.h_max, h_max);
+      atomic_max_f(&tmp->hydro.h_max_active, h_max_active);
     }
   }
 
diff --git a/src/runner_others.c b/src/runner_others.c
index b7f016136d75ccebfefdda79fbf958dfc3e7a61f..1886569708b5dbde7e95c22135aa4d428c1b868e 100644
--- a/src/runner_others.c
+++ b/src/runner_others.c
@@ -229,6 +229,11 @@ void runner_do_star_formation(struct runner *r, struct cell *c, int timer) {
         /* Update current cell using child cells */
         star_formation_logger_add(&c->stars.sfh, &cp->stars.sfh);
 
+        /* Update the h_max */
+        c->stars.h_max = max(c->stars.h_max, cp->stars.h_max);
+        c->stars.h_max_active =
+            max(c->stars.h_max_active, cp->stars.h_max_active);
+
         /* Update the dx_max */
         if (star_formation_need_update_dx_max) {
           c->hydro.dx_max_part =
@@ -324,6 +329,10 @@ void runner_do_star_formation(struct runner *r, struct cell *c, int timer) {
               /* Update the Star formation history */
               star_formation_logger_log_new_spart(sp, &c->stars.sfh);
 
+              /* Update the h_max */
+              c->stars.h_max = max(c->stars.h_max, sp->h);
+              c->stars.h_max_active = max(c->stars.h_max_active, sp->h);
+
               /* Update the displacement information */
               if (star_formation_need_update_dx_max) {
                 const float dx2_part = xp->x_diff[0] * xp->x_diff[0] +
diff --git a/src/runner_recv.c b/src/runner_recv.c
index a00e3d2dc506d45bebaee4c1682ba7942efea768..e1aec4a73897d6996c60998f7c18b06ffda980e3 100644
--- a/src/runner_recv.c
+++ b/src/runner_recv.c
@@ -49,6 +49,7 @@ void runner_do_recv_part(struct runner *r, struct cell *c, int clear_sorts,
   const struct part *restrict parts = c->hydro.parts;
   const size_t nr_parts = c->hydro.count;
   const integertime_t ti_current = r->e->ti_current;
+  const timebin_t max_active_bin = r->e->max_active_bin;
 
   TIMER_TIC;
 
@@ -57,6 +58,7 @@ void runner_do_recv_part(struct runner *r, struct cell *c, int clear_sorts,
   timebin_t time_bin_min = num_time_bins;
   timebin_t time_bin_max = 0;
   float h_max = 0.f;
+  float h_max_active = 0.f;
 
 #ifdef SWIFT_DEBUG_CHECKS
   if (c->nodeID == engine_rank) error("Updating a local cell!");
@@ -74,6 +76,8 @@ void runner_do_recv_part(struct runner *r, struct cell *c, int clear_sorts,
       time_bin_min = min(time_bin_min, parts[k].time_bin);
       time_bin_max = max(time_bin_max, parts[k].time_bin);
       h_max = max(h_max, parts[k].h);
+      if (parts[k].time_bin <= max_active_bin)
+        h_max_active = max(h_max_active, parts[k].h);
     }
 
     /* Convert into a time */
@@ -91,6 +95,7 @@ void runner_do_recv_part(struct runner *r, struct cell *c, int clear_sorts,
         ti_hydro_end_max =
             max(ti_hydro_end_max, c->progeny[k]->hydro.ti_end_max);
         h_max = max(h_max, c->progeny[k]->hydro.h_max);
+        h_max_active = max(h_max_active, c->progeny[k]->hydro.h_max_active);
       }
     }
   }
@@ -110,6 +115,7 @@ void runner_do_recv_part(struct runner *r, struct cell *c, int clear_sorts,
   // c->hydro.ti_end_max = ti_hydro_end_max;
   c->hydro.ti_old_part = ti_current;
   c->hydro.h_max = h_max;
+  c->hydro.h_max_active = h_max_active;
 
   if (timer) TIMER_TOC(timer_dorecv_part);
 
@@ -210,6 +216,7 @@ void runner_do_recv_spart(struct runner *r, struct cell *c, int clear_sorts,
   struct spart *restrict sparts = c->stars.parts;
   const size_t nr_sparts = c->stars.count;
   const integertime_t ti_current = r->e->ti_current;
+  const timebin_t max_active_bin = r->e->max_active_bin;
 
   TIMER_TIC;
 
@@ -218,6 +225,7 @@ void runner_do_recv_spart(struct runner *r, struct cell *c, int clear_sorts,
   timebin_t time_bin_min = num_time_bins;
   timebin_t time_bin_max = 0;
   float h_max = 0.f;
+  float h_max_active = 0.f;
 
 #ifdef SWIFT_DEBUG_CHECKS
   if (c->nodeID == engine_rank) error("Updating a local cell!");
@@ -238,6 +246,8 @@ void runner_do_recv_spart(struct runner *r, struct cell *c, int clear_sorts,
       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);
+      if (sparts[k].time_bin <= max_active_bin)
+        h_max_active = max(h_max_active, sparts[k].h);
     }
 
     /* Convert into a time */
@@ -255,6 +265,7 @@ void runner_do_recv_spart(struct runner *r, struct cell *c, int clear_sorts,
         ti_stars_end_max =
             max(ti_stars_end_max, c->progeny[k]->stars.ti_end_max);
         h_max = max(h_max, c->progeny[k]->stars.h_max);
+        h_max_active = max(h_max_active, c->progeny[k]->stars.h_max_active);
       }
     }
   }
@@ -273,6 +284,7 @@ void runner_do_recv_spart(struct runner *r, struct cell *c, int clear_sorts,
   // c->grav.ti_end_max = ti_gravity_end_max;
   c->stars.ti_old_part = ti_current;
   c->stars.h_max = h_max;
+  c->stars.h_max_active = h_max_active;
 
   if (timer) TIMER_TOC(timer_dorecv_spart);
 
@@ -300,6 +312,7 @@ void runner_do_recv_bpart(struct runner *r, struct cell *c, int clear_sorts,
   struct bpart *restrict bparts = c->black_holes.parts;
   const size_t nr_bparts = c->black_holes.count;
   const integertime_t ti_current = r->e->ti_current;
+  const timebin_t max_active_bin = r->e->max_active_bin;
 
   TIMER_TIC;
 
@@ -308,6 +321,7 @@ void runner_do_recv_bpart(struct runner *r, struct cell *c, int clear_sorts,
   timebin_t time_bin_min = num_time_bins;
   timebin_t time_bin_max = 0;
   float h_max = 0.f;
+  float h_max_active = 0.f;
 
 #ifdef SWIFT_DEBUG_CHECKS
   if (c->nodeID == engine_rank) error("Updating a local cell!");
@@ -322,13 +336,12 @@ void runner_do_recv_bpart(struct runner *r, struct cell *c, int clear_sorts,
       bparts[k].num_ngb_force = 0;
 #endif
 
-      /* message("Receiving bparts id=%lld time_bin=%d", */
-      /* 	      bparts[k].id, bparts[k].time_bin); */
-
       if (bparts[k].time_bin == time_bin_inhibited) continue;
       time_bin_min = min(time_bin_min, bparts[k].time_bin);
       time_bin_max = max(time_bin_max, bparts[k].time_bin);
       h_max = max(h_max, bparts[k].h);
+      if (bparts[k].time_bin <= max_active_bin)
+        h_max_active = max(h_max_active, bparts[k].h);
     }
 
     /* Convert into a time */
@@ -346,6 +359,8 @@ void runner_do_recv_bpart(struct runner *r, struct cell *c, int clear_sorts,
         ti_black_holes_end_max =
             max(ti_black_holes_end_max, c->progeny[k]->black_holes.ti_end_max);
         h_max = max(h_max, c->progeny[k]->black_holes.h_max);
+        h_max_active =
+            max(h_max_active, c->progeny[k]->black_holes.h_max_active);
       }
     }
   }
@@ -363,6 +378,7 @@ void runner_do_recv_bpart(struct runner *r, struct cell *c, int clear_sorts,
   // c->grav.ti_end_max = ti_gravity_end_max;
   c->black_holes.ti_old_part = ti_current;
   c->black_holes.h_max = h_max;
+  c->black_holes.h_max_active = h_max_active;
 
   if (timer) TIMER_TOC(timer_dorecv_bpart);
 
diff --git a/src/space_split.c b/src/space_split.c
index 38d75ad48b3e10c21d7047bb6471f41ee16bd881..af0e715a1f2e64e678e8cc8f181530e8a6ebcdf5 100644
--- a/src/space_split.c
+++ b/src/space_split.c
@@ -26,6 +26,7 @@
 #include "space.h"
 
 /* Local headers. */
+#include "active.h"
 #include "cell.h"
 #include "debug.h"
 #include "engine.h"
@@ -65,9 +66,13 @@ void space_split_recursive(struct space *s, struct cell *c,
   const int depth = c->depth;
   int maxdepth = 0;
   float h_max = 0.0f;
-  float sinks_h_max = 0.f;
+  float h_max_active = 0.0f;
   float stars_h_max = 0.f;
+  float stars_h_max_active = 0.f;
   float black_holes_h_max = 0.f;
+  float black_holes_h_max_active = 0.f;
+  float sinks_h_max = 0.f;
+  float sinks_h_max_active = 0.f;
   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,
@@ -224,14 +229,18 @@ void space_split_recursive(struct space *s, struct cell *c,
       cp->depth = c->depth + 1;
       cp->split = 0;
       cp->hydro.h_max = 0.f;
+      cp->hydro.h_max_active = 0.f;
       cp->hydro.dx_max_part = 0.f;
       cp->hydro.dx_max_sort = 0.f;
       cp->stars.h_max = 0.f;
+      cp->stars.h_max_active = 0.f;
       cp->stars.dx_max_part = 0.f;
       cp->stars.dx_max_sort = 0.f;
       cp->sinks.r_cut_max = 0.f;
+      cp->sinks.r_cut_max_active = 0.f;
       cp->sinks.dx_max_part = 0.f;
       cp->black_holes.h_max = 0.f;
+      cp->black_holes.h_max_active = 0.f;
       cp->black_holes.dx_max_part = 0.f;
       cp->nodeID = c->nodeID;
       cp->parent = c;
@@ -286,9 +295,15 @@ void space_split_recursive(struct space *s, struct cell *c,
 
         /* Update the cell-wide properties */
         h_max = max(h_max, cp->hydro.h_max);
+        h_max_active = max(h_max_active, cp->hydro.h_max_active);
         stars_h_max = max(stars_h_max, cp->stars.h_max);
+        stars_h_max_active = max(stars_h_max_active, cp->stars.h_max_active);
         black_holes_h_max = max(black_holes_h_max, cp->black_holes.h_max);
+        black_holes_h_max_active =
+            max(black_holes_h_max_active, cp->black_holes.h_max_active);
         sinks_h_max = max(sinks_h_max, cp->sinks.r_cut_max);
+        sinks_h_max_active =
+            max(sinks_h_max_active, cp->sinks.r_cut_max_active);
 
         ti_hydro_end_min = min(ti_hydro_end_min, cp->hydro.ti_end_min);
         ti_hydro_end_max = max(ti_hydro_end_max, cp->hydro.ti_end_max);
@@ -468,6 +483,9 @@ void space_split_recursive(struct space *s, struct cell *c,
 
       h_max = max(h_max, parts[k].h);
 
+      if (part_is_active(&parts[k], e))
+        h_max_active = max(h_max_active, parts[k].h);
+
       /* Collect SFR from the particles after rebuilt */
       star_formation_logger_log_inactive_part(&parts[k], &xparts[k],
                                               &c->stars.sfh);
@@ -519,6 +537,9 @@ void space_split_recursive(struct space *s, struct cell *c,
 
       stars_h_max = max(stars_h_max, sparts[k].h);
 
+      if (spart_is_active(&sparts[k], e))
+        stars_h_max_active = max(stars_h_max_active, sparts[k].h);
+
       /* Reset x_diff */
       sparts[k].x_diff[0] = 0.f;
       sparts[k].x_diff[1] = 0.f;
@@ -545,6 +566,9 @@ void space_split_recursive(struct space *s, struct cell *c,
 
       sinks_h_max = max(sinks_h_max, sinks[k].r_cut);
 
+      if (sink_is_active(&sinks[k], e))
+        sinks_h_max_active = max(sinks_h_max_active, sinks[k].r_cut);
+
       /* Reset x_diff */
       sinks[k].x_diff[0] = 0.f;
       sinks[k].x_diff[1] = 0.f;
@@ -571,6 +595,9 @@ void space_split_recursive(struct space *s, struct cell *c,
 
       black_holes_h_max = max(black_holes_h_max, bparts[k].h);
 
+      if (bpart_is_active(&bparts[k], e))
+        black_holes_h_max_active = max(black_holes_h_max_active, bparts[k].h);
+
       /* Reset x_diff */
       bparts[k].x_diff[0] = 0.f;
       bparts[k].x_diff[1] = 0.f;
@@ -611,6 +638,7 @@ void space_split_recursive(struct space *s, struct cell *c,
 
   /* Set the values for this cell. */
   c->hydro.h_max = h_max;
+  c->hydro.h_max_active = h_max_active;
   c->hydro.ti_end_min = ti_hydro_end_min;
   c->hydro.ti_end_max = ti_hydro_end_max;
   c->hydro.ti_beg_max = ti_hydro_beg_max;
@@ -621,14 +649,17 @@ void space_split_recursive(struct space *s, struct cell *c,
   c->stars.ti_end_max = ti_stars_end_max;
   c->stars.ti_beg_max = ti_stars_beg_max;
   c->stars.h_max = stars_h_max;
+  c->stars.h_max_active = stars_h_max_active;
   c->sinks.ti_end_min = ti_sinks_end_min;
   c->sinks.ti_end_max = ti_sinks_end_max;
   c->sinks.ti_beg_max = ti_sinks_beg_max;
   c->sinks.r_cut_max = sinks_h_max;
+  c->sinks.r_cut_max_active = sinks_h_max_active;
   c->black_holes.ti_end_min = ti_black_holes_end_min;
   c->black_holes.ti_end_max = ti_black_holes_end_max;
   c->black_holes.ti_beg_max = ti_black_holes_beg_max;
   c->black_holes.h_max = black_holes_h_max;
+  c->black_holes.h_max_active = black_holes_h_max_active;
   c->maxdepth = maxdepth;
 
   /* Set ownership according to the start of the parts array. */
diff --git a/src/stars_io.h b/src/stars_io.h
index 18224b1ef9a189c719d3674a037eb4b26cd14d4e..c2a095c47c8d491fe6cf97d22c367d1e9d4a7fd6 100644
--- a/src/stars_io.h
+++ b/src/stars_io.h
@@ -23,9 +23,7 @@
 #include "./const.h"
 
 /* Load the correct star type */
-#if defined(FEEDBACK_CONST)
-#include "./stars/const/stars_io.h"
-#elif defined(STARS_NONE)
+#if defined(STARS_NONE)
 #include "./stars/Default/stars_io.h"
 #elif defined(STARS_EAGLE)
 #include "./stars/EAGLE/stars_io.h"
diff --git a/src/stars_logger.h b/src/stars_logger.h
index a3d1dbb30962dcec9ad3ece094ea8e0289b5ec05..3302728d45d51cfe0a8679370acf817984053dba 100644
--- a/src/stars_logger.h
+++ b/src/stars_logger.h
@@ -30,9 +30,7 @@
 #include "timeline.h"
 
 /* Load the correct star type */
-#if defined(FEEDBACK_CONST)
-#error TODO
-#elif defined(STARS_NONE)
+#if defined(STARS_NONE)
 #include "./stars/Default/stars_logger.h"
 #elif defined(STARS_EAGLE)
 #error TODO
diff --git a/src/velociraptor_interface.c b/src/velociraptor_interface.c
index aabb58906e626f03f4ede727b7462e261d18ea19..581d1f791475578ab832577ddef12e918ecfe6c9 100644
--- a/src/velociraptor_interface.c
+++ b/src/velociraptor_interface.c
@@ -748,12 +748,6 @@ void velociraptor_invoke(struct engine *e, const int linked_with_snap) {
   const int nr_cells = s->nr_cells;
   const struct cell *cells_top = s->cells_top;
 
-  /* Start by freeing some of the unnecessary memory to give VR some breathing
-     space */
-#ifdef WITH_MPI
-  space_free_foreign_parts(e->s, /*clear_cell_pointers=*/1);
-#endif
-
   /* Allow thread to run on any core for the duration of the call to
    * VELOCIraptor so that  when OpenMP threads are spawned
    * they can run on any core on the processor. */
@@ -1117,12 +1111,6 @@ void velociraptor_invoke(struct engine *e, const int linked_with_snap) {
   /* Record we have ran stf this timestep */
   e->stf_this_timestep = 1;
 
-  /* Reallocate the memory that was freed earlier */
-#ifdef WITH_MPI
-
-  engine_allocate_foreign_particles(e);
-#endif
-
 #else
   error("SWIFT not configured to run with VELOCIraptor.");
 #endif /* HAVE_VELOCIRAPTOR */
diff --git a/tests/test125cells.c b/tests/test125cells.c
index 98a47404da2f6ac0b57a9bbdc4b4c68f759ec0e6..e1c0766fd81d9eb928a588217a16c88e46e45581 100644
--- a/tests/test125cells.c
+++ b/tests/test125cells.c
@@ -342,6 +342,7 @@ struct cell *make_cell(size_t n, const double offset[3], double size, double h,
   /* Cell properties */
   cell->split = 0;
   cell->hydro.h_max = h_max;
+  cell->hydro.h_max_active = h_max;
   cell->hydro.count = count;
   cell->grav.count = 0;
   cell->hydro.dx_max_part = 0.;
diff --git a/tests/test27cells.c b/tests/test27cells.c
index baafbc67622daffe3a3e5e40889c57952cd96abc..2361a2b1fa5d8b16dda2bbf69b9da4260afbd657 100644
--- a/tests/test27cells.c
+++ b/tests/test27cells.c
@@ -183,6 +183,7 @@ struct cell *make_cell(size_t n, double *offset, double size, double h,
   /* Cell properties */
   cell->split = 0;
   cell->hydro.h_max = h_max;
+  cell->hydro.h_max_active = h_max;
   cell->hydro.count = count;
   cell->hydro.dx_max_part = 0.;
   cell->hydro.dx_max_sort = 0.;
diff --git a/tests/test27cellsStars.c b/tests/test27cellsStars.c
index 235c0e4defbafbf538088fd56c30aed7c3e471ca..f0276075b25074a3370567b0a7dbbd6c7b62f37c 100644
--- a/tests/test27cellsStars.c
+++ b/tests/test27cellsStars.c
@@ -163,6 +163,7 @@ struct cell *make_cell(size_t n, size_t n_stars, double *offset, double size,
   cell->hydro.h_max = h_max;
   cell->hydro.count = count;
   cell->stars.h_max = stars_h_max;
+  cell->stars.h_max_active = stars_h_max;
   cell->stars.count = scount;
   cell->hydro.dx_max_part = 0.;
   cell->hydro.dx_max_sort = 0.;