Fix typo in the new f_MAC expression.

src/gravity_properties.c

@@ -67,6 +67,7 @@ void gravity_props_init(struct gravity_props *p, struct swift_params *params,
         params, "Gravity:r_cut_min", gravity_props_default_r_cut_min);
 
     p->r_s = p->a_smooth * dim[0] / p->mesh_size;
+    p->r_s_inv = 1. / p->r_s;
 
     /* Some basic checks of what we read */
     if (p->mesh_size % 2 != 0)
@@ -82,6 +83,7 @@ void gravity_props_init(struct gravity_props *p, struct swift_params *params,
     p->mesh_size = 0;
     p->a_smooth = 0.f;
     p->r_s = FLT_MAX;
+    p->r_s_inv = 0.f;
     p->r_cut_min_ratio = 0.f;
     p->r_cut_max_ratio = 0.f;
   }

src/gravity_properties.h

@@ -119,6 +119,9 @@ struct gravity_props {
   /*! Long-range gravity mesh scale. */
   float r_s;
 
+  /*! Inverse of the long-range gravity mesh scale. */
+  float r_s_inv;
+
   /*! Are we dithering the particles at every rebuild? */
   int with_dithering;
 

src/multipole_accept.h

@@ -33,22 +33,26 @@
 /**
  * @brief Compute the inverse of the force estimator entering the MAC
  *
+ * Note that in the unsofted case, the first condition is naturally
+ * never reached (as H == 0). In the non-periodic (non-truncated) case
+ * the second condition is never reached (as r_s == inf, r_s_inv == 0).
+ *
  * @param H The spline softening length.
- * @param r_s The scale of the gravity mesh.
+ * @param r_s_inv The inverse of the scale of the gravity mesh.
  * @param r2 The square of the distance between the multipoles.
  */
 __attribute__((const)) INLINE static float gravity_f_MAC_inverse(
-    const float H, const float r_s, const float r2) {
+    const float H, const float r_s_inv, const float r2) {
 
   if (r2 < (25.f / 81.f) * H * H) {
 
     /* Below softening radius */
     return (25.f / 81.f) * H * H;
 
-  } else if (r2 > (25.f / 9.f) * r_s * r_s) {
+  } else if (r_s_inv * r_s_inv * r2 > (25.f / 9.f)) {
 
     /* Above truncation radius */
-    return (25.f / 9.f) * r_s * r_s * r2 * r2;
+    return (9.f / 25.f) * r_s_inv * r_s_inv * r2 * r2;
 
   } else {
 
@@ -114,7 +118,7 @@ __attribute__((nonnull, pure)) INLINE static int gravity_M2L_accept(
 
   float f_MAC_inv;
   if (props->consider_truncation_in_MAC) {
-    f_MAC_inv = gravity_f_MAC_inverse(max_softening, props->r_s, r2);
+    f_MAC_inv = gravity_f_MAC_inverse(max_softening, props->r_s_inv, r2);
   } else {
     f_MAC_inv = r2;
   }
@@ -232,7 +236,7 @@ __attribute__((nonnull, pure)) INLINE static int gravity_M2P_accept(
 
   float f_MAC_inv;
   if (props->consider_truncation_in_MAC) {
-    f_MAC_inv = gravity_f_MAC_inverse(max_softening, props->r_s, r2);
+    f_MAC_inv = gravity_f_MAC_inverse(max_softening, props->r_s_inv, r2);
   } else {
     f_MAC_inv = r2;
   }

src/runner_doiact_grav.c

@@ -247,7 +247,7 @@ static INLINE void runner_dopair_grav_pp_truncated_no_cache(
     error("Calling truncated PP function in non-periodic setup.");
 #endif
 
-  const float r_s_inv = 1. / grav_props->r_s;
+  const float r_s_inv = grav_props->r_s_inv;
 
   /* Loop over sink particles */
   for (int i = 0; i < gcount_i; ++i) {

theory/Multipoles/fmm_mac.tex

@@ -184,7 +184,7 @@ f_{\rm MAC}(r) =
   \left(\frac{9}{5}\right)^2 H^{-2} & \mbox{if} & r <
   \frac{5}{9}H,\\
   r^{-2} & \mbox{if} & \frac{5}{9}H \leq r < \frac{5}{3}r_s, \\
-  \left(\frac{3}{5}\right)^2 r_s^{-2} r^{-4} & \mbox{if} & \frac{5}{3}r_s \leq r. \\
+  \left(\frac{5}{3}\right)^2 r_s^2 r^{-4} & \mbox{if} & \frac{5}{3}r_s \leq r. \\
 \end{array}
 \right.
 \label{eq:fmm:f_mac}
@@ -202,12 +202,12 @@ acceptance criterion instead of the $1/|\mathbf{R}|$ term:
     A}\left(|\mathbf{a}_a|\right).
   \label{eq:fmm:mac_f_mac}  
 \end{equation}
-The same change is applied to the MAC used of the M2P kernel
-(eq. \ref{eq:fmm:mac_m2p}). In the non-truncated un-softened case, this
-expression reduces to \citep{Dehnen2014} one. Using this expression
-instead of the simpler Newtonian one only makes a difference in
-simulations where a lot of particles cluster below the scale of the
-softening, which is often the case for hydrodynamical simulations
-including radiative cooling processes. The use of this term over the
-simpler $1/r^2$ estimator is a runtime parameter.
+The same change is applied to the MAC used for the M2P kernel
+(eq. \ref{eq:fmm:mac_m2p}). In the non-truncated un-softened case,
+their expressions reduce to the \citep{Dehnen2014} one. Using this
+$f_{\rm MAC}$ instead of the simpler purely-Newtonian one only makes a
+difference in simulations where a lot of particles cluster below the
+scale of the softening, which is often the case for hydrodynamical
+simulations including radiative cooling processes. The use of this
+term over the simpler $1/r^2$ estimator is a runtime parameter.