diff --git a/src/gravity/Default/gravity_iact.h b/src/gravity/Default/gravity_iact.h
index 5beba6f1b7603cb0231990974d7881942414b4eb..d4a95540de17631ad445075d672d03a1236e34e3 100644
--- a/src/gravity/Default/gravity_iact.h
+++ b/src/gravity/Default/gravity_iact.h
@@ -28,11 +28,11 @@
#include "vector.h"
/**
- * @brief Gravity forces between particles
+ * @brief Gravity forces between particles truncated by the long-range kernel
*/
-__attribute__((always_inline)) INLINE static void runner_iact_grav_pp(
- float rlr_inv, float r2, const float *dx, struct gpart *gpi,
- struct gpart *gpj, int periodic) {
+__attribute__((always_inline)) INLINE static void runner_iact_grav_pp_truncated(
+ float r2, const float *dx, struct gpart *gpi, struct gpart *gpj,
+ float rlr_inv) {
/* Apply the gravitational acceleration. */
const float r = sqrtf(r2);
@@ -41,7 +41,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pp(
const float mj = gpj->mass;
const float hi = gpi->epsilon;
const float hj = gpj->epsilon;
- const float u = r * rlr_inv;
+ const float u_lr = r * rlr_inv;
float f_lr, fi, fj, W;
#ifdef SWIFT_DEBUG_CHECKS
@@ -49,10 +49,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pp(
#endif
/* Get long-range correction */
- if (periodic)
- kernel_long_grav_eval(u, &f_lr);
- else
- f_lr = 1.f;
+ kernel_long_grav_eval(u_lr, &f_lr);
if (r >= hi) {
@@ -100,18 +97,84 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pp(
}
/**
- * @brief Gravity forces between particles (non-symmetric version)
+ * @brief Gravity forces between particles
*/
-__attribute__((always_inline)) INLINE static void runner_iact_grav_pp_nonsym(
- float rlr_inv, float r2, const float *dx, struct gpart *gpi,
- const struct gpart *gpj, int periodic) {
+__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 r = sqrtf(r2);
+ const float ir = 1.f / r;
+ const float mi = gpi->mass;
+ const float mj = gpj->mass;
+ const float hi = gpi->epsilon;
+ const float hj = gpj->epsilon;
+ float fi, fj, W;
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (r == 0.f) error("Interacting particles with 0 distance");
+#endif
+
+ if (r >= hi) {
+
+ /* Get Newtonian gravity */
+ fi = mj * ir * ir * ir;
+
+ } else {
+
+ const float hi_inv = 1.f / hi;
+ const float hi_inv3 = hi_inv * hi_inv * hi_inv;
+ const float ui = r * hi_inv;
+
+ kernel_grav_eval(ui, &W);
+
+ /* Get softened gravity */
+ fi = mj * hi_inv3 * W;
+ }
+
+ if (r >= hj) {
+
+ /* Get Newtonian gravity */
+ fj = mi * ir * ir * ir;
+
+ } else {
+
+ const float hj_inv = 1.f / hj;
+ const float hj_inv3 = hj_inv * hj_inv * hj_inv;
+ const float uj = r * hj_inv;
+
+ kernel_grav_eval(uj, &W);
+
+ /* Get softened gravity */
+ fj = mi * hj_inv3 * W;
+ }
+
+ const float fidx[3] = {fi * dx[0], fi * dx[1], fi * dx[2]};
+ gpi->a_grav[0] -= fidx[0];
+ gpi->a_grav[1] -= fidx[1];
+ gpi->a_grav[2] -= fidx[2];
+
+ const float fjdx[3] = {fj * dx[0], fj * dx[1], fj * dx[2]};
+ gpj->a_grav[0] += fjdx[0];
+ gpj->a_grav[1] += fjdx[1];
+ gpj->a_grav[2] += fjdx[2];
+}
+
+/**
+ * @brief Gravity forces between particles truncated by the long-range kernel
+ * (non-symmetric version)
+ */
+__attribute__((always_inline)) INLINE static void
+runner_iact_grav_pp_truncated_nonsym(float r2, const float *dx,
+ struct gpart *gpi, const struct gpart *gpj,
+ float rlr_inv) {
/* Apply the gravitational acceleration. */
const float r = sqrtf(r2);
const float ir = 1.f / r;
const float mj = gpj->mass;
const float hi = gpi->epsilon;
- const float u = r * rlr_inv;
+ const float u_lr = r * rlr_inv;
float f_lr, f, W;
#ifdef SWIFT_DEBUG_CHECKS
@@ -119,10 +182,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pp_nonsym(
#endif
/* Get long-range correction */
- if (periodic)
- kernel_long_grav_eval(u, &f_lr);
- else
- f_lr = 1.f;
+ kernel_long_grav_eval(u_lr, &f_lr);
if (r >= hi) {
@@ -149,13 +209,44 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pp_nonsym(
}
/**
- * @brief Gravity forces between particle and multipole
+ * @brief Gravity forces between particles (non-symmetric version)
*/
-__attribute__((always_inline)) INLINE static void runner_iact_grav_pm(
- float rlr_inv, float r2, const float *dx, struct gpart *gp,
- const struct multipole *multi) {
+__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 r = sqrtf(r2);
+ const float ir = 1.f / r;
+ const float mj = gpj->mass;
+ const float hi = gpi->epsilon;
+ float f, W;
- error("Dead function");
+#ifdef SWIFT_DEBUG_CHECKS
+ if (r == 0.f) error("Interacting particles with 0 distance");
+#endif
+
+ if (r >= hi) {
+
+ /* Get Newtonian gravity */
+ f = mj * ir * ir * ir;
+
+ } else {
+
+ const float hi_inv = 1.f / hi;
+ const float hi_inv3 = hi_inv * hi_inv * hi_inv;
+ const float ui = r * hi_inv;
+
+ kernel_grav_eval(ui, &W);
+
+ /* Get softened gravity */
+ f = mj * hi_inv3 * W;
+ }
+
+ const float fdx[3] = {f * dx[0], f * dx[1], f * dx[2]};
+
+ gpi->a_grav[0] -= fdx[0];
+ gpi->a_grav[1] -= fdx[1];
+ gpi->a_grav[2] -= fdx[2];
}
#endif /* SWIFT_DEFAULT_GRAVITY_IACT_H */
diff --git a/src/runner_doiact_grav.h b/src/runner_doiact_grav.h
index acf6478a09d6a4d9670ada8431fd288628296bf4..702274e4d1d9b2ed6163c4e2a332a1f773d5aad8 100644
--- a/src/runner_doiact_grav.h
+++ b/src/runner_doiact_grav.h
@@ -149,41 +149,21 @@ void runner_dopair_grav_mm(const struct runner *r, struct cell *restrict ci,
TIMER_TOC(timer_dopair_grav_mm);
}
-/**
- * @brief Computes the interaction of all the particles in a cell with the
- * multipole of another cell.
- *
- * @param r The #runner.
- * @param ci The #cell with particles to interct.
- * @param cj The #cell with the multipole.
- */
-void runner_dopair_grav_pm(const struct runner *r,
- const struct cell *restrict ci,
- const struct cell *restrict cj) {
-
- error("Function should not be called");
-}
-
/**
* @brief Computes the interaction of all the particles in a cell with all the
- * particles of another cell.
+ * particles of another cell using the full Newtonian potential
*
* @param r The #runner.
* @param ci The first #cell.
* @param cj The other #cell.
- *
- * @todo Use a local cache for the particles.
+ * @param shift The distance vector (periodically wrapped) between the cell
+ * centres.
*/
-void runner_dopair_grav_pp(struct runner *r, struct cell *ci, struct cell *cj) {
+void runner_dopair_grav_pp_full(struct runner *r, struct cell *ci,
+ struct cell *cj, double shift[3]) {
/* Some constants */
- const struct engine *e = r->e;
- const struct space *s = e->s;
- const int periodic = s->periodic;
- const double cell_width = s->width[0];
- const double dim[3] = {s->dim[0], s->dim[1], s->dim[2]};
- const float a_smooth = e->gravity_properties->a_smooth;
- const float rlr_inv = 1. / (a_smooth * cell_width);
+ const struct engine *const e = r->e;
/* Cell properties */
const int gcount_i = ci->gcount;
@@ -191,46 +171,121 @@ void runner_dopair_grav_pp(struct runner *r, struct cell *ci, struct cell *cj) {
struct gpart *restrict gparts_i = ci->gparts;
struct gpart *restrict gparts_j = cj->gparts;
- TIMER_TIC;
+ /* MATTHIEU: Should we use local DP accumulators ? */
- /* Anything to do here? */
- if (!cell_is_active(ci, e) && !cell_is_active(cj, e)) return;
+ /* Loop over all particles in ci... */
+ if (cell_is_active(ci, e)) {
+ for (int pid = 0; pid < gcount_i; pid++) {
- /* Let's start by drifting things */
- if (!cell_are_gpart_drifted(ci, e)) cell_drift_gpart(ci, e);
- if (!cell_are_gpart_drifted(cj, e)) cell_drift_gpart(cj, e);
+ /* Get a hold of the ith part in ci. */
+ struct gpart *restrict gpi = &gparts_i[pid];
-#if ICHECK > 0
- for (int pid = 0; pid < gcount_i; pid++) {
+ if (!gpart_is_active(gpi, e)) continue;
- /* Get a hold of the ith part in ci. */
- struct gpart *restrict gp = &gparts_i[pid];
+ /* Apply boundary condition */
+ const double pix[3] = {gpi->x[0] - shift[0], gpi->x[1] - shift[1],
+ gpi->x[2] - shift[2]};
- 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);
- }
+ /* Loop over every particle in the other cell. */
+ for (int pjd = 0; pjd < gcount_j; pjd++) {
- for (int pid = 0; pid < gcount_j; pid++) {
+ /* Get a hold of the jth part in cj. */
+ const struct gpart *restrict gpj = &gparts_j[pjd];
- /* Get a hold of the ith part in ci. */
- struct gpart *restrict gp = &gparts_j[pid];
+ /* Compute the pairwise distance. */
+ const float dx[3] = {pix[0] - gpj->x[0], // x
+ pix[1] - gpj->x[1], // y
+ pix[2] - gpj->x[2]}; // z
+ const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
- if (gp->id_or_neg_offset == ICHECK)
- message("id=%lld loc=[ %f %f %f ] size= %f count=%d",
- gp->id_or_neg_offset, ci->loc[0], ci->loc[1], ci->loc[2],
- ci->width[0], ci->gcount);
+#ifdef SWIFT_DEBUG_CHECKS
+ /* Check that particles have been drifted to the current time */
+ if (gpi->ti_drift != e->ti_current)
+ error("gpi not drifted to current time");
+ if (gpj->ti_drift != e->ti_current)
+ error("gpj not drifted to current time");
+#endif
+
+ /* Interact ! */
+ runner_iact_grav_pp_nonsym(r2, dx, gpi, gpj);
+
+#ifdef SWIFT_DEBUG_CHECKS
+ gpi->num_interacted++;
+#endif
+ }
+ }
}
+
+ /* Loop over all particles in cj... */
+ if (cell_is_active(cj, e)) {
+ for (int pjd = 0; pjd < gcount_j; pjd++) {
+
+ /* Get a hold of the ith part in ci. */
+ struct gpart *restrict gpj = &gparts_j[pjd];
+
+ if (!gpart_is_active(gpj, e)) continue;
+
+ /* Apply boundary condition */
+ const double pjx[3] = {gpj->x[0] + shift[0], gpj->x[1] + shift[1],
+ gpj->x[2] + shift[2]};
+
+ /* Loop over every particle in the other cell. */
+ for (int pid = 0; pid < gcount_i; pid++) {
+
+ /* Get a hold of the ith part in ci. */
+ const struct gpart *restrict gpi = &gparts_i[pid];
+
+ /* Compute the pairwise distance. */
+ const float dx[3] = {pjx[0] - gpi->x[0], // x
+ pjx[1] - gpi->x[1], // y
+ pjx[2] - gpi->x[2]}; // z
+ const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
+
+#ifdef SWIFT_DEBUG_CHECKS
+ /* Check that particles have been drifted to the current time */
+ if (gpi->ti_drift != e->ti_current)
+ error("gpi not drifted to current time");
+ if (gpj->ti_drift != e->ti_current)
+ error("gpj not drifted to current time");
#endif
- /* Get the relative distance between the pairs, wrapping. */
- double shift[3] = {0.0, 0.0, 0.0};
- if (periodic) {
- shift[0] = nearest(cj->loc[0] - ci->loc[0], dim[0]);
- shift[1] = nearest(cj->loc[1] - ci->loc[1], dim[1]);
- shift[2] = nearest(cj->loc[2] - ci->loc[2], dim[2]);
+ /* Interact ! */
+ runner_iact_grav_pp_nonsym(r2, dx, gpj, gpi);
+
+#ifdef SWIFT_DEBUG_CHECKS
+ gpj->num_interacted++;
+#endif
+ }
+ }
}
+}
+
+/**
+ * @brief Computes the interaction of all the particles in a cell with all the
+ * particles of another cell using the truncated Newtonian potential
+ *
+ * @param r The #runner.
+ * @param ci The first #cell.
+ * @param cj The other #cell.
+ * @param shift The distance vector (periodically wrapped) between the cell
+ * centres.
+ */
+void runner_dopair_grav_pp_truncated(struct runner *r, struct cell *ci,
+ struct cell *cj, double shift[3]) {
+
+ /* Some constants */
+ const struct engine *const e = r->e;
+ const struct space *s = e->s;
+ const double cell_width = s->width[0];
+ const double a_smooth = e->gravity_properties->a_smooth;
+ const double rlr = cell_width * a_smooth;
+ const float rlr_inv = 1. / rlr;
+
+ /* Cell properties */
+ const int gcount_i = ci->gcount;
+ const int gcount_j = cj->gcount;
+ struct gpart *restrict gparts_i = ci->gparts;
+ struct gpart *restrict gparts_j = cj->gparts;
/* MATTHIEU: Should we use local DP accumulators ? */
@@ -268,7 +323,7 @@ void runner_dopair_grav_pp(struct runner *r, struct cell *ci, struct cell *cj) {
#endif
/* Interact ! */
- runner_iact_grav_pp_nonsym(rlr_inv, r2, dx, gpi, gpj, periodic);
+ runner_iact_grav_pp_truncated_nonsym(r2, dx, gpi, gpj, rlr_inv);
#ifdef SWIFT_DEBUG_CHECKS
gpi->num_interacted++;
@@ -311,7 +366,7 @@ void runner_dopair_grav_pp(struct runner *r, struct cell *ci, struct cell *cj) {
#endif
/* Interact ! */
- runner_iact_grav_pp_nonsym(rlr_inv, r2, dx, gpj, gpi, periodic);
+ runner_iact_grav_pp_truncated_nonsym(r2, dx, gpj, gpi, rlr_inv);
#ifdef SWIFT_DEBUG_CHECKS
gpj->num_interacted++;
@@ -319,56 +374,167 @@ void runner_dopair_grav_pp(struct runner *r, struct cell *ci, struct cell *cj) {
}
}
}
+}
+
+/**
+ * @brief Computes the interaction of all the particles in a cell with all the
+ * particles of another cell (switching function between full and truncated).
+ *
+ * @param r The #runner.
+ * @param ci The first #cell.
+ * @param cj The other #cell.
+ */
+void runner_dopair_grav_pp(struct runner *r, struct cell *ci, struct cell *cj) {
+
+ /* Some properties of the space */
+ const struct engine *e = r->e;
+ const struct space *s = e->s;
+ const int periodic = s->periodic;
+ const double cell_width = s->width[0];
+ const double dim[3] = {s->dim[0], s->dim[1], s->dim[2]};
+ const double a_smooth = e->gravity_properties->a_smooth;
+ const double r_cut_min = e->gravity_properties->r_cut_min;
+ const double min_trunc = cell_width * r_cut_min * a_smooth;
+ double shift[3] = {0.0, 0.0, 0.0};
+
+ TIMER_TIC;
+
+ /* Anything to do here? */
+ if (!cell_is_active(ci, e) && !cell_is_active(cj, e)) return;
+
+ /* Let's start by drifting things */
+ if (!cell_are_gpart_drifted(ci, e)) cell_drift_gpart(ci, e);
+ if (!cell_are_gpart_drifted(cj, e)) cell_drift_gpart(cj, e);
+
+ /* Can we use the Newtonian version or do we need the truncated one ? */
+ if (!periodic) {
+ runner_dopair_grav_pp_full(r, ci, cj, shift);
+ } else {
+
+ /* Get the relative distance between the pairs, wrapping. */
+ shift[0] = nearest(cj->loc[0] - ci->loc[0], dim[0]);
+ shift[1] = nearest(cj->loc[1] - ci->loc[1], dim[1]);
+ shift[2] = nearest(cj->loc[2] - ci->loc[2], dim[2]);
+ const double r2 =
+ shift[0] * shift[0] + shift[1] * shift[1] + shift[2] * shift[2];
+
+ /* Get the maximal distance between any two particles */
+ const double max_r = sqrt(r2) + ci->multipole->r_max + cj->multipole->r_max;
+
+ /* Do we need to use the truncated interactions ? */
+ if (max_r > min_trunc)
+ runner_dopair_grav_pp_truncated(r, ci, cj, shift);
+ else
+ runner_dopair_grav_pp_full(r, ci, cj, shift);
+ }
TIMER_TOC(timer_dopair_grav_pp);
}
/**
- * @brief Computes the interaction of all the particles in a cell directly
+ * @brief Computes the interaction of all the particles in a cell using the
+ * full Newtonian potential.
*
* @param r The #runner.
* @param c The #cell.
*
* @todo Use a local cache for the particles.
*/
-void runner_doself_grav_pp(struct runner *r, struct cell *c) {
+void runner_doself_grav_pp_full(struct runner *r, struct cell *c) {
/* Some constants */
- const struct engine *e = r->e;
- const struct space *s = e->s;
- const int periodic = s->periodic;
- const double cell_width = s->width[0];
- const float a_smooth = e->gravity_properties->a_smooth;
- const float rlr_inv = 1. / (a_smooth * cell_width);
+ const struct engine *const e = r->e;
/* Cell properties */
const int gcount = c->gcount;
struct gpart *restrict gparts = c->gparts;
- TIMER_TIC;
+ /* MATTHIEU: Should we use local DP accumulators ? */
+
+ /* 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];
+
+ /* 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];
+
+ /* Compute the pairwise distance. */
+ float dx[3] = {gpi->x[0] - gpj->x[0], // x
+ gpi->x[1] - gpj->x[1], // y
+ gpi->x[2] - gpj->x[2]}; // z
+ const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
#ifdef SWIFT_DEBUG_CHECKS
- if (c->gcount == 0) error("Doing self gravity on an empty cell !");
+ /* Check that particles have been drifted to the current time */
+ if (gpi->ti_drift != e->ti_current)
+ error("gpi not drifted to current time");
+ if (gpj->ti_drift != e->ti_current)
+ error("gpj not drifted to current time");
#endif
- /* Anything to do here? */
- if (!cell_is_active(c, e)) return;
+ /* Interact ! */
+ if (gpart_is_active(gpi, e) && gpart_is_active(gpj, e)) {
- /* Do we need to start by drifting things ? */
- if (!cell_are_gpart_drifted(c, e)) cell_drift_gpart(c, e);
+ runner_iact_grav_pp(r2, dx, gpi, gpj);
-#if ICHECK > 0
- for (int pid = 0; pid < gcount; pid++) {
+#ifdef SWIFT_DEBUG_CHECKS
+ gpi->num_interacted++;
+ gpj->num_interacted++;
+#endif
- /* Get a hold of the ith part in ci. */
- struct gpart *restrict gp = &gparts[pid];
+ } else {
- if (gp->id_or_neg_offset == ICHECK)
- message("id=%lld loc=[ %f %f %f ] size= %f count= %d",
- gp->id_or_neg_offset, c->loc[0], c->loc[1], c->loc[2],
- c->width[0], c->gcount);
- }
+ if (gpart_is_active(gpi, e)) {
+
+ runner_iact_grav_pp_nonsym(r2, dx, gpi, gpj);
+
+#ifdef SWIFT_DEBUG_CHECKS
+ gpi->num_interacted++;
+#endif
+
+ } else if (gpart_is_active(gpj, e)) {
+
+ dx[0] = -dx[0];
+ dx[1] = -dx[1];
+ dx[2] = -dx[2];
+ runner_iact_grav_pp_nonsym(r2, dx, gpj, gpi);
+
+#ifdef SWIFT_DEBUG_CHECKS
+ gpj->num_interacted++;
#endif
+ }
+ }
+ }
+ }
+}
+
+/**
+ * @brief Computes the interaction of all the particles in a cell using the
+ * truncated Newtonian potential.
+ *
+ * @param r The #runner.
+ * @param c The #cell.
+ *
+ * @todo Use a local cache for the particles.
+ */
+void runner_doself_grav_pp_truncated(struct runner *r, struct cell *c) {
+
+ /* Some constants */
+ const struct engine *const e = r->e;
+ const struct space *s = e->s;
+ const double cell_width = s->width[0];
+ const double a_smooth = e->gravity_properties->a_smooth;
+ const double rlr = cell_width * a_smooth;
+ const float rlr_inv = 1. / rlr;
+
+ /* Cell properties */
+ const int gcount = c->gcount;
+ struct gpart *restrict gparts = c->gparts;
/* MATTHIEU: Should we use local DP accumulators ? */
@@ -401,7 +567,7 @@ void runner_doself_grav_pp(struct runner *r, struct cell *c) {
/* Interact ! */
if (gpart_is_active(gpi, e) && gpart_is_active(gpj, e)) {
- runner_iact_grav_pp(rlr_inv, r2, dx, gpi, gpj, periodic);
+ runner_iact_grav_pp_truncated(r2, dx, gpi, gpj, rlr_inv);
#ifdef SWIFT_DEBUG_CHECKS
gpi->num_interacted++;
@@ -412,7 +578,7 @@ void runner_doself_grav_pp(struct runner *r, struct cell *c) {
if (gpart_is_active(gpi, e)) {
- runner_iact_grav_pp_nonsym(rlr_inv, r2, dx, gpi, gpj, periodic);
+ runner_iact_grav_pp_truncated_nonsym(r2, dx, gpi, gpj, rlr_inv);
#ifdef SWIFT_DEBUG_CHECKS
gpi->num_interacted++;
@@ -423,7 +589,7 @@ void runner_doself_grav_pp(struct runner *r, struct cell *c) {
dx[0] = -dx[0];
dx[1] = -dx[1];
dx[2] = -dx[2];
- runner_iact_grav_pp_nonsym(rlr_inv, r2, dx, gpj, gpi, periodic);
+ runner_iact_grav_pp_truncated_nonsym(r2, dx, gpj, gpi, rlr_inv);
#ifdef SWIFT_DEBUG_CHECKS
gpj->num_interacted++;
@@ -432,6 +598,52 @@ void runner_doself_grav_pp(struct runner *r, struct cell *c) {
}
}
}
+}
+
+/**
+ * @brief Computes the interaction of all the particles in a cell directly
+ * (Switching function between truncated and full)
+ *
+ * @param r The #runner.
+ * @param c The #cell.
+ */
+void runner_doself_grav_pp(struct runner *r, struct cell *c) {
+
+ /* Some properties of the space */
+ const struct engine *e = r->e;
+ const struct space *s = e->s;
+ const int periodic = s->periodic;
+ const double cell_width = s->width[0];
+ const double a_smooth = e->gravity_properties->a_smooth;
+ const double r_cut_min = e->gravity_properties->r_cut_min;
+ const double min_trunc = cell_width * r_cut_min * a_smooth;
+
+ TIMER_TIC;
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (c->gcount == 0) error("Doing self gravity on an empty cell !");
+#endif
+
+ /* Anything to do here? */
+ if (!cell_is_active(c, e)) return;
+
+ /* Do we need to start by drifting things ? */
+ if (!cell_are_gpart_drifted(c, e)) cell_drift_gpart(c, e);
+
+ /* Can we use the Newtonian version or do we need the truncated one ? */
+ if (!periodic) {
+ runner_doself_grav_pp_full(r, c);
+ } else {
+
+ /* Get the maximal distance between any two particles */
+ const double max_r = 2 * c->multipole->r_max;
+
+ /* Do we need to use the truncated interactions ? */
+ if (max_r > min_trunc)
+ runner_doself_grav_pp_truncated(r, c);
+ else
+ runner_doself_grav_pp_full(r, c);
+ }
TIMER_TOC(timer_doself_grav_pp);
}
diff --git a/src/tools.c b/src/tools.c
index 0f423d5aad620edad3a430e7ea962218430da86d..68b23e1a815b510abafbd8f7c6285d3633ba0a63 100644
--- a/src/tools.c
+++ b/src/tools.c
@@ -430,7 +430,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_pp(0.f, r2, fdx, &pi, &pj, 0);
+ 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];
@@ -790,7 +790,7 @@ void gravity_n2(struct gpart *gparts, const int gcount,
if (r2 < max_d2 || 1) {
/* Apply the gravitational acceleration. */
- runner_iact_grav_pp(rlr_inv, r2, dx, gpi, gpj, 0);
+ runner_iact_grav_pp(r2, dx, gpi, gpj);
}
}
}