Improved documentation and removed debugging code.

src/gravity/Default/gravity_iact.h

@@ -44,12 +44,8 @@
  * @param pot_ij (return) The potential.
  */
 __attribute__((always_inline)) INLINE static void runner_iact_grav_pp_full(
-    float r2, float h2, float h_inv, float h_inv3, float mass, float *f_ij,
-    float *pot_ij) {
-
-  /* *f_ij = 0.f; */
-  /* *pot_ij = 0.f; */
-  /* return; */
+    const float r2, const float h2, const float h_inv, const float h_inv3,
+    const float mass, float *f_ij, float *pot_ij) {
 
   /* Get the inverse distance */
   const float r_inv = 1.f / sqrtf(r2 + FLT_MIN);
@@ -93,12 +89,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pp_full(
  * @param pot_ij (return) The potential.
  */
 __attribute__((always_inline)) INLINE static void runner_iact_grav_pp_truncated(
-    float r2, float h2, float h_inv, float h_inv3, float mass, float rlr_inv,
-    float *f_ij, float *pot_ij) {
-
-  /* *f_ij = 0.f; */
-  /* *pot_ij = 0.f; */
-  /* return; */
+    const float r2, const float h2, const float h_inv, const float h_inv3,
+    const float mass, const float rlr_inv, float *f_ij, float *pot_ij) {
 
   /* Get the inverse distance */
   const float r_inv = 1.f / sqrtf(r2 + FLT_MIN);
@@ -134,10 +126,10 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pp_truncated(
 }
 
 /**
- * @brief Computes the force at a point generated by a multipole.
+ * @brief Computes the forces at a point generated by a multipole.
  *
- * This uses the quadrupole terms only and defaults to the monopole if
- * the code is compiled with low-order gravity only.
+ * This uses the quadrupole and trace of the octupole terms only and defaults to
+ * the monopole if the code is compiled with low-order gravity only.
  *
  * @param r_x x-component of the distance vector to the multipole.
  * @param r_y y-component of the distance vector to the multipole.
@@ -152,19 +144,22 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pp_truncated(
  * @param pot (return) The potential.
  */
 __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_full(
-    float r_x, float r_y, float r_z, float r2, float h, float h_inv,
-    const struct multipole *m, float *f_x, float *f_y, float *f_z, float *pot) {
+    const float r_x, const float r_y, const float r_z, const float r2,
+    const float h, const float h_inv, const struct multipole *m, float *f_x,
+    float *f_y, float *f_z, float *pot) {
 
-  /* In the case where the order is < 3, then there is only a monopole term left. */
-  /* We can default to the normal P-P interaction with the mass of the multipole */
-  /* and its CoM as the "particle" property */
+/* In the case where the order is < 3, then there is only a monopole term left.
+ * We can default to the normal P-P interaction with the mass of the multipole
+ * and its CoM as the "particle" property */
 #if SELF_GRAVITY_MULTIPOLE_ORDER < 3
+
   float f_ij;
   runner_iact_grav_pp_full(r2, h * h, h_inv, h_inv * h_inv * h_inv, m->M_000,
                            &f_ij, pot);
   *f_x = f_ij * r_x;
   *f_y = f_ij * r_y;
   *f_z = f_ij * r_z;
+
 #else
 
   /* Get the inverse distance */
@@ -172,7 +167,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_full(
 
   /* Compute the derivatives of the potential */
   struct potential_derivatives_M2P d;
-  compute_potential_derivatives_M2P(r_x, r_y, r_z, r2, r_inv, h, h_inv, 0, FLT_MAX, &d);
+  compute_potential_derivatives_M2P(r_x, r_y, r_z, r2, r_inv, h, h_inv, 0,
+                                    FLT_MAX, &d);
 
   /* 1st order terms (monopole) */
   *f_x = m->M_000 * d.D_100;
@@ -199,7 +195,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_full(
   *f_z += m->M_300 * d.D_301 + m->M_030 * d.D_031 + m->M_003 * d.D_004;
   *pot -= m->M_300 * d.D_100 + m->M_030 * d.D_030 + m->M_003 * d.D_003;
 #endif
-  
+
   /* Take care of the the sign convention */
   *f_x *= -1.f;
   *f_y *= -1.f;
@@ -208,15 +204,36 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_full(
 #endif
 }
 
+/**
+ * @brief Computes the forces at a point generated by a multipole, truncated for
+ * long-range periodicity.
+ *
+ * This uses the quadrupole term and trace of the octupole terms only and
+ * defaults to the monopole if the code is compiled with low-order gravity only.
+ *
+ * @param r_x x-component of the distance vector to the multipole.
+ * @param r_y y-component of the distance vector to the multipole.
+ * @param r_z z-component of the distance vector to the multipole.
+ * @param r2 Square of the distance vector to the multipole.
+ * @param h The softening length.
+ * @param h_inv Inverse of the softening length.
+ * @param rlr_inv The inverse of the gravity mesh-smoothing scale.
+ * @param m The multipole.
+ * @param f_x (return) The x-component of the acceleration.
+ * @param f_y (return) The y-component of the acceleration.
+ * @param f_z (return) The z-component of the acceleration.
+ * @param pot (return) The potential.
+ */
 __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_truncated(
-    float r_x, float r_y, float r_z, float r2, float h, float h_inv,
-    float rlr_inv, const struct multipole *m, float *f_x, float *f_y,
-    float *f_z, float *pot) {
+    const float r_x, const float r_y, const float r_z, const float r2,
+    const float h, const float h_inv, const float rlr_inv,
+    const struct multipole *m, float *f_x, float *f_y, float *f_z, float *pot) {
 
-  /* In the case where the order is < 3, then there is only a monopole term left. */
-  /* We can default to the normal P-P interaction with the mass of the multipole */
-  /* and its CoM as the "particle" property */
+/* In the case where the order is < 3, then there is only a monopole term left.
+ * We can default to the normal P-P interaction with the mass of the multipole
+ * and its CoM as the "particle" property */
 #if SELF_GRAVITY_MULTIPOLE_ORDER < 3
+
   float f_ij;
   runner_iact_grav_pp_truncated(r2, h * h, h_inv, h_inv * h_inv * h_inv,
                                 m->M_000, rlr_inv, &f_ij, pot);
@@ -231,7 +248,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_truncated(
 
   /* Compute the derivatives of the potential */
   struct potential_derivatives_M2P d;
-  compute_potential_derivatives_M2P(r_x, r_y, r_z, r2, r_inv, h, h_inv, 1, rlr_inv, &d);
+  compute_potential_derivatives_M2P(r_x, r_y, r_z, r2, r_inv, h, h_inv, 1,
+                                    rlr_inv, &d);
 
   /* 1st order terms (monopole) */
   *f_x = m->M_000 * d.D_100;
@@ -258,7 +276,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_truncated(
   *f_z += m->M_300 * d.D_301 + m->M_030 * d.D_031 + m->M_003 * d.D_004;
   *pot -= m->M_300 * d.D_100 + m->M_030 * d.D_030 + m->M_003 * d.D_003;
 #endif
-  
+
   /* Take care of the the sign convention */
   *f_x *= -1.f;
   *f_y *= -1.f;

src/gravity_cache.h

@@ -142,6 +142,8 @@ static INLINE void gravity_cache_init(struct gravity_cache *c,
  * more expensive P2P.
  *
  * @param max_active_bin The largest active bin in the current time-step.
+ * @param dim The size of the simulation volume along each dimension.
+ * @param periodic Are we using periodic BCs ?
  * @param c The #gravity_cache to fill.
  * @param gparts The #gpart array to read from.
  * @param gcount The number of particles to read.
@@ -154,10 +156,10 @@ static INLINE void gravity_cache_init(struct gravity_cache *c,
  * @param grav_props The global gravity properties.
  */
 __attribute__((always_inline)) INLINE static void gravity_cache_populate(
-									 const timebin_t max_active_bin, const float dim[3], const int periodic, 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 float dim[3], const int periodic,
+    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;
@@ -184,18 +186,19 @@ __attribute__((always_inline)) INLINE static void gravity_cache_populate(
     m[i] = gparts[i].mass;
     active[i] = (int)(gparts[i].time_bin <= max_active_bin);
 
-    /* Check whether we can use the multipole instead of P-P */
+    /* Distance to the CoM of the other cell. */
     float dx = x[i] - CoM[0];
     float dy = y[i] - CoM[1];
     float dz = z[i] - CoM[2];
 
-    if(periodic) {
+    if (periodic) {
       dx = nearestf(dx, dim[0]);
       dy = nearestf(dy, dim[1]);
       dz = nearestf(dz, dim[2]);
     }
-      
-    float r2 = dx * dx + dy * dy + dz * dz;
+    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);
   }
 
@@ -292,7 +295,8 @@ gravity_cache_populate_no_mpole(const timebin_t max_active_bin,
  * @param gcount The number of particles to write.
  */
 __attribute__((always_inline)) INLINE void gravity_cache_write_back(
-    const struct gravity_cache *c, struct gpart *restrict gparts, int gcount) {
+    const struct gravity_cache *c, struct gpart *restrict gparts,
+    const int gcount) {
 
   /* Make the compiler understand we are in happy vectorization land */
   swift_declare_aligned_ptr(float, a_x, c->a_x, SWIFT_CACHE_ALIGNMENT);

src/gravity_derivatives.h

@@ -136,9 +136,11 @@ struct potential_derivatives_M2P {
  * @param pot (return) The structure containing all the derivatives.
  */
 __attribute__((always_inline)) INLINE static void
-compute_potential_derivatives_M2L(const float r_x, const float r_y, const float r_z, const float r2,
-                                  const float r_inv, const float eps, const float eps_inv,
-                                  const int periodic, const float rs_inv,
+compute_potential_derivatives_M2L(const float r_x, const float r_y,
+                                  const float r_z, const float r2,
+                                  const float r_inv, const float eps,
+                                  const float eps_inv, const int periodic,
+                                  const float rs_inv,
                                   struct potential_derivatives_M2L *pot) {
 
   float Dt_1;
@@ -390,9 +392,11 @@ compute_potential_derivatives_M2L(const float r_x, const float r_y, const float
  * @param pot (return) The structure containing all the derivatives.
  */
 __attribute__((always_inline)) INLINE static void
-compute_potential_derivatives_M2P(const float r_x, const float r_y, const float r_z, const float r2,
-                                  const float r_inv, const float eps, const float eps_inv,
-				  const int periodic, const float rs_inv,
+compute_potential_derivatives_M2P(const float r_x, const float r_y,
+                                  const float r_z, const float r2,
+                                  const float r_inv, const float eps,
+                                  const float eps_inv, const int periodic,
+                                  const float rs_inv,
                                   struct potential_derivatives_M2P *pot) {
 
   float Dt_1;
@@ -400,7 +404,7 @@ compute_potential_derivatives_M2P(const float r_x, const float r_y, const float
   float Dt_5;
   float Dt_7;
   float Dt_9;
-  
+
   /* Un-softened un-truncated case (Newtonian potential) */
   if (!periodic && r2 > eps * eps) {
 
@@ -429,8 +433,11 @@ compute_potential_derivatives_M2P(const float r_x, const float r_y, const float
     Dt_1 = d.chi_0 * r_inv;
     Dt_3 = (r * d.chi_1 - d.chi_0) * r_inv3;
     Dt_5 = (r * r * d.chi_2 - 3.f * r * d.chi_1 + 3.f * d.chi_0) * r_inv5;
-    Dt_7 = (r * r * r * d.chi_3 - 6.f * r * r * d.chi_2 + 15.f * r * d.chi_1 - 15.f * d.chi_0) * r_inv7;
-    Dt_9 = (r * r * r * r * d.chi_4 - 10.f * r * r * r * d.chi_3 +  45.f * r * r * d.chi_2 - 105.f * r * d.chi_1 + 105.f * d.chi_0) *
+    Dt_7 = (r * r * r * d.chi_3 - 6.f * r * r * d.chi_2 + 15.f * r * d.chi_1 -
+            15.f * d.chi_0) *
+           r_inv7;
+    Dt_9 = (r * r * r * r * d.chi_4 - 10.f * r * r * r * d.chi_3 +
+            45.f * r * r * d.chi_2 - 105.f * r * d.chi_1 + 105.f * d.chi_0) *
            r_inv9;
 
     /* Softened case */

src/multipole.h

@@ -1516,6 +1516,7 @@ INLINE static void gravity_M2M(struct multipole *m_a,
  * @param props The #gravity_props of this calculation.
  * @param periodic Is the calculation periodic ?
  * @param dim The size of the simulation box.
+ * @param rs_inv The inverse of the gravity mesh-smoothing scale.
  */
 INLINE static void gravity_M2L(struct grav_tensor *l_b,
                                const struct multipole *m_a,

src/runner_doiact_grav.h

@@ -403,7 +403,7 @@ static INLINE void runner_dopair_grav_pp_truncated(
       /* Interact! */
       float f_ij, pot_ij;
       runner_iact_grav_pp_truncated(r2, h2_i, h_inv_i, h_inv3_i, mass_j,
-				    rlr_inv, &f_ij, &pot_ij);
+                                    rlr_inv, &f_ij, &pot_ij);
 
       /* if (id_i == 1000 && id_j == 900 && pjd < gcount_j) { */
       /*   message("--- Interacting part"); */
@@ -504,10 +504,10 @@ static INLINE void runner_dopair_grav_pm_full(
     const float r_max2 = r_max_j * r_max_j;
     const float theta_crit2 = e->gravity_properties->theta_crit2;
 
-    if(!gravity_M2P_accept(r_max2, theta_crit2, r2))
+    if (!gravity_M2P_accept(r_max2, theta_crit2, r2))
       error("use_mpole[i] set when M2P accept fails");
 #endif
-    
+
     /* Interact! */
     float f_x, f_y, f_z, pot_ij;
     runner_iact_grav_pm_full(dx, dy, dz, r2, h_i, h_inv_i, multi_j, &f_x, &f_y,
@@ -594,10 +594,10 @@ static INLINE void runner_dopair_grav_pm_truncated(
     const float r_max2 = r_max_j * r_max_j;
     const float theta_crit2 = e->gravity_properties->theta_crit2;
 
-    if(!gravity_M2P_accept(r_max2, theta_crit2, r2))
+    if (!gravity_M2P_accept(r_max2, theta_crit2, r2))
       error("use_mpole[i] set when M2P accept fails");
 #endif
-    
+
     /* Interact! */
     float f_x, f_y, f_z, pot_ij;
     runner_iact_grav_pm_truncated(dx, dy, dz, r2, h_i, h_inv_i, rlr_inv,
@@ -713,12 +713,12 @@ static INLINE void runner_dopair_grav_pp(struct runner *r, struct cell *ci,
 
   /* Fill the caches */
   const float dim[3] = {e->s->dim[0], e->s->dim[1], e->s->dim[2]};
-  gravity_cache_populate(e->max_active_bin, dim, s->periodic, 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, dim, s->periodic, 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, dim, s->periodic, 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, dim, s->periodic, 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) {
@@ -1061,7 +1061,7 @@ static INLINE void runner_doself_grav_pp_truncated(struct runner *r,
       /* Interact! */
       float f_ij, pot_ij;
       runner_iact_grav_pp_truncated(r2, h2_i, h_inv_i, h_inv3_i, mass_j,
-				    rlr_inv, &f_ij, &pot_ij);
+                                    rlr_inv, &f_ij, &pot_ij);
 
       /* if (id_i == 1000 && id_j == 901) { */
       /*   message("--- Interacting part"); */