Only drift the particles that need to be.

src/cell.c

@@ -262,7 +262,7 @@ int cell_pack_ti_ends(struct cell *c, int *ti_ends) {
  * @brief Unpack the time information of a given cell and its sub-cells.
  *
  * @param c The #cell
- * @param s ti_ends The time information to unpack
+ * @param ti_ends The time information to unpack
  *
  * @return The number of cells created.
  */
@@ -794,9 +794,9 @@ void cell_check_multipole(struct cell *c, void *data) {
 }
 
 /**
- * @brief Frees up the memory allocated for this #cell
+ * @brief Frees up the memory allocated for this #cell.
  *
- * @param c The #cell
+ * @param c The #cell.
  */
 void cell_clean(struct cell *c) {
 
@@ -806,3 +806,34 @@ void cell_clean(struct cell *c) {
   for (int k = 0; k < 8; k++)
     if (c->progeny[k]) cell_clean(c->progeny[k]);
 }
+
+/**
+ * @brief Checks whether a given cell needs drifting or not.
+ *
+ * @param c the #cell.
+ * @param ti_current The current time on the integer time-line.
+ *
+ * @return 1 If the cell needs drifting, 0 otherwise.
+ */
+int cell_is_drift_needed(struct cell *c, int ti_current) {
+
+  /* Do we have at least one active particle in the cell ?*/
+  if (c->ti_end_min == ti_current) return 1;
+
+  /* Loop over the pair tasks that involve this cell */
+  for (struct link *l = c->density; l != NULL; l = l->next) {
+
+    if (l->t->type != task_type_pair && l->t->type != task_type_sub_pair)
+      continue;
+
+    /* Left or right? */
+    if (l->t->ci == c) {
+      if (l->t->cj->ti_end_min == ti_current) return 1;
+    } else if (l->t->cj == c) {
+      if (l->t->ci->ti_end_min == ti_current) return 1;
+    }
+  }
+
+  /* No neighbouring cell has active particles. Drift not necessary */
+  return 0;
+}

src/cell.h

@@ -47,6 +47,16 @@ struct space;
 /* Global variables. */
 extern int cell_next_tag;
 
+/* Mini struct to link cells to tasks. Used as a linked list. */
+struct link {
+
+  /* The task pointer. */
+  struct task *t;
+
+  /* The next pointer. */
+  struct link *next;
+};
+
 /* Packed cell. */
 struct pcell {
 
@@ -79,6 +89,9 @@ struct cell {
   /* Minimum and maximum end of time step in this cell. */
   int ti_end_min, ti_end_max;
 
+  /* Last time the cell's content was drifted forward in time. */
+  int ti_old;
+
   /* Minimum dimension, i.e. smallest edge of this cell. */
   float dmin;
 
@@ -208,5 +221,6 @@ 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);
+int cell_is_drift_needed(struct cell *c, int ti_current);
 
 #endif /* SWIFT_CELL_H */

src/engine.c

@@ -2088,6 +2088,14 @@ void engine_prepare(struct engine *e) {
 
   /* Did this not go through? */
   if (rebuild) {
+
+    /* First drift all particles to the current time */
+    e->drift_all = 1;
+    threadpool_map(&e->threadpool, runner_do_drift_mapper, e->s->cells,
+                   e->s->nr_cells, sizeof(struct cell), 1, e);
+    e->drift_all = 0;
+
+    /* And now rebuild */
     engine_rebuild(e);
   }
 
@@ -2468,8 +2476,10 @@ void engine_step(struct engine *e) {
     snapshot_drift_time = e->timeStep;
 
     /* Drift everybody to the snapshot position */
+    e->drift_all = 1;
     threadpool_map(&e->threadpool, runner_do_drift_mapper, e->s->cells,
                    e->s->nr_cells, sizeof(struct cell), 1, e);
+    e->drift_all = 0;
 
     /* Dump... */
     engine_dump_snapshot(e);
@@ -2486,7 +2496,7 @@ void engine_step(struct engine *e) {
   e->timeOld = e->ti_old * e->timeBase + e->timeBegin;
   e->timeStep = (e->ti_current - e->ti_old) * e->timeBase + snapshot_drift_time;
 
-  /* Drift everybody */
+  /* Drift only the necessary particles */
   threadpool_map(&e->threadpool, runner_do_drift_mapper, e->s->cells,
                  e->s->nr_cells, sizeof(struct cell), 1, e);
 
@@ -2915,6 +2925,7 @@ void engine_init(struct engine *e, struct space *s,
   e->timeStep = 0.;
   e->timeBase = 0.;
   e->timeBase_inv = 0.;
+  e->drift_all = 0;
   e->internalUnits = internal_units;
   e->timeFirstSnapshot =
       parser_get_param_double(params, "Snapshots:time_first");

src/engine.h

@@ -81,16 +81,6 @@ extern int engine_rank;
 /* The maximal number of timesteps in a simulation */
 #define max_nr_timesteps (1 << 28)
 
-/* Mini struct to link cells to density/force tasks. */
-struct link {
-
-  /* The task pointer. */
-  struct task *t;
-
-  /* The next pointer. */
-  struct link *next;
-};
-
 /* Data structure for the engine. */
 struct engine {
 
@@ -139,6 +129,9 @@ struct engine {
   /* Minimal ti_end for the next time-step */
   int ti_end_min;
 
+  /* Are we drifting all particles now ? */
+  int drift_all;
+
   /* Number of particles updated */
   size_t updates, g_updates;
 

src/runner.c

@@ -592,15 +592,19 @@ void runner_do_ghost(struct runner *r, struct cell *c) {
 static void runner_do_drift(struct cell *c, struct engine *e) {
 
   const double timeBase = e->timeBase;
-  const double dt = (e->ti_current - e->ti_old) * timeBase;
-  const int ti_old = e->ti_old;
+  const int ti_old = c->ti_old;
   const int ti_current = e->ti_current;
-
   struct part *const parts = c->parts;
   struct xpart *const xparts = c->xparts;
   struct gpart *const gparts = c->gparts;
-  float dx_max = 0.f, dx2_max = 0.f, h_max = 0.f;
 
+  /* Do we need to drift ? */
+  if (!e->drift_all && !cell_is_drift_needed(c, ti_current)) return;
+
+  /* Drift from the last time the cell was drifted to the current time */
+  const double dt = (ti_current - ti_old) * timeBase;
+
+  float dx_max = 0.f, dx2_max = 0.f, h_max = 0.f;
   double e_kin = 0.0, e_int = 0.0, e_pot = 0.0, entropy = 0.0, mass = 0.0;
   double mom[3] = {0.0, 0.0, 0.0};
   double ang_mom[3] = {0.0, 0.0, 0.0};
@@ -722,6 +726,9 @@ static void runner_do_drift(struct cell *c, struct engine *e) {
   c->ang_mom[0] = ang_mom[0];
   c->ang_mom[1] = ang_mom[1];
   c->ang_mom[2] = ang_mom[2];
+
+  /* Update the time of the last drift */
+  c->ti_old = ti_current;
 }
 
 /**

src/space.c

@@ -164,6 +164,7 @@ void space_regrid(struct space *s, double cell_max, int verbose) {
   const size_t nr_parts = s->nr_parts;
   struct cell *restrict c;
   ticks tic = getticks();
+  const int ti_current = (s->e != NULL) ? s->e->ti_current : 0;
 
   /* Run through the cells and get the current h_max. */
   // tic = getticks();
@@ -294,6 +295,7 @@ void space_regrid(struct space *s, double cell_max, int verbose) {
           c->count = 0;
           c->gcount = 0;
           c->super = c;
+          c->ti_old = ti_current;
           lock_init(&c->lock);
         }
 
@@ -392,6 +394,7 @@ void space_rebuild(struct space *s, double cell_max, int verbose) {
   size_t nr_parts = s->nr_parts;
   size_t nr_gparts = s->nr_gparts;
   struct cell *restrict cells = s->cells;
+  const int ti_current = (s->e != NULL) ? s->e->ti_current : 0;
 
   const double ih[3] = {s->iwidth[0], s->iwidth[1], s->iwidth[2]};
   const double dim[3] = {s->dim[0], s->dim[1], s->dim[2]};
@@ -643,6 +646,7 @@ void space_rebuild(struct space *s, double cell_max, int verbose) {
   struct gpart *gfinger = s->gparts;
   for (int k = 0; k < s->nr_cells; k++) {
     struct cell *restrict c = &cells[k];
+    c->ti_old = ti_current;
     c->parts = finger;
     c->xparts = xfinger;
     c->gparts = gfinger;
@@ -1234,6 +1238,7 @@ void space_split_mapper(void *map_data, int num_elements, void *extra_data) {
   /* Unpack the inputs. */
   struct space *s = (struct space *)extra_data;
   struct cell *cells = (struct cell *)map_data;
+  struct engine *e = s->e;
 
   for (int ind = 0; ind < num_elements; ind++) {
 
@@ -1263,6 +1268,7 @@ void space_split_mapper(void *map_data, int num_elements, void *extra_data) {
         temp = space_getcell(s);
         temp->count = 0;
         temp->gcount = 0;
+        temp->ti_old = e->ti_current;
         temp->loc[0] = c->loc[0];
         temp->loc[1] = c->loc[1];
         temp->loc[2] = c->loc[2];