diff --git a/src/multipole.h b/src/multipole.h
index 2e828baf6bb8ae05ab58ae03324cd419a968d7ad..8a35557264dc7ad2995cfb95376e87b5fbd24cc2 100644
--- a/src/multipole.h
+++ b/src/multipole.h
@@ -906,55 +906,45 @@ INLINE static void gravity_M2M(struct multipole *m_a,
m_a->M_001 = m_b->M_001 + X_001(dx) * m_b->M_000;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
- m_a->M_200 = m_b->M_200 + X_100(dx) * m_b->M_100 + X_200(dx) * m_b->M_000;
+ /* Shift 2nd order term */
+ m_a->M_200 = m_b->M_200 + X_100(dx) * m_b->M_100 + X_200(dx) * m_b->M_000;
m_a->M_020 = m_b->M_020 + X_010(dx) * m_b->M_010 + X_020(dx) * m_b->M_000;
-
m_a->M_002 = m_b->M_002 + X_001(dx) * m_b->M_001 + X_002(dx) * m_b->M_000;
-
m_a->M_110 = m_b->M_110 + X_100(dx) * m_b->M_010 + X_010(dx) * m_b->M_100 +
X_110(dx) * m_b->M_000;
-
m_a->M_101 = m_b->M_101 + X_100(dx) * m_b->M_001 + X_001(dx) * m_b->M_100 +
X_101(dx) * m_b->M_000;
-
m_a->M_011 = m_b->M_011 + X_010(dx) * m_b->M_001 + X_001(dx) * m_b->M_010 +
X_011(dx) * m_b->M_000;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
+
+ /* Shift 3rd order term */
m_a->M_300 = m_b->M_300 + X_100(dx) * m_b->M_200 + X_200(dx) * m_b->M_100 +
X_300(dx) * m_b->M_000;
-
m_a->M_030 = m_b->M_030 + X_010(dx) * m_b->M_020 + X_020(dx) * m_b->M_010 +
X_030(dx) * m_b->M_000;
-
m_a->M_003 = m_b->M_003 + X_001(dx) * m_b->M_002 + X_002(dx) * m_b->M_001 +
X_003(dx) * m_b->M_000;
-
m_a->M_210 = m_b->M_210 + X_100(dx) * m_b->M_110 + X_010(dx) * m_b->M_200 +
X_200(dx) * m_b->M_010 + X_110(dx) * m_b->M_100 +
X_210(dx) * m_b->M_000;
-
m_a->M_201 = m_b->M_201 + X_100(dx) * m_b->M_101 + X_001(dx) * m_b->M_200 +
X_200(dx) * m_b->M_001 + X_101(dx) * m_b->M_100 +
X_201(dx) * m_b->M_000;
-
m_a->M_120 = m_b->M_120 + X_010(dx) * m_b->M_110 + X_100(dx) * m_b->M_020 +
X_020(dx) * m_b->M_100 + X_110(dx) * m_b->M_010 +
X_120(dx) * m_b->M_000;
-
m_a->M_021 = m_b->M_021 + X_010(dx) * m_b->M_011 + X_001(dx) * m_b->M_020 +
X_020(dx) * m_b->M_001 + X_011(dx) * m_b->M_010 +
X_021(dx) * m_b->M_000;
-
m_a->M_102 = m_b->M_102 + X_001(dx) * m_b->M_101 + X_100(dx) * m_b->M_002 +
X_002(dx) * m_b->M_100 + X_101(dx) * m_b->M_001 +
X_102(dx) * m_b->M_000;
-
m_a->M_012 = m_b->M_012 + X_001(dx) * m_b->M_011 + X_010(dx) * m_b->M_002 +
X_002(dx) * m_b->M_010 + X_011(dx) * m_b->M_001 +
X_012(dx) * m_b->M_000;
-
m_a->M_111 = m_b->M_111 + X_100(dx) * m_b->M_011 + X_010(dx) * m_b->M_101 +
X_001(dx) * m_b->M_110 + X_110(dx) * m_b->M_001 +
X_101(dx) * m_b->M_010 + X_011(dx) * m_b->M_100 +
@@ -1013,7 +1003,6 @@ INLINE static void gravity_M2L(struct grav_tensor *l_b,
l_b->F_010 += m_a->M_000 * D_010(dx, dy, dz, r_inv);
l_b->F_001 += m_a->M_000 * D_001(dx, dy, dz, r_inv);
#endif
-
#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
/* 2nd order multipole term (addition to rank 0)*/
@@ -1043,7 +1032,6 @@ INLINE static void gravity_M2L(struct grav_tensor *l_b,
l_b->F_101 += m_a->M_000 * D_101(dx, dy, dz, r_inv);
l_b->F_011 += m_a->M_000 * D_011(dx, dy, dz, r_inv);
#endif
-
#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
/* 3rd order multipole term (addition to rank 0)*/
@@ -1156,7 +1144,6 @@ INLINE static void gravity_L2L(struct grav_tensor *la,
la->F_010 += X_000(dx) * lb->F_010;
la->F_001 += X_000(dx) * lb->F_001;
#endif
-
#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
/* Shift 2nd order multipole term (addition to rank 0)*/
@@ -1181,7 +1168,6 @@ INLINE static void gravity_L2L(struct grav_tensor *la,
la->F_101 += X_000(dx) * lb->F_101;
la->F_011 += X_000(dx) * lb->F_011;
#endif
-
#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
/* Shift 3rd order multipole term (addition to rank 0)*/
@@ -1271,6 +1257,7 @@ INLINE static void gravity_L2P(const struct grav_tensor *lb,
gp->a_grav[2] += X_000(dx) * lb->F_001;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
+
/* 1st order interaction */
gp->a_grav[0] +=
X_100(dx) * lb->F_200 + X_010(dx) * lb->F_110 + X_001(dx) * lb->F_101;
@@ -1280,6 +1267,7 @@ INLINE static void gravity_L2P(const struct grav_tensor *lb,
X_100(dx) * lb->F_101 + X_010(dx) * lb->F_011 + X_001(dx) * lb->F_002;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
+
/* 2nd order interaction */
gp->a_grav[0] +=
X_200(dx) * lb->F_300 + X_020(dx) * lb->F_120 + X_002(dx) * lb->F_102;
diff --git a/src/runner_doiact_grav.h b/src/runner_doiact_grav.h
index 4ff2abcfb095965518c185446b2703816dccbbc2..c3a7f45b5ac4c3bbfeed04db6203b120f2e96e41 100644
--- a/src/runner_doiact_grav.h
+++ b/src/runner_doiact_grav.h
@@ -25,6 +25,14 @@
#include "gravity.h"
#include "part.h"
+/**
+ * @brief Recursively propagate the multipoles down the tree by applying the
+ * L2L and L2P kernels.
+ *
+ * @param r The #runner.
+ * @param c The #cell we are working on.
+ * @param timer Are we timing this ?
+ */
void runner_do_grav_down(struct runner *r, struct cell *c, int timer) {
const struct engine *e = r->e;
@@ -114,62 +122,6 @@ void runner_dopair_grav_pm(const struct runner *r,
const struct cell *restrict cj) {
error("Function should not be called");
-
- /* const struct engine *e = r->e; */
- /* const int gcount = ci->gcount; */
- /* struct gpart *restrict gparts = ci->gparts; */
- /* const struct gravity_tensors *multi = cj->multipole; */
- /* const float a_smooth = e->gravity_properties->a_smooth; */
- /* const float rlr_inv = 1. / (a_smooth * ci->super->width[0]); */
-
- /* TIMER_TIC; */
-
- /* #ifdef SWIFT_DEBUG_CHECKS */
- /* if (gcount == 0) error("Empty cell!"); */
-
- /* if (multi->m_pole.mass == 0.0) */
- /* error("Multipole does not seem to have been set."); */
- /* #endif */
-
- /* /\* Anything to do here? *\/ */
- /* if (!cell_is_active(ci, e)) return; */
-
- /* #if ICHECK > 0 */
- /* for (int pid = 0; pid < gcount; pid++) { */
-
- /* /\* Get a hold of the ith part in ci. *\/ */
- /* struct gpart *restrict gp = &gparts[pid]; */
-
- /* if (gp->id_or_neg_offset == ICHECK) */
- /* message("id=%lld loc=[ %f %f %f ] size= %f count= %d", */
- /* gp->id_or_neg_offset, cj->loc[0], cj->loc[1], cj->loc[2], */
- /* cj->width[0], cj->gcount); */
- /* } */
- /* #endif */
-
- /* /\* Loop over every particle in leaf. *\/ */
- /* for (int pid = 0; pid < gcount; pid++) { */
-
- /* /\* Get a hold of the ith part in ci. *\/ */
- /* struct gpart *restrict gp = &gparts[pid]; */
-
- /* if (!gpart_is_active(gp, e)) continue; */
-
- /* /\* Compute the pairwise distance. *\/ */
- /* const float dx[3] = {multi->CoM[0] - gp->x[0], // x */
- /* multi->CoM[1] - gp->x[1], // y */
- /* multi->CoM[2] - gp->x[2]}; // z */
- /* const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2]; */
-
- /* /\* Interact !*\/ */
- /* runner_iact_grav_pm(rlr_inv, r2, dx, gp, &multi->m_pole); */
-
- /* #ifdef SWIFT_DEBUG_CHECKS */
- /* gp->mass_interacted += multi->m_pole.mass; */
- /* #endif */
- /* } */
-
- /* TIMER_TOC(timer_dopair_grav_pm); */
}
/**
@@ -307,8 +259,7 @@ void runner_doself_grav_pp(struct runner *r, struct cell *c) {
TIMER_TIC;
#ifdef SWIFT_DEBUG_CHECKS
- if (c->gcount == 0) // MATTHIEU sanity check
- error("Empty cell !");
+ if (c->gcount == 0) error("Doing self gravity on an empty cell !");
#endif
/* Anything to do here? */
@@ -403,7 +354,8 @@ void runner_dopair_grav(struct runner *r, struct cell *ci, struct cell *cj,
const int gcount_j = cj->gcount;
/* Early abort? */
- if (gcount_i == 0 || gcount_j == 0) error("Empty cell !");
+ if (gcount_i == 0 || gcount_j == 0)
+ error("Doing pair gravity on an empty cell !");
/* Bad stuff will happen if cell sizes are different */
if (ci->width[0] != cj->width[0])
@@ -500,7 +452,7 @@ void runner_doself_grav(struct runner *r, struct cell *c, int gettimer) {
#ifdef SWIFT_DEBUG_CHECKS
/* Early abort? */
- if (c->gcount == 0) error("Empty cell !");
+ if (c->gcount == 0) error("Doing self gravity on an empty cell !");
#endif
TIMER_TIC;
@@ -594,8 +546,7 @@ void runner_do_grav_long_range(struct runner *r, struct cell *ci, int timer) {
/* cj->loc[1] - pos_i[1], // y */
/* cj->loc[2] - pos_i[2]}; // z */
/* const double r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2]; */
-
- // if (r2 > max_d2) continue;
+ /* if (r2 > max_d2) continue; */
if (!cell_are_neighbours(ci, cj)) runner_dopair_grav_mm(r, ci, cj);
}