Merge branch 'gravity_order2_mac' into 'master'

Add Gadget-4 MAC

See merge request !1679

doc/RTD/source/ParameterFiles/parameter_description.rst

@@ -259,14 +259,16 @@ The accuracy of the gravity calculation is governed by the following four parame
 * The accuracy criterion used in the adaptive MAC:  :math:`\epsilon_{\rm fmm}`: ``epsilon_fmm``,
 * The time-step size pre-factor :math:`\eta`: ``eta``,
 
-The first three parameters govern the way the Fast-Multipole method
-tree-walk is done (see the theory documents for full details).  The ``MAC``
-parameter can take two values: ``adaptive`` or ``geometric``. In the first
-case, the tree recursion decision is based on the estimated accelerations
-that a given tree node will produce, trying to recurse to levels where the
-fractional contribution of the accelerations to the cell is less than
-:math:`\epsilon_{\rm fmm}`. In the second case, a fixed Barnes-Hut-like
-opening angle :math:`\theta_{\rm cr}` is used.
+The first three parameters govern the way the Fast-Multipole method tree-walk is
+done (see the theory documents for full details).  The ``MAC`` parameter can
+take three values: ``adaptive``, ``geometric``, or ``gadget``. In the first
+case, the tree recursion decision is based on the estimated accelerations that a
+given tree node will produce, trying to recurse to levels where the fractional
+contribution of the accelerations to the cell is less than :math:`\epsilon_{\rm
+fmm}`. In the second case, a fixed Barnes-Hut-like opening angle
+:math:`\theta_{\rm cr}` is used. The final case corresponds to the choice made
+in the Gadget-4 code. It is an implementation using eq. 36 of `Springel et
+al. (2021) <https://adsabs.harvard.edu/abs/2021MNRAS.506.2871S>`_.
 
 The time-step of a given particle is given by :math:`\Delta t =
 \sqrt{2\eta\epsilon_i/|\overrightarrow{a}_i|}`, where

src/gravity_properties.c

@@ -37,7 +37,7 @@
 
 #define gravity_props_default_a_smooth 1.25f
 #define gravity_props_default_r_cut_max 4.5f
-#define gravity_props_default_r_cut_min 0.1f
+#define gravity_props_default_r_cut_min 0.0f
 #define gravity_props_default_rebuild_frequency 0.01f
 #define gravity_props_default_rebuild_active_fraction 1.01f  // > 1 means never
 #define gravity_props_default_distributed_mesh 0
@@ -124,12 +124,16 @@ void gravity_props_init(struct gravity_props *p, struct swift_params *params,
 
   if (strcmp(buffer, "adaptive") == 0) {
     p->use_adaptive_tolerance = 1;
+    p->use_gadget_tolerance = 0;
+  } else if (strcmp(buffer, "gadget") == 0) {
+    p->use_adaptive_tolerance = 1;
+    p->use_gadget_tolerance = 1;
   } else if (strcmp(buffer, "geometric") == 0) {
     p->use_adaptive_tolerance = 0;
   } else {
     error(
         "Invalid choice of multipole acceptance criterion: '%s'. Should be "
-        "'adaptive' or 'geometric'",
+        "'adaptive', 'gadget', or 'geometric'",
         buffer);
   }
 
@@ -295,9 +299,15 @@ void gravity_props_print(const struct gravity_props *p) {
   message("Self-gravity time integration: eta=%.4f", p->eta);
 
   if (p->use_adaptive_tolerance) {
-    message("Self-gravity opening angle scheme:  adaptive");
-    message("Self-gravity opening angle:  epsilon_fmm=%.6f",
-            p->adaptive_tolerance);
+    if (p->use_gadget_tolerance) {
+      message("Self-gravity opening angle scheme:  Gadget");
+      message("Self-gravity opening angle:  epsilon_fmm=%.6f",
+              p->adaptive_tolerance);
+    } else {
+      message("Self-gravity opening angle scheme:  adaptive");
+      message("Self-gravity opening angle:  epsilon_fmm=%.6f",
+              p->adaptive_tolerance);
+    }
   } else {
     message("Self-gravity opening angle scheme:  fixed");
     message("Self-gravity opening angle:  theta_cr=%.4f", p->theta_crit);

src/gravity_properties.h

@@ -66,6 +66,10 @@ struct gravity_props {
   /*! Are we using the adaptive opening angle? (as read from param file) */
   int use_adaptive_tolerance;
 
+  /*! Are we using the Gadget adaptive opening angle? (as read from param file)
+   */
+  int use_gadget_tolerance;
+
   /*! Accuracy parameter of the advanced MAC */
   float adaptive_tolerance;
 

src/multipole_accept.h

@@ -84,12 +84,15 @@ __attribute__((nonnull, pure)) INLINE static int gravity_M2L_accept(
     const int use_rebuild_sizes, const int periodic) {
 
   /* Order of the expansion */
-  const int p = SELF_GRAVITY_MULTIPOLE_ORDER;
+  const int p = 2;
 
   /* Sizes of the multipoles */
   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;
 
+  /* Max size of both multipoles */
+  const float rho_max = max(rho_A, rho_B);
+
   /* Get the softening */
   const float max_softening =
       max(A->m_pole.max_softening, B->m_pole.max_softening);
@@ -102,19 +105,16 @@ __attribute__((nonnull, pure)) INLINE static int gravity_M2L_accept(
   }
   E_BA_term *= 8.f;
   if (rho_A + rho_B > 0.f) {
-    E_BA_term *= max(rho_A, rho_B);
+    E_BA_term *= rho_max;
+    ;
     E_BA_term /= (rho_A + rho_B);
   }
 
-  /* Compute r^p */
-#if SELF_GRAVITY_MULTIPOLE_ORDER % 2 == 1
-  const float r_to_p = integer_powf(sqrtf(r2), p);
-#else
+  /* Compute r^p = (r^2)^(p/2) */
   const float r_to_p = integer_powf(r2, (p / 2));
-#endif
 
   float f_MAC_inv;
-  if (props->consider_truncation_in_MAC) {
+  if (periodic && props->consider_truncation_in_MAC) {
     f_MAC_inv = gravity_f_MAC_inverse(max_softening, props->r_s_inv, r2);
   } else {
     f_MAC_inv = r2;
@@ -123,6 +123,9 @@ __attribute__((nonnull, pure)) INLINE static int gravity_M2L_accept(
   /* Get the mimimal acceleration in A */
   const float min_a_grav = A->m_pole.min_old_a_grav_norm;
 
+  /* Maximal mass */
+  const float M_max = max(A->m_pole.M_000, B->m_pole.M_000);
+
   /* Get the relative tolerance */
   const float eps = props->adaptive_tolerance;
 
@@ -133,7 +136,19 @@ __attribute__((nonnull, pure)) INLINE static int gravity_M2L_accept(
   /* Get the sum of the multipole sizes */
   const float rho_sum = rho_A + rho_B;
 
-  if (props->use_advanced_MAC) {
+  if (props->use_advanced_MAC && props->use_gadget_tolerance) {
+
+    /* Gadget 4 paper -- eq. 36 */
+    const int power = SELF_GRAVITY_MULTIPOLE_ORDER - 1;
+    const float ratio = integer_powf(rho_max / sqrtf(r2), power);
+    const int cond_1 = M_max * ratio < eps * min_a_grav * f_MAC_inv;
+
+    const int cond_2 =
+        props->use_tree_below_softening || max_softening * max_softening < r2;
+
+    return cond_1 && cond_2;
+
+  } else if (props->use_advanced_MAC && !props->use_gadget_tolerance) {
 
 #ifdef SWIFT_DEBUG_CHECKS
     if (min_a_grav == 0.) error("Acceleration is 0");
@@ -277,7 +292,7 @@ __attribute__((nonnull, pure)) INLINE static int gravity_M2P_accept(
     const struct gravity_tensors *B, const float r2, const int periodic) {
 
   /* Order of the expansion */
-  const int p = SELF_GRAVITY_MULTIPOLE_ORDER;
+  const int p = 2;
 
   /* Sizes of the multipoles */
   const float rho_B = B->r_max;
@@ -293,15 +308,11 @@ __attribute__((nonnull, pure)) INLINE static int gravity_M2P_accept(
   /* 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];
 
-  /* Compute r^p */
-#if SELF_GRAVITY_MULTIPOLE_ORDER % 2 == 1
-  const float r_to_p = integer_powf(sqrtf(r2), p);
-#else
+  /* Compute r^p = (r^2)^(p/2) */
   const float r_to_p = integer_powf(r2, (p / 2));
-#endif
 
   float f_MAC_inv;
-  if (props->consider_truncation_in_MAC) {
+  if (periodic && props->consider_truncation_in_MAC) {
     f_MAC_inv = gravity_f_MAC_inverse(max_softening, props->r_s_inv, r2);
   } else {
     f_MAC_inv = r2;
@@ -317,7 +328,19 @@ __attribute__((nonnull, pure)) INLINE static int gravity_M2P_accept(
   const float theta_crit = props->theta_crit;
   const float theta_crit2 = theta_crit * theta_crit;
 
-  if (props->use_advanced_MAC) {
+  if (props->use_advanced_MAC && props->use_gadget_tolerance) {
+
+    /* Gadget 4 paper -- eq. 12 */
+    const int power = SELF_GRAVITY_MULTIPOLE_ORDER;
+    const float ratio = integer_powf(rho_B / sqrtf(r2), power);
+    const int cond_1 = B->m_pole.M_000 * ratio < eps * old_a_grav * f_MAC_inv;
+
+    const int cond_2 =
+        props->use_tree_below_softening || max_softening * max_softening < r2;
+
+    return cond_1 && cond_2;
+
+  } else if (props->use_advanced_MAC && !props->use_gadget_tolerance) {
 
 #ifdef SWIFT_DEBUG_CHECKS
     if (old_a_grav == 0.) error("Acceleration is 0");