Tighter condition for active particles/cells. Also inline the call to cell_is_drift_needed()

src/active.h

@@ -24,6 +24,7 @@
 
 /* Local includes. */
 #include "cell.h"
+#include "const.h"
 #include "engine.h"
 #include "part.h"
 
@@ -36,7 +37,13 @@
 __attribute__((always_inline)) INLINE static int cell_is_active(
     const struct cell *c, const struct engine *e) {
 
-  return (c->ti_end_min <= e->ti_current);
+#ifdef SWIFT_DEBUG_CHECKS
+  if (c->ti_end_min < e->ti_current)
+    error("cell in an impossible time-zone! c->ti_end_min=%d e->ti_current=%d",
+          c->ti_end_min, e->ti_current);
+#endif
+
+  return (c->ti_end_min == e->ti_current);
 }
 
 /**
@@ -48,9 +55,46 @@ __attribute__((always_inline)) INLINE static int cell_is_active(
 __attribute__((always_inline)) INLINE static int cell_is_all_active(
     const struct cell *c, const struct engine *e) {
 
-  return (c->ti_end_max <= e->ti_current);
+#ifdef SWIFT_DEBUG_CHECKS
+  if (c->ti_end_max < e->ti_current)
+    error("cell in an impossible time-zone! c->ti_end_max=%d e->ti_current=%d",
+          c->ti_end_max, e->ti_current);
+#endif
+
+  return (c->ti_end_max == e->ti_current);
 }
 
+/**
+ * @brief Checks whether a given cell needs drifting or not.
+ *
+ * @param c the #cell.
+ * @param e The #engine (holding current time information).
+ *
+ * @return 1 If the cell needs drifting, 0 otherwise.
+ */
+INLINE static int cell_is_drift_needed(struct cell *c, const struct engine *e) {
+
+  /* Do we have at least one active particle in the cell ?*/
+  if (cell_is_active(c, e)) 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;
+
+    /* Is the other cell in the pair active ? */
+    if ((l->t->ci == c && cell_is_active(l->t->cj, e)) ||
+        (l->t->cj == c && cell_is_active(l->t->ci, e)))
+      return 1;
+  }
+
+  /* No neighbouring cell has active particles. Drift not necessary */
+  return 0;
+}
+
+
+
 /**
  * @brief Is this particle active ?
  *
@@ -60,7 +104,13 @@ __attribute__((always_inline)) INLINE static int cell_is_all_active(
 __attribute__((always_inline)) INLINE static int part_is_active(
     const struct part *p, const struct engine *e) {
 
-  return (p->ti_end <= e->ti_current);
+#ifdef SWIFT_DEBUG_CHECKS
+  if (p->ti_end < e->ti_current)
+    error("particle in an impossible time-zone! p->ti_end=%d e->ti_current=%d",
+          p->ti_end, e->ti_current);
+#endif
+
+  return (p->ti_end == e->ti_current);
 }
 
 /**
@@ -72,7 +122,14 @@ __attribute__((always_inline)) INLINE static int part_is_active(
 __attribute__((always_inline)) INLINE static int gpart_is_active(
     const struct gpart *gp, const struct engine *e) {
 
-  return (gp->ti_end <= e->ti_current);
+#ifdef SWIFT_DEBUG_CHECKS
+  if (gp->ti_end < e->ti_current)
+    error(
+        "g-particle in an impossible time-zone! gp->ti_end=%d e->ti_current=%d",
+        gp->ti_end, e->ti_current);
+#endif
+
+  return (gp->ti_end == e->ti_current);
 }
 
 #endif /* SWIFT_ACTIVE_H */

src/cell.c

@@ -863,35 +863,6 @@ void cell_clean(struct cell *c) {
     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;
-
-    /* Does the other cell in the pair have an active particle ? */
-    if ((l->t->ci == c && l->t->cj->ti_end_min == ti_current) ||
-        (l->t->cj == c && l->t->ci->ti_end_min == ti_current))
-      return 1;
-  }
-
-  /* No neighbouring cell has active particles. Drift not necessary */
-  return 0;
-}
-
 /**
  * @brief Un-skips all the tasks associated with a given cell and checks
  * if the space needs to be rebuilt.

src/cell.h

@@ -290,7 +290,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);
 int cell_unskip_tasks(struct cell *c, struct scheduler *s);
 void cell_set_super(struct cell *c, struct cell *super);
 

src/runner.c

@@ -758,20 +758,18 @@ void runner_do_ghost(struct runner *r, struct cell *c) {
  */
 static void runner_do_drift(struct cell *c, struct engine *e, int drift) {
 
-  const int ti_current = e->ti_current;
-
   /* Unskip any active tasks. */
-  if (c->ti_end_min == e->ti_current) {
+  if (cell_is_active(c, e)) {
     const int forcerebuild = cell_unskip_tasks(c, &e->sched);
     if (forcerebuild) atomic_inc(&e->forcerebuild);
   }
 
   /* Do we really need to drift? */
   if (drift) {
-    if (!e->drift_all && !cell_is_drift_needed(c, ti_current)) return;
+    if (!e->drift_all && !cell_is_drift_needed(c, e)) return;
   } else {
 
-    /* Not drifting, but may still need to recurse for task skipping. */
+    /* Not drifting, but may still need to recurse for task un-skipping. */
     if (c->split) {
       for (int k = 0; k < 8; k++) {
         if (c->progeny[k] != NULL) {
@@ -783,8 +781,12 @@ static void runner_do_drift(struct cell *c, struct engine *e, int drift) {
     return;
   }
 
+  /* Now, we can drift */
+
+  /* Get some information first */
   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;
@@ -797,7 +799,7 @@ static void runner_do_drift(struct cell *c, struct engine *e, int drift) {
   if (!c->split) {
 
     /* Check that we are actually going to move forward. */
-    if (ti_current >= ti_old) {
+    if (ti_current > ti_old) {
 
       /* Loop over all the g-particles in the cell */
       const size_t nr_gparts = c->gcount;

src/statistics.c

@@ -157,7 +157,7 @@ void stats_collect_part_mapper(void *map_data, int nr_parts, void *extra_data) {
     stats.E_pot_self += 0.f;
     if (gp != NULL)
       stats.E_pot_ext +=
-        m * external_gravity_get_potential_energy(potential, phys_const, gp);
+          m * external_gravity_get_potential_energy(potential, phys_const, gp);
     stats.E_int += m * hydro_get_internal_energy(p, dt);
     stats.E_rad += cooling_get_radiated_energy(xp);