diff --git a/src/cell.c b/src/cell.c
index dfae80d5a2cf4a46fe8bc42052220be559efe31b..0b949f3a27edf6e35e4c7679a9daf06f3829397b 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -5822,7 +5822,11 @@ int cell_can_use_pair_mm(const struct cell *ci, const struct cell *cj,
}
const double r2 = dx * dx + dy * dy + dz * dz;
- return gravity_M2L_accept(multi_i->r_max, multi_j->r_max, theta_crit2, r2);
+ const double epsilon_i = multi_i->m_pole.max_softening;
+ const double epsilon_j = multi_j->m_pole.max_softening;
+
+ return gravity_M2L_accept(multi_i->r_max, multi_j->r_max, theta_crit2, r2,
+ epsilon_i, epsilon_j);
}
/**
@@ -5884,6 +5888,9 @@ int cell_can_use_pair_mm_rebuild(const struct cell *ci, const struct cell *cj,
}
const double r2 = dx * dx + dy * dy + dz * dz;
+ const double epsilon_i = multi_i->m_pole.max_softening;
+ const double epsilon_j = multi_j->m_pole.max_softening;
+
return gravity_M2L_accept(multi_i->r_max_rebuild, multi_j->r_max_rebuild,
- theta_crit2, r2);
+ theta_crit2, r2, epsilon_i, epsilon_j);
}
diff --git a/src/engine.c b/src/engine.c
index 4afbeb9b2f891d339d1c4d96d072e0b2c5a525ec..e9c250939c20908392086cb811351343aab1c26c 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -4405,19 +4405,26 @@ void engine_makeproxies(struct engine *e) {
sqrt(min_dist_centres2) - 2. * delta_CoM;
const double min_dist_CoM2 = min_dist_CoM * min_dist_CoM;
+ /* We also assume that the softening is negligible compared
+ to the cell size */
+ const double epsilon_i = 0.;
+ const double epsilon_j = 0.;
+
/* Are we beyond the distance where the truncated forces are 0
* but not too far such that M2L can be used? */
if (periodic) {
if ((min_dist_CoM2 < max_mesh_dist2) &&
(!gravity_M2L_accept(r_max, r_max, theta_crit2,
- min_dist_CoM2)))
+ min_dist_CoM2, epsilon_i,
+ epsilon_j)))
proxy_type |= (int)proxy_cell_type_gravity;
} else {
if (!gravity_M2L_accept(r_max, r_max, theta_crit2,
- min_dist_CoM2))
+ min_dist_CoM2, epsilon_i,
+ epsilon_j))
proxy_type |= (int)proxy_cell_type_gravity;
}
}
diff --git a/src/multipole.h b/src/multipole.h
index 0c699ea4d782aeb744b78b472a52c358480184b7..956ec3e47a5c9b031c6cc396581ae21613d6f577 100644
--- a/src/multipole.h
+++ b/src/multipole.h
@@ -2575,23 +2575,31 @@ INLINE static void gravity_L2P(const struct grav_tensor *lb,
* We use the multipole acceptance criterion of Dehnen, 2002, JCoPh, Volume 179,
* Issue 1, pp.27-42, equation 10.
*
+ * We also additionally check that the distance between the multipoles
+ * is larger than the softening lengths (here the distance at which
+ * the gravity becomes Newtonian again, not the Plummer-equivalent quantity).
+ *
* @param r_crit_a The size of the multipole A.
* @param r_crit_b The size of the multipole B.
* @param theta_crit2 The square of the critical opening angle.
* @param r2 Square of the distance (periodically wrapped) between the
* multipoles.
+ * @param epsilon_a The maximal softening length of any particle in A.
+ * @param epsilon_a The maximal softening length of any particle in B.
*/
__attribute__((always_inline, const)) INLINE static int gravity_M2L_accept(
const double r_crit_a, const double r_crit_b, const double theta_crit2,
- const double r2) {
+ const double r2, const double epsilon_a, const double epsilon_b) {
const double size = r_crit_a + r_crit_b;
const double size2 = size * size;
+ const double epsilon_a2 = epsilon_a * epsilon_a;
+ const double epsilon_b2 = epsilon_b * epsilon_b;
// MATTHIEU: Make this mass-dependent ?
/* Multipole acceptance criterion (Dehnen 2002, eq.10) */
- return (r2 * theta_crit2 > size2);
+ return (r2 * theta_crit2 > size2) && (r2 > epsilon_a2) && (r2 > epsilon_b2);
}
/**
@@ -2602,7 +2610,7 @@ __attribute__((always_inline, const)) INLINE static int gravity_M2L_accept(
* Issue 1, pp.27-42, equation 10.
*
* We also additionally check that the distance between the particle and the
- * multipole is larger than then softening length (here the distance at which
+ * multipole is larger than the softening length (here the distance at which
* the gravity becomes Newtonian again, not the Plummer-equivalent quantity).
*
* @param r_max2 The square of the size of the multipole.
diff --git a/src/runner_doiact_grav.h b/src/runner_doiact_grav.h
index 5609b295a9c416c3f9bd69bec006baff0134eac1..bf61a70457200bec4686d9eddfca8f676ddfa961 100644
--- a/src/runner_doiact_grav.h
+++ b/src/runner_doiact_grav.h
@@ -1585,7 +1585,9 @@ static INLINE void runner_dopair_recursive_grav(struct runner *r,
* option... */
/* Can we use M-M interactions ? */
- if (gravity_M2L_accept(multi_i->r_max, multi_j->r_max, theta_crit2, r2)) {
+ if (gravity_M2L_accept(multi_i->r_max, multi_j->r_max, theta_crit2, r2,
+ multi_i->m_pole.max_softening,
+ multi_j->m_pole.max_softening)) {
/* Go M-M */
runner_dopair_grav_mm(r, ci, cj);
@@ -1807,7 +1809,9 @@ static INLINE void runner_do_grav_long_range(struct runner *r, struct cell *ci,
/* Are we in charge of this cell pair? */
if (gravity_M2L_accept(multi_top->r_max_rebuild, multi_j->r_max_rebuild,
- theta_crit2, r2_rebuild)) {
+ theta_crit2, r2_rebuild,
+ multi_top->m_pole.max_softening,
+ multi_j->m_pole.max_softening)) {
/* Call the PM interaction fucntion on the active sub-cells of ci */
runner_dopair_grav_mm_nonsym(r, ci, cj);