Make the drift task recurse and updated all its children's ti_old

src/cell.c

@@ -49,6 +49,7 @@
 /* Local headers. */
 #include "active.h"
 #include "atomic.h"
+#include "drift.h"
 #include "error.h"
 #include "gravity.h"
 #include "hydro.h"
@@ -985,3 +986,109 @@ void cell_set_super(struct cell *c, struct cell *super) {
     for (int k = 0; k < 8; k++)
       if (c->progeny[k] != NULL) cell_set_super(c->progeny[k], super);
 }
+
+void cell_drift(struct cell *c, struct engine *e) {
+
+  const double timeBase = e->timeBase;
+  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;
+
+  /* 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;
+  
+  /* Check that we are actually going to move forward. */
+  if (ti_current <= ti_old) return;
+
+  /* Are we in a leaf ? */
+  if (c->split) {
+
+    /* Loop over the progeny and collect their data. */
+    for (int k = 0; k < 8; k++)
+      if (c->progeny[k] != NULL) {
+        struct cell *cp = c->progeny[k];
+	cell_drift(cp, e);
+	dx_max = max(dx_max, cp->dx_max);
+        h_max = max(h_max, cp->h_max);
+      }
+    
+  } else {
+      
+    /* Loop over all the g-particles in the cell */
+    const size_t nr_gparts = c->gcount;
+    for (size_t k = 0; k < nr_gparts; k++) {
+      
+      /* Get a handle on the gpart. */
+      struct gpart *const gp = &gparts[k];
+      
+      /* Drift... */
+      drift_gpart(gp, dt, timeBase, ti_old, ti_current);
+      
+      /* Compute (square of) motion since last cell construction */
+      const float dx2 = gp->x_diff[0] * gp->x_diff[0] +
+	gp->x_diff[1] * gp->x_diff[1] +
+	gp->x_diff[2] * gp->x_diff[2];
+      dx2_max = (dx2_max > dx2) ? dx2_max : dx2;
+    }
+    
+    /* Loop over all the particles in the cell */
+    const size_t nr_parts = c->count;
+    for (size_t k = 0; k < nr_parts; k++) {
+
+      /* Get a handle on the part. */
+      struct part *const p = &parts[k];
+      struct xpart *const xp = &xparts[k];
+
+      /* Drift... */
+      drift_part(p, xp, dt, timeBase, ti_old, ti_current);
+
+      /* Compute (square of) motion since last cell construction */
+      const float dx2 = xp->x_diff[0] * xp->x_diff[0] +
+                        xp->x_diff[1] * xp->x_diff[1] +
+                        xp->x_diff[2] * xp->x_diff[2];
+      dx2_max = (dx2_max > dx2) ? dx2_max : dx2;
+
+      /* Maximal smoothing length */
+      h_max = (h_max > p->h) ? h_max : p->h;
+    }
+
+    /* Now, get the maximal particle motion from its square */
+    dx_max = sqrtf(dx2_max);
+
+    /* Set ti_old on the sub-cells */
+    cell_set_ti_old(c, e->ti_current);
+
+  } /* Check that we are actually going to move forward. */
+
+  /* Store the values */
+  c->h_max = h_max;
+  c->dx_max = dx_max;
+
+  /* Update the time of the last drift */
+  c->ti_old = ti_current;
+}
+
+/**
+ * Set ti_old of a #cell and all its progenies to a new value.
+ *
+ * @param c The #cell.
+ * @param ti_current The new value of ti_old.
+ */
+void cell_set_ti_old(struct cell *c, int ti_current) {
+
+  /* Set this cell */
+  c->ti_old = ti_current;
+
+  /* Recurse */
+  if (c->split) {
+    for (int k = 0; k < 8; ++k) {
+      if (c->progeny[k] != NULL) {
+        struct cell *cp = c->progeny[k];
+        cell_set_ti_old(cp, ti_current);
+      }
+    }
+  }
+}

src/cell.h

@@ -298,5 +298,6 @@ void cell_check_drift_point(struct cell *c, void *data);
 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(struct cell *c, struct engine *e);
+void cell_set_ti_old(struct cell *c, int ti_current);
 #endif /* SWIFT_CELL_H */

src/engine.c

@@ -2702,7 +2702,7 @@ void engine_drift_all(struct engine *e) {
   e->drift_all = 1;
 
   const ticks tic = getticks();
-  threadpool_map(&e->threadpool, runner_do_unskip_mapper, e->s->cells_top,
+  threadpool_map(&e->threadpool, runner_do_drift_mapper, e->s->cells_top,
                  e->s->nr_cells, sizeof(struct cell), 1, e);
 
   e->drift_all = e->policy & engine_policy_drift_all;

src/runner.c

@@ -779,23 +779,23 @@ static void runner_do_unskip(struct cell *c, struct engine *e, int drift) {
 
   /* Now, we can drift */
   /* Get some information first */
-  const double timeBase = e->timeBase;
+  // const double timeBase = e->timeBase;
   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;
+  // struct part *const parts = c->parts;
+  // struct xpart *const xparts = c->xparts;
+  // struct gpart *const gparts = c->gparts;
 
   /* 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;
+  // const double dt = (ti_current - ti_old) * timeBase;
+  //  float dx_max = 0.f, dx2_max = 0.f, h_max = 0.f;
 
   /* No children? */
   if (!c->split) {
 
     /* Check that we are actually going to move forward. */
     if (ti_current > ti_old) {
-
+#if 0
       /* Loop over all the g-particles in the cell */
       const size_t nr_gparts = c->gcount;
       for (size_t k = 0; k < nr_gparts; k++) {
@@ -836,13 +836,15 @@ static void runner_do_unskip(struct cell *c, struct engine *e, int drift) {
 
       /* Now, get the maximal particle motion from its square */
       dx_max = sqrtf(dx2_max);
-
+#endif
     } /* Check that we are actually going to move forward. */
 
     else {
+#if 0
       /* ti_old == ti_current, just keep the current cell values. */
       h_max = c->h_max;
       dx_max = c->dx_max;
+#endif
     }
   }
 
@@ -853,20 +855,23 @@ static void runner_do_unskip(struct cell *c, struct engine *e, int drift) {
     for (int k = 0; k < 8; k++)
       if (c->progeny[k] != NULL) {
         struct cell *cp = c->progeny[k];
-
+	message("aaa");
         /* Recurse. */
         runner_do_unskip(cp, e, drift);
+#if 0
         dx_max = max(dx_max, cp->dx_max);
         h_max = max(h_max, cp->h_max);
+#endif
       }
   }
-
+#if 0
   /* Store the values */
   c->h_max = h_max;
   c->dx_max = dx_max;
 
   /* Update the time of the last drift */
   c->ti_old = ti_current;
+#endif
 }
 
 /**
@@ -876,7 +881,6 @@ static void runner_do_unskip(struct cell *c, struct engine *e, int drift) {
  * @param num_elements Chunk size.
  * @param extra_data Pointer to an #engine.
  */
-
 void runner_do_unskip_mapper(void *map_data, int num_elements,
                              void *extra_data) {
 
@@ -893,7 +897,22 @@ void runner_do_unskip_mapper(void *map_data, int num_elements,
   }
 }
 
-void runner_do_drift(struct runner *r, struct cell *c, int timer) {}
+void runner_do_drift(struct runner *r, struct cell *c, int timer) {
+
+  cell_drift(c, r->e);
+}
+
+void runner_do_drift_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(c, e);
+  }
+}
 
 /**
  * @brief Kick particles in momentum space and collect statistics (floating

src/runner.h

@@ -58,5 +58,6 @@ 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);
 
 #endif /* SWIFT_RUNNER_H */