Commit dab5f545 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Check whether particles or cells are active via an inlined function rather...

Check whether particles or cells are active via an inlined function rather than via a direct read of engine->ti_current
parent 06b50e49
......@@ -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 \
......
/*******************************************************************************
* 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 */
......@@ -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;
......
......@@ -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;
......
......@@ -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;
......
......@@ -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;
......
......@@ -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;
......
......@@ -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",
......
......@@ -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)) {