Call the new spart drifting function when needed.

src/cell.c

@@ -390,6 +390,7 @@ int cell_pack(struct cell *restrict c, struct pcell *restrict pc,
   pc->hydro.ti_old_part = c->hydro.ti_old_part;
   pc->grav.ti_old_part = c->grav.ti_old_part;
   pc->grav.ti_old_multipole = c->grav.ti_old_multipole;
+  pc->stars.ti_old_part = c->stars.ti_old_part;
   pc->hydro.count = c->hydro.count;
   pc->grav.count = c->grav.count;
   pc->stars.count = c->stars.count;
@@ -495,6 +496,7 @@ int cell_unpack(struct pcell *restrict pc, struct cell *restrict c,
   c->hydro.ti_old_part = pc->hydro.ti_old_part;
   c->grav.ti_old_part = pc->grav.ti_old_part;
   c->grav.ti_old_multipole = pc->grav.ti_old_multipole;
+  c->stars.ti_old_part = pc->stars.ti_old_part;
   c->hydro.count = pc->hydro.count;
   c->grav.count = pc->grav.count;
   c->stars.count = pc->stars.count;
@@ -1477,7 +1479,7 @@ void cell_check_part_drift_point(struct cell *c, void *data) {
 }
 
 /**
- * @brief Checks that the #gpart and #spart in a cell are at the
+ * @brief Checks that the #gpart in a cell are at the
  * current point in time
  *
  * Calls error() if the cell is not at the current time.
@@ -1508,11 +1510,42 @@ void cell_check_gpart_drift_point(struct cell *c, void *data) {
         c->grav.parts[i].time_bin != time_bin_inhibited)
       error("g-part in an incorrect time-zone! gp->ti_drift=%lld ti_drift=%lld",
             c->grav.parts[i].ti_drift, ti_drift);
+#else
+  error("Calling debugging code without debugging flag activated.");
+#endif
+}
+
+/**
+ * @brief Checks that the #spart in a cell are at the
+ * current point in time
+ *
+ * Calls error() if the cell is not at the current time.
+ *
+ * @param c Cell to act upon
+ * @param data The current time on the integer time-line
+ */
+void cell_check_spart_drift_point(struct cell *c, void *data) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+
+  const integertime_t ti_drift = *(integertime_t *)data;
+
+  /* Only check local cells */
+  if (c->nodeID != engine_rank) return;
+
+  /* Only check cells with content */
+  if (c->stars.count == 0) return;
+
+  if (c->stars.ti_old_part != ti_drift)
+    error(
+        "Cell in an incorrect time-zone! c->stars.ti_old_part=%lld "
+        "ti_drift=%lld",
+        c->stars.ti_old_part, ti_drift);
 
   for (int i = 0; i < c->stars.count; ++i)
     if (c->stars.parts[i].ti_drift != ti_drift &&
         c->stars.parts[i].time_bin != time_bin_inhibited)
-      error("s-part in an incorrect time-zone! sp->ti_drift=%lld ti_drift=%lld",
+      error("g-part in an incorrect time-zone! gp->ti_drift=%lld ti_drift=%lld",
             c->stars.parts[i].ti_drift, ti_drift);
 #else
   error("Calling debugging code without debugging flag activated.");
@@ -1802,6 +1835,8 @@ void cell_clear_drift_flags(struct cell *c, void *data) {
   c->hydro.do_sub_drift = 0;
   c->grav.do_drift = 0;
   c->grav.do_sub_drift = 0;
+  c->stars.do_drift = 0;
+  c->stars.do_sub_drift = 0;
 }
 
 /**
@@ -1898,8 +1933,39 @@ void cell_activate_drift_gpart(struct cell *c, struct scheduler *s) {
  * @brief Activate the #spart drifts on the given cell.
  */
 void cell_activate_drift_spart(struct cell *c, struct scheduler *s) {
-  // MATTHIEU: This will need changing
-  cell_activate_drift_gpart(c, s);
+
+  /* If this cell is already marked for drift, quit early. */
+  if (c->stars.do_drift) return;
+
+  /* Mark this cell for drifting. */
+  c->stars.do_drift = 1;
+
+  /* Set the do_stars_sub_drifts all the way up and activate the super drift
+     if this has not yet been done. */
+  if (c == c->hydro.super) {
+#ifdef SWIFT_DEBUG_CHECKS
+    if (c->stars.drift == NULL)
+      error("Trying to activate un-existing c->stars.drift");
+#endif
+    scheduler_activate(s, c->stars.drift);
+  } else {
+    for (struct cell *parent = c->parent;
+         parent != NULL && !parent->stars.do_sub_drift;
+         parent = parent->parent) {
+
+      /* Mark this cell for drifting */
+      parent->stars.do_sub_drift = 1;
+
+      if (parent == c->hydro.super) {
+#ifdef SWIFT_DEBUG_CHECKS
+        if (parent->stars.drift == NULL)
+          error("Trying to activate un-existing parent->stars.drift");
+#endif
+        scheduler_activate(s, parent->stars.drift);
+        break;
+      }
+    }
+  }
 }
 
 /**
@@ -3995,7 +4061,7 @@ void cell_drift_spart(struct cell *c, const struct engine *e, int force) {
         struct cell *cp = c->progeny[k];
 
         /* Recurse */
-        cell_drift_gpart(cp, e, force);
+        cell_drift_spart(cp, e, force);
 
         /* Update */
         dx_max = max(dx_max, cp->stars.dx_max_part);

src/cell.h

@@ -150,6 +150,9 @@ struct pcell {
     /*! Minimal integer end-of-timestep in this cell for stars tasks */
     integertime_t ti_end_min;
 
+    /*! Integer time of the last drift of the #spart in this cell */
+    integertime_t ti_old_part;
+
   } stars;
 
   /*! Maximal depth in that part of the tree */
@@ -721,6 +724,7 @@ void cell_check_foreign_multipole(const struct cell *c);
 void cell_clean(struct cell *c);
 void cell_check_part_drift_point(struct cell *c, void *data);
 void cell_check_gpart_drift_point(struct cell *c, void *data);
+void cell_check_spart_drift_point(struct cell *c, void *data);
 void cell_check_multipole_drift_point(struct cell *c, void *data);
 void cell_reset_task_counters(struct cell *c);
 int cell_unskip_hydro_tasks(struct cell *c, struct scheduler *s);

src/engine_drift.c

@@ -124,6 +124,54 @@ void engine_do_drift_all_gpart_mapper(void *map_data, int num_elements,
   }
 }
 
+/**
+ * @brief Mapper function to drift *all* the #spart to the current time.
+ *
+ * @param map_data An array of #cell%s.
+ * @param num_elements Chunk size.
+ * @param extra_data Pointer to an #engine.
+ */
+void engine_do_drift_all_spart_mapper(void *map_data, int num_elements,
+                                      void *extra_data) {
+
+  const struct engine *e = (const struct engine *)extra_data;
+  const int restarting = e->restarting;
+  struct space *s = e->s;
+  struct cell *cells_top;
+  int *local_cells_top;
+
+  if (restarting) {
+
+    /* When restarting, we loop over all top-level cells */
+    cells_top = (struct cell *)map_data;
+    local_cells_top = NULL;
+
+  } else {
+
+    /* In any other case, we use the list of local cells with tasks */
+    cells_top = s->cells_top;
+    local_cells_top = (int *)map_data;
+  }
+
+  for (int ind = 0; ind < num_elements; ind++) {
+
+    struct cell *c;
+
+    /* When restarting, the list of local cells does not
+       yet exist. We use the raw list of top-level cells instead */
+    if (restarting)
+      c = &cells_top[ind];
+    else
+      c = &cells_top[local_cells_top[ind]];
+
+    if (c->nodeID == e->nodeID) {
+
+      /* Drift all the particles */
+      cell_drift_spart(c, e, /* force the drift=*/1);
+    }
+  }
+}
+
 /**
  * @brief Mapper function to drift *all* the multipoles to the current time.
  *
@@ -204,6 +252,11 @@ void engine_drift_all(struct engine *e, const int drift_mpoles) {
                      e->s->local_cells_top, e->s->nr_local_cells, sizeof(int),
                      /* default chunk */ 0, e);
     }
+    if (e->s->nr_sparts > 0) {
+      threadpool_map(&e->threadpool, engine_do_drift_all_spart_mapper,
+                     e->s->local_cells_top, e->s->nr_local_cells, sizeof(int),
+                     /* default chunk */ 0, e);
+    }
     if (drift_mpoles && (e->policy & engine_policy_self_gravity)) {
       threadpool_map(&e->threadpool, engine_do_drift_all_multipole_mapper,
                      e->s->local_cells_with_tasks_top,

src/space.c

@@ -242,8 +242,9 @@ void space_rebuild_recycle_mapper(void *map_data, int num_elements,
     c->stars.parts = NULL;
     c->hydro.do_sub_sort = 0;
     c->stars.do_sub_sort = 0;
-    c->grav.do_sub_drift = 0;
     c->hydro.do_sub_drift = 0;
+    c->grav.do_sub_drift = 0;
+    c->stars.do_sub_drift = 0;
     c->hydro.do_sub_limiter = 0;
     c->hydro.do_limiter = 0;
     c->hydro.ti_end_min = -1;
@@ -540,6 +541,7 @@ void space_regrid(struct space *s, int verbose) {
           c->grav.super = c;
           c->hydro.ti_old_part = ti_current;
           c->grav.ti_old_part = ti_current;
+          c->stars.ti_old_part = ti_current;
           c->grav.ti_old_multipole = ti_current;
 #ifdef WITH_MPI
           c->mpi.tag = -1;
@@ -1529,6 +1531,7 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
     c->hydro.ti_old_part = ti_current;
     c->grav.ti_old_part = ti_current;
     c->grav.ti_old_multipole = ti_current;
+    c->stars.ti_old_part = ti_current;
 
 #ifdef SWIFT_DEBUG_CHECKS
     c->cellID = -last_cell_id;
@@ -2688,6 +2691,7 @@ void space_split_recursive(struct space *s, struct cell *c,
       cp->hydro.ti_old_part = c->hydro.ti_old_part;
       cp->grav.ti_old_part = c->grav.ti_old_part;
       cp->grav.ti_old_multipole = c->grav.ti_old_multipole;
+      cp->stars.ti_old_part = c->stars.ti_old_part;
       cp->loc[0] = c->loc[0];
       cp->loc[1] = c->loc[1];
       cp->loc[2] = c->loc[2];
@@ -2715,6 +2719,7 @@ 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->stars.do_sub_drift = 0;
       cp->hydro.do_sub_limiter = 0;
       cp->hydro.do_limiter = 0;
 #ifdef WITH_MPI
@@ -4292,6 +4297,7 @@ void space_check_drift_point(struct space *s, integertime_t ti_drift,
   /* Recursively check all cells */
   space_map_cells_pre(s, 1, cell_check_part_drift_point, &ti_drift);
   space_map_cells_pre(s, 1, cell_check_gpart_drift_point, &ti_drift);
+  space_map_cells_pre(s, 1, cell_check_spart_drift_point, &ti_drift);
   if (multipole)
     space_map_cells_pre(s, 1, cell_check_multipole_drift_point, &ti_drift);
 #else