From c3b6422320f5aff96d53ffc99661cd2bb1ff2c65 Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <schaller@strw.leidenuniv.nl>
Date: Fri, 31 Aug 2018 17:17:11 +0200
Subject: [PATCH] First version of a symmetric M2L kernel. Benchmark with vTune
needed now.
---
src/gravity/Default/gravity_iact.h | 4 +-
src/gravity/Potential/gravity_iact.h | 4 +-
src/gravity_derivatives.h | 35 ++++++-
src/multipole.h | 141 ++++++++++++++++++---------
src/runner_doiact_grav.h | 76 +++++++++++++--
5 files changed, 203 insertions(+), 57 deletions(-)
diff --git a/src/gravity/Default/gravity_iact.h b/src/gravity/Default/gravity_iact.h
index 71e5007a49..6fce3ddd51 100644
--- a/src/gravity/Default/gravity_iact.h
+++ b/src/gravity/Default/gravity_iact.h
@@ -166,7 +166,7 @@ __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, 0.f,
+ potential_derivatives_compute_M2P(r_x, r_y, r_z, r2, r_inv, h, h_inv, 0, 0.f,
&d);
/* 0th order contributions */
@@ -271,7 +271,7 @@ __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,
+ potential_derivatives_compute_M2P(r_x, r_y, r_z, r2, r_inv, h, h_inv, 1,
r_s_inv, &d);
/* 0th order contributions */
diff --git a/src/gravity/Potential/gravity_iact.h b/src/gravity/Potential/gravity_iact.h
index fdc8c17da1..f2094f6ecd 100644
--- a/src/gravity/Potential/gravity_iact.h
+++ b/src/gravity/Potential/gravity_iact.h
@@ -169,7 +169,7 @@ __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, 0.f,
+ potential_derivatives_compute_M2P(r_x, r_y, r_z, r2, r_inv, h, h_inv, 0, 0.f,
&d);
/* 0th order contributions */
@@ -281,7 +281,7 @@ __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,
+ potential_derivatives_compute_M2P(r_x, r_y, r_z, r2, r_inv, h, h_inv, 1,
r_s_inv, &d);
/* 0th order contributions */
diff --git a/src/gravity_derivatives.h b/src/gravity_derivatives.h
index 756fb7af66..2516d29d17 100644
--- a/src/gravity_derivatives.h
+++ b/src/gravity_derivatives.h
@@ -125,6 +125,37 @@ struct potential_derivatives_M2P {
#endif
};
+__attribute__((always_inline)) INLINE static void
+potential_derivatives_flip_signs(struct potential_derivatives_M2L *pot) {
+
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
+ /* 1st order terms */
+ pot->D_100 = -pot->D_100;
+ pot->D_010 = -pot->D_010;
+ pot->D_001 = -pot->D_001;
+#endif
+
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
+ /* 3rd order terms */
+ pot->D_300 = -pot->D_300;
+ pot->D_030 = -pot->D_030;
+ pot->D_003 = -pot->D_003;
+ pot->D_210 = -pot->D_210;
+ pot->D_201 = -pot->D_201;
+ pot->D_021 = -pot->D_021;
+ pot->D_120 = -pot->D_120;
+ pot->D_012 = -pot->D_012;
+ pot->D_102 = -pot->D_102;
+ pot->D_111 = -pot->D_111;
+#endif
+
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 4
+ //MATTHIEU
+ error("oo");
+#endif
+
+}
+
/**
* @brief Compute all the relevent derivatives of the softened and truncated
* gravitational potential for the M2L kernel.
@@ -141,7 +172,7 @@ 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,
+potential_derivatives_compute_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,
@@ -397,7 +428,7 @@ compute_potential_derivatives_M2L(const float r_x, const float r_y,
* @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,
+potential_derivatives_compute_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,
diff --git a/src/multipole.h b/src/multipole.h
index c3d5d84252..994131dbe8 100644
--- a/src/multipole.h
+++ b/src/multipole.h
@@ -1551,51 +1551,11 @@ INLINE static void gravity_M2M(struct multipole *m_a,
#endif
}
-/**
- * @brief Compute the field tensors due to a multipole.
- *
- * Corresponds to equation (28b).
- *
- * @param l_b The field tensor to compute.
- * @param m_a The multipole creating the field.
- * @param pos_b The position of the field tensor.
- * @param pos_a The position of the multipole.
- * @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,
- const double pos_b[3], const double pos_a[3],
- const struct gravity_props *props, int periodic,
- const double dim[3], float rs_inv) {
-
- /* Recover some constants */
- const float eps = props->epsilon_cur;
- const float eps_inv = props->epsilon_cur_inv;
-
- /* Compute distance vector */
- float dx = (float)(pos_b[0] - pos_a[0]);
- float dy = (float)(pos_b[1] - pos_a[1]);
- float dz = (float)(pos_b[2] - pos_a[2]);
-
- /* Apply BC */
- if (periodic) {
- dx = nearest(dx, dim[0]);
- dy = nearest(dy, dim[1]);
- dz = nearest(dz, dim[2]);
- }
-
- /* Compute distance */
- const float r2 = dx * dx + dy * dy + dz * dz;
- const float r_inv = 1. / sqrtf(r2);
-
- /* Compute all derivatives */
- struct potential_derivatives_M2L pot;
- compute_potential_derivatives_M2L(dx, dy, dz, r2, r_inv, eps, eps_inv,
- periodic, rs_inv, &pot);
+INLINE static void gravity_M2L_apply(struct grav_tensor *l_b,
+ const struct multipole *m_a,
+ struct potential_derivatives_M2L pot){
+
#ifdef SWIFT_DEBUG_CHECKS
/* Count interactions */
l_b->num_interacted += m_a->num_gpart;
@@ -1926,6 +1886,99 @@ INLINE static void gravity_M2L(struct grav_tensor *l_b,
#endif
}
+/**
+ * @brief Compute the field tensors due to a multipole.
+ *
+ * Corresponds to equation (28b).
+ *
+ * @param l_b The field tensor to compute.
+ * @param m_a The multipole creating the field.
+ * @param pos_b The position of the field tensor.
+ * @param pos_a The position of the multipole.
+ * @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_nonsym(struct grav_tensor *l_b,
+ const struct multipole *m_a,
+ const double pos_b[3], const double pos_a[3],
+ const struct gravity_props *props, int periodic,
+ const double dim[3], float rs_inv) {
+
+ /* Recover some constants */
+ const float eps = props->epsilon_cur;
+ const float eps_inv = props->epsilon_cur_inv;
+
+ /* Compute distance vector */
+ float dx = (float)(pos_b[0] - pos_a[0]);
+ float dy = (float)(pos_b[1] - pos_a[1]);
+ float dz = (float)(pos_b[2] - pos_a[2]);
+
+ /* Apply BC */
+ if (periodic) {
+ dx = nearest(dx, dim[0]);
+ dy = nearest(dy, dim[1]);
+ dz = nearest(dz, dim[2]);
+ }
+
+ /* Compute distance */
+ const float r2 = dx * dx + dy * dy + dz * dz;
+ const float r_inv = 1. / sqrtf(r2);
+
+ /* Compute all derivatives */
+ struct potential_derivatives_M2L pot;
+ potential_derivatives_compute_M2L(dx, dy, dz, r2, r_inv, eps, eps_inv,
+ periodic, rs_inv, &pot);
+
+ /* Do the M2L tensor multiplication */
+ gravity_M2L_apply(l_b, m_a, pot);
+}
+
+
+INLINE static void gravity_M2L_symmetric(struct grav_tensor *l_a,
+ struct grav_tensor *l_b,
+ const struct multipole *m_a,
+ const struct multipole *m_b,
+ const double pos_b[3], const double pos_a[3],
+ const struct gravity_props *props, int periodic,
+ const double dim[3], float rs_inv) {
+
+ /* Recover some constants */
+ const float eps = props->epsilon_cur;
+ const float eps_inv = props->epsilon_cur_inv;
+
+ /* Compute distance vector */
+ float dx = (float)(pos_b[0] - pos_a[0]);
+ float dy = (float)(pos_b[1] - pos_a[1]);
+ float dz = (float)(pos_b[2] - pos_a[2]);
+
+ /* Apply BC */
+ if (periodic) {
+ dx = nearest(dx, dim[0]);
+ dy = nearest(dy, dim[1]);
+ dz = nearest(dz, dim[2]);
+ }
+
+ /* Compute distance */
+ const float r2 = dx * dx + dy * dy + dz * dz;
+ const float r_inv = 1. / sqrtf(r2);
+
+ /* Compute all derivatives */
+ struct potential_derivatives_M2L pot;
+ potential_derivatives_compute_M2L(dx, dy, dz, r2, r_inv, eps, eps_inv,
+ periodic, rs_inv, &pot);
+
+ /* Do the first M2L tensor multiplication */
+ gravity_M2L_apply(l_b, m_a, pot);
+
+ /* Flip the signs */
+ potential_derivatives_flip_signs(&pot);
+
+ /* Do the first M2L tensor multiplication */
+ gravity_M2L_apply(l_a, m_b, pot);
+}
+
/**
* @brief Creates a copy of #grav_tensor shifted to a new location.
*
diff --git a/src/runner_doiact_grav.h b/src/runner_doiact_grav.h
index a4bc0509ce..bc45000965 100644
--- a/src/runner_doiact_grav.h
+++ b/src/runner_doiact_grav.h
@@ -642,8 +642,8 @@ static INLINE void runner_dopair_grav_pp(struct runner *r, struct cell *ci,
TIMER_TIC;
/* Record activity status */
- const int ci_active = cell_is_active_gravity(ci, e);
- const int cj_active = cell_is_active_gravity(cj, e);
+ const int ci_active = cell_is_active_gravity(ci, e) && (ci->nodeID == e->nodeID);
+ const int cj_active = cell_is_active_gravity(cj, e) && (cj->nodeID == e->nodeID);
/* Anything to do here? */
if (!ci_active && !cj_active) return;
@@ -1122,9 +1122,71 @@ static INLINE void runner_doself_grav_pp(struct runner *r, struct cell *c) {
TIMER_TOC(timer_doself_grav_pp);
}
+/**
+ * @brief Computes the interaction of the field tensor and multipole
+ * of two cells symmetrically.
+ *
+ * @param r The #runner.
+ * @param ci The first #cell.
+ * @param cj The second #cell.
+ */
static INLINE void runner_dopair_grav_mm_symmetric(struct runner *r,
struct cell *restrict ci,
struct cell *restrict cj) {
+
+ /* Some constants */
+ const struct engine *e = r->e;
+ const struct gravity_props *props = e->gravity_properties;
+ const int periodic = e->mesh->periodic;
+ const double dim[3] = {e->mesh->dim[0], e->mesh->dim[1], e->mesh->dim[2]};
+ const float r_s_inv = e->mesh->r_s_inv;
+
+ TIMER_TIC;
+
+ /* Anything to do here? */
+ if ((!cell_is_active_gravity_mm(ci, e) || ci->nodeID != engine_rank) ||
+ (!cell_is_active_gravity_mm(cj, e) || cj->nodeID != engine_rank))
+ error("Invalid state in symmetric M-M calculation!");
+
+ /* Short-cut to the multipole */
+ const struct multipole *multi_i = &ci->multipole->m_pole;
+ const struct multipole *multi_j = &cj->multipole->m_pole;
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (ci == cj) error("Interacting a cell with itself using M2L");
+
+ if (multi_i->num_gpart == 0)
+ error("Multipole i does not seem to have been set.");
+
+ if (multi_j->num_gpart == 0)
+ error("Multipole j does not seem to have been set.");
+
+ if (ci->multipole->pot.ti_init != e->ti_current)
+ error("ci->grav tensor not initialised.");
+
+ if (ci->multipole->pot.ti_init != e->ti_current)
+ error("cj->grav tensor not initialised.");
+
+ if (ci->ti_old_multipole != e->ti_current)
+ error(
+ "Undrifted multipole ci->ti_old_multipole=%lld ci->nodeID=%d "
+ "cj->nodeID=%d e->ti_current=%lld",
+ ci->ti_old_multipole, ci->nodeID, cj->nodeID, e->ti_current);
+
+ if (cj->ti_old_multipole != e->ti_current)
+ error(
+ "Undrifted multipole cj->ti_old_multipole=%lld cj->nodeID=%d "
+ "ci->nodeID=%d e->ti_current=%lld",
+ cj->ti_old_multipole, cj->nodeID, ci->nodeID, e->ti_current);
+#endif
+
+
+ /* Let's interact at this level */
+ gravity_M2L_symmetric(&ci->multipole->pot, &cj->multipole->pot, multi_i, multi_j, ci->multipole->CoM,
+ cj->multipole->CoM, props, periodic, dim, r_s_inv);
+
+
+ TIMER_TOC(timer_dopair_grav_mm);
}
/**
@@ -1171,8 +1233,8 @@ static INLINE void runner_dopair_grav_mm_nonsym(struct runner *r,
#endif
/* Let's interact at this level */
- gravity_M2L(&ci->multipole->pot, multi_j, ci->multipole->CoM,
- cj->multipole->CoM, props, periodic, dim, r_s_inv);
+ gravity_M2L_nonsym(&ci->multipole->pot, multi_j, ci->multipole->CoM,
+ cj->multipole->CoM, props, periodic, dim, r_s_inv);
TIMER_TOC(timer_dopair_grav_mm);
}
@@ -1185,8 +1247,8 @@ static INLINE void runner_dopair_grav_mm_nonsym(struct runner *r,
* @param cj The second #cell.
*/
static INLINE void runner_dopair_grav_mm(struct runner *r,
- struct cell *ci,
- struct cell *cj) {
+ struct cell *ci,
+ struct cell *cj) {
const struct engine *e = r->e;
@@ -1593,7 +1655,7 @@ static INLINE void runner_do_grav_long_range(struct runner *r, struct cell *ci,
theta_crit2, r2_rebuild)) {
/* Call the PM interaction fucntion on the active sub-cells of ci */
- runner_dopair_grav_mm(r, ci, cj);
+ runner_dopair_grav_mm_nonsym(r, ci, cj);
// runner_dopair_recursive_grav_pm(r, ci, cj);
/* Record that this multipole received a contribution */
--
GitLab