Add a parameter to the gravity pair P-P interaction function to forbid the use of P2M interactions

examples/ZeldovichPancake_3D/zeldovichPancake.yml

@@ -46,7 +46,7 @@ Cosmology:
 Gravity:
   mesh_side_length:   16
   eta: 0.025
-  theta: 0.5
+  theta: 0.7
   r_cut_max: 6.
   comoving_softening: 0.0001
   max_physical_softening: 0.0001

src/approx_math.h

@@ -32,7 +32,7 @@
  *
  * @param x The number to take the exponential of.
  */
-__attribute__((always_inline)) INLINE static float approx_expf(float x) {
+__attribute__((always_inline, const)) INLINE static float approx_expf(float x) {
   return 1.f + x * (1.f + x * (0.5f + x * (1.f / 6.f + 1.f / 24.f * x)));
 }
 
@@ -40,7 +40,8 @@ __attribute__((always_inline)) INLINE static float approx_expf(float x) {
  * @brief Approximate version of expf(x) using a 6th order Taylor expansion
  *
  */
-__attribute__((always_inline)) INLINE static float good_approx_expf(float x) {
+__attribute__((always_inline, const)) INLINE static float good_approx_expf(
+    float x) {
   return 1.f +
          x * (1.f +
               x * (0.5f +
@@ -52,7 +53,8 @@ __attribute__((always_inline)) INLINE static float good_approx_expf(float x) {
 /**
  * @brief Approximate version of exp(x) using a 6th order Taylor expansion
  */
-__attribute__((always_inline)) INLINE static double good_approx_exp(double x) {
+__attribute__((always_inline, const)) INLINE static double good_approx_exp(
+    double x) {
   return 1. +
          x * (1. +
               x * (0.5 +

src/gravity_cache.h

@@ -171,6 +171,7 @@ __attribute__((always_inline)) INLINE static void gravity_cache_zero_output(
  * more expensive P2P.
  *
  * @param max_active_bin The largest active bin in the current time-step.
+ * @param allow_mpole Are we allowing the use of multipoles?
  * @param periodic Are we using periodic BCs ?
  * @param dim The size of the simulation volume along each dimension.
  * @param c The #gravity_cache to fill.
@@ -185,10 +186,11 @@ __attribute__((always_inline)) INLINE static void gravity_cache_zero_output(
  * @param grav_props The global gravity properties.
  */
 __attribute__((always_inline)) INLINE static void gravity_cache_populate(
-    const timebin_t max_active_bin, const int periodic, const float dim[3],
-    struct gravity_cache *c, const struct gpart *restrict gparts,
-    const int gcount, const int gcount_padded, const double shift[3],
-    const float CoM[3], const float r_max2, const struct cell *cell,
+    const timebin_t max_active_bin, const int allow_mpole, const int periodic,
+    const float dim[3], struct gravity_cache *c,
+    const struct gpart *restrict gparts, const int gcount,
+    const int gcount_padded, const double shift[3], const float CoM[3],
+    const float r_max2, const struct cell *cell,
     const struct gravity_props *grav_props) {
 
   const float theta_crit2 = grav_props->theta_crit2;
@@ -227,7 +229,7 @@ __attribute__((always_inline)) INLINE static void gravity_cache_populate(
     const float r2 = dx * dx + dy * dy + dz * dz;
 
     /* Check whether we can use the multipole instead of P-P */
-    use_mpole[i] = gravity_M2P_accept(r_max2, theta_crit2, r2);
+    use_mpole[i] = allow_mpole && gravity_M2P_accept(r_max2, theta_crit2, r2);
   }
 
 #ifdef SWIFT_DEBUG_CHECKS

src/gravity_properties.c

@@ -113,7 +113,8 @@ void gravity_props_print(const struct gravity_props *p) {
 
   message("Self-gravity opening angle:  theta=%.4f", p->theta_crit);
 
-  message("Self-gravity softening: %s", kernel_gravity_softening_name);
+  message("Self-gravity softening functional form: %s",
+          kernel_gravity_softening_name);
 
   message(
       "Self-gravity comoving softening:    epsilon=%.4f (Plummer equivalent: "

src/kernel_long_gravity.h

@@ -31,7 +31,7 @@
 #include <float.h>
 #include <math.h>
 
-//#define GADGET2_LONG_RANGE_CORRECTION
+#define GADGET2_LONG_RANGE_CORRECTION
 
 #ifdef GADGET2_LONG_RANGE_CORRECTION
 #define kernel_long_gravity_truncation_name "Gadget-like (using erfc())"

src/multipole.h

@@ -2396,7 +2396,7 @@ INLINE static void gravity_L2P(const struct grav_tensor *lb,
  * @param r2 Square of the distance (periodically wrapped) between the
  * multipoles.
  */
-__attribute__((always_inline)) INLINE static int gravity_M2L_accept(
+__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) {
 
@@ -2421,7 +2421,7 @@ __attribute__((always_inline)) INLINE static int gravity_M2L_accept(
  * @param r2 Square of the distance (periodically wrapped) between the
  * particle and the multipole.
  */
-__attribute__((always_inline)) INLINE static int gravity_M2P_accept(
+__attribute__((always_inline, const)) INLINE static int gravity_M2P_accept(
     const float r_max2, const float theta_crit2, const float r2) {
 
   // MATTHIEU: Make this mass-dependent ?

src/runner_doiact_grav.h

@@ -621,9 +621,11 @@ static INLINE void runner_dopair_grav_pm_truncated(
  * @param ci The first #cell.
  * @param cj The other #cell.
  * @param symmetric Are we updating both cells (1) or just ci (0) ?
+ * @param Are we allowing the use of P2M interactions ?
  */
 static INLINE void runner_dopair_grav_pp(struct runner *r, struct cell *ci,
-                                         struct cell *cj, const int symmetric) {
+                                         struct cell *cj, const int symmetric,
+                                         const int allow_mpole) {
 
   /* Recover some useful constants */
   const struct engine *e = r->e;
@@ -691,12 +693,12 @@ static INLINE void runner_dopair_grav_pp(struct runner *r, struct cell *ci,
 #endif
 
   /* Fill the caches */
-  gravity_cache_populate(e->max_active_bin, periodic, dim, ci_cache, ci->gparts,
-                         gcount_i, gcount_padded_i, shift_i, CoM_j, rmax2_j, ci,
-                         e->gravity_properties);
-  gravity_cache_populate(e->max_active_bin, periodic, dim, cj_cache, cj->gparts,
-                         gcount_j, gcount_padded_j, shift_j, CoM_i, rmax2_i, cj,
-                         e->gravity_properties);
+  gravity_cache_populate(e->max_active_bin, allow_mpole, periodic, dim,
+                         ci_cache, ci->gparts, gcount_i, gcount_padded_i,
+                         shift_i, CoM_j, rmax2_j, ci, e->gravity_properties);
+  gravity_cache_populate(e->max_active_bin, allow_mpole, periodic, dim,
+                         cj_cache, cj->gparts, gcount_j, gcount_padded_j,
+                         shift_j, CoM_i, rmax2_i, cj, e->gravity_properties);
 
   /* Can we use the Newtonian version or do we need the truncated one ? */
   if (!periodic) {
@@ -712,8 +714,9 @@ static INLINE void runner_dopair_grav_pp(struct runner *r, struct cell *ci,
                                  cj->gparts);
 
       /* Then the M2P */
-      runner_dopair_grav_pm_full(ci_cache, gcount_padded_i, CoM_j, multi_j,
-                                 periodic, dim, e, ci->gparts, gcount_i, cj);
+      if (allow_mpole)
+        runner_dopair_grav_pm_full(ci_cache, gcount_padded_i, CoM_j, multi_j,
+                                   periodic, dim, e, ci->gparts, gcount_i, cj);
     }
     if (cj_active && symmetric) {
 
@@ -723,8 +726,9 @@ static INLINE void runner_dopair_grav_pp(struct runner *r, struct cell *ci,
                                  ci->gparts);
 
       /* Then the M2P */
-      runner_dopair_grav_pm_full(cj_cache, gcount_padded_j, CoM_i, multi_i,
-                                 periodic, dim, e, cj->gparts, gcount_j, ci);
+      if (allow_mpole)
+        runner_dopair_grav_pm_full(cj_cache, gcount_padded_j, CoM_i, multi_i,
+                                   periodic, dim, e, cj->gparts, gcount_j, ci);
     }
 
   } else { /* Periodic BC */
@@ -751,9 +755,10 @@ static INLINE void runner_dopair_grav_pp(struct runner *r, struct cell *ci,
                                         ci->gparts, cj->gparts);
 
         /* Then the M2P */
-        runner_dopair_grav_pm_truncated(ci_cache, gcount_padded_i, CoM_j,
-                                        multi_j, dim, r_s_inv, e, ci->gparts,
-                                        gcount_i, cj);
+        if (allow_mpole)
+          runner_dopair_grav_pm_truncated(ci_cache, gcount_padded_i, CoM_j,
+                                          multi_j, dim, r_s_inv, e, ci->gparts,
+                                          gcount_i, cj);
       }
       if (cj_active && symmetric) {
 
@@ -763,9 +768,10 @@ static INLINE void runner_dopair_grav_pp(struct runner *r, struct cell *ci,
                                         cj->gparts, ci->gparts);
 
         /* Then the M2P */
-        runner_dopair_grav_pm_truncated(cj_cache, gcount_padded_j, CoM_i,
-                                        multi_i, dim, r_s_inv, e, cj->gparts,
-                                        gcount_j, ci);
+        if (allow_mpole)
+          runner_dopair_grav_pm_truncated(cj_cache, gcount_padded_j, CoM_i,
+                                          multi_i, dim, r_s_inv, e, cj->gparts,
+                                          gcount_j, ci);
       }
 
     } else {
@@ -781,8 +787,10 @@ static INLINE void runner_dopair_grav_pp(struct runner *r, struct cell *ci,
                                    ci->gparts, cj->gparts);
 
         /* Then the M2P */
-        runner_dopair_grav_pm_full(ci_cache, gcount_padded_i, CoM_j, multi_j,
-                                   periodic, dim, e, ci->gparts, gcount_i, cj);
+        if (allow_mpole)
+          runner_dopair_grav_pm_full(ci_cache, gcount_padded_i, CoM_j, multi_j,
+                                     periodic, dim, e, ci->gparts, gcount_i,
+                                     cj);
       }
       if (cj_active && symmetric) {
 
@@ -792,8 +800,10 @@ static INLINE void runner_dopair_grav_pp(struct runner *r, struct cell *ci,
                                    cj->gparts, ci->gparts);
 
         /* Then the M2P */
-        runner_dopair_grav_pm_full(cj_cache, gcount_padded_j, CoM_i, multi_i,
-                                   periodic, dim, e, cj->gparts, gcount_j, ci);
+        if (allow_mpole)
+          runner_dopair_grav_pm_full(cj_cache, gcount_padded_j, CoM_i, multi_i,
+                                     periodic, dim, e, cj->gparts, gcount_j,
+                                     ci);
       }
     }
   }
@@ -1337,7 +1347,7 @@ static INLINE void runner_dopair_recursive_grav(struct runner *r,
   } else if (!ci->split && !cj->split) {
 
     /* We have two leaves. Go P-P. */
-    runner_dopair_grav_pp(r, ci, cj, 1);
+    runner_dopair_grav_pp(r, ci, cj, /*symmetric*/ 1, /*allow_mpoles*/ 0);
 
   } else {
 
@@ -1546,13 +1556,13 @@ static INLINE void runner_do_grav_long_range(struct runner *r, struct cell *ci,
       /* Call the PM interaction fucntion on the active sub-cells of ci */
 
       // runner_dopair_recursive_grav_pm(r, ci, cj);
-      runner_dopair_grav_mm(r, ci, cj);
-      // runner_dopair_grav_pp(r, ci, cj, 0);
+      // runner_dopair_grav_mm(r, ci, cj);
+      runner_dopair_grav_pp(r, ci, cj, 0, 1);
 
     } /* We are in charge of this pair */
   }   /* Loop over top-level cells */
 
-  message("cc=%d", cc);
+  // message("cc=%d", cc);
 
   if (timer) TIMER_TOC(timer_dograv_long_range);
 }