From cbfefe63d568c647bff5fc377cc8a8fbac117a9f Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <schaller@strw.leidenuniv.nl>
Date: Sat, 10 Aug 2019 17:33:02 +0200
Subject: [PATCH] Make the M2L MAC check that the distance between the two
multipoles is larger tan the softenings
---
src/cell.c | 11 +++++++++--
src/engine.c | 11 +++++++++--
src/multipole.h | 14 +++++++++++---
src/runner_doiact_grav.h | 8 ++++++--
4 files changed, 35 insertions(+), 9 deletions(-)
diff --git a/src/cell.c b/src/cell.c
index dfae80d5a2..0b949f3a27 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 4afbeb9b2f..e9c250939c 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 0c699ea4d7..956ec3e47a 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 5609b295a9..bf61a70457 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);
--
GitLab