Allow M2L on cells with size 1. Symmetries correctly the long-range task and...

Allow M2L on cells with size 1. Symmetries correctly the long-range task and the task creation loop.

src/cell.c

@@ -6396,10 +6396,6 @@ int cell_can_use_pair_mm(const struct cell *restrict ci,
   const int periodic = s->periodic;
   const double dim[3] = {s->dim[0], s->dim[1], s->dim[2]};
 
-  /* Check for trivial cases */
-  if (ci->grav.count <= 1) return 0;
-  if (cj->grav.count <= 1) return 0;
-
   /* Recover the multipole information */
   const struct gravity_tensors *restrict multi_i = ci->grav.multipole;
   const struct gravity_tensors *restrict multi_j = cj->grav.multipole;

src/multipole_accept.h

@@ -94,11 +94,6 @@ __attribute__((nonnull, pure)) INLINE static int gravity_M2L_accept(
   const float max_softening =
       max(A->m_pole.max_softening, B->m_pole.max_softening);
 
-#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) {
@@ -106,10 +101,12 @@ __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;
-  E_BA_term *= max(rho_A, rho_B);
-  E_BA_term /= (rho_A + rho_B);
+  if (rho_A + rho_B > 0.f) {
+    E_BA_term *= max(rho_A, rho_B);
+    E_BA_term /= (rho_A + rho_B);
+  }
 
-  /* Compute r^p */
+    /* Compute r^p */
 #if SELF_GRAVITY_MULTIPOLE_ORDER % 2 == 1
   const float r_to_p = integer_powf(sqrtf(r2), p);
 #else

src/runner_doiact_grav.c

@@ -2357,12 +2357,6 @@ void runner_do_grav_long_range(struct runner *r, struct cell *ci,
   struct cell *top = ci;
   while (top->parent != NULL) top = top->parent;
 
-  /* Recover the top-level multipole (for distance checks) */
-  struct gravity_tensors *const multi_top = top->grav.multipole;
-  const double CoM_rebuild_top[3] = {multi_top->CoM_rebuild[0],
-                                     multi_top->CoM_rebuild[1],
-                                     multi_top->CoM_rebuild[2]};
-
   /* Loop over all the top-level cells and go for a M-M interaction if
    * well-separated */
   for (int n = 0; n < nr_cells_with_particles; ++n) {
@@ -2407,23 +2401,7 @@ void runner_do_grav_long_range(struct runner *r, struct cell *ci,
       }
     }
 
-    /* Get the distance between the CoMs at the last rebuild*/
-    double dx_r = CoM_rebuild_top[0] - multi_j->CoM_rebuild[0];
-    double dy_r = CoM_rebuild_top[1] - multi_j->CoM_rebuild[1];
-    double dz_r = CoM_rebuild_top[2] - multi_j->CoM_rebuild[2];
-
-    /* Apply BC */
-    if (periodic) {
-      dx_r = nearest(dx_r, dim[0]);
-      dy_r = nearest(dy_r, dim[1]);
-      dz_r = nearest(dz_r, dim[2]);
-    }
-    const double r2_rebuild = dx_r * dx_r + dy_r * dy_r + dz_r * dz_r;
-
-    /* Are we in charge of this cell pair? */
-    if (gravity_M2L_accept_symmetric(e->gravity_properties, multi_top, multi_j,
-                                     r2_rebuild, /*use_rebuild_sizes=*/1,
-                                     periodic)) {
+    if (cell_can_use_pair_mm(top, cj, e, e->s, /*use_rebuild_data=*/1)) {
 
       /* Call the PM interaction fucntion on the active sub-cells of ci */
       runner_dopair_grav_mm_nonsym(r, ci, cj);