Use the long-range truncation in the M/r^2 term of the MAC

examples/main.c

@@ -1076,7 +1076,7 @@ int main(int argc, char *argv[]) {
       gravity_props_init(&gravity_properties, params, &prog_const, &cosmo,
                          with_cosmology, with_external_gravity,
                          with_baryon_particles, with_DM_particles,
-                         with_DM_background_particles, periodic);
+                         with_DM_background_particles, periodic, s.dim);
 
     /* Initialise the external potential properties */
     bzero(&potential, sizeof(struct external_potential));

examples/main_fof.c

@@ -526,7 +526,7 @@ int main(int argc, char *argv[]) {
   gravity_props_init(&gravity_properties, params, &prog_const, &cosmo,
                      /*with_cosmology=*/1, /*with_external_gravity=*/0,
                      with_baryon_particles, with_DM_particles,
-                     with_DM_background_particles, periodic);
+                     with_DM_background_particles, periodic, s.dim);
 
   /* Initialise the long-range gravity mesh */
   if (periodic) {

src/gravity_properties.c

@@ -45,7 +45,8 @@ void gravity_props_init(struct gravity_props *p, struct swift_params *params,
                         const struct cosmology *cosmo, const int with_cosmology,
                         const int with_external_potential,
                         const int has_baryons, const int has_DM,
-                        const int is_zoom_simulation, const int periodic) {
+                        const int is_zoom_simulation, const int periodic,
+                        const double dim[3]) {
 
   /* Tree updates */
   p->rebuild_frequency =
@@ -65,6 +66,9 @@ void gravity_props_init(struct gravity_props *p, struct swift_params *params,
     p->r_cut_min_ratio = parser_get_opt_param_float(
         params, "Gravity:r_cut_min", gravity_props_default_r_cut_min);
 
+    const double r_s = p->a_smooth * dim[0] / p->mesh_size;
+    p->r_s_inv = 1. / r_s;
+
     /* Some basic checks of what we read */
     if (p->mesh_size % 2 != 0)
       error("The mesh side-length must be an even number.");
@@ -101,7 +105,7 @@ void gravity_props_init(struct gravity_props *p, struct swift_params *params,
   }
 
   /* We always start with the geometric MAC */
-  p->use_advanced_mac = 0;
+  p->use_advanced_MAC = 0;
 
   /* Geometric opening angle */
   p->theta_crit = parser_get_param_double(params, "Gravity:theta_cr");
@@ -112,6 +116,11 @@ void gravity_props_init(struct gravity_props *p, struct swift_params *params,
     p->adaptive_tolerance =
         parser_get_param_float(params, "Gravity:epsilon_fmm");
 
+  /* Consider truncated forces in the MAC? */
+  if (p->use_adaptive_tolerance)
+    p->consider_truncation_in_MAC =
+        parser_get_param_int(params, "Gravity:allow_truncation_in_MAC");
+
   /* Are we allowing tree use below softening? */
   p->use_tree_below_softening =
       parser_get_opt_param_int(params, "Gravity:use_tree_below_softening", 1);
@@ -206,7 +215,7 @@ void gravity_props_update_MAC_choice(struct gravity_props *p) {
 
   /* Now that we have run initial accelerations,
    * switch to the better MAC */
-  if (p->use_adaptive_tolerance) p->use_advanced_mac = 1;
+  if (p->use_adaptive_tolerance) p->use_advanced_MAC = 1;
 }
 
 void gravity_props_update(struct gravity_props *p,

src/gravity_properties.h

@@ -55,7 +55,7 @@ struct gravity_props {
   /* -------------- Properties of the FFM gravity ---------------------- */
 
   /*! What MAC are we currently using? */
-  int use_advanced_mac;
+  int use_advanced_MAC;
 
   /*! Are we using the adaptive opening angle? (as read from param file) */
   int use_adaptive_tolerance;
@@ -69,6 +69,9 @@ struct gravity_props {
   /*! Are we allowing tree gravity below softening? */
   int use_tree_below_softening;
 
+  /*! Are we applying long-range truncation to the forces in the MAC? */
+  int consider_truncation_in_MAC;
+
   /* ------------- Properties of the softened gravity ------------------ */
 
   /*! Co-moving softening length for for high-res. DM particles */
@@ -113,12 +116,17 @@ struct gravity_props {
    * a_smooth */
   float r_cut_max_ratio;
 
+  /*! Inverse of the long-range gravity mesh scale. */
+  float r_s_inv;
+
   /*! Are we dithering the particles at every rebuild? */
   int with_dithering;
 
   /*! Fraction of the top-level cell size used to normalize the dithering */
   double dithering_ratio;
 
+  /* ------------- Physical constants ---------------------------------- */
+
   /*! Gravitational constant (in internal units, copied from the physical
    * constants) */
   float G_Newton;
@@ -130,7 +138,8 @@ void gravity_props_init(struct gravity_props *p, struct swift_params *params,
                         const struct cosmology *cosmo, const int with_cosmology,
                         const int with_external_potential,
                         const int has_baryons, const int has_DM,
-                        const int is_zoom_simulation, const int periodic);
+                        const int is_zoom_simulation, const int periodic,
+                        const double dim[3]);
 void gravity_props_update(struct gravity_props *p,
                           const struct cosmology *cosmo);
 void gravity_props_update_MAC_choice(struct gravity_props *p);

src/mesh_gravity.c

@@ -726,7 +726,7 @@ void pm_mesh_free(struct pm_mesh* mesh) {
  * @param nr_threads The number of threads on this MPI rank.
  */
 void pm_mesh_init(struct pm_mesh* mesh, const struct gravity_props* props,
-                  double dim[3], int nr_threads) {
+                  const double dim[3], int nr_threads) {
 
 #ifdef HAVE_FFTW
 

src/mesh_gravity.h

@@ -68,7 +68,7 @@ struct pm_mesh {
 };
 
 void pm_mesh_init(struct pm_mesh *mesh, const struct gravity_props *props,
-                  double dim[3], int nr_threads);
+                  const double dim[3], int nr_threads);
 void pm_mesh_init_no_mesh(struct pm_mesh *mesh, double dim[3]);
 void pm_mesh_compute_potential(struct pm_mesh *mesh, const struct space *s,
                                struct threadpool *tp, int verbose);

src/multipole_accept.h

@@ -25,6 +25,7 @@
 /* Local includes */
 #include "binomial.h"
 #include "integer_power.h"
+#include "kernel_long_gravity.h"
 #include "minmax.h"
 #include "multipole_struct.h"
 
@@ -74,12 +75,24 @@ __attribute__((nonnull, pure)) INLINE static int gravity_M2L_accept(
     E_BA_term /= (rho_A + rho_B);
   }
 
-    /* Compute r^(p+2) */
+  float r_to_p_plus2, W;
+  if (periodic && props->consider_truncation_in_MAC) {
+
+    /* Compute r^(p+2) and the long-range correction */
+    const float r = sqrtf(r2);
+    r_to_p_plus2 = integer_powf(r, (p + 2));
+    W = kernel_long_grav_force_eval(r * props->r_s_inv);
+
+  } else {
+
+  /* Compute r^(p+2) */
 #if SELF_GRAVITY_MULTIPOLE_ORDER % 2 == 1
-  const float r_to_p_plus2 = integer_powf(sqrtf(r2), (p + 2));
+    r_to_p_plus2 = integer_powf(sqrtf(r2), (p + 2));
 #else
-  const float r_to_p_plus2 = integer_powf(r2, ((p / 2) + 1));
+    r_to_p_plus2 = integer_powf(r2, ((p / 2) + 1));
 #endif
+    W = 1.f;
+  }
 
   /* Get the mimimal acceleration in A */
   const float min_a_grav = A->m_pole.min_old_a_grav_norm;
@@ -98,7 +111,7 @@ __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) {
 
 #ifdef SWIFT_DEBUG_CHECKS
     if (min_a_grav == 0.) error("Acceleration is 0");
@@ -114,8 +127,8 @@ __attribute__((nonnull, pure)) INLINE static int gravity_M2L_accept(
         props->use_tree_below_softening || 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;
+     * (E_BA * (1 / r^(p)) * ((1 / r^2) * W) < eps * a_min) */
+    const int cond_3 = E_BA_term * W < eps * min_a_grav * r_to_p_plus2;
 
     return cond_1 && cond_2 && cond_3;
 
@@ -185,12 +198,24 @@ __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];
 
+  float r_to_p_plus2, W;
+  if (periodic && props->consider_truncation_in_MAC) {
+
+    /* Compute r^(p+2) and the long-range correction */
+    const float r = sqrtf(r2);
+    r_to_p_plus2 = integer_powf(r, (p + 2));
+    W = kernel_long_grav_force_eval(r * props->r_s_inv);
+
+  } else {
+
   /* Compute r^(p+2) */
 #if SELF_GRAVITY_MULTIPOLE_ORDER % 2 == 1
-  const float r_to_p_plus2 = integer_powf(sqrtf(r2), (p + 2));
+    r_to_p_plus2 = integer_powf(sqrtf(r2), (p + 2));
 #else
-  const float r_to_p_plus2 = integer_powf(r2, ((p / 2) + 1));
+    r_to_p_plus2 = integer_powf(r2, ((p / 2) + 1));
 #endif
+    W = 1.f;
+  }
 
   /* Get the estimate of the acceleration */
   const float old_a_grav = pa->old_a_grav_norm;
@@ -206,7 +231,7 @@ __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) {
 
 #ifdef SWIFT_DEBUG_CHECKS
     if (old_a_grav == 0.) error("Acceleration is 0");
@@ -222,8 +247,8 @@ __attribute__((nonnull, pure)) INLINE static int gravity_M2P_accept(
         props->use_tree_below_softening || 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;
+     * (E_BA * (1 / r^(p)) * ((1 / r^2) * W) < eps * a_min) */
+    const int cond_3 = E_BA_term * W < eps * old_a_grav * r_to_p_plus2;
 
     return cond_1 && cond_2 && cond_3;