diff --git a/src/gravity_derivatives.h b/src/gravity_derivatives.h
index 2516d29d179ae2273c8737b12882cd8f56622e4b..e6ab3729e979caaf3b6ec05dec92411636bff35f 100644
--- a/src/gravity_derivatives.h
+++ b/src/gravity_derivatives.h
@@ -125,6 +125,15 @@ struct potential_derivatives_M2P {
#endif
};
+/**
+ * @brief Converts the derivatives from a distance vector to its opposite.
+ *
+ * From a series of tensors D_xxx(\vec{r}), compute D_xxx(-\vec{r}).
+ * This can be computed efficiently by flipping the sign of all the odd
+ * derivative terms.
+ *
+ * @param pot The derivatives of the potential.
+ */
__attribute__((always_inline)) INLINE static void
potential_derivatives_flip_signs(struct potential_derivatives_M2L *pot) {
@@ -150,10 +159,29 @@ potential_derivatives_flip_signs(struct potential_derivatives_M2L *pot) {
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 4
- //MATTHIEU
- error("oo");
+ /* 5th order terms */
+ pot->D_500 = -pot->D_500;
+ pot->D_050 = -pot->D_050;
+ pot->D_005 = -pot->D_005;
+ pot->D_410 = -pot->D_410;
+ pot->D_401 = -pot->D_401;
+ pot->D_041 = -pot->D_041;
+ pot->D_140 = -pot->D_140;
+ pot->D_014 = -pot->D_014;
+ pot->D_104 = -pot->D_104;
+ pot->D_320 = -pot->D_320;
+ pot->D_302 = -pot->D_302;
+ pot->D_032 = -pot->D_032;
+ pot->D_230 = -pot->D_230;
+ pot->D_023 = -pot->D_023;
+ pot->D_203 = -pot->D_203;
+ pot->D_311 = -pot->D_311;
+ pot->D_131 = -pot->D_131;
+ pot->D_113 = -pot->D_113;
+ pot->D_122 = -pot->D_122;
+ pot->D_212 = -pot->D_212;
+ pot->D_221 = -pot->D_221;
#endif
-
}
/**
diff --git a/src/multipole.h b/src/multipole.h
index 994131dbe85151a4ac2ab1f3672ec9c2bd9eaf72..681b9dfaed14d9cedef2af39ce53212fb7d76f80 100644
--- a/src/multipole.h
+++ b/src/multipole.h
@@ -1551,11 +1551,19 @@ 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 pot The derivatives of the potential.
+ */
INLINE static void gravity_M2L_apply(struct grav_tensor *l_b,
- const struct multipole *m_a,
- struct potential_derivatives_M2L pot){
+ const struct multipole *m_a,
+ struct potential_derivatives_M2L pot) {
-
#ifdef SWIFT_DEBUG_CHECKS
/* Count interactions */
l_b->num_interacted += m_a->num_gpart;
@@ -1887,12 +1895,10 @@ INLINE static void gravity_M2L_apply(struct grav_tensor *l_b,
}
/**
- * @brief Compute the field tensors due to a multipole.
- *
- * Corresponds to equation (28b).
+ * @brief Compute the field tensor due to a multipole.
*
* @param l_b The field tensor to compute.
- * @param m_a The multipole creating the field.
+ * @param m_a The multipole.
* @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.
@@ -1900,11 +1906,10 @@ INLINE static void gravity_M2L_apply(struct grav_tensor *l_b,
* @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) {
+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;
@@ -1935,14 +1940,27 @@ INLINE static void gravity_M2L_nonsym(struct grav_tensor *l_b,
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) {
+/**
+ * @brief Compute the field tensor due to a multipole and the symmetric
+ * equivalent.
+ *
+ * @param l_a The first field tensor to compute.
+ * @param l_b The second field tensor to compute.
+ * @param m_a The first multipole.
+ * @param m_b The second multipole.
+ * @param pos_b The position of the first m-pole and field tensor.
+ * @param pos_a The position of the second m-pole and field tensor.
+ * @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_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;
@@ -1972,10 +1990,10 @@ INLINE static void gravity_M2L_symmetric(struct grav_tensor *l_a,
/* Do the first M2L tensor multiplication */
gravity_M2L_apply(l_b, m_a, pot);
- /* Flip the signs */
+ /* Flip the signs of odd derivatives */
potential_derivatives_flip_signs(&pot);
-
- /* Do the first M2L tensor multiplication */
+
+ /* Do the second M2L tensor multiplication */
gravity_M2L_apply(l_a, m_b, pot);
}
diff --git a/src/runner_doiact_grav.h b/src/runner_doiact_grav.h
index bc4500096553676e0bf9d8e7f0c0dcb36c1e522f..29a00c7833cb9d0c39a492e90ea3dc3f816e1c87 100644
--- a/src/runner_doiact_grav.h
+++ b/src/runner_doiact_grav.h
@@ -642,8 +642,10 @@ 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) && (ci->nodeID == e->nodeID);
- const int cj_active = cell_is_active_gravity(cj, e) && (cj->nodeID == e->nodeID);
+ 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;
@@ -1131,8 +1133,8 @@ static INLINE void runner_doself_grav_pp(struct runner *r, struct cell *c) {
* @param cj The second #cell.
*/
static INLINE void runner_dopair_grav_mm_symmetric(struct runner *r,
- struct cell *restrict ci,
- struct cell *restrict cj) {
+ struct cell *restrict ci,
+ struct cell *restrict cj) {
/* Some constants */
const struct engine *e = r->e;
@@ -1160,7 +1162,7 @@ static INLINE void runner_dopair_grav_mm_symmetric(struct runner *r,
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.");
@@ -1180,13 +1182,12 @@ static INLINE void runner_dopair_grav_mm_symmetric(struct runner *r,
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);
-
+ 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);
+ TIMER_TOC(timer_dopair_grav_mm);
}
/**
@@ -1198,8 +1199,8 @@ static INLINE void runner_dopair_grav_mm_symmetric(struct runner *r,
* @param cj The #cell with the multipole.
*/
static INLINE void runner_dopair_grav_mm_nonsym(struct runner *r,
- struct cell *restrict ci,
- struct cell *restrict cj) {
+ struct cell *restrict ci,
+ struct cell *restrict cj) {
/* Some constants */
const struct engine *e = r->e;
@@ -1234,7 +1235,7 @@ static INLINE void runner_dopair_grav_mm_nonsym(struct runner *r,
/* Let's interact at this level */
gravity_M2L_nonsym(&ci->multipole->pot, multi_j, ci->multipole->CoM,
- cj->multipole->CoM, props, periodic, dim, r_s_inv);
+ cj->multipole->CoM, props, periodic, dim, r_s_inv);
TIMER_TOC(timer_dopair_grav_mm);
}
@@ -1246,9 +1247,8 @@ static INLINE void runner_dopair_grav_mm_nonsym(struct runner *r,
* @param ci The first #cell.
* @param cj The second #cell.
*/
-static INLINE void runner_dopair_grav_mm(struct runner *r,
- struct cell *ci,
- struct cell *cj) {
+static INLINE void runner_dopair_grav_mm(struct runner *r, struct cell *ci,
+ struct cell *cj) {
const struct engine *e = r->e;