Recursively compute an upper-limit to the radius containing all particles in a given cell

src/cell.c

@@ -1165,6 +1165,16 @@ void cell_check_multipole(struct cell *c, void *data) {
       gravity_multipole_print(&c->multipole->m_pole);
       error("Aborting");
     }
+
+    /* Check that the upper limit of r_max is good enough */
+    if (!(c->multipole->r_max >= ma.r_max)) {
+      error("Upper-limit r_max=%e too small. Should be >=%e.",
+            c->multipole->r_max, ma.r_max);
+    } else if (c->multipole->r_max * c->multipole->r_max >
+               3. * c->width[0] * c->width[0]) {
+      error("r_max=%e larger than cell diagonal %e.", c->multipole->r_max,
+            sqrt(3. * c->width[0] * c->width[0]));
+    }
   }
 #else
   error("Calling debugging code without debugging flag activated.");

src/multipole.h

@@ -176,6 +176,9 @@ struct gravity_tensors {
       /*! Centre of mass of the matter dsitribution */
       double CoM[3];
 
+      /*! Upper limit of the CoM<->gpart distance */
+      double r_max;
+
       /*! Multipole mass */
       struct multipole m_pole;
 
@@ -967,9 +970,9 @@ INLINE static void gravity_P2M(struct gravity_tensors *m,
   double com[3] = {0.0, 0.0, 0.0};
   double vel[3] = {0.f, 0.f, 0.f};
 
-  /* Collect the particle data. */
+  /* Collect the particle data for CoM. */
   for (int k = 0; k < gcount; k++) {
-    const float m = gparts[k].mass;
+    const double m = gparts[k].mass;
 
     mass += m;
     com[0] += gparts[k].x[0] * m;
@@ -989,7 +992,8 @@ INLINE static void gravity_P2M(struct gravity_tensors *m,
   vel[1] *= imass;
   vel[2] *= imass;
 
-/* Prepare some local counters */
+  /* Prepare some local counters */
+  double r_max2 = 0.;
 #if SELF_GRAVITY_MULTIPOLE_ORDER > 0
   double M_100 = 0., M_010 = 0., M_001 = 0.;
 #endif
@@ -1026,11 +1030,15 @@ INLINE static void gravity_P2M(struct gravity_tensors *m,
   /* Construce the higher order terms */
   for (int k = 0; k < gcount; k++) {
 
-#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
-    const double m = gparts[k].mass;
     const double dx[3] = {gparts[k].x[0] - com[0], gparts[k].x[1] - com[1],
                           gparts[k].x[2] - com[2]};
 
+    /* Maximal distance CoM<->gpart */
+    r_max2 = max(r_max2, dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2]);
+
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
+    const double m = gparts[k].mass;
+
     /* 1st order terms */
     M_100 += -m * X_100(dx);
     M_010 += -m * X_010(dx);
@@ -1115,6 +1123,7 @@ INLINE static void gravity_P2M(struct gravity_tensors *m,
 
   /* Store the data on the multipole. */
   m->m_pole.M_000 = mass;
+  m->r_max = sqrt(r_max2);
   m->CoM[0] = com[0];
   m->CoM[1] = com[1];
   m->CoM[2] = com[2];

src/space.c

@@ -2083,15 +2083,38 @@ void space_split_recursive(struct space *s, struct cell *c,
 
       /* Now shift progeny multipoles and add them up */
       struct multipole temp;
+      double r_max = 0.;
       for (int k = 0; k < 8; ++k) {
         if (c->progeny[k] != NULL) {
           const struct cell *cp = c->progeny[k];
           const struct multipole *m = &cp->multipole->m_pole;
+
+          /* Contribution to multipole */
           gravity_M2M(&temp, m, c->multipole->CoM, cp->multipole->CoM,
                       s->periodic);
           gravity_multipole_add(&c->multipole->m_pole, &temp);
+
+          /* Upper limit of max CoM<->gpart distance */
+          const float dx = c->multipole->CoM[0] - cp->multipole->CoM[0];
+          const float dy = c->multipole->CoM[1] - cp->multipole->CoM[1];
+          const float dz = c->multipole->CoM[2] - cp->multipole->CoM[2];
+          const float r2 = dx * dx + dy * dy + dz * dz;
+          r_max = max(r_max, cp->multipole->r_max + sqrtf(r2));
         }
       }
+      /* Alternative upper limit of max CoM<->gpart distance */
+      const double dx = c->multipole->CoM[0] > c->loc[0] + c->width[0] / 2.
+                            ? c->multipole->CoM[0] - c->loc[0]
+                            : c->loc[0] + c->width[0] - c->multipole->CoM[0];
+      const double dy = c->multipole->CoM[1] > c->loc[1] + c->width[1] / 2.
+                            ? c->multipole->CoM[1] - c->loc[1]
+                            : c->loc[1] + c->width[1] - c->multipole->CoM[1];
+      const double dz = c->multipole->CoM[2] > c->loc[2] + c->width[2] / 2.
+                            ? c->multipole->CoM[2] - c->loc[2]
+                            : c->loc[2] + c->width[2] - c->multipole->CoM[2];
+
+      /* Take minimum of both limits */
+      c->multipole->r_max = min(r_max, sqrt(dx * dx + dy * dy + dz * dz));
     }
   }
 
@@ -2151,13 +2174,24 @@ void space_split_recursive(struct space *s, struct cell *c,
 
     /* Construct the multipole and the centre of mass*/
     if (s->gravity) {
-      if (gcount > 0)
+      if (gcount > 0) {
         gravity_P2M(c->multipole, c->gparts, c->gcount);
-      else {
+        const double dx = c->multipole->CoM[0] > c->loc[0] + c->width[0] / 2.
+                              ? c->multipole->CoM[0] - c->loc[0]
+                              : c->loc[0] + c->width[0] - c->multipole->CoM[0];
+        const double dy = c->multipole->CoM[1] > c->loc[1] + c->width[1] / 2.
+                              ? c->multipole->CoM[1] - c->loc[1]
+                              : c->loc[1] + c->width[1] - c->multipole->CoM[1];
+        const double dz = c->multipole->CoM[2] > c->loc[2] + c->width[2] / 2.
+                              ? c->multipole->CoM[2] - c->loc[2]
+                              : c->loc[2] + c->width[2] - c->multipole->CoM[2];
+        c->multipole->r_max = sqrt(dx * dx + dy * dy + dz * dz);
+      } else {
         gravity_multipole_init(&c->multipole->m_pole);
         c->multipole->CoM[0] = c->loc[0] + c->width[0] / 2.;
         c->multipole->CoM[1] = c->loc[1] + c->width[1] / 2.;
         c->multipole->CoM[2] = c->loc[2] + c->width[2] / 2.;
+        c->multipole->r_max = 0.;
       }
     }
   }