diff --git a/src/Makefile.am b/src/Makefile.am index 30afe6cf16960fe5b13f7a0bf6cb0cae86eaa43b..8b8dce69877b1c561d96b6d5720fdab1d068436b 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -60,7 +60,7 @@ AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \ nobase_noinst_HEADERS = align.h approx_math.h atomic.h cycle.h error.h inline.h kernel_hydro.h kernel_gravity.h \ kernel_long_gravity.h vector.h runner_doiact.h runner_doiact_grav.h runner_doiact_fft.h \ units.h intrinsics.h minmax.h kick.h timestep.h drift.h adiabatic_index.h io_properties.h \ - dimension.h equation_of_state.h \ + dimension.h equation_of_state.h active.h \ gravity.h gravity_io.h \ gravity/Default/gravity.h gravity/Default/gravity_iact.h gravity/Default/gravity_io.h \ gravity/Default/gravity_debug.h gravity/Default/gravity_part.h \ diff --git a/src/active.h b/src/active.h new file mode 100644 index 0000000000000000000000000000000000000000..0ba79f1dec2e37c45c268b24f4c938866500d5f7 --- /dev/null +++ b/src/active.h @@ -0,0 +1,78 @@ +/******************************************************************************* + * This file is part of SWIFT. + * Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published + * by the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + * + ******************************************************************************/ +#ifndef SWIFT_ACTIVE_H +#define SWIFT_ACTIVE_H + +/* Config parameters. */ +#include "../config.h" + +/* Local includes. */ +#include "cell.h" +#include "engine.h" +#include "part.h" + +/** + * @brief Does a cell contain any active particle ? + * + * @param c The #cell. + * @param e The #engine containing information about the current time. + */ +__attribute__((always_inline)) INLINE static int cell_is_active( + const struct cell *c, const struct engine *e) { + + return (c->ti_end_min <= e->ti_current); +} + +/** + * @brief Are *all* particles in a cell active ? + * + * @param c The #cell. + * @param e The #engine containing information about the current time. + */ +__attribute__((always_inline)) INLINE static int cell_is_all_active( + const struct cell *c, const struct engine *e) { + + return (c->ti_end_max <= e->ti_current); +} + +/** + * @brief Is this particle active ? + * + * @param p The #part. + * @param e The #engine containing information about the current time. + */ +__attribute__((always_inline)) INLINE static int part_is_active( + const struct part *p, const struct engine *e) { + + return (p->ti_end <= e->ti_current); +} + +/** + * @brief Is this g-particle active ? + * + * @param gp The #gpart. + * @param e The #engine containing information about the current time. + */ +__attribute__((always_inline)) INLINE static int gpart_is_active( + const struct gpart *gp, const struct engine *e) { + + return (gp->ti_end <= e->ti_current); +} + +#endif /* SWIFT_ACTIVE_H */ diff --git a/src/hydro/Default/hydro.h b/src/hydro/Default/hydro.h index ccdd0cee32b9386eff54da655b75285b8e08a598..3fd357a2d8778f5ca8b014935d538350eccb99c6 100644 --- a/src/hydro/Default/hydro.h +++ b/src/hydro/Default/hydro.h @@ -199,10 +199,9 @@ __attribute__((always_inline)) INLINE static void hydro_init_part( * and add the self-contribution term. * * @param p The particle to act upon - * @param time The current time */ __attribute__((always_inline)) INLINE static void hydro_end_density( - struct part *restrict p, float time) { + struct part *restrict p) { /* Some smoothing length multiples. */ const float h = p->h; diff --git a/src/hydro/Gadget2/hydro.h b/src/hydro/Gadget2/hydro.h index c999e20d401570ac6518291df8cf315569fe78bd..157893bc9e27806d2b97ac5f5a81d0f6fbb1c589 100644 --- a/src/hydro/Gadget2/hydro.h +++ b/src/hydro/Gadget2/hydro.h @@ -246,10 +246,9 @@ __attribute__((always_inline)) INLINE static void hydro_init_part( * and add the self-contribution term. * * @param p The particle to act upon - * @param ti_current The current time (on the integer timeline) */ __attribute__((always_inline)) INLINE static void hydro_end_density( - struct part *restrict p, int ti_current) { + struct part *restrict p) { /* Some smoothing length multiples. */ const float h = p->h; diff --git a/src/hydro/Gizmo/hydro.h b/src/hydro/Gizmo/hydro.h index bd970795bdf070a7bd7915cc4f493218dbf319d1..1c64291ee64dd770b1f1a76371f67a34230365c7 100644 --- a/src/hydro/Gizmo/hydro.h +++ b/src/hydro/Gizmo/hydro.h @@ -109,10 +109,9 @@ __attribute__((always_inline)) INLINE static void hydro_init_part( * passive particles. * * @param p The particle to act upon. - * @param The current physical time. */ __attribute__((always_inline)) INLINE static void hydro_end_density( - struct part* restrict p, float time) { + struct part* restrict p) { /* Some smoothing length multiples. */ const float h = p->h; diff --git a/src/hydro/Minimal/hydro.h b/src/hydro/Minimal/hydro.h index 3015f26c6bad7f006e8cda16695c750ff40d74df..3b3454f1bb348b178ac57899da4f7611802a69cd 100644 --- a/src/hydro/Minimal/hydro.h +++ b/src/hydro/Minimal/hydro.h @@ -256,10 +256,9 @@ __attribute__((always_inline)) INLINE static void hydro_init_part( * added to them here. * * @param p The particle to act upon - * @param time The current time */ __attribute__((always_inline)) INLINE static void hydro_end_density( - struct part *restrict p, float time) { + struct part *restrict p) { /* Some smoothing length multiples. */ const float h = p->h; diff --git a/src/hydro/PressureEntropy/hydro.h b/src/hydro/PressureEntropy/hydro.h index 0c69728a9ce6317a8d7ddab945e7aab6e94ee247..8c063596efd3be97ebb4da6b6879ac06122bd357 100644 --- a/src/hydro/PressureEntropy/hydro.h +++ b/src/hydro/PressureEntropy/hydro.h @@ -254,10 +254,9 @@ __attribute__((always_inline)) INLINE static void hydro_init_part( * and add the self-contribution term. * * @param p The particle to act upon - * @param ti_current The current time (on the integer timeline) */ __attribute__((always_inline)) INLINE static void hydro_end_density( - struct part *restrict p, int ti_current) { + struct part *restrict p) { /* Some smoothing length multiples. */ const float h = p->h; diff --git a/src/runner.c b/src/runner.c index 472f11c067d3e83018d58d200cff2ca0271eb3d2..1d8335aaf09de770fb96b781acdb5fbbdc3dcf08 100644 --- a/src/runner.c +++ b/src/runner.c @@ -38,6 +38,7 @@ #include "runner.h" /* Local headers. */ +#include "active.h" #include "approx_math.h" #include "atomic.h" #include "cell.h" @@ -159,15 +160,15 @@ void runner_do_grav_external(struct runner *r, struct cell *c, int timer) { struct gpart *restrict gparts = c->gparts; const int gcount = c->gcount; - const int ti_current = r->e->ti_current; - const struct external_potential *potential = r->e->external_potential; - const struct phys_const *constants = r->e->physical_constants; + const struct engine *e = r->e; + const struct external_potential *potential = e->external_potential; + const struct phys_const *constants = e->physical_constants; const double time = r->e->time; TIMER_TIC; /* Anything to do here? */ - if (c->ti_end_min > ti_current) return; + if (!cell_is_active(c, e)) return; /* Recurse? */ if (c->split) { @@ -187,7 +188,7 @@ void runner_do_grav_external(struct runner *r, struct cell *c, int timer) { struct gpart *restrict gp = &gparts[i]; /* Is this part within the time step? */ - if (gp->ti_end <= ti_current) { + if (gpart_is_active(gp, e)) { external_gravity_acceleration(time, potential, constants, gp); } @@ -492,12 +493,12 @@ void runner_do_init(struct runner *r, struct cell *c, int timer) { struct gpart *restrict gparts = c->gparts; const int count = c->count; const int gcount = c->gcount; - const int ti_current = r->e->ti_current; + const struct engine *e = r->e; TIMER_TIC; /* Anything to do here? */ - if (c->ti_end_min > ti_current) return; + if (!cell_is_active(c, e)) return; /* Recurse? */ if (c->split) { @@ -512,7 +513,7 @@ void runner_do_init(struct runner *r, struct cell *c, int timer) { /* Get a direct pointer on the part. */ struct part *restrict p = &parts[i]; - if (p->ti_end <= ti_current) { + if (part_is_active(p, e)) { /* Get ready for a density calculation */ hydro_init_part(p); @@ -525,7 +526,7 @@ void runner_do_init(struct runner *r, struct cell *c, int timer) { /* Get a direct pointer on the part. */ struct gpart *restrict gp = &gparts[i]; - if (gp->ti_end <= ti_current) { + if (gpart_is_active(gp, e)) { /* Get ready for a density calculation */ gravity_init_gpart(gp); @@ -549,10 +550,10 @@ void runner_do_extra_ghost(struct runner *r, struct cell *c) { struct part *restrict parts = c->parts; const int count = c->count; - const int ti_current = r->e->ti_current; + const struct engine *e = r->e; /* Anything to do here? */ - if (c->ti_end_min > ti_current) return; + if (!cell_is_active(c, e)) return; /* Recurse? */ if (c->split) { @@ -567,7 +568,7 @@ void runner_do_extra_ghost(struct runner *r, struct cell *c) { /* Get a direct pointer on the part. */ struct part *restrict p = &parts[i]; - if (p->ti_end <= ti_current) { + if (part_is_active(p, e)) { /* Get ready for a force calculation */ hydro_end_gradient(p); @@ -592,20 +593,20 @@ void runner_do_ghost(struct runner *r, struct cell *c) { struct part *restrict parts = c->parts; struct xpart *restrict xparts = c->xparts; int redo, count = c->count; - const int ti_current = r->e->ti_current; - const double timeBase = r->e->timeBase; - const float target_wcount = r->e->hydro_properties->target_neighbours; + const struct engine *e = r->e; + const int ti_current = e->ti_current; + const double timeBase = e->timeBase; + const float target_wcount = e->hydro_properties->target_neighbours; const float max_wcount = - target_wcount + r->e->hydro_properties->delta_neighbours; + target_wcount + e->hydro_properties->delta_neighbours; const float min_wcount = - target_wcount - r->e->hydro_properties->delta_neighbours; - const int max_smoothing_iter = - r->e->hydro_properties->max_smoothing_iterations; + target_wcount - e->hydro_properties->delta_neighbours; + const int max_smoothing_iter = e->hydro_properties->max_smoothing_iterations; TIMER_TIC; /* Anything to do here? */ - if (c->ti_end_min > ti_current) return; + if (!cell_is_active(c, e)) return; /* Recurse? */ if (c->split) { @@ -635,10 +636,10 @@ void runner_do_ghost(struct runner *r, struct cell *c) { struct xpart *restrict xp = &xparts[pid[i]]; /* Is this part within the timestep? */ - if (p->ti_end <= ti_current) { + if (part_is_active(p, e)) { /* Finish the density calculation */ - hydro_end_density(p, ti_current); + hydro_end_density(p); float h_corr = 0.f; @@ -901,18 +902,17 @@ void runner_do_kick(struct runner *r, struct cell *c, int timer) { const struct engine *e = r->e; const double timeBase = e->timeBase; - const int ti_current = r->e->ti_current; const int count = c->count; const int gcount = c->gcount; struct part *restrict parts = c->parts; struct xpart *restrict xparts = c->xparts; struct gpart *restrict gparts = c->gparts; - const double const_G = r->e->physical_constants->const_newton_G; + const double const_G = e->physical_constants->const_newton_G; TIMER_TIC; /* Anything to do here? */ - if (c->ti_end_min > ti_current) { + if (!cell_is_active(c, e)) { c->updated = 0; c->g_updated = 0; return; @@ -937,7 +937,7 @@ void runner_do_kick(struct runner *r, struct cell *c, int timer) { /* If the g-particle has no counterpart and needs to be kicked */ if (gp->id_or_neg_offset > 0) { - if (gp->ti_end <= ti_current) { + if (gpart_is_active(gp, e)) { /* First, finish the force calculation */ gravity_end_force(gp, const_G); @@ -968,7 +968,7 @@ void runner_do_kick(struct runner *r, struct cell *c, int timer) { struct xpart *restrict xp = &xparts[k]; /* If particle needs to be kicked */ - if (p->ti_end <= ti_current) { + if (part_is_active(p, e)) { /* First, finish the force loop */ hydro_end_force(p); @@ -1126,7 +1126,7 @@ void *runner_main(void *data) { /* Check that we haven't scheduled an inactive task */ #ifdef SWIFT_DEBUG_CHECKS if (cj == NULL) { /* self */ - if (ci->ti_end_min > e->ti_current && t->type != task_type_sort) + if (!cell_is_active(ci, e) && t->type != task_type_sort) error( "Task (type='%s/%s') should have been skipped ti_current=%d " "c->ti_end_min=%d", @@ -1134,7 +1134,7 @@ void *runner_main(void *data) { ci->ti_end_min); /* Special case for sorts */ - if (ci->ti_end_min > e->ti_current && t->type == task_type_sort && + if (!cell_is_active(ci, e) && t->type == task_type_sort && t->flags == 0) error( "Task (type='%s/%s') should have been skipped ti_current=%d " @@ -1143,7 +1143,7 @@ void *runner_main(void *data) { ci->ti_end_min, t->flags); } else { /* pair */ - if (ci->ti_end_min > e->ti_current && cj->ti_end_min > e->ti_current) + if (!cell_is_active(ci, e) && !cell_is_active(cj, e)) error( "Task (type='%s/%s') should have been skipped ti_current=%d " "ci->ti_end_min=%d cj->ti_end_min=%d", diff --git a/src/runner_doiact.h b/src/runner_doiact.h index 3c968cbf7d955198ad6bb44ab70e93af17735e99..969afe33292f874cb3cd511496fd50d225c947df 100644 --- a/src/runner_doiact.h +++ b/src/runner_doiact.h @@ -109,7 +109,6 @@ void DOPAIR_NAIVE(struct runner *r, struct cell *restrict ci, struct cell *restrict cj) { const struct engine *e = r->e; - const int ti_current = e->ti_current; error("Don't use in actual runs ! Slow code !"); @@ -124,7 +123,7 @@ void DOPAIR_NAIVE(struct runner *r, struct cell *restrict ci, TIMER_TIC; /* Anything to do here? */ - if (ci->ti_end_min > ti_current && cj->ti_end_min > ti_current) return; + if (!cell_is_active(ci, e) && !cell_is_active(cj, e)) return; const int count_i = ci->count; const int count_j = cj->count; @@ -210,7 +209,7 @@ void DOPAIR_NAIVE(struct runner *r, struct cell *restrict ci, void DOSELF_NAIVE(struct runner *r, struct cell *restrict c) { - const int ti_current = r->e->ti_current; + const struct engine *e = r->e; error("Don't use in actual runs ! Slow code !"); @@ -226,7 +225,7 @@ void DOSELF_NAIVE(struct runner *r, struct cell *restrict c) { TIMER_TIC; /* Anything to do here? */ - if (c->ti_end_min > ti_current) return; + if (!cell_is_active(c, e)) return; const int count = c->count; struct part *restrict parts = c->parts; @@ -706,8 +705,7 @@ void DOSELF_SUBSET(struct runner *r, struct cell *restrict ci, */ void DOPAIR1(struct runner *r, struct cell *ci, struct cell *cj) { - struct engine *restrict e = r->e; - const int ti_current = e->ti_current; + const struct engine *restrict e = r->e; #ifdef WITH_VECTORIZATION int icount = 0; @@ -721,7 +719,7 @@ void DOPAIR1(struct runner *r, struct cell *ci, struct cell *cj) { TIMER_TIC; /* Anything to do here? */ - if (ci->ti_end_min > ti_current && cj->ti_end_min > ti_current) return; + if (!cell_is_active(ci, e) && !cell_is_active(cj, e)) return; /* Get the sort ID. */ double shift[3] = {0.0, 0.0, 0.0}; @@ -756,7 +754,7 @@ void DOPAIR1(struct runner *r, struct cell *ci, struct cell *cj) { /* Get a hold of the ith part in ci. */ struct part *restrict pi = &parts_i[sort_i[pid].i]; - if (pi->ti_end > ti_current) continue; + if (!part_is_active(pi, e)) continue; const float hi = pi->h; const double di = sort_i[pid].d + hi * kernel_gamma + dx_max - rshift; if (di < dj_min) continue; @@ -818,7 +816,7 @@ void DOPAIR1(struct runner *r, struct cell *ci, struct cell *cj) { /* Get a hold of the jth part in cj. */ struct part *restrict pj = &parts_j[sort_j[pjd].i]; - if (pj->ti_end > ti_current) continue; + if (!part_is_active(pj, e)) continue; const float hj = pj->h; const double dj = sort_j[pjd].d - hj * kernel_gamma - dx_max - rshift; if (dj > di_max) continue; @@ -894,7 +892,6 @@ void DOPAIR1(struct runner *r, struct cell *ci, struct cell *cj) { void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj) { struct engine *restrict e = r->e; - const int ti_current = e->ti_current; #ifdef WITH_VECTORIZATION int icount1 = 0; @@ -914,7 +911,7 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj) { TIMER_TIC; /* Anything to do here? */ - if (ci->ti_end_min > ti_current && cj->ti_end_min > ti_current) return; + if (!cell_is_active(ci, e) && !cell_is_active(cj, e)) return; /* Get the shift ID. */ double shift[3] = {0.0, 0.0, 0.0}; @@ -946,28 +943,28 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj) { /* Collect the number of parts left and right below dt. */ int countdt_i = 0, countdt_j = 0; struct entry *restrict sortdt_i = NULL, *restrict sortdt_j = NULL; - if (ci->ti_end_max <= ti_current) { + if (cell_is_all_active(ci, e)) { sortdt_i = sort_i; countdt_i = count_i; - } else if (ci->ti_end_min <= ti_current) { + } else if (cell_is_active(ci, e)) { if (posix_memalign((void *)&sortdt_i, VEC_SIZE * sizeof(float), sizeof(struct entry) * count_i) != 0) error("Failed to allocate dt sortlists."); for (int k = 0; k < count_i; k++) - if (parts_i[sort_i[k].i].ti_end <= ti_current) { + if (part_is_active(&parts_i[sort_i[k].i], e)) { sortdt_i[countdt_i] = sort_i[k]; countdt_i += 1; } } - if (cj->ti_end_max <= ti_current) { + if (cell_is_all_active(cj, e)) { sortdt_j = sort_j; countdt_j = count_j; - } else if (cj->ti_end_min <= ti_current) { + } else if (cell_is_active(cj, e)) { if (posix_memalign((void *)&sortdt_j, VEC_SIZE * sizeof(float), sizeof(struct entry) * count_j) != 0) error("Failed to allocate dt sortlists."); for (int k = 0; k < count_j; k++) - if (parts_j[sort_j[k].i].ti_end <= ti_current) { + if (part_is_active(&parts_j[sort_j[k].i], e)) { sortdt_j[countdt_j] = sort_j[k]; countdt_j += 1; } @@ -988,7 +985,7 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj) { const float hig2 = hi * hi * kernel_gamma2; /* Look at valid dt parts only? */ - if (pi->ti_end > ti_current) { + if (!part_is_active(pi, e)) { /* Loop over the parts in cj within dt. */ for (int pjd = 0; pjd < countdt_j && sortdt_j[pjd].d < di; pjd++) { @@ -1062,7 +1059,7 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj) { #ifndef WITH_VECTORIZATION /* Does pj need to be updated too? */ - if (pj->ti_end <= ti_current) + if (part_is_active(pj, e)) IACT(r2, dx, hi, hj, pi, pj); else IACT_NONSYM(r2, dx, hi, hj, pi, pj); @@ -1070,7 +1067,7 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj) { #else /* Does pj need to be updated too? */ - if (pj->ti_end <= ti_current) { + if (part_is_active(pj, e)) { /* Add this interaction to the symmetric queue. */ r2q2[icount2] = r2; @@ -1132,7 +1129,7 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj) { const float hjg2 = hj * hj * kernel_gamma2; /* Is this particle outside the dt? */ - if (pj->ti_end > ti_current) { + if (!part_is_active(pj, e)) { /* Loop over the parts in ci. */ for (int pid = countdt_i - 1; pid >= 0 && sortdt_i[pid].d > dj; pid--) { @@ -1205,7 +1202,7 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj) { #ifndef WITH_VECTORIZATION /* Does pi need to be updated too? */ - if (pi->ti_end <= ti_current) + if (part_is_active(pi, e)) IACT(r2, dx, hj, hi, pj, pi); else IACT_NONSYM(r2, dx, hj, hi, pj, pi); @@ -1213,7 +1210,7 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj) { #else /* Does pi need to be updated too? */ - if (pi->ti_end <= ti_current) { + if (part_is_active(pi, e)) { /* Add this interaction to the symmetric queue. */ r2q2[icount2] = r2; @@ -1270,10 +1267,9 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj) { IACT(r2q2[k], &dxq2[3 * k], hiq2[k], hjq2[k], piq2[k], pjq2[k]); #endif - if (ci->ti_end_max > ti_current && ci->ti_end_min <= ti_current) - free(sortdt_i); - if (cj->ti_end_max > ti_current && cj->ti_end_min <= ti_current) - free(sortdt_j); + /* Clean-up if necessary */ + if (cell_is_active(ci, e) && !cell_is_all_active(ci, e)) free(sortdt_i); + if (cell_is_active(cj, e) && !cell_is_all_active(cj, e)) free(sortdt_j); TIMER_TOC(TIMER_DOPAIR); } @@ -1286,7 +1282,7 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj) { */ void DOSELF1(struct runner *r, struct cell *restrict c) { - const int ti_current = r->e->ti_current; + const struct engine *e = r->e; #ifdef WITH_VECTORIZATION int icount1 = 0; @@ -1305,8 +1301,7 @@ void DOSELF1(struct runner *r, struct cell *restrict c) { TIMER_TIC; - if (c->ti_end_min > ti_current) return; - if (c->ti_end_max < ti_current) error("Cell in an impossible time-zone"); + if (!cell_is_active(c, e)) return; struct part *restrict parts = c->parts; const int count = c->count; @@ -1318,7 +1313,7 @@ void DOSELF1(struct runner *r, struct cell *restrict c) { count * sizeof(int)) != 0) error("Failed to allocate indt."); for (int k = 0; k < count; k++) - if (parts[k].ti_end <= ti_current) { + if (part_is_active(&parts[k], e)) { indt[countdt] = k; countdt += 1; } @@ -1336,7 +1331,7 @@ void DOSELF1(struct runner *r, struct cell *restrict c) { const float hig2 = hi * hi * kernel_gamma2; /* Is the ith particle inactive? */ - if (pi->ti_end > ti_current) { + if (!part_is_active(pi, e)) { /* Loop over the other particles .*/ for (int pjd = firstdt; pjd < countdt; pjd++) { @@ -1407,7 +1402,7 @@ void DOSELF1(struct runner *r, struct cell *restrict c) { r2 += dx[k] * dx[k]; } const int doj = - (pj->ti_end <= ti_current) && (r2 < hj * hj * kernel_gamma2); + (part_is_active(pj, e)) && (r2 < hj * hj * kernel_gamma2); /* Hit or miss? */ if (r2 < hig2 || doj) { @@ -1518,7 +1513,7 @@ void DOSELF1(struct runner *r, struct cell *restrict c) { */ void DOSELF2(struct runner *r, struct cell *restrict c) { - const int ti_current = r->e->ti_current; + const struct engine *e = r->e; #ifdef WITH_VECTORIZATION int icount1 = 0; @@ -1537,8 +1532,7 @@ void DOSELF2(struct runner *r, struct cell *restrict c) { TIMER_TIC; - if (c->ti_end_min > ti_current) return; - if (c->ti_end_max < ti_current) error("Cell in an impossible time-zone"); + if (!cell_is_active(c, e)) return; struct part *restrict parts = c->parts; const int count = c->count; @@ -1550,7 +1544,7 @@ void DOSELF2(struct runner *r, struct cell *restrict c) { count * sizeof(int)) != 0) error("Failed to allocate indt."); for (int k = 0; k < count; k++) - if (parts[k].ti_end <= ti_current) { + if (part_is_active(&parts[k], e)) { indt[countdt] = k; countdt += 1; } @@ -1568,7 +1562,7 @@ void DOSELF2(struct runner *r, struct cell *restrict c) { const float hig2 = hi * hi * kernel_gamma2; /* Is the ith particle not active? */ - if (pi->ti_end > ti_current) { + if (!part_is_active(pi, e)) { /* Loop over the other particles .*/ for (int pjd = firstdt; pjd < countdt; pjd++) { @@ -1645,7 +1639,7 @@ void DOSELF2(struct runner *r, struct cell *restrict c) { #ifndef WITH_VECTORIZATION /* Does pj need to be updated too? */ - if (pj->ti_end <= ti_current) + if (part_is_active(pj, e)) IACT(r2, dx, hi, hj, pi, pj); else IACT_NONSYM(r2, dx, hi, hj, pi, pj); @@ -1653,7 +1647,7 @@ void DOSELF2(struct runner *r, struct cell *restrict c) { #else /* Does pj need to be updated too? */ - if (pj->ti_end <= ti_current) { + if (part_is_active(pj, e)) { /* Add this interaction to the symmetric queue. */ r2q2[icount2] = r2; @@ -1731,12 +1725,12 @@ void DOSUB_PAIR1(struct runner *r, struct cell *ci, struct cell *cj, int sid, int gettimer) { struct space *s = r->e->s; - const int ti_current = r->e->ti_current; + const struct engine *e = r->e; TIMER_TIC; /* Should we even bother? */ - if (ci->ti_end_min > ti_current && cj->ti_end_min > ti_current) return; + if (!cell_is_active(ci, e) && !cell_is_active(cj, e)) return; /* Get the cell dimensions. */ const float h = min(ci->width[0], min(ci->width[1], ci->width[2])); @@ -1950,7 +1944,7 @@ void DOSUB_PAIR1(struct runner *r, struct cell *ci, struct cell *cj, int sid, } /* Otherwise, compute the pair directly. */ - else if (ci->ti_end_min <= ti_current || cj->ti_end_min <= ti_current) { + else if (cell_is_active(ci, e) || cell_is_active(cj, e)) { /* Do any of the cells need to be sorted first? */ if (!(ci->sorted & (1 << sid))) runner_do_sort(r, ci, (1 << sid), 1); @@ -1972,12 +1966,10 @@ void DOSUB_PAIR1(struct runner *r, struct cell *ci, struct cell *cj, int sid, */ void DOSUB_SELF1(struct runner *r, struct cell *ci, int gettimer) { - const int ti_current = r->e->ti_current; - TIMER_TIC; /* Should we even bother? */ - if (ci->ti_end_min > ti_current) return; + if (!cell_is_active(ci, r->e)) return; /* Recurse? */ if (ci->split) { @@ -2014,13 +2006,13 @@ void DOSUB_SELF1(struct runner *r, struct cell *ci, int gettimer) { void DOSUB_PAIR2(struct runner *r, struct cell *ci, struct cell *cj, int sid, int gettimer) { - struct space *s = r->e->s; - const int ti_current = r->e->ti_current; + const struct engine *e = r->e; + struct space *s = e->s; TIMER_TIC; /* Should we even bother? */ - if (ci->ti_end_min > ti_current && cj->ti_end_min > ti_current) return; + if (!cell_is_active(ci, e) && cell_is_active(cj, e)) return; /* Get the cell dimensions. */ const float h = min(ci->width[0], min(ci->width[1], ci->width[2])); @@ -2234,7 +2226,7 @@ void DOSUB_PAIR2(struct runner *r, struct cell *ci, struct cell *cj, int sid, } /* Otherwise, compute the pair directly. */ - else if (ci->ti_end_min <= ti_current || cj->ti_end_min <= ti_current) { + else if (cell_is_active(ci, e) || cell_is_active(cj, e)) { /* Do any of the cells need to be sorted first? */ if (!(ci->sorted & (1 << sid))) runner_do_sort(r, ci, (1 << sid), 1); @@ -2256,12 +2248,10 @@ void DOSUB_PAIR2(struct runner *r, struct cell *ci, struct cell *cj, int sid, */ void DOSUB_SELF2(struct runner *r, struct cell *ci, int gettimer) { - const int ti_current = r->e->ti_current; - TIMER_TIC; /* Should we even bother? */ - if (ci->ti_end_min > ti_current) return; + if (!cell_is_active(ci, r->e)) return; /* Recurse? */ if (ci->split) { @@ -2287,8 +2277,8 @@ void DOSUB_SELF2(struct runner *r, struct cell *ci, int gettimer) { void DOSUB_SUBSET(struct runner *r, struct cell *ci, struct part *parts, int *ind, int count, struct cell *cj, int sid, int gettimer) { - struct space *s = r->e->s; - const int ti_current = r->e->ti_current; + const struct engine *e = r->e; + struct space *s = e->s; TIMER_TIC; @@ -2850,7 +2840,7 @@ void DOSUB_SUBSET(struct runner *r, struct cell *ci, struct part *parts, } /* Otherwise, compute the pair directly. */ - else if (ci->ti_end_min <= ti_current || cj->ti_end_min <= ti_current) { + else if (cell_is_active(ci, e) || cell_is_active(cj, e)) { /* Get the relative distance between the pairs, wrapping. */ double shift[3] = {0.0, 0.0, 0.0}; diff --git a/src/runner_doiact_grav.h b/src/runner_doiact_grav.h index 0fcd2d2e80a72b92588acd5b8275b9dafc68df45..59a5ae496680390c23458bde65b4bba321ffe7a1 100644 --- a/src/runner_doiact_grav.h +++ b/src/runner_doiact_grav.h @@ -71,7 +71,6 @@ __attribute__((always_inline)) INLINE static void runner_dopair_grav_pm( const int gcount = ci->gcount; struct gpart *restrict gparts = ci->gparts; const struct multipole multi = cj->multipole; - const int ti_current = e->ti_current; const float rlr_inv = 1. / (const_gravity_a_smooth * ci->super->width[0]); TIMER_TIC; @@ -84,7 +83,7 @@ __attribute__((always_inline)) INLINE static void runner_dopair_grav_pm( #endif /* Anything to do here? */ - if (ci->ti_end_min > ti_current) return; + if (!cell_is_active(ci, e)) return; #if ICHECK > 0 for (int pid = 0; pid < gcount; pid++) { @@ -105,7 +104,7 @@ __attribute__((always_inline)) INLINE static void runner_dopair_grav_pm( /* Get a hold of the ith part in ci. */ struct gpart *restrict gp = &gparts[pid]; - if (gp->ti_end > ti_current) continue; + if (!gpart_is_active(gp, e)) continue; /* Compute the pairwise distance. */ const float dx[3] = {multi.CoM[0] - gp->x[0], // x @@ -138,7 +137,6 @@ __attribute__((always_inline)) INLINE static void runner_dopair_grav_pp( const int gcount_j = cj->gcount; struct gpart *restrict gparts_i = ci->gparts; struct gpart *restrict gparts_j = cj->gparts; - const int ti_current = e->ti_current; const float rlr_inv = 1. / (const_gravity_a_smooth * ci->super->width[0]); TIMER_TIC; @@ -150,7 +148,7 @@ __attribute__((always_inline)) INLINE static void runner_dopair_grav_pp( #endif /* Anything to do here? */ - if (ci->ti_end_min > ti_current && cj->ti_end_min > ti_current) return; + if (!cell_is_active(ci, e) && !cell_is_active(cj, e)) return; #if ICHECK > 0 for (int pid = 0; pid < gcount_i; pid++) { @@ -195,17 +193,17 @@ __attribute__((always_inline)) INLINE static void runner_dopair_grav_pp( const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2]; /* Interact ! */ - if (gpi->ti_end <= ti_current && gpj->ti_end <= ti_current) { + if (gpart_is_active(gpi, e) && gpart_is_active(gpj, e)) { runner_iact_grav_pp(rlr_inv, r2, dx, gpi, gpj); } else { - if (gpi->ti_end <= ti_current) { + if (gpart_is_active(gpi, e)) { runner_iact_grav_pp_nonsym(rlr_inv, r2, dx, gpi, gpj); - } else if (gpj->ti_end <= ti_current) { + } else if (gpart_is_active(gpj, e)) { dx[0] = -dx[0]; dx[1] = -dx[1]; @@ -233,7 +231,6 @@ __attribute__((always_inline)) INLINE static void runner_doself_grav_pp( const struct engine *e = r->e; const int gcount = c->gcount; struct gpart *restrict gparts = c->gparts; - const int ti_current = e->ti_current; const float rlr_inv = 1. / (const_gravity_a_smooth * c->super->width[0]); TIMER_TIC; @@ -244,7 +241,7 @@ __attribute__((always_inline)) INLINE static void runner_doself_grav_pp( #endif /* Anything to do here? */ - if (c->ti_end_min > ti_current) return; + if (!cell_is_active(c, e)) return; #if ICHECK > 0 for (int pid = 0; pid < gcount; pid++) { @@ -278,17 +275,17 @@ __attribute__((always_inline)) INLINE static void runner_doself_grav_pp( const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2]; /* Interact ! */ - if (gpi->ti_end <= ti_current && gpj->ti_end <= ti_current) { + if (gpart_is_active(gpi, e) && gpart_is_active(gpj, e)) { runner_iact_grav_pp(rlr_inv, r2, dx, gpi, gpj); } else { - if (gpi->ti_end <= ti_current) { + if (gpart_is_active(gpi, e)) { runner_iact_grav_pp_nonsym(rlr_inv, r2, dx, gpi, gpj); - } else if (gpj->ti_end <= ti_current) { + } else if (gpart_is_active(gpj, e)) { dx[0] = -dx[0]; dx[1] = -dx[1]; @@ -490,14 +487,13 @@ static void runner_do_grav_mm(struct runner *r, struct cell *ci, int timer) { const struct engine *e = r->e; struct cell *cells = e->s->cells_top; const int nr_cells = e->s->nr_cells; - const int ti_current = e->ti_current; const double max_d = const_gravity_a_smooth * const_gravity_r_cut * ci->width[0]; const double max_d2 = max_d * max_d; const double pos_i[3] = {ci->loc[0], ci->loc[1], ci->loc[2]}; /* Anything to do here? */ - if (ci->ti_end_min > ti_current) return; + if (!cell_is_active(ci, e)) return; /* Loop over all the cells and go for a p-m interaction if far enough but not * too far */ diff --git a/src/tools.c b/src/tools.c index 060bf1439f30dc6237938c060bc4ddc8d9be822b..e526bb1b838f6d97b72eadb4070f3f2a94938c04 100644 --- a/src/tools.c +++ b/src/tools.c @@ -481,7 +481,7 @@ void engine_single_density(double *dim, long long int pid, } /* Dump the result. */ - hydro_end_density(&p, 0); + hydro_end_density(&p); message("part %lli (h=%e) has wcount=%e, rho=%e.", p.id, p.h, p.density.wcount, hydro_get_density(&p)); fflush(stdout);