Use the multipole acceptance criterion to build the long-range interactions...

Use the multipole acceptance criterion to build the long-range interactions and not just the geometry of the mesh.

src/engine.c

@@ -1647,9 +1647,9 @@ void engine_exchange_strays(struct engine *e, size_t offset_parts,
  * @brief Constructs the top-level tasks for the short-range gravity
  * interactions.
  *
- * All top-cells get a self task.
- * All neighbouring pairs get a pair task.
- * All non-neighbouring pairs within a range of 6 cells get a M-M task.
+ * - All top-cells get a self task.
+ * - All pairs within range according to the multipole acceptance
+ *   criterion get a pair task.
  *
  * @param e The #engine.
  */
@@ -1658,6 +1658,7 @@ void engine_make_self_gravity_tasks(struct engine *e) {
   struct space *s = e->s;
   struct scheduler *sched = &e->sched;
   const int nodeID = e->nodeID;
+  const double theta_crit_inv = e->gravity_properties->theta_crit_inv;
   struct cell *cells = s->cells_top;
   const int nr_cells = s->nr_cells;
 
@@ -1668,13 +1669,14 @@ void engine_make_self_gravity_tasks(struct engine *e) {
     /* Skip cells without gravity particles */
     if (ci->gcount == 0) continue;
 
-    /* Is that neighbour local ? */
+    /* Is that cell local ? */
     if (ci->nodeID != nodeID) continue;
 
     /* If the cells is local build a self-interaction */
     scheduler_addtask(sched, task_type_self, task_subtype_grav, 0, 0, ci, NULL,
                       0);
 
+    /* Loop over every other cell */
     for (int cjd = cid + 1; cjd < nr_cells; ++cjd) {
 
       struct cell *cj = &cells[cjd];
@@ -1683,9 +1685,11 @@ void engine_make_self_gravity_tasks(struct engine *e) {
       if (cj->gcount == 0) continue;
 
       /* Is that neighbour local ? */
-      if (cj->nodeID != nodeID) continue;
+      if (cj->nodeID != nodeID) continue;  // MATTHIEU
 
-      if (cell_are_neighbours(ci, cj))
+      /* Are the cells to close for a MM interaction ? */
+      if (!gravity_multipole_accept(ci->multipole, cj->multipole,
+                                    theta_crit_inv))
         scheduler_addtask(sched, task_type_pair, task_subtype_grav, 0, 0, ci,
                           cj, 1);
     }

src/gravity_properties.c

@@ -48,8 +48,8 @@ void gravity_props_init(struct gravity_props *p,
   p->eta = parser_get_param_float(params, "Gravity:eta");
 
   /* Opening angle */
-  p->theta = parser_get_param_double(params, "Gravity:theta");
-  p->theta_inv = 1. / p->theta;
+  p->theta_crit = parser_get_param_double(params, "Gravity:theta");
+  p->theta_crit_inv = 1. / p->theta_crit;
 
   /* Softening lengths */
   p->epsilon = parser_get_param_double(params, "Gravity:epsilon");
@@ -63,7 +63,7 @@ void gravity_props_print(const struct gravity_props *p) {
 
   message("Self-gravity time integration: eta=%.4f", p->eta);
 
-  message("Self-gravity opening angle:  theta=%.4f", p->theta);
+  message("Self-gravity opening angle:  theta=%.4f", p->theta_crit);
 
   message("Self-gravity softening:    epsilon=%.4f", p->epsilon);
 
@@ -80,7 +80,7 @@ void gravity_props_print_snapshot(hid_t h_grpgrav,
 
   io_write_attribute_f(h_grpgrav, "Time integration eta", p->eta);
   io_write_attribute_f(h_grpgrav, "Softening length", p->epsilon);
-  io_write_attribute_f(h_grpgrav, "Opening angle", p->theta);
+  io_write_attribute_f(h_grpgrav, "Opening angle", p->theta_crit);
   io_write_attribute_f(h_grpgrav, "MM a_smooth", p->a_smooth);
   io_write_attribute_f(h_grpgrav, "MM r_cut", p->r_cut);
 }

src/gravity_properties.h

@@ -42,10 +42,10 @@ struct gravity_props {
   float eta;
 
   /*! Tree opening angle (Multipole acceptance criterion) */
-  double theta;
+  double theta_crit;
 
   /*! Inverse of opening angle */
-  double theta_inv;
+  double theta_crit_inv;
 
   /*! Softening length */
   double epsilon;

src/multipole.h

@@ -208,9 +208,9 @@ INLINE static void gravity_reset(struct gravity_tensors *m) {
 INLINE static void gravity_drift(struct gravity_tensors *m, double dt) {
 
   /* Move the whole thing according to bulk motion */
-  m->CoM[0] += m->m_pole.vel[0] * dt;
-  m->CoM[1] += m->m_pole.vel[1] * dt;
-  m->CoM[2] += m->m_pole.vel[2] * dt;
+  m->CoM[0] += m->m_pole.vel[0] * dt * 0;  // MATTHIEU
+  m->CoM[1] += m->m_pole.vel[1] * dt * 0;
+  m->CoM[2] += m->m_pole.vel[2] * dt * 0;
 }
 
 /**
@@ -2508,4 +2508,29 @@ INLINE static void gravity_L2P(const struct grav_tensor *lb,
   gp->a_grav[2] += a_grav[2];
 }
 
+/**
+ * @brief Checks whether a cell-cell interaction can be appromixated by a M-M
+ * interaction.
+ *
+ * @param ma The #multipole of the first #cell.
+ * @param mb The #multipole of the second #cell.
+ * @param theta_crit_inv The inverse of the critical opening angle.
+ */
+__attribute__((always_inline)) INLINE static int gravity_multipole_accept(
+    const struct gravity_tensors *ma, const struct gravity_tensors *mb,
+    double theta_crit_inv) {
+
+  const double r_crit_a = ma->r_max * theta_crit_inv;
+  const double r_crit_b = mb->r_max * theta_crit_inv;
+
+  const double dx = ma->CoM[0] - mb->CoM[0];
+  const double dy = ma->CoM[1] - mb->CoM[1];
+  const double dz = ma->CoM[2] - mb->CoM[2];
+
+  const double r2 = dx * dx + dy * dy + dz * dz;
+
+  /* Multipole acceptance criterion (Dehnen 2002, eq.10) */
+  return (r2 > (r_crit_a + r_crit_b) * (r_crit_a + r_crit_b));
+}
+
 #endif /* SWIFT_MULTIPOLE_H */

src/runner_doiact_grav.h

@@ -23,6 +23,7 @@
 /* Includes. */
 #include "cell.h"
 #include "gravity.h"
+#include "inline.h"
 #include "part.h"
 
 /**
@@ -379,14 +380,15 @@ void runner_dopair_grav(struct runner *r, struct cell *ci, struct cell *cj,
         "The impossible has happened: pair interaction between a cell and "
         "itself.");
 
-  /* Are the cells direct neighbours? */
-  if (!cell_are_neighbours(ci, cj))
-    error(
-        "Non-neighbouring cells ! ci->x=[%f %f %f] ci->width=%f cj->loc=[%f %f "
-        "%f] "
-        "cj->width=%f",
-        ci->loc[0], ci->loc[1], ci->loc[2], ci->width[0], cj->loc[0],
-        cj->loc[1], cj->loc[2], cj->width[0]);
+/* Are the cells direct neighbours? */
+/* if (!cell_are_neighbours(ci, cj)) */
+/*   error( */
+/*       "Non-neighbouring cells ! ci->x=[%f %f %f] ci->width=%f cj->loc=[%f %f
+ * " */
+/*       "%f] " */
+/*       "cj->width=%f", */
+/*       ci->loc[0], ci->loc[1], ci->loc[2], ci->width[0], cj->loc[0], */
+/*       cj->loc[1], cj->loc[2], cj->width[0]); */
 
 #endif
 
@@ -515,6 +517,14 @@ void runner_dosub_grav(struct runner *r, struct cell *ci, struct cell *cj,
   }
 }
 
+/**
+ * @brief Performs all M-M interactions between a given top-level cell and all
+ * the other top-levels that are far enough.
+ *
+ * @param r The thread #runner.
+ * @param ci The #cell of interest.
+ * @param timer Are we timing this ?
+ */
 void runner_do_grav_long_range(struct runner *r, struct cell *ci, int timer) {
 
 #if ICHECK > 0
@@ -530,39 +540,40 @@ void runner_do_grav_long_range(struct runner *r, struct cell *ci, int timer) {
   }
 #endif
 
+  /* Some constants */
+  const struct engine *e = r->e;
+  const struct gravity_props *props = e->gravity_properties;
+  const double theta_crit_inv = props->theta_crit_inv;
+
   TIMER_TIC;
 
   /* Recover the list of top-level cells */
-  const struct engine *e = r->e;
   struct cell *cells = e->s->cells_top;
   const int nr_cells = e->s->nr_cells;
-  /* const double max_d = */
-  /*     const_gravity_a_smooth * const_gravity_r_cut * ci->width[0]; */
-  /* const double max_d2 = max_d * max_d; */
-  // const double pos_i[3] = {ci->loc[0], ci->loc[1], ci->loc[2]};
 
   /* Anything to do here? */
-  if (!cell_is_active(ci, e)) return;
+  if (!cell_is_active(ci, e)) return;  // MATTHIEU (should never happen)
 
-  /* Drift our own multipole if need be */
-  if (ci->ti_old_multipole != e->ti_current) cell_drift_multipole(ci, e);
+  /* Check multipole has been drifted */
+  if (ci->ti_old_multipole != e->ti_current)
+    error("Interacting un-drifted multipole");
 
-  /* Loop over all the cells and go for a p-m interaction if far enough but not
-   * too far */
+  /* Loop over all the top-level cells and go for a M-M interaction if
+   * well-separated */
   for (int i = 0; i < nr_cells; ++i) {
 
+    /* Handle on the top-level cell */
     struct cell *cj = &cells[i];
 
-    if (ci == cj) continue;
-    if (cj->gcount == 0) continue;
-
-    /* const double dx[3] = {cj->loc[0] - pos_i[0],   // x */
-    /*                       cj->loc[1] - pos_i[1],   // y */
-    /*                       cj->loc[2] - pos_i[2]};  // z */
-    /* const double r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2]; */
-    /* if (r2 > max_d2) continue; */
+    /* Avoid stupid cases */
+    if (ci == cj || cj->gcount == 0) continue;
 
-    if (!cell_are_neighbours(ci, cj)) runner_dopair_grav_mm(r, ci, cj);
+    /* Check the multipole acceptance criterion */
+    if (gravity_multipole_accept(ci->multipole, cj->multipole,
+                                 theta_crit_inv)) {
+      /* Go for a (non-symmetric) M-M calculation */
+      runner_dopair_grav_mm(r, ci, cj);
+    }
   }
 
   if (timer) TIMER_TOC(timer_dograv_long_range);