diff --git a/src/runner_doiact.h b/src/runner_doiact.h index 473c0604ee60d755f8aef87a96a37bef30ca6279..cd52702b619843ac11ebc666e91214f65ac8624d 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) @@ -171,6 +174,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) { @@ -255,6 +266,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,14 +299,14 @@ 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; @@ -308,10 +327,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 +339,125 @@ 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; + +#ifndef SWIFT_DEBUG_CHECKS + error("Don't use in actual runs ! Slow code !"); +#endif + + 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); } @@ -523,6 +629,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 +676,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 +744,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) {