From 4c992de23f5f77ecdcf8f8386167b782cd6e7bd3 Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <matthieu.schaller@durham.ac.uk>
Date: Mon, 11 Jun 2018 20:59:13 +0200
Subject: [PATCH] Improved documentation and removed debugging code.
---
src/gravity/Default/gravity_iact.h | 78 ++++++++++++++++++------------
src/gravity_cache.h | 22 +++++----
src/gravity_derivatives.h | 25 ++++++----
src/multipole.h | 1 +
src/runner_doiact_grav.h | 24 ++++-----
5 files changed, 90 insertions(+), 60 deletions(-)
diff --git a/src/gravity/Default/gravity_iact.h b/src/gravity/Default/gravity_iact.h
index f4da4bffcd..99d636a307 100644
--- a/src/gravity/Default/gravity_iact.h
+++ b/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;
diff --git a/src/gravity_cache.h b/src/gravity_cache.h
index addb4bb12f..62b8a5f648 100644
--- a/src/gravity_cache.h
+++ b/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);
diff --git a/src/gravity_derivatives.h b/src/gravity_derivatives.h
index 73f7d39c39..6b6806a7da 100644
--- a/src/gravity_derivatives.h
+++ b/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 */
diff --git a/src/multipole.h b/src/multipole.h
index 7831e2532a..3e0533a6a6 100644
--- a/src/multipole.h
+++ b/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,
diff --git a/src/runner_doiact_grav.h b/src/runner_doiact_grav.h
index b9782bd65d..0b6292ebb6 100644
--- a/src/runner_doiact_grav.h
+++ b/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"); */
--
GitLab