Merge branch 'gravity_speedup' into 'master'

Gravity speedup



See merge request !420

configure.ac

@@ -842,10 +842,10 @@ esac
 #  Gravity multipole order
 AC_ARG_WITH([multipole-order],
    [AS_HELP_STRING([--with-multipole-order=<order>],
-      [order of the multipole and gravitational field expansion @<:@ default: 4@:>@]
+      [order of the multipole and gravitational field expansion @<:@ default: 5@:>@]
    )],
    [with_multipole_order="$withval"],
-   [with_multipole_order="4"]
+   [with_multipole_order="5"]
 )
 AC_DEFINE_UNQUOTED([SELF_GRAVITY_MULTIPOLE_ORDER], [$with_multipole_order], [Multipole order])
 

examples/EAGLE_6/eagle_6.yml

@@ -30,7 +30,7 @@ Statistics:
 Gravity:
   eta:                   0.025    # Constant dimensionless multiplier for time integration.
   theta:                 0.7      # Opening angle (Multipole acceptance criterion)
-  epsilon:               0.0001   # Softening length (in internal units).
+  epsilon:               0.001   # Softening length (in internal units).
   
 # Parameters for the hydrodynamics scheme
 SPH:

src/cell.c

@@ -1290,45 +1290,18 @@ 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 e The #engine (holding current time information).
- *
- * @return 1 If the cell needs drifting, 0 otherwise.
- */
-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 Clear the drift flags on the given cell.
  */
 void cell_clear_drift_flags(struct cell *c, void *data) {
   c->do_drift = 0;
   c->do_sub_drift = 0;
+  c->do_grav_drift = 0;
+  c->do_grav_sub_drift = 0;
 }
 
 /**
- * @brief Activate the drifts on the given cell.
+ * @brief Activate the #part drifts on the given cell.
  */
 void cell_activate_drift_part(struct cell *c, struct scheduler *s) {
 
@@ -1355,9 +1328,36 @@ void cell_activate_drift_part(struct cell *c, struct scheduler *s) {
 }
 
 /**
- * @brief Activate the sorts up a cell hierarchy.
+ * @brief Activate the #gpart drifts on the given cell.
  */
+void cell_activate_drift_gpart(struct cell *c, struct scheduler *s) {
+
+  /* If this cell is already marked for drift, quit early. */
+  if (c->do_grav_drift) return;
 
+  /* Mark this cell for drifting. */
+  c->do_grav_drift = 1;
+
+  /* Set the do_grav_sub_drifts all the way up and activate the super drift
+     if this has not yet been done. */
+  if (c == c->super) {
+    scheduler_activate(s, c->drift_gpart);
+  } else {
+    for (struct cell *parent = c->parent;
+         parent != NULL && !parent->do_grav_sub_drift;
+         parent = parent->parent) {
+      parent->do_grav_sub_drift = 1;
+      if (parent == c->super) {
+        scheduler_activate(s, parent->drift_gpart);
+        break;
+      }
+    }
+  }
+}
+
+/**
+ * @brief Activate the sorts up a cell hierarchy.
+ */
 void cell_activate_sorts_up(struct cell *c, struct scheduler *s) {
   if (c == c->super) {
     scheduler_activate(s, c->sorts);
@@ -1401,7 +1401,13 @@ void cell_activate_sorts(struct cell *c, int sid, struct scheduler *s) {
 }
 
 /**
- * @brief Traverse a sub-cell task and activate the sort tasks along the way.
+ * @brief Traverse a sub-cell task and activate the hydro drift tasks that are
+ * required
+ * by a hydro task
+ *
+ * @param ci The first #cell we recurse in.
+ * @param cj The second #cell we recurse in.
+ * @param s The task #scheduler.
  */
 void cell_activate_subcell_tasks(struct cell *ci, struct cell *cj,
                                  struct scheduler *s) {
@@ -1668,6 +1674,172 @@ void cell_activate_subcell_tasks(struct cell *ci, struct cell *cj,
   }
 }
 
+/**
+ * @brief Traverse a sub-cell task and activate the gravity drift tasks that are
+ * required
+ * by a self gravity task.
+ *
+ * @param ci The first #cell we recurse in.
+ * @param cj The second #cell we recurse in.
+ * @param s The task #scheduler.
+ */
+void cell_activate_subcell_grav_tasks(struct cell *ci, struct cell *cj,
+                                      struct scheduler *s) {
+  /* Some constants */
+  const struct space *sp = s->space;
+  const struct engine *e = sp->e;
+  const int periodic = sp->periodic;
+  const double dim[3] = {sp->dim[0], sp->dim[1], sp->dim[2]};
+  const double theta_crit2 = e->gravity_properties->theta_crit2;
+
+  /* Self interaction? */
+  if (cj == NULL) {
+
+    /* Do anything? */
+    if (!cell_is_active(ci, e)) return;
+
+    /* Recurse? */
+    if (ci->split) {
+
+      /* Loop over all progenies and pairs of progenies */
+      for (int j = 0; j < 8; j++) {
+        if (ci->progeny[j] != NULL) {
+          cell_activate_subcell_grav_tasks(ci->progeny[j], NULL, s);
+          for (int k = j + 1; k < 8; k++)
+            if (ci->progeny[k] != NULL)
+              cell_activate_subcell_grav_tasks(ci->progeny[j], ci->progeny[k],
+                                               s);
+        }
+      }
+    } else {
+
+      /* We have reached the bottom of the tree: activate gpart drift */
+      cell_activate_drift_gpart(ci, s);
+    }
+  }
+
+  /* Pair interaction */
+  else {
+
+    /* Anything to do here? */
+    if (!cell_is_active(ci, e) && !cell_is_active(cj, e)) return;
+
+    /* Recover the multipole information */
+    struct gravity_tensors *const multi_i = ci->multipole;
+    struct gravity_tensors *const multi_j = cj->multipole;
+    const double ri_max = multi_i->r_max;
+    const double rj_max = multi_j->r_max;
+
+    /* Get the distance between the CoMs */
+    double dx = multi_i->CoM[0] - multi_j->CoM[0];
+    double dy = multi_i->CoM[1] - multi_j->CoM[1];
+    double dz = multi_i->CoM[2] - multi_j->CoM[2];
+
+    /* Apply BC */
+    if (periodic) {
+      dx = nearest(dx, dim[0]);
+      dy = nearest(dy, dim[1]);
+      dz = nearest(dz, dim[2]);
+    }
+    const double r2 = dx * dx + dy * dy + dz * dz;
+
+    /* Can we use multipoles ? */
+    if (gravity_M2L_accept(multi_i->r_max, multi_j->r_max, theta_crit2, r2)) {
+
+      /* Ok, no need to drift anything */
+      return;
+    }
+    /* Otherwise, activate the gpart drifts if we are at the bottom. */
+    else if (!ci->split && !cj->split) {
+
+      /* Activate the drifts if the cells are local. */
+      if (cell_is_active(ci, e) || cell_is_active(cj, e)) {
+        if (ci->nodeID == engine_rank) cell_activate_drift_gpart(ci, s);
+        if (cj->nodeID == engine_rank) cell_activate_drift_gpart(cj, s);
+      }
+    }
+    /* Ok, we can still recurse */
+    else {
+
+      if (ri_max > rj_max) {
+        if (ci->split) {
+
+          /* Loop over ci's children */
+          for (int k = 0; k < 8; k++) {
+            if (ci->progeny[k] != NULL)
+              cell_activate_subcell_grav_tasks(ci->progeny[k], cj, s);
+          }
+
+        } else if (cj->split) {
+
+          /* Loop over cj's children */
+          for (int k = 0; k < 8; k++) {
+            if (cj->progeny[k] != NULL)
+              cell_activate_subcell_grav_tasks(ci, cj->progeny[k], s);
+          }
+
+        } else {
+          error("Fundamental error in the logic");
+        }
+      } else if (rj_max >= ri_max) {
+        if (cj->split) {
+
+          /* Loop over cj's children */
+          for (int k = 0; k < 8; k++) {
+            if (cj->progeny[k] != NULL)
+              cell_activate_subcell_grav_tasks(ci, cj->progeny[k], s);
+          }
+
+        } else if (ci->split) {
+
+          /* Loop over ci's children */
+          for (int k = 0; k < 8; k++) {
+            if (ci->progeny[k] != NULL)
+              cell_activate_subcell_grav_tasks(ci->progeny[k], cj, s);
+          }
+
+        } else {
+          error("Fundamental error in the logic");
+        }
+      }
+    }
+  }
+}
+
+/**
+ * @brief Traverse a sub-cell task and activate the gravity drift tasks that are
+ * required
+ * by an external gravity task.
+ *
+ * @param ci The #cell we recurse in.
+ * @param s The task #scheduler.
+ */
+void cell_activate_subcell_external_grav_tasks(struct cell *ci,
+                                               struct scheduler *s) {
+
+  /* Some constants */
+  const struct space *sp = s->space;
+  const struct engine *e = sp->e;
+
+  /* Do anything? */
+  if (!cell_is_active(ci, e)) return;
+
+  /* Recurse? */
+  if (ci->split) {
+
+    /* Loop over all progenies (no need for pairs for self-gravity) */
+    for (int j = 0; j < 8; j++) {
+      if (ci->progeny[j] != NULL) {
+        cell_activate_subcell_external_grav_tasks(ci->progeny[j], s);
+      }
+    }
+  } else {
+
+    /* We have reached the bottom of the tree: activate gpart drift */
+    cell_activate_drift_gpart(ci, s);
+  }
+}
+
 /**
  * @brief Un-skips all the tasks associated with a given cell and checks
  * if the space needs to be rebuilt.
@@ -1693,8 +1865,13 @@ int cell_unskip_tasks(struct cell *c, struct scheduler *s) {
         (cj != NULL && cell_is_active(cj, e) && cj->nodeID == engine_rank)) {
       scheduler_activate(s, t);
 
-      /* Set the correct sorting flags */
-      if (t->type == task_type_pair) {
+      /* Activate hydro drift */
+      if (t->type == task_type_self) {
+        if (ci->nodeID == engine_rank) cell_activate_drift_part(ci, s);
+      }
+
+      /* Set the correct sorting flags and activate hydro drifts */
+      else if (t->type == task_type_pair) {
         /* Store some values. */
         atomic_or(&ci->requires_sorts, 1 << t->flags);
         atomic_or(&cj->requires_sorts, 1 << t->flags);
@@ -1843,6 +2020,29 @@ int cell_unskip_tasks(struct cell *c, struct scheduler *s) {
     }
   }
 
+  /* Un-skip the gravity tasks involved with this cell. */
+  for (struct link *l = c->grav; l != NULL; l = l->next) {
+    struct task *t = l->t;
+    struct cell *ci = t->ci;
+    struct cell *cj = t->cj;
+
+    /* Only activate tasks that involve a local active cell. */
+    if ((cell_is_active(ci, e) && ci->nodeID == engine_rank) ||
+        (cj != NULL && cell_is_active(cj, e) && cj->nodeID == engine_rank)) {
+      scheduler_activate(s, t);
+
+      /* Set the drifting flags */
+      if (t->type == task_type_self &&
+          t->subtype == task_subtype_external_grav) {
+        cell_activate_subcell_external_grav_tasks(t->ci, s);
+      } else if (t->type == task_type_self && t->subtype == task_subtype_grav) {
+        cell_activate_subcell_grav_tasks(t->ci, NULL, s);
+      } else if (t->type == task_type_pair) {
+        cell_activate_subcell_grav_tasks(t->ci, t->cj, s);
+      }
+    }
+  }
+
   /* Unskip all the other task types. */
   if (c->nodeID == engine_rank && cell_is_active(c, e)) {
 
@@ -1850,15 +2050,12 @@ int cell_unskip_tasks(struct cell *c, struct scheduler *s) {
       scheduler_activate(s, l->t);
     for (struct link *l = c->force; l != NULL; l = l->next)
       scheduler_activate(s, l->t);
-    for (struct link *l = c->grav; l != NULL; l = l->next)
-      scheduler_activate(s, l->t);
 
     if (c->extra_ghost != NULL) scheduler_activate(s, c->extra_ghost);
     if (c->ghost_in != NULL) scheduler_activate(s, c->ghost_in);
     if (c->ghost_out != NULL) scheduler_activate(s, c->ghost_out);
     if (c->ghost != NULL) scheduler_activate(s, c->ghost);
     if (c->init_grav != NULL) scheduler_activate(s, c->init_grav);
-    if (c->drift_gpart != NULL) scheduler_activate(s, c->drift_gpart);
     if (c->kick1 != NULL) scheduler_activate(s, c->kick1);
     if (c->kick2 != NULL) scheduler_activate(s, c->kick2);
     if (c->timestep != NULL) scheduler_activate(s, c->timestep);
@@ -1931,7 +2128,7 @@ void cell_drift_part(struct cell *c, const struct engine *e, int force) {
 
   /* Check that we are actually going to move forward. */
   if (ti_current < ti_old_part) error("Attempt to drift to the past");
-#endif  // SWIFT_DEBUG_CHECKS
+#endif
 
   /* Are we not in a leaf ? */
   if (c->split && (force || c->do_sub_drift)) {
@@ -2016,8 +2213,9 @@ void cell_drift_part(struct cell *c, const struct engine *e, int force) {
  *
  * @param c The #cell.
  * @param e The #engine (to get ti_current).
+ * @param force Drift the particles irrespective of the #cell flags.
  */
-void cell_drift_gpart(struct cell *c, const struct engine *e) {
+void cell_drift_gpart(struct cell *c, const struct engine *e, int force) {
 
   const double timeBase = e->timeBase;
   const integertime_t ti_old_gpart = c->ti_old_gpart;
@@ -2029,11 +2227,19 @@ void cell_drift_gpart(struct cell *c, const struct engine *e) {
   const double dt = (ti_current - ti_old_gpart) * timeBase;
   float dx_max = 0.f, dx2_max = 0.f;
 
+  /* Drift irrespective of cell flags? */
+  force |= c->do_grav_drift;
+
+#ifdef SWIFT_DEBUG_CHECKS
+  /* Check that we only drift local cells. */
+  if (c->nodeID != engine_rank) error("Drifting a foreign cell is nope.");
+
   /* Check that we are actually going to move forward. */
   if (ti_current < ti_old_gpart) error("Attempt to drift to the past");
+#endif
 
   /* Are we not in a leaf ? */
-  if (c->split) {
+  if (c->split && (force || c->do_grav_sub_drift)) {
 
     /* Loop over the progeny and collect their data. */
     for (int k = 0; k < 8; k++)
@@ -2041,13 +2247,19 @@ void cell_drift_gpart(struct cell *c, const struct engine *e) {
         struct cell *cp = c->progeny[k];
 
         /* Recurse */
-        cell_drift_gpart(cp, e);
+        cell_drift_gpart(cp, e, force);
 
         /* Update */
         dx_max = max(dx_max, cp->dx_max_gpart);
       }
 
-  } else if (ti_current > ti_old_gpart) {
+    /* Store the values */
+    c->dx_max_gpart = dx_max;
+
+    /* Update the time of the last drift */
+    c->ti_old_gpart = ti_current;
+
+  } else if (!c->split && force && ti_current > ti_old_gpart) {
 
     /* Loop over all the g-particles in the cell */
     const size_t nr_gparts = c->gcount;
@@ -2087,16 +2299,16 @@ void cell_drift_gpart(struct cell *c, const struct engine *e) {
     /* Now, get the maximal particle motion from its square */
     dx_max = sqrtf(dx2_max);
 
-  } else {
+    /* Store the values */
+    c->dx_max_gpart = dx_max;
 
-    dx_max = c->dx_max_gpart;
+    /* Update the time of the last drift */
+    c->ti_old_gpart = ti_current;
   }
 
-  /* Store the values */
-  c->dx_max_gpart = dx_max;
-
-  /* Update the time of the last drift */
-  c->ti_old_gpart = ti_current;
+  /* Clear the drift flags. */
+  c->do_grav_drift = 0;
+  c->do_grav_sub_drift = 0;
 }
 
 /**
@@ -2118,7 +2330,8 @@ void cell_drift_all_multipoles(struct cell *c, const struct engine *e) {
   if (ti_current < ti_old_multipole) error("Attempt to drift to the past");
 
   /* Drift the multipole */
-  if (ti_current > ti_old_multipole) gravity_drift(c->multipole, dt);
+  if (ti_current > ti_old_multipole)
+    gravity_drift(c->multipole, dt, c->dx_max_gpart);
 
   /* Are we not in a leaf ? */
   if (c->split) {
@@ -2153,7 +2366,8 @@ void cell_drift_multipole(struct cell *c, const struct engine *e) {
   /* Check that we are actually going to move forward. */
   if (ti_current < ti_old_multipole) error("Attempt to drift to the past");
 
-  if (ti_current > ti_old_multipole) gravity_drift(c->multipole, dt);
+  if (ti_current > ti_old_multipole)
+    gravity_drift(c->multipole, dt, c->dx_max_gpart);
 
   /* Update the time of the last drift */
   c->ti_old_multipole = ti_current;

src/cell.h

@@ -152,9 +152,13 @@ struct cell {
   /*! The multipole initialistation task */
   struct task *init_grav;
 
-  /*! The ghost tasks */
+  /*! Dependency implicit task for the ghost  (in->ghost->out)*/
   struct task *ghost_in;
+
+  /*! Dependency implicit task for the ghost  (in->ghost->out)*/
   struct task *ghost_out;
+
+  /*! The ghost task itself */
   struct task *ghost;
 
   /*! The extra ghost task for complex hydro schemes */
@@ -311,6 +315,21 @@ struct cell {
   /*! Is the #spart data of this cell being used in a sub-cell? */
   int shold;
 
+  /*! Values of dx_max before the drifts, used for sub-cell tasks. */
+  float dx_max_old;
+
+  /*! Values of h_max before the drifts, used for sub-cell tasks. */
+  float h_max_old;
+
+  /*! Values of dx_max_sort before the drifts, used for sub-cell tasks. */
+  float dx_max_sort_old;
+
+  /*! Bit mask of sort directions that will be needed in the next timestep. */
+  unsigned int requires_sorts;
+
+  /*! Bit mask of sorts that need to be computed for this cell. */
+  unsigned int do_sort;
+
   /*! Number of tasks that are associated with this cell. */
   short int nr_tasks;
 
@@ -323,22 +342,17 @@ struct cell {
   /*! The maximal depth of this cell and its progenies */
   char maxdepth;
 
-  /*! Values of dx_max and h_max before the drifts, used for sub-cell tasks. */
-  float dx_max_old;
-  float h_max_old;
-  float dx_max_sort_old;
-
-  /* Bit mask of sort directions that will be needed in the next timestep. */
-  unsigned int requires_sorts;
-
-  /*! Does this cell need to be drifted? */
+  /*! Does this cell need to be drifted (hydro)? */
   char do_drift;