diff --git a/configure.ac b/configure.ac
index 9ec78c8b8b8168191df002a01ed5de6fb1580243..5851324ec76a466aa52f0cbecf972cf32d4a256f 100644
--- a/configure.ac
+++ b/configure.ac
@@ -226,6 +226,18 @@ if test "$enable_debugging_checks" = "yes"; then
AC_DEFINE([SWIFT_DEBUG_CHECKS],1,[Enable expensive debugging])
fi
+# Check whether we want to default to naive cell interactions
+AC_ARG_ENABLE([naive-interactions],
+ [AS_HELP_STRING([--enable-naive-interactions],
+ [Activate use of naive cell interaction functions @<:@yes/no@:>@]
+ )],
+ [enable_naive_interactions="$enableval"],
+ [enable_naive_interactions="no"]
+)
+if test "$enable_naive_interactions" = "yes"; then
+ AC_DEFINE([SWIFT_USE_NAIVE_INTERACTIONS],1,[Enable use of naive cell interaction functions])
+fi
+
# Check if gravity force checks are on for some particles.
AC_ARG_ENABLE([gravity-force-checks],
[AS_HELP_STRING([--enable-gravity-force-checks],
@@ -919,6 +931,7 @@ AC_MSG_RESULT([
Task debugging : $enable_task_debugging
Threadpool debugging : $enable_threadpool_debugging
Debugging checks : $enable_debugging_checks
+ Naive interactions : $enable_naive_interactions
Gravity checks : $gravity_force_checks
------------------------])
diff --git a/examples/main.c b/examples/main.c
index 56f463940d51cc1c30904c2472334768fa3d58d1..0d87ea350793e0d089b41d3bd05a68cd97753ef9 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -361,6 +361,13 @@ int main(int argc, char *argv[]) {
message("WARNING: Debugging checks activated. Code will be slower !");
#endif
+/* Do we have debugging checks ? */
+#ifdef SWIFT_USE_NAIVE_INTERACTIONS
+ if (myrank == 0)
+ message(
+ "WARNING: Naive cell interactions activated. Code will be slower !");
+#endif
+
/* Do we have gravity accuracy checks ? */
#ifdef SWIFT_GRAVITY_FORCE_CHECKS
if (myrank == 0)
diff --git a/src/runner_doiact.h b/src/runner_doiact.h
index 473c0604ee60d755f8aef87a96a37bef30ca6279..bcbd0b7425747ea5e9cf0bcccb1d653d8b8d8126 100644
--- a/src/runner_doiact.h
+++ b/src/runner_doiact.h
@@ -50,6 +50,9 @@
#define _DOPAIR2_NAIVE(f) PASTE(runner_dopair2_naive, f)
#define DOPAIR2_NAIVE _DOPAIR2_NAIVE(FUNCTION)
+#define _DOSELF1_NAIVE(f) PASTE(runner_doself1_naive, f)
+#define DOSELF1_NAIVE _DOSELF1_NAIVE(FUNCTION)
+
#define _DOSELF2_NAIVE(f) PASTE(runner_doself2_naive, f)
#define DOSELF2_NAIVE _DOSELF2_NAIVE(FUNCTION)
@@ -121,10 +124,6 @@ void DOPAIR1_NAIVE(struct runner *r, struct cell *restrict ci,
const struct engine *e = r->e;
-#ifndef SWIFT_DEBUG_CHECKS
- error("Don't use in actual runs ! Slow code !");
-#endif
-
TIMER_TIC;
/* Anything to do here? */
@@ -171,6 +170,14 @@ void DOPAIR1_NAIVE(struct runner *r, struct cell *restrict ci,
float dx[3] = {pix[0] - pjx[0], pix[1] - pjx[1], pix[2] - pjx[2]};
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 (pi->ti_drift != e->ti_current)
+ error("Particle pi not drifted to current time");
+ if (pj->ti_drift != e->ti_current)
+ error("Particle pj not drifted to current time");
+#endif
+
/* Hit or miss? */
if (r2 < hig2 && pi_active) {
@@ -205,10 +212,6 @@ void DOPAIR2_NAIVE(struct runner *r, struct cell *restrict ci,
const struct engine *e = r->e;
-#ifndef SWIFT_DEBUG_CHECKS
- error("Don't use in actual runs ! Slow code !");
-#endif
-
TIMER_TIC;
/* Anything to do here? */
@@ -255,6 +258,14 @@ void DOPAIR2_NAIVE(struct runner *r, struct cell *restrict ci,
float dx[3] = {pix[0] - pjx[0], pix[1] - pjx[1], pix[2] - pjx[2]};
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 (pi->ti_drift != e->ti_current)
+ error("Particle pi not drifted to current time");
+ if (pj->ti_drift != e->ti_current)
+ error("Particle pj not drifted to current time");
+#endif
+
/* Hit or miss? */
if (r2 < hig2 || r2 < hjg2) {
@@ -280,21 +291,17 @@ void DOPAIR2_NAIVE(struct runner *r, struct cell *restrict ci,
}
/**
- * @brief Compute the interactions within a cell (symmetric case).
+ * @brief Compute the interactions within a cell (non-symmetric case).
*
* Inefficient version using a brute-force algorithm.
*
* @param r The #runner.
* @param c The #cell.
*/
-void DOSELF2_NAIVE(struct runner *r, struct cell *restrict c) {
+void DOSELF1_NAIVE(struct runner *r, struct cell *restrict c) {
const struct engine *e = r->e;
-#ifndef SWIFT_DEBUG_CHECKS
- error("Don't use in actual runs ! Slow code !");
-#endif
-
TIMER_TIC;
/* Anything to do here? */
@@ -308,10 +315,11 @@ void DOSELF2_NAIVE(struct runner *r, struct cell *restrict c) {
/* Get a hold of the ith part in ci. */
struct part *restrict pi = &parts[pid];
- const double pix[3] = {pi->x[0], pi->x[1], pi->x[2]};
+ const int pi_active = part_is_active(pi, e);
const float hi = pi->h;
const float hig2 = hi * hi * kernel_gamma2;
- const int pi_active = part_is_active(pi, e);
+ const float pix[3] = {pi->x[0] - c->loc[0], pi->x[1] - c->loc[1],
+ pi->x[2] - c->loc[2]};
/* Loop over the parts in cj. */
for (int pjd = pid + 1; pjd < count; pjd++) {
@@ -319,39 +327,121 @@ void DOSELF2_NAIVE(struct runner *r, struct cell *restrict c) {
/* Get a pointer to the jth particle. */
struct part *restrict pj = &parts[pjd];
const float hj = pj->h;
+ const float hjg2 = hj * hj * kernel_gamma2;
const int pj_active = part_is_active(pj, e);
/* Compute the pairwise distance. */
- float r2 = 0.0f;
- float dx[3];
- for (int k = 0; k < 3; k++) {
- dx[k] = pix[k] - pj->x[k];
- r2 += dx[k] * dx[k];
+ const float pjx[3] = {pj->x[0] - c->loc[0], pj->x[1] - c->loc[1],
+ pj->x[2] - c->loc[2]};
+ float dx[3] = {pix[0] - pjx[0], pix[1] - pjx[1], pix[2] - pjx[2]};
+ const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
+
+ const int doi = pi_active && (r2 < hig2);
+ const int doj = pj_active && (r2 < hjg2);
+
+#ifdef SWIFT_DEBUG_CHECKS
+ /* Check that particles have been drifted to the current time */
+ if (pi->ti_drift != e->ti_current)
+ error("Particle pi not drifted to current time");
+ if (pj->ti_drift != e->ti_current)
+ error("Particle pj not drifted to current time");
+#endif
+
+ /* Hit or miss? */
+ if (doi && doj) {
+
+ IACT(r2, dx, hi, hj, pi, pj);
+ } else if (doi) {
+
+ IACT_NONSYM(r2, dx, hi, hj, pi, pj);
+ } else if (doj) {
+
+ dx[0] = -dx[0];
+ dx[1] = -dx[1];
+ dx[2] = -dx[2];
+
+ IACT_NONSYM(r2, dx, hj, hi, pj, pi);
}
+ } /* loop over the parts in cj. */
+ } /* loop over the parts in ci. */
+
+ TIMER_TOC(TIMER_DOSELF);
+}
+
+/**
+ * @brief Compute the interactions within a cell (symmetric case).
+ *
+ * Inefficient version using a brute-force algorithm.
+ *
+ * @param r The #runner.
+ * @param c The #cell.
+ */
+void DOSELF2_NAIVE(struct runner *r, struct cell *restrict c) {
+
+ const struct engine *e = r->e;
+
+ TIMER_TIC;
+
+ /* Anything to do here? */
+ if (!cell_is_active(c, e)) return;
+
+ const int count = c->count;
+ struct part *restrict parts = c->parts;
+
+ /* Loop over the parts in ci. */
+ for (int pid = 0; pid < count; pid++) {
+
+ /* Get a hold of the ith part in ci. */
+ struct part *restrict pi = &parts[pid];
+ const int pi_active = part_is_active(pi, e);
+ const float hi = pi->h;
+ const float hig2 = hi * hi * kernel_gamma2;
+ const float pix[3] = {pi->x[0] - c->loc[0], pi->x[1] - c->loc[1],
+ pi->x[2] - c->loc[2]};
+
+ /* Loop over the parts in cj. */
+ for (int pjd = pid + 1; pjd < count; pjd++) {
+
+ /* Get a pointer to the jth particle. */
+ struct part *restrict pj = &parts[pjd];
+ const float hj = pj->h;
const float hjg2 = hj * hj * kernel_gamma2;
+ const int pj_active = part_is_active(pj, e);
- /* Hit or miss? */
- if (r2 < hig2 || r2 < hjg2) {
+ /* Compute the pairwise distance. */
+ const float pjx[3] = {pj->x[0] - c->loc[0], pj->x[1] - c->loc[1],
+ pj->x[2] - c->loc[2]};
+ float dx[3] = {pix[0] - pjx[0], pix[1] - pjx[1], pix[2] - pjx[2]};
+ const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
- if (pi_active && pj_active) {
+ const int doi = pi_active && ((r2 < hig2) || (r2 < hjg2));
+ const int doj = pj_active && ((r2 < hig2) || (r2 < hjg2));
- IACT(r2, dx, hi, hj, pi, pj);
- } else if (pi_active) {
+#ifdef SWIFT_DEBUG_CHECKS
+ /* Check that particles have been drifted to the current time */
+ if (pi->ti_drift != e->ti_current)
+ error("Particle pi not drifted to current time");
+ if (pj->ti_drift != e->ti_current)
+ error("Particle pj not drifted to current time");
+#endif
- IACT_NONSYM(r2, dx, hi, hj, pi, pj);
- } else if (pj_active) {
+ /* Hit or miss? */
+ if (doi && doj) {
- dx[0] = -dx[0];
- dx[1] = -dx[1];
- dx[2] = -dx[2];
+ IACT(r2, dx, hi, hj, pi, pj);
+ } else if (doi) {
- IACT_NONSYM(r2, dx, hj, hi, pj, pi);
- }
- }
+ IACT_NONSYM(r2, dx, hi, hj, pi, pj);
+ } else if (doj) {
- } /* loop over the parts in cj. */
+ dx[0] = -dx[0];
+ dx[1] = -dx[1];
+ dx[2] = -dx[2];
- } /* loop over the parts in ci. */
+ IACT_NONSYM(r2, dx, hj, hi, pj, pi);
+ }
+ } /* loop over the parts in cj. */
+ } /* loop over the parts in ci. */
TIMER_TOC(TIMER_DOSELF);
}
@@ -375,10 +465,6 @@ void DOPAIR_SUBSET_NAIVE(struct runner *r, struct cell *restrict ci,
const struct engine *e = r->e;
-#ifndef SWIFT_DEBUG_CHECKS
- error("Don't use in actual runs ! Slow code !");
-#endif
-
TIMER_TIC;
const int count_j = cj->count;
@@ -459,6 +545,11 @@ void DOPAIR_SUBSET(struct runner *r, struct cell *restrict ci,
const struct engine *e = r->e;
+#ifdef SWIFT_USE_NAIVE_INTERACTIONS
+ DOPAIR_SUBSET_NAIVE(r, ci, parts_i, ind, count, cj);
+ return;
+#endif
+
TIMER_TIC;
const int count_j = cj->count;
@@ -523,6 +614,14 @@ void DOPAIR_SUBSET(struct runner *r, struct cell *restrict ci,
float dx[3] = {pix - pjx, piy - pjy, piz - pjz};
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 (pi->ti_drift != e->ti_current)
+ error("Particle pi not drifted to current time");
+ if (pj->ti_drift != e->ti_current)
+ error("Particle pj not drifted to current time");
+#endif
+
/* Hit or miss? */
if (r2 < hig2) {
@@ -562,6 +661,14 @@ void DOPAIR_SUBSET(struct runner *r, struct cell *restrict ci,
float dx[3] = {pix - pjx, piy - pjy, piz - pjz};
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 (pi->ti_drift != e->ti_current)
+ error("Particle pi not drifted to current time");
+ if (pj->ti_drift != e->ti_current)
+ error("Particle pj not drifted to current time");
+#endif
+
/* Hit or miss? */
if (r2 < hig2) {
@@ -622,6 +729,14 @@ void DOSELF_SUBSET(struct runner *r, struct cell *restrict ci,
float dx[3] = {pix[0] - pjx[0], pix[1] - pjx[1], pix[2] - pjx[2]};
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 (pi->ti_drift != e->ti_current)
+ error("Particle pi not drifted to current time");
+ if (pj->ti_drift != e->ti_current)
+ error("Particle pj not drifted to current time");
+#endif
+
/* Hit or miss? */
if (r2 > 0.f && r2 < hig2) {
@@ -647,8 +762,10 @@ void DOPAIR1(struct runner *r, struct cell *ci, struct cell *cj, const int sid,
const struct engine *restrict e = r->e;
- // DOPAIR1_NAIVE(r, ci, cj);
- // return;
+#ifdef SWIFT_USE_NAIVE_INTERACTIONS
+ DOPAIR1_NAIVE(r, ci, cj);
+ return;
+#endif
TIMER_TIC;
@@ -869,8 +986,10 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj) {
struct engine *restrict e = r->e;
- // DOPAIR2_NAIVE(r, ci, cj);
- // return;
+#ifdef SWIFT_USE_NAIVE_INTERACTIONS
+ DOPAIR2_NAIVE(r, ci, cj);
+ return;
+#endif
TIMER_TIC;
@@ -1207,6 +1326,11 @@ void DOSELF1(struct runner *r, struct cell *restrict c) {
const struct engine *e = r->e;
+#ifdef SWIFT_USE_NAIVE_INTERACTIONS
+ DOSELF1_NAIVE(r, c);
+ return;
+#endif
+
TIMER_TIC;
if (!cell_is_active(c, e)) return;
@@ -1339,6 +1463,11 @@ void DOSELF2(struct runner *r, struct cell *restrict c) {
const struct engine *e = r->e;
+#ifdef SWIFT_USE_NAIVE_INTERACTIONS
+ DOSELF2_NAIVE(r, c);
+ return;
+#endif
+
TIMER_TIC;
if (!cell_is_active(c, e)) return;