Drift multipoles independently of the g-particles.

src/cell.c

@@ -1326,7 +1326,7 @@ void cell_set_super(struct cell *c, struct cell *super) {
 }
 
 /**
- * @brief Recursively drifts all particles and g-particles in a cell hierarchy.
+ * @brief Recursively drifts particles of all kinds in a cell hierarchy.
  *
  * @param c The #cell.
  * @param e The #engine (to get ti_current).
@@ -1430,6 +1430,41 @@ void cell_drift_particles(struct cell *c, const struct engine *e) {
   c->ti_old = ti_current;
 }
 
+/**
+ * @brief Recursively drifts all multipoles in a cell hierarchy.
+ *
+ * @param c The #cell.
+ * @param e The #engine (to get ti_current).
+ */
+void cell_drift_multipole(struct cell *c, const struct engine *e) {
+
+  const double timeBase = e->timeBase;
+  const integertime_t ti_old_multipole = c->ti_old_multipole;
+  const integertime_t ti_current = e->ti_current;
+
+  /* Drift from the last time the cell was drifted to the current time */
+  const double dt = (ti_current - ti_old_multipole) * timeBase;
+
+  /* Check that we are actually going to move forward. */
+  if (ti_current < ti_old_multipole) error("Attempt to drift to the past");
+
+  /* Are we not in a leaf ? */
+  if (c->split) {
+
+    /* Loop over the progeny and drift the multipoles. */
+    for (int k = 0; k < 8; k++)
+      if (c->progeny[k] != NULL) cell_drift_particles(c->progeny[k], e);
+
+  } else if (ti_current > ti_old_multipole) {
+
+    /* Drift the multipole */
+    multipole_drift(c->multipole, dt);
+  }
+
+  /* Update the time of the last drift */
+  c->ti_old_multipole = ti_current;
+}
+
 /**
  * @brief Recursively checks that all particles in a cell have a time-step
  */

src/cell.h

@@ -230,9 +230,12 @@ struct cell {
   /*! Maximum beginning of (integer) time step in this cell. */
   integertime_t ti_beg_max;
 
-  /*! Last (integer) time the cell's content was drifted forward in time. */
+  /*! Last (integer) time the cell's particle was drifted forward in time. */
   integertime_t ti_old;
 
+  /*! Last (integer) time the cell's multipole was drifted forward in time. */
+  integertime_t ti_old_multipole;
+
   /*! Minimum dimension, i.e. smallest edge of this cell (min(width)). */
   float dmin;
 
@@ -343,6 +346,7 @@ int cell_is_drift_needed(struct cell *c, const struct engine *e);
 int cell_unskip_tasks(struct cell *c, struct scheduler *s);
 void cell_set_super(struct cell *c, struct cell *super);
 void cell_drift_particles(struct cell *c, const struct engine *e);
+void cell_drift_multipole(struct cell *c, const struct engine *e);
 void cell_check_timesteps(struct cell *c);
 
 #endif /* SWIFT_CELL_H */

src/engine.c

@@ -3144,7 +3144,7 @@ void engine_unskip(struct engine *e) {
 void engine_drift_all(struct engine *e) {
 
   const ticks tic = getticks();
-  threadpool_map(&e->threadpool, runner_do_drift_mapper, e->s->cells_top,
+  threadpool_map(&e->threadpool, runner_do_drift_all_mapper, e->s->cells_top,
                  e->s->nr_cells, sizeof(struct cell), 1, e);
 
 #ifdef SWIFT_DEBUG_CHECKS

src/multipole.h

@@ -185,6 +185,20 @@ INLINE static int multipole_equal(const struct multipole *ma,
   return 1;
 }
 
+/**
+ * @brief Drifts a #multipole forward in time.
+ *
+ * @param m The #multipole.
+ * @param dt The drift time-step.
+ */
+INLINE static void multipole_drift(struct multipole *m, double dt) {
+
+  /* Move the whole thing according to bulk motion */
+  m->CoM[0] += m->vel[0];
+  m->CoM[1] += m->vel[1];
+  m->CoM[2] += m->vel[2];
+}
+
 /**
  * @brief Applies the forces due to particles j onto particles i directly.
  *

src/runner.c

@@ -826,15 +826,21 @@ void runner_do_drift_particles(struct runner *r, struct cell *c, int timer) {
  * @param num_elements Chunk size.
  * @param extra_data Pointer to an #engine.
  */
-void runner_do_drift_mapper(void *map_data, int num_elements,
-                            void *extra_data) {
+void runner_do_drift_all_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) cell_drift_particles(c, e);
+    if (c != NULL && c->nodeID == e->nodeID) {
+      /* Drift all the particles */
+      cell_drift_particles(c, e);
+
+      /* Drift the multipole */
+      if (e->policy & engine_policy_self_gravity) cell_drift_multipole(c, e);
+    }
   }
 }
 

src/runner.h

@@ -66,6 +66,7 @@ void runner_do_grav_external(struct runner *r, struct cell *c, int timer);
 void *runner_main(void *data);
 void runner_do_unskip_mapper(void *map_data, int num_elements,
                              void *extra_data);
-void runner_do_drift_mapper(void *map_data, int num_elements, void *extra_data);
+void runner_do_drift_all_mapper(void *map_data, int num_elements,
+                                void *extra_data);
 
 #endif /* SWIFT_RUNNER_H */

src/space.c

@@ -441,6 +441,7 @@ void space_regrid(struct space *s, int verbose) {
           c->scount = 0;
           c->super = c;
           c->ti_old = ti_old;
+          c->ti_old_multipole = ti_old;
           if (s->gravity) c->multipole = &s->multipoles_top[cid];
         }
 
@@ -922,6 +923,7 @@ void space_rebuild(struct space *s, int verbose) {
   for (int k = 0; k < s->nr_cells; k++) {
     struct cell *restrict c = &cells_top[k];
     c->ti_old = ti_old;
+    c->ti_old_multipole = ti_old;
     c->parts = finger;
     c->xparts = xfinger;
     c->gparts = gfinger;
@@ -2030,6 +2032,7 @@ void space_split_recursive(struct space *s, struct cell *c,
       cp->gcount = 0;
       cp->scount = 0;
       cp->ti_old = c->ti_old;
+      cp->ti_old_multipole = c->ti_old_multipole;
       cp->loc[0] = c->loc[0];
       cp->loc[1] = c->loc[1];
       cp->loc[2] = c->loc[2];