diff --git a/src/gravity/Default/gravity_iact.h b/src/gravity/Default/gravity_iact.h
index 6d8aa3261fe30b765a58755856f29c98b0c10fd8..671ca6a343237aab5262338ae522d50cf7072f3f 100644
--- a/src/gravity/Default/gravity_iact.h
+++ b/src/gravity/Default/gravity_iact.h
@@ -23,18 +23,97 @@
/* Includes. */
#include "const.h"
#include "kernel_gravity.h"
+#include "multipole.h"
#include "vector.h"
/**
- * @brief Gravity potential
+ * @brief Gravity forces between particles
*/
-__attribute__((always_inline)) INLINE static void runner_iact_grav(
- float r2, float *dx, struct gpart *pi, struct gpart *pj) {}
+__attribute__((always_inline)) INLINE static void runner_iact_grav_pp(
+ float r2, const float *dx, struct gpart *gpi, struct gpart *gpj) {
+
+ /* Apply the gravitational acceleration. */
+ const float ir = 1.0f / sqrtf(r2);
+ const float w = ir * ir * ir;
+ const float wdx[3] = {w * dx[0], w * dx[1], w * dx[2]};
+ const float mi = gpi->mass;
+ const float mj = gpj->mass;
+
+ gpi->a_grav[0] -= wdx[0] * mj;
+ gpi->a_grav[1] -= wdx[1] * mj;
+ gpi->a_grav[2] -= wdx[2] * mj;
+ gpi->mass_interacted += mj;
+
+ gpj->a_grav[0] += wdx[0] * mi;
+ gpj->a_grav[1] += wdx[1] * mi;
+ gpj->a_grav[2] += wdx[2] * mi;
+ gpj->mass_interacted += mi;
+}
+
+__attribute__((always_inline)) INLINE static void runner_iact_grav_pp_nonsym(
+ float r2, const float *dx, struct gpart *gpi, const struct gpart *gpj) {
+
+ /* Apply the gravitational acceleration. */
+ const float ir = 1.0f / sqrtf(r2);
+ const float w = ir * ir * ir;
+ const float wdx[3] = {w * dx[0], w * dx[1], w * dx[2]};
+ const float mj = gpj->mass;
+
+ gpi->a_grav[0] -= wdx[0] * mj;
+ gpi->a_grav[1] -= wdx[1] * mj;
+ gpi->a_grav[2] -= wdx[2] * mj;
+ gpi->mass_interacted += mj;
+}
/**
- * @brief Gravity potential (Vectorized version)
+ * @brief Gravity forces between particles
*/
-__attribute__((always_inline)) INLINE static void runner_iact_vec_grav(
- float *R2, float *Dx, struct gpart **pi, struct gpart **pj) {}
+__attribute__((always_inline)) INLINE static void runner_iact_grav_pm(
+ float r2, const float *dx, struct gpart *gp,
+ const struct multipole *multi) {
+
+ /* Apply the gravitational acceleration. */
+ const float r = sqrtf(r2);
+ const float ir = 1.f / r;
+ const float mrinv3 = multi->mass * ir * ir * ir;
+
+#if multipole_order < 2
+
+ /* 0th and 1st order terms */
+ gp->a_grav[0] += mrinv3 * dx[0];
+ gp->a_grav[1] += mrinv3 * dx[1];
+ gp->a_grav[2] += mrinv3 * dx[2];
+
+ gp->mass_interacted += multi->mass;
+#elif multipole_order == 2
+ /* Terms up to 2nd order (quadrupole) */
+
+ /* Follows the notation in Bonsai */
+ const float mrinv5 = mrinv3 * ir * ir;
+ const float mrinv7 = mrinv5 * ir * ir;
+
+ const float D1 = -mrinv3;
+ const float D2 = 3.f * mrinv5;
+ const float D3 = -15.f * mrinv7;
+
+ const float q = multi->I_xx + multi->I_yy + multi->I_zz;
+ const float qRx =
+ multi->I_xx * dx[0] + multi->I_xy * dx[1] + multi->I_xz * dx[2];
+ const float qRy =
+ multi->I_xy * dx[0] + multi->I_yy * dx[1] + multi->I_yz * dx[2];
+ const float qRz =
+ multi->I_xz * dx[0] + multi->I_yz * dx[1] + multi->I_zz * dx[2];
+ const float qRR = qRx * dx[0] + qRy * dx[1] + qRz * dx[2];
+ const float C = D1 + 0.5f * D2 * q + 0.5f * D3 * qRR;
+
+ gp->a_grav[0] -= C * dx[0] + D2 * qRx;
+ gp->a_grav[1] -= C * dx[1] + D2 * qRy;
+ gp->a_grav[2] -= C * dx[2] + D2 * qRz;
+
+ gp->mass_interacted += multi->mass;
+#else
+#error "Multipoles of order >2 not yet implemented."
+#endif
+}
#endif /* SWIFT_RUNNER_IACT_GRAV_H */
diff --git a/src/kernel_gravity.h b/src/kernel_gravity.h
index bef650347649915416028156213aa7c90c0eb14f..1e39f96854fd5b3c325edaa2e597a1a6f3b81f60 100644
--- a/src/kernel_gravity.h
+++ b/src/kernel_gravity.h
@@ -31,7 +31,8 @@
r. The resulting value should be post-multiplied with r^-3, resulting
in a polynomial with terms ranging from r^-3 to r^3, which are
sufficient to model both the direct potential as well as the splines
- near the origin. */
+ near the origin.
+ As in the hydro case, the 1/h^3 needs to be multiplied in afterwards */
/* Coefficients for the kernel. */
#define kernel_grav_name "Gadget-2 softening kernel"
diff --git a/src/runner_doiact_grav.h b/src/runner_doiact_grav.h
index 0b17d6e2bd883c51b7048957567cd2b5c7f272f5..9dc8432a42c18dca6c93c5a665630f89c758c431 100644
--- a/src/runner_doiact_grav.h
+++ b/src/runner_doiact_grav.h
@@ -22,7 +22,7 @@
/* Includes. */
#include "cell.h"
-#include "clocks.h"
+#include "gravity.h"
#include "part.h"
#define ICHECK -1
@@ -113,19 +113,20 @@ __attribute__((always_inline)) INLINE static void runner_dopair_grav_pm(
error("Multipole does not seem to have been set.");
#endif
- /* Anything to do here? */
- // if (ci->ti_end_min > ti_current) return;
- /* Loop over every particle in leaf. */
- /* for (int pid = 0; pid < gcount; pid++) { */
+/* Anything to do here? */
+// if (ci->ti_end_min > ti_current) 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]; */
+ /* Get a hold of the ith part in ci. */
+ struct gpart *restrict gp = &gparts[pid];
- /* if (gp->id == -ICHECK) */
- /* message("id=%lld loc=[ %f %f %f ] size= %f count= %d", gp->id,
- * cj->loc[0], */
- /* cj->loc[1], cj->loc[2], cj->h[0], cj->gcount); */
- /* } */
+ if (gp->id == -ICHECK)
+ message("id=%lld loc=[ %f %f %f ] size= %f count= %d", gp->id, cj->loc[0],
+ cj->loc[1], cj->loc[2], cj->h[0], cj->gcount);
+ }
+#endif
/* Loop over every particle in leaf. */
for (int pid = 0; pid < gcount; pid++) {
@@ -143,46 +144,8 @@ __attribute__((always_inline)) INLINE static void runner_dopair_grav_pm(
multi.CoM[2] - gp->x[2]}; // z
const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
- /* Apply the gravitational acceleration. */
- const float ir = 1.0f / sqrtf(r2);
- const float mrinv3 = multi.mass * ir * ir * ir;
-
-#if multipole_order < 2
-
- /* 0th and 1st order terms */
- gp->a_grav[0] += mrinv3 * dx[0];
- gp->a_grav[1] += mrinv3 * dx[1];
- gp->a_grav[2] += mrinv3 * dx[2];
-
-#elif multipole_order == 2
- /* Terms up to 2nd order (quadrupole) */
-
- /* Follows the notation in Bonsai */
- const float mrinv5 = mrinv3 * ir * ir;
- const float mrinv7 = mrinv5 * ir * ir;
-
- const float D1 = -mrinv3;
- const float D2 = 3.f * mrinv5;
- const float D3 = -15.f * mrinv7;
-
- const float q = multi.I_xx + multi.I_yy + multi.I_zz;
- const float qRx =
- multi.I_xx * dx[0] + multi.I_xy * dx[1] + multi.I_xz * dx[2];
- const float qRy =
- multi.I_xy * dx[0] + multi.I_yy * dx[1] + multi.I_yz * dx[2];
- const float qRz =
- multi.I_xz * dx[0] + multi.I_yz * dx[1] + multi.I_zz * dx[2];
- const float qRR = qRx * dx[0] + qRy * dx[1] + qRz * dx[2];
- const float C = D1 + 0.5f * D2 * q + 0.5f * D3 * qRR;
-
- gp->a_grav[0] -= C * dx[0] + D2 * qRx;
- gp->a_grav[1] -= C * dx[1] + D2 * qRy;
- gp->a_grav[2] -= C * dx[2] + D2 * qRz;
-
- gp->mass_interacted += multi.mass;
-#else
-#error "Multipoles of order >2 not yet implemented."
-#endif
+ /* Interact !*/
+ runner_iact_grav_pm(r2, dx, gp, &multi);
}
TIMER_TOC(TIMER_DOPAIR); // MATTHIEU
@@ -215,44 +178,42 @@ __attribute__((always_inline)) INLINE static void runner_dopair_grav_pp(
error("Non matching cell sizes !! h_i=%f h_j=%f", ci->h[0], cj->h[0]);
#endif
- /* Anything to do here? */
- // if (ci->ti_end_min > ti_current && cj->ti_end_min > ti_current) return;
+/* Anything to do here? */
+// if (ci->ti_end_min > ti_current && cj->ti_end_min > ti_current) return;
- /* for (int pid = 0; pid < gcount_i; pid++) { */
+#if ICHECK > 0
+ for (int pid = 0; pid < gcount_i; pid++) {
- /* /\* Get a hold of the ith part in ci. *\/ */
- /* struct gpart *restrict gp = &gparts_i[pid]; */
+ /* Get a hold of the ith part in ci. */
+ struct gpart *restrict gp = &gparts_i[pid];
- /* if (gp->id == -ICHECK) */
- /* message("id=%lld loc=[ %f %f %f ] size= %f count= %d", gp->id,
- * cj->loc[0], */
- /* cj->loc[1], cj->loc[2], cj->h[0], cj->gcount); */
- /* } */
+ if (gp->id == -ICHECK)
+ message("id=%lld loc=[ %f %f %f ] size= %f count= %d", gp->id, cj->loc[0],
+ cj->loc[1], cj->loc[2], cj->h[0], cj->gcount);
+ }
- /* for (int pid = 0; pid < gcount_j; pid++) { */
+ for (int pid = 0; pid < gcount_j; pid++) {
- /* /\* Get a hold of the ith part in ci. *\/ */
- /* struct gpart *restrict gp = &gparts_j[pid]; */
+ /* Get a hold of the ith part in ci. */
+ struct gpart *restrict gp = &gparts_j[pid];
- /* if (gp->id == -ICHECK) */
- /* message("id=%lld loc=[ %f %f %f ] size= %f count=%d", gp->id,
- * ci->loc[0], */
- /* ci->loc[1], ci->loc[2], ci->h[0], ci->gcount); */
- /* } */
+ if (gp->id == -ICHECK)
+ message("id=%lld loc=[ %f %f %f ] size= %f count=%d", gp->id, ci->loc[0],
+ ci->loc[1], ci->loc[2], ci->h[0], ci->gcount);
+ }
+#endif
/* Loop over all particles in ci... */
for (int pid = 0; pid < gcount_i; pid++) {
/* Get a hold of the ith part in ci. */
struct gpart *restrict gpi = &gparts_i[pid];
- const float mi = gpi->mass;
/* Loop over every particle in the other cell. */
for (int pjd = 0; pjd < gcount_j; pjd++) {
/* Get a hold of the jth part in cj. */
struct gpart *restrict gpj = &gparts_j[pjd];
- const float mj = gpj->mass;
/* Compute the pairwise distance. */
const float dx[3] = {gpi->x[0] - gpj->x[0], // x
@@ -260,23 +221,8 @@ __attribute__((always_inline)) INLINE static void runner_dopair_grav_pp(
gpi->x[2] - gpj->x[2]}; // z
const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
- /* Apply the gravitational acceleration. */
- const float ir = 1.0f / sqrtf(r2);
- const float w = ir * ir * ir;
- const float wdx[3] = {w * dx[0], w * dx[1], w * dx[2]};
-
- // if (gpi->ti_end <= ti_current) {
- gpi->a_grav[0] -= wdx[0] * mj;
- gpi->a_grav[1] -= wdx[1] * mj;
- gpi->a_grav[2] -= wdx[2] * mj;
- gpi->mass_interacted += mj;
- //}
- // if (gpj->ti_end <= ti_current) {
- gpj->a_grav[0] += wdx[0] * mi;
- gpj->a_grav[1] += wdx[1] * mi;
- gpj->a_grav[2] += wdx[2] * mi;
- gpj->mass_interacted += mi;
- // }
+ /* Interact ! */
+ runner_iact_grav_pp(r2, dx, gpi, gpj);
}
}
@@ -306,33 +252,32 @@ __attribute__((always_inline)) INLINE static void runner_doself_grav_pp(
error("Empty cell !");
#endif
- /* Anything to do here? */
- // if (c->ti_end_min > ti_current) return;
+/* Anything to do here? */
+// if (c->ti_end_min > ti_current) return;
- /* for (int pid = 0; pid < gcount; pid++) { */
+#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]; */
+ /* Get a hold of the ith part in ci. */
+ struct gpart *restrict gp = &gparts[pid];
- /* if (gp->id == -ICHECK) */
- /* message("id=%lld loc=[ %f %f %f ] size= %f count= %d", gp->id,
- * c->loc[0], */
- /* c->loc[1], c->loc[2], c->h[0], c->gcount); */
- /* } */
+ if (gp->id == -ICHECK)
+ message("id=%lld loc=[ %f %f %f ] size= %f count= %d", gp->id, c->loc[0],
+ c->loc[1], c->loc[2], c->h[0], c->gcount);
+ }
+#endif
/* Loop over all particles in ci... */
for (int pid = 0; pid < gcount; pid++) {
/* Get a hold of the ith part in ci. */
struct gpart *restrict gpi = &gparts[pid];
- const float mi = gpi->mass;
/* Loop over every particle in the other cell. */
for (int pjd = pid + 1; pjd < gcount; pjd++) {
/* Get a hold of the jth part in ci. */
struct gpart *restrict gpj = &gparts[pjd];
- const float mj = gpj->mass;
/* Compute the pairwise distance. */
const float dx[3] = {gpi->x[0] - gpj->x[0], // x
@@ -340,23 +285,8 @@ __attribute__((always_inline)) INLINE static void runner_doself_grav_pp(
gpi->x[2] - gpj->x[2]}; // z
const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
- /* Apply the gravitational acceleration. */
- const float ir = 1.0f / sqrtf(r2);
- const float w = ir * ir * ir;
- const float wdx[3] = {w * dx[0], w * dx[1], w * dx[2]};
-
- // if (gpi->ti_end <= ti_current) {
- gpi->a_grav[0] -= wdx[0] * mj;
- gpi->a_grav[1] -= wdx[1] * mj;
- gpi->a_grav[2] -= wdx[2] * mj;
- gpi->mass_interacted += mj;
- //}
- // if (gpj->ti_end <= ti_current) {
- gpj->a_grav[0] += wdx[0] * mi;
- gpj->a_grav[1] += wdx[1] * mi;
- gpj->a_grav[2] += wdx[2] * mi;
- gpj->mass_interacted += mi;
- //}
+ /* Interact ! */
+ runner_iact_grav_pp(r2, dx, gpi, gpj);
}
}
@@ -404,27 +334,27 @@ static void runner_dopair_grav(struct runner *r, struct cell *ci,
#endif
- /* for (int pid = 0; pid < gcount_i; pid++) { */
+#if ICHECK > 0
+ for (int pid = 0; pid < gcount_i; pid++) {
- /* /\* Get a hold of the ith part in ci. *\/ */
- /* struct gpart *restrict gp = &ci->gparts[pid]; */
+ /* Get a hold of the ith part in ci. */
+ struct gpart *restrict gp = &ci->gparts[pid];
- /* if (gp->id == -ICHECK) */
- /* message("id=%lld loc=[ %f %f %f ] size= %f count= %d", gp->id,
- * cj->loc[0], */
- /* cj->loc[1], cj->loc[2], cj->h[0], cj->gcount); */
- /* } */
+ if (gp->id == -ICHECK)
+ message("id=%lld loc=[ %f %f %f ] size= %f count= %d", gp->id, cj->loc[0],
+ cj->loc[1], cj->loc[2], cj->h[0], cj->gcount);
+ }
- /* for (int pid = 0; pid < gcount_j; pid++) { */
+ for (int pid = 0; pid < gcount_j; pid++) {
- /* /\* Get a hold of the ith part in ci. *\/ */
- /* struct gpart *restrict gp = &cj->gparts[pid]; */
+ /* Get a hold of the ith part in ci. */
+ struct gpart *restrict gp = &cj->gparts[pid];
- /* if (gp->id == -ICHECK) */
- /* message("id=%lld loc=[ %f %f %f ] size= %f count= %d", gp->id,
- * ci->loc[0], */
- /* ci->loc[1], ci->loc[2], ci->h[0], ci->gcount); */
- /* } */
+ if (gp->id == -ICHECK)
+ message("id=%lld loc=[ %f %f %f ] size= %f count= %d", gp->id, ci->loc[0],
+ ci->loc[1], ci->loc[2], ci->h[0], ci->gcount);
+ }
+#endif
/* Are both cells split ? */
if (ci->split && cj->split) {
@@ -501,7 +431,9 @@ static void runner_dosub_grav(struct runner *r, struct cell *ci,
} else {
+#ifdef SANITY_CHECKS
if (!are_neighbours(ci, cj)) error("Non-neighbouring cells in pair task !");
+#endif
runner_dopair_grav(r, ci, cj);
}
@@ -510,35 +442,37 @@ static void runner_dosub_grav(struct runner *r, struct cell *ci,
static void runner_do_grav_mm(struct runner *r, struct cell *ci,
struct cell *cj) {
- /* for (int pid = 0; pid < ci->gcount; pid++) { */
+#ifdef SANITY_CHECKS
+ if (are_neighbours(ci, cj)) {
+
+ error("Non-neighbouring cells in mm task !");
- /* /\* Get a hold of the ith part in ci. *\/ */
- /* struct gpart *restrict gp = &ci->gparts[pid]; */
+#endif
- /* if (gp->id == -ICHECK) */
- /* message("id=%lld loc=[ %f %f %f ] size= %f count= %d", gp->id,
- * cj->loc[0], */
- /* cj->loc[1], cj->loc[2], cj->h[0], cj->gcount); */
- /* } */
+#if ICHECK > 0
+ for (int pid = 0; pid < ci->gcount; pid++) {
- /* for (int pid = 0; pid < cj->gcount; pid++) { */
+ /* Get a hold of the ith part in ci. */
+ struct gpart *restrict gp = &ci->gparts[pid];
- /* /\* Get a hold of the ith part in ci. *\/ */
- /* struct gpart *restrict gp = &cj->gparts[pid]; */
+ if (gp->id == -ICHECK)
+ message("id=%lld loc=[ %f %f %f ] size= %f count= %d", gp->id,
+ cj->loc[0], cj->loc[1], cj->loc[2], cj->h[0], cj->gcount);
+ }
- /* if (gp->id == -ICHECK) */
- /* message("id=%lld loc=[ %f %f %f ] size= %f count= %d", gp->id,
- * ci->loc[0], */
- /* ci->loc[1], ci->loc[2], ci->h[0], ci->gcount); */
- /* } */
+ for (int pid = 0; pid < cj->gcount; pid++) {
- if (are_neighbours(ci, cj)) {
+ /* Get a hold of the ith part in ci. */
+ struct gpart *restrict gp = &cj->gparts[pid];
- error("Non-neighbouring cells in mm task !");
- }
+ if (gp->id == -ICHECK)
+ message("id=%lld loc=[ %f %f %f ] size= %f count= %d", gp->id,
+ ci->loc[0], ci->loc[1], ci->loc[2], ci->h[0], ci->gcount);
+ }
+#endif
- runner_dopair_grav_pm(r, ci, cj);
- runner_dopair_grav_pm(r, cj, ci);
-}
+ runner_dopair_grav_pm(r, ci, cj);
+ runner_dopair_grav_pm(r, cj, ci);
+ }
#endif /* SWIFT_RUNNER_DOIACT_GRAV_H */
diff --git a/src/tools.c b/src/tools.c
index e4a4f44268081aa819752b7bbfb15ad17f86f600..30950381d0463a3e2fb7079ab01f4c1d27d1daac 100644
--- a/src/tools.c
+++ b/src/tools.c
@@ -305,7 +305,7 @@ void pairs_single_grav(double *dim, long long int pid,
fdx[i] = dx[i];
}
r2 = fdx[0] * fdx[0] + fdx[1] * fdx[1] + fdx[2] * fdx[2];
- runner_iact_grav(r2, fdx, &pi, &pj);
+ runner_iact_grav_pp(r2, fdx, &pi, &pj);
a[0] += pi.a_grav[0];
a[1] += pi.a_grav[1];
a[2] += pi.a_grav[2];