Check for multipoles of size 0 in M2L

src/multipole_accept.h

@@ -56,6 +56,11 @@ __attribute__((nonnull, pure)) INLINE static int gravity_M2L_accept(
   const float rho_A = use_rebuild_sizes ? A->r_max_rebuild : A->r_max;
   const float rho_B = use_rebuild_sizes ? B->r_max_rebuild : B->r_max;
 
+#ifdef SWIFT_DEBUG_CHECKS
+  if (rho_A == 0.) error("Size of multipole A is 0!");
+  if (rho_B == 0.) error("Size of multipole B is 0!");
+#endif
+
   /* Compute the error estimator (without the 1/M_B term that cancels out) */
   float E_BA_term = 0.f;
   for (int n = 0; n <= p; ++n) {
@@ -63,12 +68,10 @@ __attribute__((nonnull, pure)) INLINE static int gravity_M2L_accept(
         binomial(p, n) * B->m_pole.power[n] * integer_powf(rho_A, p - n);
   }
   E_BA_term *= 8.f;
-  if (rho_A != 0. && rho_B != 0.) {
-    E_BA_term *= max(rho_A, rho_B);
-    E_BA_term /= (rho_A + rho_B);
-  }
+  E_BA_term *= max(rho_A, rho_B);
+  E_BA_term /= (rho_A + rho_B);
 
-    /* Compute r^(p+2) */
+  /* Compute r^(p+2) */
 #if SELF_GRAVITY_MULTIPOLE_ORDER % 2 == 1
   const float r_to_p_plus2 = integer_powf(sqrtf(r2), (p + 2));
 #else
@@ -164,6 +167,13 @@ __attribute__((nonnull, pure)) INLINE static int gravity_M2P_accept(
   /* Order of the expansion */
   const int p = SELF_GRAVITY_MULTIPOLE_ORDER;
 
+  /* Sizes of the multipoles */
+  const float rho_B = B->r_max;
+
+#ifdef SWIFT_DEBUG_CHECKS
+  if (rho_B == 0.) error("Size of multipole B is 0!");
+#endif
+
   /* Compute the error estimator (without the 1/M_B term that cancels out) */
   const float E_BA_term = 8.f * B->m_pole.power[p];
 
@@ -196,7 +206,7 @@ __attribute__((nonnull, pure)) INLINE static int gravity_M2P_accept(
     /* Test the different conditions */
 
     /* Condition 1: We are in the converging part of the Taylor expansion */
-    const int cond_1 = (B->r_max) * (B->r_max) < r2;
+    const int cond_1 = rho_B * rho_B < r2;
 
     /* Condition 2: We are not below softening */
     const int cond_2 = max_softening * max_softening < r2;
@@ -210,7 +220,7 @@ __attribute__((nonnull, pure)) INLINE static int gravity_M2P_accept(
   } else {
 
     /* Condition 1: We are obeying the purely geometric criterion */
-    const int cond_1 = (B->r_max) * (B->r_max) < theta_crit2 * r2;
+    const int cond_1 = rho_B * rho_B < theta_crit2 * r2;
 
     /* Condition 2: We are not below softening */
     const int cond_2 = max_softening * max_softening < r2;