diff --git a/configure.ac b/configure.ac
index d56d01ecc42fcf642c6b30437e09ef2f6557940f..fbf097976a7e67b8b35e20f1d66d231738fd80f6 100644
--- a/configure.ac
+++ b/configure.ac
@@ -306,6 +306,18 @@ if test "$enable_naive_interactions" = "yes"; then
AC_DEFINE([SWIFT_USE_NAIVE_INTERACTIONS],1,[Enable use of naive cell interaction functions])
fi
+# Check whether we want to default to naive cell interactions (stars)
+AC_ARG_ENABLE([naive-interactions-stars],
+ [AS_HELP_STRING([--enable-naive-interactions-stars],
+ [Activate use of naive cell interaction functions for stars @<:@yes/no@:>@]
+ )],
+ [enable_naive_interactions_stars="$enableval"],
+ [enable_naive_interactions_stars="no"]
+)
+if test "$enable_naive_interactions_stars" = "yes"; then
+ AC_DEFINE([SWIFT_USE_NAIVE_INTERACTIONS_STARS],1,[Enable use of naive cell interaction functions for stars])
+fi
+
# Check if gravity force checks are on for some particles.
AC_ARG_ENABLE([gravity-force-checks],
[AS_HELP_STRING([--enable-gravity-force-checks],
@@ -1943,6 +1955,7 @@ AC_MSG_RESULT([
Interaction debugging : $enable_debug_interactions
Stars interaction debugging : $enable_debug_interactions_stars
Naive interactions : $enable_naive_interactions
+ Naive stars interactions : $enable_naive_interactions_stars
Gravity checks : $gravity_force_checks
Custom icbrtf : $enable_custom_icbrtf
diff --git a/src/runner.c b/src/runner.c
index a44858659069b3ceea2e600e7e21eda798fff876..719c193e3d739bd929b008abff35eb6ad102a701 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -74,6 +74,7 @@
#define TASK_LOOP_GRADIENT 1
#define TASK_LOOP_FORCE 2
#define TASK_LOOP_LIMITER 3
+#define TASK_LOOP_FEEDBACK 4
/* Import the density loop functions. */
#define FUNCTION density
@@ -110,14 +111,16 @@
/* Import the stars density loop functions. */
#define FUNCTION density
-#define UPDATE_STARS 1
+#define FUNCTION_TASK_LOOP TASK_LOOP_DENSITY
#include "runner_doiact_stars.h"
-#undef UPDATE_STARS
+#undef FUNCTION_TASK_LOOP
#undef FUNCTION
/* Import the stars feedback loop functions. */
#define FUNCTION feedback
+#define FUNCTION_TASK_LOOP TASK_LOOP_FEEDBACK
#include "runner_doiact_stars.h"
+#undef FUNCTION_TASK_LOOP
#undef FUNCTION
/**
diff --git a/src/runner_doiact_stars.h b/src/runner_doiact_stars.h
index f208f14ac98a31a55df6741a2cc5f9cb0a829762..1ce0c370f31d860dd36e5af54020c8feb0da6c4f 100644
--- a/src/runner_doiact_stars.h
+++ b/src/runner_doiact_stars.h
@@ -28,15 +28,23 @@
#define _DOSELF1_STARS(f) PASTE(runner_doself_stars, f)
#define DOSELF1_STARS _DOSELF1_STARS(FUNCTION)
-#define _DO_NONSYM_PAIR1_STARS(f) PASTE(runner_do_nonsym_pair_stars, f)
-#define DO_NONSYM_PAIR1_STARS _DO_NONSYM_PAIR1_STARS(FUNCTION)
+#define _DO_SYM_PAIR1_STARS(f) PASTE(runner_do_sym_pair_stars, f)
+#define DO_SYM_PAIR1_STARS _DO_SYM_PAIR1_STARS(FUNCTION)
-#define _DOPAIR1_STARS(f) PASTE(runner_dopair_stars, f)
-#define DOPAIR1_STARS _DOPAIR1_STARS(FUNCTION)
+#define _DO_NONSYM_PAIR1_STARS_NAIVE(f) \
+ PASTE(runner_do_nonsym_pair_stars_naive, f)
+#define DO_NONSYM_PAIR1_STARS_NAIVE _DO_NONSYM_PAIR1_STARS_NAIVE(FUNCTION)
+
+#define _DOPAIR1_STARS_NAIVE(f) PASTE(runner_dopair_stars_naive, f)
+#define DOPAIR1_STARS_NAIVE _DOPAIR1_STARS_NAIVE(FUNCTION)
#define _DOPAIR1_SUBSET_STARS(f) PASTE(runner_dopair_subset_stars, f)
#define DOPAIR1_SUBSET_STARS _DOPAIR1_SUBSET_STARS(FUNCTION)
+#define _DOPAIR1_SUBSET_STARS_NAIVE(f) \
+ PASTE(runner_dopair_subset_stars_naive, f)
+#define DOPAIR1_SUBSET_STARS_NAIVE _DOPAIR1_SUBSET_STARS_NAIVE(FUNCTION)
+
#define _DOSELF1_SUBSET_STARS(f) PASTE(runner_doself_subset_stars, f)
#define DOSELF1_SUBSET_STARS _DOSELF1_SUBSET_STARS(FUNCTION)
@@ -83,8 +91,8 @@ void DOSELF1_STARS(struct runner *r, struct cell *c, int timer) {
const struct cosmology *cosmo = e->cosmology;
/* Anything to do here? */
+ if (c->hydro.count == 0 || c->stars.count == 0) return;
if (!cell_is_active_stars(c, e)) return;
- if (c->hydro.count == 0 && c->stars.count == 0) return;
/* Cosmological terms */
const float a = cosmo->a;
@@ -100,7 +108,7 @@ void DOSELF1_STARS(struct runner *r, struct cell *c, int timer) {
/* Get a hold of the ith spart in ci. */
struct spart *restrict si = &sparts[sid];
- if (!spart_is_active(si, e) || spart_is_inhibited(si, e)) continue;
+ if (!spart_is_active(si, e)) continue;
const float hi = si->h;
const float hig2 = hi * hi * kernel_gamma2;
@@ -113,7 +121,9 @@ void DOSELF1_STARS(struct runner *r, struct cell *c, int timer) {
/* Get a pointer to the jth particle. */
struct part *restrict pj = &parts[pjd];
- const float hj = pj->h;
+
+ /* Early abort? */
+ if (part_is_inhibited(pj, e)) continue;
/* Compute the pairwise distance. */
const float pjx[3] = {(float)(pj->x[0] - c->loc[0]),
@@ -128,8 +138,8 @@ void DOSELF1_STARS(struct runner *r, struct cell *c, int timer) {
error("Particle pj not drifted to current time");
#endif
- if (r2 > 0.f && r2 < hig2) {
- IACT_STARS(r2, dx, hi, hj, si, pj, a, H);
+ if (r2 < hig2) {
+ IACT_STARS(r2, dx, hi, pj->h, si, pj, a, H);
}
} /* loop over the parts in ci. */
} /* loop over the sparts in ci. */
@@ -142,11 +152,11 @@ void DOSELF1_STARS(struct runner *r, struct cell *c, int timer) {
* @param ci The first #cell
* @param cj The second #cell
*/
-void DO_NONSYM_PAIR1_STARS(struct runner *r, struct cell *restrict ci,
- struct cell *restrict cj) {
+void DO_NONSYM_PAIR1_STARS_NAIVE(struct runner *r, struct cell *restrict ci,
+ struct cell *restrict cj) {
#ifdef SWIFT_DEBUG_CHECKS
-#ifdef UPDATE_STARS
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
if (ci->nodeID != engine_rank) error("Should be run on a different node");
#else
if (cj->nodeID != engine_rank) error("Should be run on a different node");
@@ -157,6 +167,7 @@ void DO_NONSYM_PAIR1_STARS(struct runner *r, struct cell *restrict ci,
const struct cosmology *cosmo = e->cosmology;
/* Anything to do here? */
+ if (cj->hydro.count == 0 || ci->stars.count == 0) return;
if (!cell_is_active_stars(ci, e)) return;
const int scount_i = ci->stars.count;
@@ -182,7 +193,8 @@ void DO_NONSYM_PAIR1_STARS(struct runner *r, struct cell *restrict ci,
/* Get a hold of the ith spart in ci. */
struct spart *restrict si = &sparts_i[sid];
- if (!spart_is_active(si, e) || spart_is_inhibited(si, e)) continue;
+ if (!spart_is_active(si, e)) continue;
+
const float hi = si->h;
const float hig2 = hi * hi * kernel_gamma2;
const float six[3] = {(float)(si->x[0] - (cj->loc[0] + shift[0])),
@@ -194,7 +206,9 @@ void DO_NONSYM_PAIR1_STARS(struct runner *r, struct cell *restrict ci,
/* Get a pointer to the jth particle. */
struct part *restrict pj = &parts_j[pjd];
- const float hj = pj->h;
+
+ /* Skip inhibited particles. */
+ if (part_is_inhibited(pj, e)) continue;
/* Compute the pairwise distance. */
const float pjx[3] = {(float)(pj->x[0] - cj->loc[0]),
@@ -209,27 +223,297 @@ void DO_NONSYM_PAIR1_STARS(struct runner *r, struct cell *restrict ci,
error("Particle pj not drifted to current time");
#endif
- if (r2 < hig2) IACT_STARS(r2, dx, hi, hj, si, pj, a, H);
+ if (r2 < hig2) IACT_STARS(r2, dx, hi, pj->h, si, pj, a, H);
} /* loop over the parts in cj. */
} /* loop over the parts in ci. */
}
-void DOPAIR1_STARS(struct runner *r, struct cell *restrict ci,
- struct cell *restrict cj, int timer) {
+/**
+ * @brief Compute the interactions between a cell pair.
+ *
+ * @param r The #runner.
+ * @param ci The first #cell.
+ * @param cj The second #cell.
+ * @param sid The direction of the pair.
+ * @param shift The shift vector to apply to the particles in ci.
+ */
+void DO_SYM_PAIR1_STARS(struct runner *r, struct cell *ci, struct cell *cj,
+ const int sid, const double *shift) {
+
+ const struct engine *restrict e = r->e;
+ const struct cosmology *restrict cosmo = e->cosmology;
+
+ /* Get the cutoff shift. */
+ double rshift = 0.0;
+ for (int k = 0; k < 3; k++) rshift += shift[k] * runner_shift[sid][k];
+
+ /* Cosmological terms */
+ const float a = cosmo->a;
+ const float H = cosmo->H;
+
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
+ const int do_ci_stars = (ci->nodeID == e->nodeID) && (ci->stars.count != 0) &&
+ (cj->hydro.count != 0) && cell_is_active_stars(ci, e);
+ const int do_cj_stars = (cj->nodeID == e->nodeID) && (cj->stars.count != 0) &&
+ (ci->hydro.count != 0) && cell_is_active_stars(cj, e);
+#else
+ /* here we are updating the hydro -> switch ci, cj for local */
+ const int do_ci_stars = (cj->nodeID == e->nodeID) && (ci->stars.count != 0) &&
+ (cj->hydro.count != 0) && cell_is_active_stars(ci, e);
+ const int do_cj_stars = (ci->nodeID == e->nodeID) && (cj->stars.count != 0) &&
+ (ci->hydro.count != 0) && cell_is_active_stars(cj, e);
+#endif
+
+ if (do_ci_stars) {
+
+ /* Pick-out the sorted lists. */
+ const struct entry *restrict sort_j = cj->hydro.sort[sid];
+ const struct entry *restrict sort_i = ci->stars.sort[sid];
+
+#ifdef SWIFT_DEBUG_CHECKS
+ /* Some constants used to checks that the parts are in the right frame */
+ const float shift_threshold_x =
+ 2. * ci->width[0] +
+ 2. * max(ci->stars.dx_max_part, cj->hydro.dx_max_part);
+ const float shift_threshold_y =
+ 2. * ci->width[1] +
+ 2. * max(ci->stars.dx_max_part, cj->hydro.dx_max_part);
+ const float shift_threshold_z =
+ 2. * ci->width[2] +
+ 2. * max(ci->stars.dx_max_part, cj->hydro.dx_max_part);
+#endif /* SWIFT_DEBUG_CHECKS */
+
+ /* Get some other useful values. */
+ const double hi_max = ci->stars.h_max * kernel_gamma - rshift;
+ const int count_i = ci->stars.count;
+ const int count_j = cj->hydro.count;
+ struct spart *restrict sparts_i = ci->stars.parts;
+ struct part *restrict parts_j = cj->hydro.parts;
+ const double dj_min = sort_j[0].d;
+ const float dx_max_rshift =
+ (ci->stars.dx_max_sort + cj->hydro.dx_max_sort) - rshift;
+ const float dx_max = (ci->stars.dx_max_sort + cj->hydro.dx_max_sort);
+
+ /* Loop over the sparts in ci. */
+ for (int pid = count_i - 1;
+ pid >= 0 && sort_i[pid].d + hi_max + dx_max > dj_min; pid--) {
+
+ /* Get a hold of the ith part in ci. */
+ struct spart *restrict spi = &sparts_i[pid];
+ const float hi = spi->h;
+
+ /* Skip inactive particles */
+ if (!spart_is_active(spi, e)) continue;
+
+ /* Compute distance from the other cell. */
+ const double px[3] = {spi->x[0], spi->x[1], spi->x[2]};
+ float dist = px[0] * runner_shift[sid][0] + px[1] * runner_shift[sid][1] +
+ px[2] * runner_shift[sid][2];
+
+ /* Is there anything we need to interact with ? */
+ const double di = dist + hi * kernel_gamma + dx_max_rshift;
+ if (di < dj_min) continue;
+
+ /* Get some additional information about pi */
+ const float hig2 = hi * hi * kernel_gamma2;
+ const float pix = spi->x[0] - (cj->loc[0] + shift[0]);
+ const float piy = spi->x[1] - (cj->loc[1] + shift[1]);
+ const float piz = spi->x[2] - (cj->loc[2] + shift[2]);
+
+ /* Loop over the parts in cj. */
+ for (int pjd = 0; pjd < count_j && sort_j[pjd].d < di; pjd++) {
+
+ /* Recover pj */
+ struct part *pj = &parts_j[sort_j[pjd].i];
+
+ /* Skip inhibited particles. */
+ if (part_is_inhibited(pj, e)) continue;
+
+ const float hj = pj->h;
+ const float pjx = pj->x[0] - cj->loc[0];
+ const float pjy = pj->x[1] - cj->loc[1];
+ const float pjz = pj->x[2] - cj->loc[2];
+
+ /* Compute the pairwise distance. */
+ 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 are in the correct frame after the shifts */
+ if (pix > shift_threshold_x || pix < -shift_threshold_x)
+ error(
+ "Invalid particle position in X for pi (pix=%e ci->width[0]=%e)",
+ pix, ci->width[0]);
+ if (piy > shift_threshold_y || piy < -shift_threshold_y)
+ error(
+ "Invalid particle position in Y for pi (piy=%e ci->width[1]=%e)",
+ piy, ci->width[1]);
+ if (piz > shift_threshold_z || piz < -shift_threshold_z)
+ error(
+ "Invalid particle position in Z for pi (piz=%e ci->width[2]=%e)",
+ piz, ci->width[2]);
+ if (pjx > shift_threshold_x || pjx < -shift_threshold_x)
+ error(
+ "Invalid particle position in X for pj (pjx=%e ci->width[0]=%e)",
+ pjx, ci->width[0]);
+ if (pjy > shift_threshold_y || pjy < -shift_threshold_y)
+ error(
+ "Invalid particle position in Y for pj (pjy=%e ci->width[1]=%e)",
+ pjy, ci->width[1]);
+ if (pjz > shift_threshold_z || pjz < -shift_threshold_z)
+ error(
+ "Invalid particle position in Z for pj (pjz=%e ci->width[2]=%e)",
+ pjz, ci->width[2]);
+
+ /* Check that particles have been drifted to the current time */
+ if (spi->ti_drift != e->ti_current)
+ error("Particle spi 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) {
+
+ IACT_STARS(r2, dx, hi, hj, spi, pj, a, H);
+ }
+ } /* loop over the parts in cj. */
+ } /* loop over the parts in ci. */
+ } /* do_ci_stars */
+
+ if (do_cj_stars) {
+ /* Pick-out the sorted lists. */
+ const struct entry *restrict sort_i = ci->hydro.sort[sid];
+ const struct entry *restrict sort_j = cj->stars.sort[sid];
+
+#ifdef SWIFT_DEBUG_CHECKS
+ /* Some constants used to checks that the parts are in the right frame */
+ const float shift_threshold_x =
+ 2. * ci->width[0] +
+ 2. * max(ci->hydro.dx_max_part, cj->stars.dx_max_part);
+ const float shift_threshold_y =
+ 2. * ci->width[1] +
+ 2. * max(ci->hydro.dx_max_part, cj->stars.dx_max_part);
+ const float shift_threshold_z =
+ 2. * ci->width[2] +
+ 2. * max(ci->hydro.dx_max_part, cj->stars.dx_max_part);
+#endif /* SWIFT_DEBUG_CHECKS */
+
+ /* Get some other useful values. */
+ const double hj_max = cj->hydro.h_max * kernel_gamma;
+ const int count_i = ci->hydro.count;
+ const int count_j = cj->stars.count;
+ struct part *restrict parts_i = ci->hydro.parts;
+ struct spart *restrict sparts_j = cj->stars.parts;
+ const double di_max = sort_i[count_i - 1].d - rshift;
+ const float dx_max_rshift =
+ (ci->hydro.dx_max_sort + cj->stars.dx_max_sort) + rshift;
+ const float dx_max = (ci->hydro.dx_max_sort + cj->stars.dx_max_sort);
+
+ /* Loop over the parts in cj. */
+ for (int pjd = 0; pjd < count_j && sort_j[pjd].d - hj_max - dx_max < di_max;
+ pjd++) {
+
+ /* Get a hold of the jth part in cj. */
+ struct spart *spj = &sparts_j[pjd];
+ const float hj = spj->h;
+
+ /* Skip inactive particles */
+ if (!spart_is_active(spj, e)) continue;
+
+ /* Compute distance from the other cell. */
+ const double px[3] = {spj->x[0], spj->x[1], spj->x[2]};
+ float dist = px[0] * runner_shift[sid][0] + px[1] * runner_shift[sid][1] +
+ px[2] * runner_shift[sid][2];
+
+ /* Is there anything we need to interact with ? */
+ const double dj = dist - hj * kernel_gamma - dx_max_rshift;
+ if (dj - rshift > di_max) continue;
+
+ /* Get some additional information about pj */
+ const float hjg2 = hj * hj * kernel_gamma2;
+ const float pjx = spj->x[0] - cj->loc[0];
+ const float pjy = spj->x[1] - cj->loc[1];
+ const float pjz = spj->x[2] - cj->loc[2];
+
+ /* Loop over the parts in ci. */
+ for (int pid = count_i - 1; pid >= 0 && sort_i[pid].d > dj; pid--) {
+
+ /* Recover pi */
+ struct part *pi = &parts_i[sort_i[pid].i];
+
+ /* Skip inhibited particles. */
+ if (part_is_inhibited(pi, e)) continue;
+
+ const float hi = pi->h;
+ const float pix = pi->x[0] - (cj->loc[0] + shift[0]);
+ const float piy = pi->x[1] - (cj->loc[1] + shift[1]);
+ const float piz = pi->x[2] - (cj->loc[2] + shift[2]);
+
+ /* Compute the pairwise distance. */
+ float dx[3] = {pjx - pix, pjy - piy, pjz - piz};
+ const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
+
+#ifdef SWIFT_DEBUG_CHECKS
+ /* Check that particles are in the correct frame after the shifts */
+ if (pix > shift_threshold_x || pix < -shift_threshold_x)
+ error(
+ "Invalid particle position in X for pi (pix=%e ci->width[0]=%e)",
+ pix, ci->width[0]);
+ if (piy > shift_threshold_y || piy < -shift_threshold_y)
+ error(
+ "Invalid particle position in Y for pi (piy=%e ci->width[1]=%e)",
+ piy, ci->width[1]);
+ if (piz > shift_threshold_z || piz < -shift_threshold_z)
+ error(
+ "Invalid particle position in Z for pi (piz=%e ci->width[2]=%e)",
+ piz, ci->width[2]);
+ if (pjx > shift_threshold_x || pjx < -shift_threshold_x)
+ error(
+ "Invalid particle position in X for pj (pjx=%e ci->width[0]=%e)",
+ pjx, ci->width[0]);
+ if (pjy > shift_threshold_y || pjy < -shift_threshold_y)
+ error(
+ "Invalid particle position in Y for pj (pjy=%e ci->width[1]=%e)",
+ pjy, ci->width[1]);
+ if (pjz > shift_threshold_z || pjz < -shift_threshold_z)
+ error(
+ "Invalid particle position in Z for pj (pjz=%e ci->width[2]=%e)",
+ pjz, ci->width[2]);
+
+ /* 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 (spj->ti_drift != e->ti_current)
+ error("Particle spj not drifted to current time");
+#endif
+
+ /* Hit or miss? */
+ if (r2 < hjg2) {
+
+ IACT_STARS(r2, dx, hj, hi, spj, pi, a, H);
+ }
+ } /* loop over the parts in ci. */
+ } /* loop over the parts in cj. */
+ } /* Cell cj is active */
+}
+
+void DOPAIR1_STARS_NAIVE(struct runner *r, struct cell *restrict ci,
+ struct cell *restrict cj, int timer) {
-#ifdef UPDATE_STARS
- const int ci_local = ci->nodeID == engine_rank;
- const int cj_local = cj->nodeID == engine_rank;
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
+ const int do_ci_stars = ci->nodeID == r->e->nodeID;
+ const int do_cj_stars = cj->nodeID == r->e->nodeID;
#else
/* here we are updating the hydro -> switch ci, cj */
- const int ci_local = cj->nodeID == engine_rank;
- const int cj_local = ci->nodeID == engine_rank;
+ const int do_ci_stars = cj->nodeID == r->e->nodeID;
+ const int do_cj_stars = ci->nodeID == r->e->nodeID;
#endif
- if (ci_local && ci->stars.count != 0 && cj->hydro.count != 0)
- DO_NONSYM_PAIR1_STARS(r, ci, cj);
- if (cj_local && cj->stars.count != 0 && ci->hydro.count != 0)
- DO_NONSYM_PAIR1_STARS(r, cj, ci);
+ if (do_ci_stars && ci->stars.count != 0 && cj->hydro.count != 0)
+ DO_NONSYM_PAIR1_STARS_NAIVE(r, ci, cj);
+ if (do_cj_stars && cj->stars.count != 0 && ci->hydro.count != 0)
+ DO_NONSYM_PAIR1_STARS_NAIVE(r, cj, ci);
}
/**
@@ -244,12 +528,151 @@ void DOPAIR1_STARS(struct runner *r, struct cell *restrict ci,
* @param ind The list of indices of particles in @c ci to interact with.
* @param scount The number of particles in @c ind.
* @param cj The second #cell.
+ * @param sid The direction of the pair.
+ * @param flipped Flag to check whether the cells have been flipped or not.
* @param shift The shift vector to apply to the particles in ci.
*/
void DOPAIR1_SUBSET_STARS(struct runner *r, struct cell *restrict ci,
struct spart *restrict sparts_i, int *restrict ind,
- int scount, struct cell *restrict cj,
- const double *shift) {
+ int scount, struct cell *restrict cj, const int sid,
+ const int flipped, const double *shift) {
+
+ const struct engine *e = r->e;
+ const struct cosmology *cosmo = e->cosmology;
+
+ const int count_j = cj->hydro.count;
+ struct part *restrict parts_j = cj->hydro.parts;
+
+ /* Early abort? */
+ if (count_j == 0) return;
+
+ /* Cosmological terms */
+ const float a = cosmo->a;
+ const float H = cosmo->H;
+
+ /* Pick-out the sorted lists. */
+ const struct entry *restrict sort_j = cj->hydro.sort[sid];
+ const float dxj = cj->hydro.dx_max_sort;
+
+ /* Sparts are on the left? */
+ if (!flipped) {
+
+ /* Loop over the sparts_i. */
+ for (int pid = 0; pid < scount; pid++) {
+
+ /* Get a hold of the ith spart in ci. */
+ struct spart *restrict spi = &sparts_i[ind[pid]];
+ const double pix = spi->x[0] - (shift[0]);
+ const double piy = spi->x[1] - (shift[1]);
+ const double piz = spi->x[2] - (shift[2]);
+ const float hi = spi->h;
+ const float hig2 = hi * hi * kernel_gamma2;
+ const double di = hi * kernel_gamma + dxj + pix * runner_shift[sid][0] +
+ piy * runner_shift[sid][1] + piz * runner_shift[sid][2];
+
+ /* Loop over the parts in cj. */
+ for (int pjd = 0; pjd < count_j && sort_j[pjd].d < di; pjd++) {
+
+ /* Get a pointer to the jth particle. */
+ struct part *restrict pj = &parts_j[sort_j[pjd].i];
+
+ /* Skip inhibited particles. */
+ if (part_is_inhibited(pj, e)) continue;
+
+ const double pjx = pj->x[0];
+ const double pjy = pj->x[1];
+ const double pjz = pj->x[2];
+
+ /* Compute the pairwise distance. */
+ float dx[3] = {(float)(pix - pjx), (float)(piy - pjy),
+ (float)(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 (spi->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) {
+ IACT_STARS(r2, dx, hi, pj->h, spi, pj, a, H);
+ }
+ } /* loop over the parts in cj. */
+ } /* loop over the sparts in ci. */
+ }
+
+ /* Sparts are on the right. */
+ else {
+
+ /* Loop over the sparts_i. */
+ for (int pid = 0; pid < scount; pid++) {
+
+ /* Get a hold of the ith spart in ci. */
+ struct spart *restrict spi = &sparts_i[ind[pid]];
+ const double pix = spi->x[0] - (shift[0]);
+ const double piy = spi->x[1] - (shift[1]);
+ const double piz = spi->x[2] - (shift[2]);
+ const float hi = spi->h;
+ const float hig2 = hi * hi * kernel_gamma2;
+ const double di = -hi * kernel_gamma - dxj + pix * runner_shift[sid][0] +
+ piy * runner_shift[sid][1] + piz * runner_shift[sid][2];
+
+ /* Loop over the parts in cj. */
+ for (int pjd = count_j - 1; pjd >= 0 && di < sort_j[pjd].d; pjd--) {
+
+ /* Get a pointer to the jth particle. */
+ struct part *restrict pj = &parts_j[sort_j[pjd].i];
+
+ /* Skip inhibited particles. */
+ if (part_is_inhibited(pj, e)) continue;
+
+ const double pjx = pj->x[0];
+ const double pjy = pj->x[1];
+ const double pjz = pj->x[2];
+
+ /* Compute the pairwise distance. */
+ float dx[3] = {(float)(pix - pjx), (float)(piy - pjy),
+ (float)(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 (spi->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) {
+ IACT_STARS(r2, dx, hi, pj->h, spi, pj, a, H);
+ }
+ } /* loop over the parts in cj. */
+ } /* loop over the sparts in ci. */
+ }
+}
+
+/**
+ * @brief Compute the interactions between a cell pair, but only for the
+ * given indices in ci.
+ *
+ * Version using a brute-force algorithm.
+ *
+ * @param r The #runner.
+ * @param ci The first #cell.
+ * @param sparts_i The #part to interact with @c cj.
+ * @param ind The list of indices of particles in @c ci to interact with.
+ * @param scount The number of particles in @c ind.
+ * @param cj The second #cell.
+ * @param shift The shift vector to apply to the particles in ci.
+ */
+void DOPAIR1_SUBSET_STARS_NAIVE(struct runner *r, struct cell *restrict ci,
+ struct spart *restrict sparts_i,
+ int *restrict ind, int scount,
+ struct cell *restrict cj, const double *shift) {
#ifdef SWIFT_DEBUG_CHECKS
if (ci->nodeID != engine_rank) error("Should be run on a different node");
@@ -260,6 +683,9 @@ void DOPAIR1_SUBSET_STARS(struct runner *r, struct cell *restrict ci,
const int count_j = cj->hydro.count;
struct part *restrict parts_j = cj->hydro.parts;
+ /* Early abort? */
+ if (count_j == 0) return;
+
/* Cosmological terms */
const float a = cosmo->a;
const float H = cosmo->H;
@@ -285,6 +711,9 @@ void DOPAIR1_SUBSET_STARS(struct runner *r, struct cell *restrict ci,
/* Get a pointer to the jth particle. */
struct part *restrict pj = &parts_j[pjd];
+ /* Skip inhibited particles */
+ if (part_is_inhibited(pj, e)) continue;
+
/* Compute the pairwise distance. */
float r2 = 0.0f;
float dx[3];
@@ -334,6 +763,9 @@ void DOSELF1_SUBSET_STARS(struct runner *r, struct cell *restrict ci,
const int count_i = ci->hydro.count;
struct part *restrict parts_j = ci->hydro.parts;
+ /* Early abort? */
+ if (count_i == 0) return;
+
/* Loop over the parts in ci. */
for (int spid = 0; spid < scount; spid++) {
@@ -355,7 +787,9 @@ void DOSELF1_SUBSET_STARS(struct runner *r, struct cell *restrict ci,
/* Get a pointer to the jth particle. */
struct part *restrict pj = &parts_j[pjd];
- const float hj = pj->h;
+
+ /* Early abort? */
+ if (part_is_inhibited(pj, e)) continue;
/* Compute the pairwise distance. */
const float pjx[3] = {(float)(pj->x[0] - ci->loc[0]),
@@ -371,8 +805,8 @@ void DOSELF1_SUBSET_STARS(struct runner *r, struct cell *restrict ci,
#endif
/* Hit or miss? */
- if (r2 > 0.f && r2 < hig2) {
- IACT_STARS(r2, dx, hi, hj, spi, pj, a, H);
+ if (r2 < hig2) {
+ IACT_STARS(r2, dx, hi, pj->h, spi, pj, a, H);
}
} /* loop over the parts in cj. */
} /* loop over the parts in ci. */
@@ -414,6 +848,9 @@ void DOPAIR1_SUBSET_BRANCH_STARS(struct runner *r, struct cell *restrict ci,
const struct engine *e = r->e;
+ /* Anything to do here? */
+ if (cj->hydro.count == 0) return;
+
/* Get the relative distance between the pairs, wrapping. */
double shift[3] = {0.0, 0.0, 0.0};
for (int k = 0; k < 3; k++) {
@@ -423,7 +860,27 @@ void DOPAIR1_SUBSET_BRANCH_STARS(struct runner *r, struct cell *restrict ci,
shift[k] = -e->s->dim[k];
}
- DOPAIR1_SUBSET_STARS(r, ci, sparts_i, ind, scount, cj, shift);
+#ifdef SWIFT_USE_NAIVE_INTERACTIONS_STARS
+ DOPAIR1_SUBSET_STARS_NAIVE(r, ci, sparts_i, ind, scount, cj, shift);
+#else
+ /* Get the sorting index. */
+ int sid = 0;
+ for (int k = 0; k < 3; k++)
+ sid = 3 * sid + ((cj->loc[k] - ci->loc[k] + shift[k] < 0)
+ ? 0
+ : (cj->loc[k] - ci->loc[k] + shift[k] > 0) ? 2 : 1);
+
+ /* Switch the cells around? */
+ const int flipped = runner_flip[sid];
+ sid = sortlistID[sid];
+
+ /* Has the cell cj been sorted? */
+ if (!(cj->hydro.sorted & (1 << sid)) ||
+ cj->hydro.dx_max_sort_old > space_maxreldx * cj->dmin)
+ error("Interacting unsorted cells.");
+
+ DOPAIR1_SUBSET_STARS(r, ci, sparts_i, ind, scount, cj, sid, flipped, shift);
+#endif
}
void DOSUB_SUBSET_STARS(struct runner *r, struct cell *ci, struct spart *sparts,
@@ -1075,18 +1532,18 @@ void DOPAIR1_BRANCH_STARS(struct runner *r, struct cell *ci, struct cell *cj) {
const int ci_active = cell_is_active_stars(ci, e);
const int cj_active = cell_is_active_stars(cj, e);
-#ifdef UPDATE_STARS
- const int ci_local = ci->nodeID == engine_rank;
- const int cj_local = cj->nodeID == engine_rank;
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
+ const int do_ci_stars = ci->nodeID == e->nodeID;
+ const int do_cj_stars = cj->nodeID == e->nodeID;
#else
/* here we are updating the hydro -> switch ci, cj */
- const int ci_local = cj->nodeID == engine_rank;
- const int cj_local = ci->nodeID == engine_rank;
+ const int do_ci_stars = cj->nodeID == e->nodeID;
+ const int do_cj_stars = ci->nodeID == e->nodeID;
#endif
- const int do_ci =
- (ci->stars.count != 0 && cj->hydro.count != 0 && ci_active && ci_local);
- const int do_cj =
- (cj->stars.count != 0 && ci->hydro.count != 0 && cj_active && cj_local);
+ const int do_ci = (ci->stars.count != 0 && cj->hydro.count != 0 &&
+ ci_active && do_ci_stars);
+ const int do_cj = (cj->stars.count != 0 && ci->hydro.count != 0 &&
+ cj_active && do_cj_stars);
/* Anything to do here? */
if (!do_ci && !do_cj) return;
@@ -1130,7 +1587,11 @@ void DOPAIR1_BRANCH_STARS(struct runner *r, struct cell *ci, struct cell *cj) {
}
#endif /* SWIFT_DEBUG_CHECKS */
- DOPAIR1_STARS(r, ci, cj, 1);
+#ifdef SWIFT_USE_NAIVE_INTERACTIONS_STARS
+ DOPAIR1_STARS_NAIVE(r, ci, cj, 1);
+#else
+ DO_SYM_PAIR1_STARS(r, ci, cj, sid, shift);
+#endif
}
/**
@@ -1367,18 +1828,18 @@ void DOSUB_PAIR1_STARS(struct runner *r, struct cell *ci, struct cell *cj,
/* Otherwise, compute the pair directly. */
else {
-#ifdef UPDATE_STARS
- const int ci_local = ci->nodeID == engine_rank;
- const int cj_local = cj->nodeID == engine_rank;
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
+ const int do_ci_stars = ci->nodeID == e->nodeID;
+ const int do_cj_stars = cj->nodeID == e->nodeID;
#else
/* here we are updating the hydro -> switch ci, cj */
- const int ci_local = cj->nodeID == engine_rank;
- const int cj_local = ci->nodeID == engine_rank;
+ const int do_ci_stars = cj->nodeID == e->nodeID;
+ const int do_cj_stars = ci->nodeID == e->nodeID;
#endif
const int do_ci = ci->stars.count != 0 && cj->hydro.count != 0 &&
- cell_is_active_stars(ci, e) && ci_local;
+ cell_is_active_stars(ci, e) && do_ci_stars;
const int do_cj = cj->stars.count != 0 && ci->hydro.count != 0 &&
- cell_is_active_stars(cj, e) && cj_local;
+ cell_is_active_stars(cj, e) && do_cj_stars;
if (do_ci) {