From a651fa110bdb9dc7df89c93aad8ce0e9f58b0727 Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <schaller@strw.leidenuniv.nl>
Date: Thu, 23 Apr 2020 11:52:41 +0200
Subject: [PATCH] Implemented the new updated MAC
---
src/gravity_properties.c | 3 +-
src/gravity_properties.h | 16 ++--
src/multipole_accept.h | 179 +++++++++++++++++++++++++++++----------
3 files changed, 143 insertions(+), 55 deletions(-)
diff --git a/src/gravity_properties.c b/src/gravity_properties.c
index cd7a890917..4297ca1b45 100644
--- a/src/gravity_properties.c
+++ b/src/gravity_properties.c
@@ -32,6 +32,7 @@
#include "gravity.h"
#include "kernel_gravity.h"
#include "kernel_long_gravity.h"
+#include "restart.h"
#define gravity_props_default_a_smooth 1.25f
#define gravity_props_default_r_cut_max 4.5f
@@ -86,8 +87,6 @@ void gravity_props_init(struct gravity_props *p, struct swift_params *params,
/* Opening angle */
p->theta_crit = parser_get_param_double(params, "Gravity:theta");
if (p->theta_crit >= 1.) error("Theta too large. FMM won't converge.");
- p->theta_crit2 = p->theta_crit * p->theta_crit;
- p->theta_crit_inv = 1. / p->theta_crit;
/* Mesh dithering */
if (periodic && !with_external_potential) {
diff --git a/src/gravity_properties.h b/src/gravity_properties.h
index 8a4abe4bd3..56f69bd253 100644
--- a/src/gravity_properties.h
+++ b/src/gravity_properties.h
@@ -26,10 +26,6 @@
#include <hdf5.h>
#endif
-/* Local includes. */
-#include "kernel_gravity.h"
-#include "restart.h"
-
/* Forward declarations */
struct cosmology;
struct phys_const;
@@ -58,14 +54,14 @@ struct gravity_props {
/* -------------- Properties of the FFM gravity ---------------------- */
- /*! Tree opening angle (Multipole acceptance criterion) */
- double theta_crit;
+ /*! Are we using the adaptive opening angle? */
+ int use_adaptive_tolerance;
- /*! Square of opening angle */
- double theta_crit2;
+ /*! Accuracy parameter of the advanced MAC */
+ float adaptive_tolerance;
- /*! Inverse of opening angle */
- double theta_crit_inv;
+ /*! Tree opening angle (Multipole acceptance criterion) */
+ double theta_crit;
/* ------------- Properties of the softened gravity ------------------ */
diff --git a/src/multipole_accept.h b/src/multipole_accept.h
index a673098639..272b470bd4 100644
--- a/src/multipole_accept.h
+++ b/src/multipole_accept.h
@@ -23,67 +23,160 @@
#include "../config.h"
/* Local includes */
+#include "binomial.h"
+#include "integer_power.h"
+#include "minmax.h"
#include "multipole_struct.h"
/**
- * @brief Checks whether a cell-cell interaction can be appromixated by a M-M
- * interaction using the distance and cell radius.
+ * @brief Checks whether The multipole in B can be used to update the field
+ * tensor in A.
*
- * We use the multipole acceptance criterion of Dehnen, 2002, JCoPh, Volume 179,
- * Issue 1, pp.27-42, equation 10.
+ * We use the MAC of Dehnen 2014 eq. 16.
*
- * 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).
+ * Note: this is *not* symmetric in A<->B unless the purely geometric criterion
+ * is used.
*
- * @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_b The maximal softening length of any particle in B.
+ * @param props The properties of the gravity scheme.
+ * @param A The gravity tensors that we want to update (sink).
+ * @param B The gravity tensors that act as a source.
+ * @param r2 The square of the distance between the centres of mass of A and 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 epsilon_a, const double epsilon_b) {
+__attribute__((nonnull, pure)) INLINE static int gravity_M2L_accept(
+ const struct gravity_props *props, const struct gravity_tensors *restrict A,
+ const struct gravity_tensors *restrict B, const float r2) {
- 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;
+ /* Order of the expansion */
+ const int p = SELF_GRAVITY_MULTIPOLE_ORDER;
- // MATTHIEU: Make this mass-dependent ?
+ /* 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) {
+ E_BA_term +=
+ binomial(p, n) * B->m_pole.power[n] * integer_powf(A->r_max, p - n);
+ }
+ E_BA_term *= 8.f;
+ E_BA_term *= max(A->r_max, B->r_max);
+ E_BA_term /= (A->r_max + B->r_max);
- /* Multipole acceptance criterion (Dehnen 2002, eq.10) */
- return (r2 * theta_crit2 > size2) && (r2 > epsilon_a2) && (r2 > epsilon_b2);
+ /* Compute r^(p+2) */
+#if SELF_GRAVITY_MULTIPOLE_ORDER % 2 == 1
+ const float r_to_p_plus2 = integer_powf(sqrtf(r2), (p + 2));
+#else
+ const float r_to_p_plus2 = integer_powf(r2, ((p / 2) + 1));
+#endif
+
+ /* Get the mimimal acceleration in A */
+ const float min_a_grav = A->m_pole.min_old_a_grav_norm;
+
+ /* Get the maximal softening length in B */
+ const float max_softening = B->m_pole.max_softening;
+
+ /* Get the relative tolerance */
+ const float eps = props->adaptive_tolerance;
+
+ /* Get the basic geometric critical angle */
+ const float theta_crit = props->theta_crit;
+ const float theta_crit2 = theta_crit * theta_crit;
+
+ /* Get the sum of the multipole sizes */
+ const float rho_sum = A->r_max + B->r_max;
+
+ if (props->use_adaptive_tolerance) {
+
+ /* Test the different conditions */
+
+ /* Condition 1: We are in the converging part of the Taylor expansion */
+ const int cond_1 = rho_sum * rho_sum < r2;
+
+ /* Condition 2: We are not below softening */
+ const int cond_2 = max_softening * max_softening < r2;
+
+ /* Condition 3: The contribution is accurate enough
+ * (E_BA / r^(p+2) < eps a_min) */
+ const int cond_3 = E_BA_term < eps * min_a_grav * r_to_p_plus2;
+
+ return cond_1 && cond_2 && cond_3;
+
+ } else {
+
+ /* Condition 1: We are in the converging part of the Taylor expansion */
+ const int cond_1 = rho_sum * rho_sum < theta_crit2 * r2;
+
+ /* Condition 2: We are not below softening */
+ const int cond_2 = max_softening * max_softening < r2;
+
+ return cond_1 && cond_2;
+ }
}
/**
- * @brief Checks whether a particle-cell interaction can be appromixated by a
- * M2P interaction using the distance and cell radius.
- *
- * We use the multipole acceptance criterion of Dehnen, 2002, JCoPh, Volume 179,
- * Issue 1, pp.27-42, equation 10.
+ * @brief Checks whether The multipole in B can be used to update the particle
+ * pa
*
- * We also additionally check that the distance between the particle and the
- * multipole is larger than the softening length (here the distance at which
- * the gravity becomes Newtonian again, not the Plummer-equivalent quantity).
+ * We use the MAC of Dehnen 2014 eq. 16.
*
- * @param r_max2 The square of the size of the multipole.
- * @param theta_crit2 The square of the critical opening angle.
- * @param r2 Square of the distance (periodically wrapped) between the
- * particle and the multipole.
- * @param epsilon The softening length of the particle.
+ * @param props The properties of the gravity scheme.
+ * @param pa The particle we want to compute forces for (sink)
+ * @param B The gravity tensors that act as a source.
+ * @param r2 The square of the distance between pa and the centres of mass of B.
*/
-__attribute__((always_inline, const)) INLINE static int gravity_M2P_accept(
- const float r_max2, const float theta_crit2, const float r2,
- const float epsilon) {
+__attribute__((nonnull, pure)) INLINE static int gravity_M2P_accept(
+ const struct gravity_props *props, const struct gpart *pa,
+ const struct gravity_tensors *B, const float r2) {
+
+ /* Order of the expansion */
+ const int p = SELF_GRAVITY_MULTIPOLE_ORDER;
+
+ /* Compute the error estimator (without the 1/M_B term that cancels out) */
+ float E_BA_term = 8.f * B->m_pole.power[p];
+
+ /* Compute r^(p+2) */
+#if SELF_GRAVITY_MULTIPOLE_ORDER % 2 == 1
+ const float r_to_p_plus2 = integer_powf(sqrtf(r2), (p + 2));
+#else
+ const float r_to_p_plus2 = integer_powf(r2, ((p / 2) + 1));
+#endif
+
+ /* Get the estimate of the acceleration */
+ const float old_a_grav = pa->old_a_grav_norm;
+
+ /* Get the maximal softening length in B */
+ const float max_softening = B->m_pole.max_softening;
+
+ /* Get the relative tolerance */
+ const float eps = props->adaptive_tolerance;
+
+ /* Get the basic geometric critical angle */
+ const float theta_crit = props->theta_crit;
+ const float theta_crit2 = theta_crit * theta_crit;
+
+ if (props->use_adaptive_tolerance) {
+
+ /* 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;
+
+ /* Condition 2: We are not below softening */
+ const int cond_2 = max_softening * max_softening < r2;
+
+ /* Condition 3: The contribution is accurate enough
+ * (E_BA / r^(p+2) < eps * a) */
+ const int cond_3 = E_BA_term < eps * old_a_grav * r_to_p_plus2;
+
+ return cond_1 && cond_2 && cond_3;
+
+ } else {
+
+ /* Condition 1: We are in the converging part of the Taylor expansion */
+ const int cond_1 = (B->r_max) * (B->r_max) < theta_crit2 * r2;
- // MATTHIEU: Make this mass-dependent ?
+ /* Condition 2: We are not below softening */
+ const int cond_2 = max_softening * max_softening < r2;
- /* Multipole acceptance criterion (Dehnen 2002, eq.10) */
- return (r2 * theta_crit2 > r_max2) && (r2 > epsilon * epsilon);
+ return cond_1 && cond_2;
+ }
}
#endif /* SWIFT_MULTIPOLE_ACCEPT_H */
--
GitLab