Drift all the top-level multipoles to current time before launching the tasks.

src/cell.c

@@ -1017,14 +1017,15 @@ void cell_clean_links(struct cell *c, void *data) {
 }
 
 /**
- * @brief Checks that a cell is at the current point in time
+ * @brief Checks that the particles 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_drift_point(struct cell *c, void *data) {
+void cell_check_particle_drift_point(struct cell *c, void *data) {
 
 #ifdef SWIFT_DEBUG_CHECKS
 
@@ -1056,6 +1057,34 @@ void cell_check_drift_point(struct cell *c, void *data) {
 #endif
 }
 
+/**
+ * @brief Checks that the multipole of a cell is 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_multipole_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;
+
+  if (c->ti_old_multipole != ti_drift)
+    error(
+        "Cell multipole in an incorrect time-zone! c->ti_old_multipole=%lld "
+        "ti_drift=%lld",
+        c->ti_old_multipole, ti_drift);
+
+#else
+  error("Calling debugging code without debugging flag activated.");
+#endif
+}
+
 /**
  * @brief Resets all the individual cell task counters to 0.
  *

src/cell.h

@@ -352,7 +352,8 @@ int cell_are_neighbours(const struct cell *restrict ci,
                         const struct cell *restrict cj);
 void cell_check_multipole(struct cell *c, void *data);
 void cell_clean(struct cell *c);
-void cell_check_drift_point(struct cell *c, void *data);
+void cell_check_particle_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_is_drift_needed(struct cell *c, const struct engine *e);
 int cell_unskip_tasks(struct cell *c, struct scheduler *s);

src/engine.c

@@ -846,7 +846,8 @@ void engine_repartition(struct engine *e) {
   fflush(stdout);
 
   /* Check that all cells have been drifted to the current time */
-  space_check_drift_point(e->s, e->ti_old);
+  space_check_drift_point(e->s, e->ti_old,
+                          e->policy & engine_policy_self_gravity);
 #endif
 
   /* Clear the repartition flag. */
@@ -2664,7 +2665,8 @@ void engine_rebuild(struct engine *e) {
   /* Check that all cells have been drifted to the current time.
    * That can include cells that have not
    * previously been active on this rank. */
-  space_check_drift_point(e->s, e->ti_old);
+  space_check_drift_point(e->s, e->ti_old,
+                          e->policy & engine_policy_self_gravity);
 #endif
 
   if (e->verbose)
@@ -2686,7 +2688,8 @@ void engine_prepare(struct engine *e) {
     /* Check that all cells have been drifted to the current time.
      * That can include cells that have not
      * previously been active on this rank. */
-    space_check_drift_point(e->s, e->ti_old);
+    space_check_drift_point(e->s, e->ti_old,
+                            e->policy & engine_policy_self_gravity);
   }
 #endif
 
@@ -2887,7 +2890,8 @@ void engine_print_stats(struct engine *e) {
   /* Check that all cells have been drifted to the current time.
    * That can include cells that have not
    * previously been active on this rank. */
-  space_check_drift_point(e->s, e->ti_current);
+  space_check_drift_point(e->s, e->ti_current,
+                          e->policy & engine_policy_self_gravity);
 
   /* Be verbose about this */
   message("Saving statistics at t=%e.", e->time);
@@ -3195,6 +3199,9 @@ void engine_step(struct engine *e) {
   gravity_exact_force_compute(e->s, e);
 #endif
 
+  /* Do we need to drift the top-level multipoles ? */
+  if (e->policy & engine_policy_self_gravity) engine_drift_top_multipoles(e);
+
   /* Start all the tasks. */
   TIMER_TIC;
   engine_launch(e, e->nr_threads);
@@ -3286,8 +3293,8 @@ void engine_unskip(struct engine *e) {
 }
 
 /**
- * @brief Mapper function to drift ALL particle types and multipoles forward in
- * time.
+ * @brief Mapper function to drift *all* particle types and multipoles forward
+ * in time.
  *
  * @param map_data An array of #cell%s.
  * @param num_elements Chunk size.
@@ -3313,7 +3320,8 @@ void engine_do_drift_all_mapper(void *map_data, int num_elements,
 }
 
 /**
- * @brief Drift *all* particles forward to the current time.
+ * @brief Drift *all* particles and multipoles at all levels
+ * forward to the current time.
  *
  * @param e The #engine.
  */
@@ -3330,7 +3338,54 @@ void engine_drift_all(struct engine *e) {
 
 #ifdef SWIFT_DEBUG_CHECKS
   /* Check that all cells have been drifted to the current time. */
-  space_check_drift_point(e->s, e->ti_current);
+  space_check_drift_point(e->s, e->ti_current,
+                          e->policy & engine_policy_self_gravity);
+#endif
+
+  if (e->verbose)
+    message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
+            clocks_getunit());
+}
+
+/**
+ * @brief Mapper function to drift *all* top-level multipoles forward in
+ * 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_top_multipoles_mapper(void *map_data, int num_elements,
+                                           void *extra_data) {
+
+  struct engine *e = (struct engine *)extra_data;
+  struct cell *cells = (struct cell *)map_data;
+
+  for (int ind = 0; ind < num_elements; ind++) {
+    struct cell *c = &cells[ind];
+    if (c != NULL && c->nodeID == e->nodeID) {
+
+      /* Drift the multipole at this level only */
+      cell_drift_multipole(c, e);
+    }
+  }
+}
+
+/**
+ * @brief Drift *all* top-level multipoles forward to the current time.
+ *
+ * @param e The #engine.
+ */
+void engine_drift_top_multipoles(struct engine *e) {
+
+  const ticks tic = getticks();
+
+  threadpool_map(&e->threadpool, engine_do_drift_top_multipoles_mapper,
+                 e->s->cells_top, e->s->nr_cells, sizeof(struct cell), 10, e);
+
+#ifdef SWIFT_DEBUG_CHECKS
+  /* Check that all cells have been drifted to the current time. */
+  space_check_top_multipoles_drift_point(e->s, e->ti_current);
 #endif
 
   if (e->verbose)
@@ -3546,7 +3601,8 @@ void engine_dump_snapshot(struct engine *e) {
   /* Check that all cells have been drifted to the current time.
    * That can include cells that have not
    * previously been active on this rank. */
-  space_check_drift_point(e->s, e->ti_current);
+  space_check_drift_point(e->s, e->ti_current,
+                          e->policy & engine_policy_self_gravity);
 
   /* Be verbose about this */
   message("writing snapshot at t=%e.", e->time);

src/engine.h

@@ -250,6 +250,7 @@ void engine_barrier(struct engine *e, int tid);
 void engine_compute_next_snapshot_time(struct engine *e);
 void engine_unskip(struct engine *e);
 void engine_drift_all(struct engine *e);
+void engine_drift_top_multipoles(struct engine *e);
 void engine_dump_snapshot(struct engine *e);
 void engine_init(struct engine *e, struct space *s,
                  const struct swift_params *params, int nr_nodes, int nodeID,

src/space.c

@@ -2819,11 +2819,26 @@ void space_link_cleanup(struct space *s) {
  *
  * @param s The #space to check.
  * @param ti_drift The (integer) time.
+ * @param multipole Are we also checking the multipoles ?
  */
-void space_check_drift_point(struct space *s, integertime_t ti_drift) {
+void space_check_drift_point(struct space *s, integertime_t ti_drift,
+                             int multipole) {
 #ifdef SWIFT_DEBUG_CHECKS
   /* Recursively check all cells */
-  space_map_cells_pre(s, 1, cell_check_drift_point, &ti_drift);
+  space_map_cells_pre(s, 1, cell_check_particle_drift_point, &ti_drift);
+  if (multipole)
+    space_map_cells_pre(s, 1, cell_check_multipole_drift_point, &ti_drift);
+#else
+  error("Calling debugging code without debugging flag activated.");
+#endif
+}
+
+void space_check_top_multipoles_drift_point(struct space *s,
+                                            integertime_t ti_drift) {
+#ifdef SWIFT_DEBUG_CHECKS
+  for (int i = 0; i < s->nr_cells; ++i) {
+    cell_check_multipole_drift_point(&s->cells_top[i], &ti_drift);
+  }
 #else
   error("Calling debugging code without debugging flag activated.");
 #endif

src/space.h

@@ -213,7 +213,10 @@ void space_init_parts(struct space *s);
 void space_init_gparts(struct space *s);
 void space_init_sparts(struct space *s);
 void space_link_cleanup(struct space *s);
-void space_check_drift_point(struct space *s, integertime_t ti_drift);
+void space_check_drift_point(struct space *s, integertime_t ti_drift,
+                             int multipole);
+void space_check_top_multipoles_drift_point(struct space *s,
+                                            integertime_t ti_drift);
 void space_check_timesteps(struct space *s);
 void space_replicate(struct space *s, int replicate, int verbose);
 void space_reset_task_counters(struct space *s);