diff --git a/src/runner_doiact_functions_hydro.h b/src/runner_doiact_functions_hydro.h index 8559d80ee51b54633045a8593ba8144891f85293..79a3ee7a8b84a6880d6efc52c176a1831074a670 100644 --- a/src/runner_doiact_functions_hydro.h +++ b/src/runner_doiact_functions_hydro.h @@ -67,9 +67,15 @@ void DOPAIR1_NAIVE(struct runner *r, const struct cell *restrict ci, if (ci->dmin != cj->dmin) error("Cells of different size!"); #endif + /* Get the depth limits (if any) */ + const char min_depth = limit_max_h ? ci->depth : 0; + const char max_depth = limit_min_h ? ci->depth : CHAR_MAX; + +#ifdef SWIFT_DEBUG_CHECKS /* Get the limits in h (if any) */ const float h_min = limit_min_h ? ci->h_min_allowed : 0.; const float h_max = limit_max_h ? ci->h_max_allowed : FLT_MAX; +#endif /* Get the relative distance between the pairs, wrapping. */ double shift[3] = {0.0, 0.0, 0.0}; @@ -90,6 +96,7 @@ void DOPAIR1_NAIVE(struct runner *r, const struct cell *restrict ci, if (part_is_inhibited(pi, e)) continue; const int pi_active = part_is_active(pi, e); + const char depth_i = pi->depth_h; const float hi = pi->h; const float hig2 = hi * hi * kernel_gamma2; const float pix[3] = {(float)(pi->x[0] - (cj->loc[0] + shift[0])), @@ -105,9 +112,10 @@ void DOPAIR1_NAIVE(struct runner *r, const struct cell *restrict ci, /* Skip inhibited particles. */ if (part_is_inhibited(pj, e)) continue; + const int pj_active = part_is_active(pj, e); + const char depth_j = pj->depth_h; 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. */ const float pjx[3] = {(float)(pj->x[0] - cj->loc[0]), @@ -124,12 +132,18 @@ void DOPAIR1_NAIVE(struct runner *r, const struct cell *restrict ci, error("Particle pj not drifted to current time"); #endif - const int doi = pi_active && (r2 < hig2) && (hi >= h_min) && (hi < h_max); - const int doj = pj_active && (r2 < hjg2) && (hj >= h_min) && (hj < h_max); + const int doi = pi_active && (r2 < hig2) && (depth_i >= min_depth) && + (depth_i <= max_depth); + const int doj = pj_active && (r2 < hjg2) && (depth_j >= min_depth) && + (depth_j <= max_depth); /* Hit or miss? */ if (doi) { +#ifdef SWIFT_DEBUG_CHECKS + if (hi < h_min || hi >= h_max) error("Inappropriate h for this level!"); +#endif + IACT_NONSYM(r2, dx, hi, hj, pi, pj, a, H); #if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY) runner_iact_nonsym_chemistry(r2, dx, hi, hj, pi, pj, a, H); @@ -144,6 +158,10 @@ void DOPAIR1_NAIVE(struct runner *r, const struct cell *restrict ci, } if (doj) { +#ifdef SWIFT_DEBUG_CHECKS + if (hj < h_min || hj >= h_max) error("Inappropriate h for this level!"); +#endif + dx[0] = -dx[0]; dx[1] = -dx[1]; dx[2] = -dx[2]; @@ -337,9 +355,15 @@ void DOSELF1_NAIVE(struct runner *r, const struct cell *c, const int count = c->hydro.count; struct part *parts = c->hydro.parts; + /* Get the depth limits (if any) */ + const char min_depth = limit_max_h ? c->depth : 0; + const char max_depth = limit_min_h ? c->depth : CHAR_MAX; + +#ifdef SWIFT_DEBUG_CHECKS /* Get the limits in h (if any) */ - const float h_min = limit_min_h ? c->dmin * 0.5 * (1. / kernel_gamma) : 0.; - const float h_max = limit_max_h ? c->dmin * (1. / kernel_gamma) : FLT_MAX; + const float h_min = limit_min_h ? c->h_min_allowed : 0.; + const float h_max = limit_max_h ? c->h_max_allowed : FLT_MAX; +#endif /* Loop over the parts in ci. */ for (int pid = 0; pid < count; pid++) { @@ -351,6 +375,7 @@ void DOSELF1_NAIVE(struct runner *r, const struct cell *c, if (part_is_inhibited(pi, e)) continue; const int pi_active = part_is_active(pi, e); + const char depth_i = pi->depth_h; const float hi = pi->h; const float hig2 = hi * hi * kernel_gamma2; const float pix[3] = {(float)(pi->x[0] - c->loc[0]), @@ -369,6 +394,7 @@ void DOSELF1_NAIVE(struct runner *r, const struct cell *c, const float hj = pj->h; const float hjg2 = hj * hj * kernel_gamma2; const int pj_active = part_is_active(pj, e); + const char depth_j = pj->depth_h; /* Compute the pairwise distance. */ const float pjx[3] = {(float)(pj->x[0] - c->loc[0]), @@ -377,8 +403,10 @@ void DOSELF1_NAIVE(struct runner *r, const struct cell *c, 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) && (hi >= h_min) && (hi < h_max); - const int doj = pj_active && (r2 < hjg2) && (hj >= h_min) && (hj < h_max); + const int doi = pi_active && (r2 < hig2) && (depth_i >= min_depth) && + (depth_i <= max_depth); + const int doj = pj_active && (r2 < hjg2) && (depth_j >= min_depth) && + (depth_j <= max_depth); #ifdef SWIFT_DEBUG_CHECKS /* Check that particles have been drifted to the current time */ @@ -391,6 +419,11 @@ void DOSELF1_NAIVE(struct runner *r, const struct cell *c, /* Hit or miss? */ if (doi && doj) { +#ifdef SWIFT_DEBUG_CHECKS + if (hi < h_min || hi >= h_max) error("Inappropriate h for this level!"); + if (hj < h_min || hj >= h_max) error("Inappropriate h for this level!"); +#endif + IACT(r2, dx, hi, hj, pi, pj, a, H); #if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY) runner_iact_chemistry(r2, dx, hi, hj, pi, pj, a, H); @@ -404,6 +437,10 @@ void DOSELF1_NAIVE(struct runner *r, const struct cell *c, #endif } else if (doi) { +#ifdef SWIFT_DEBUG_CHECKS + if (hi < h_min || hi >= h_max) error("Inappropriate h for this level!"); +#endif + IACT_NONSYM(r2, dx, hi, hj, pi, pj, a, H); #if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY) runner_iact_nonsym_chemistry(r2, dx, hi, hj, pi, pj, a, H); @@ -417,6 +454,10 @@ void DOSELF1_NAIVE(struct runner *r, const struct cell *c, #endif } else if (doj) { +#ifdef SWIFT_DEBUG_CHECKS + if (hj < h_min || hj >= h_max) error("Inappropriate h for this level!"); +#endif + dx[0] = -dx[0]; dx[1] = -dx[1]; dx[2] = -dx[2]; @@ -473,8 +514,8 @@ void DOSELF2_NAIVE(struct runner *r, const struct cell *c, struct part *parts = c->hydro.parts; /* Get the limits in h (if any) */ - const float h_min = limit_min_h ? c->dmin * 0.5 * (1. / kernel_gamma) : 0.; - const float h_max = limit_max_h ? c->dmin * (1. / kernel_gamma) : FLT_MAX; + const float h_min = limit_min_h ? c->h_min_allowed : 0.; + const float h_max = limit_max_h ? c->h_max_allowed : FLT_MAX; /* Loop over the parts in ci. */ for (int pid = 0; pid < count; pid++) { @@ -1983,8 +2024,8 @@ void DOSELF1(struct runner *r, const struct cell *c, const int limit_min_h, const int count = c->hydro.count; /* Get the limits in h (if any) */ - const float h_min = limit_min_h ? c->dmin * 0.5 * (1. / kernel_gamma) : 0.; - const float h_max = limit_max_h ? c->dmin * (1. / kernel_gamma) : FLT_MAX; + const float h_min = limit_min_h ? c->h_min_allowed : 0.; + const float h_max = limit_max_h ? c->h_max_allowed : FLT_MAX; /* Set up a list of the particles for which we want to compute interactions */ int *indt = NULL; @@ -2249,8 +2290,8 @@ void DOSELF2(struct runner *r, const struct cell *c, const int limit_min_h, const int count = c->hydro.count; /* Get the limits in h (if any) */ - const float h_min = limit_min_h ? c->dmin * 0.5 * (1. / kernel_gamma) : 0.; - const float h_max = limit_max_h ? c->dmin * (1. / kernel_gamma) : FLT_MAX; + const float h_min = limit_min_h ? c->h_min_allowed : 0.; + const float h_max = limit_max_h ? c->h_max_allowed : FLT_MAX; /* Set up a list of the particles for which we want to compute interactions */ int *indt = NULL;