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..17adfd07d2bfe0519f0b432217d5253de13c4b78
--- /dev/null
+++ b/src/active.h
@@ -0,0 +1,104 @@
+/*******************************************************************************
+ * 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 "const.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) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (c->ti_end_min < e->ti_current)
+ error("cell in an impossible time-zone! c->ti_end_min=%d e->ti_current=%d",
+ c->ti_end_min, e->ti_current);
+#endif
+
+ 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) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (c->ti_end_max < e->ti_current)
+ error("cell in an impossible time-zone! c->ti_end_max=%d e->ti_current=%d",
+ c->ti_end_max, e->ti_current);
+#endif
+
+ 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) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (p->ti_end < e->ti_current)
+ error("particle in an impossible time-zone! p->ti_end=%d e->ti_current=%d",
+ p->ti_end, e->ti_current);
+#endif
+
+ 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) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (gp->ti_end < e->ti_current)
+ error(
+ "g-particle in an impossible time-zone! gp->ti_end=%d e->ti_current=%d",
+ gp->ti_end, e->ti_current);
+#endif
+
+ return (gp->ti_end == e->ti_current);
+}
+
+#endif /* SWIFT_ACTIVE_H */
diff --git a/src/cell.c b/src/cell.c
index 573272d05839d6d082dac61c97f6abd18d8eb41a..18c2279b76cdd54a626dfe91f655c9e4e01d6b1b 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -47,6 +47,7 @@
#include "cell.h"
/* Local headers. */
+#include "active.h"
#include "atomic.h"
#include "error.h"
#include "gravity.h"
@@ -867,14 +868,14 @@ void cell_clean(struct cell *c) {
* @brief Checks whether a given cell needs drifting or not.
*
* @param c the #cell.
- * @param ti_current The current time on the integer time-line.
+ * @param e The #engine (holding current time information).
*
* @return 1 If the cell needs drifting, 0 otherwise.
*/
-int cell_is_drift_needed(struct cell *c, int ti_current) {
+int cell_is_drift_needed(struct cell *c, const struct engine *e) {
/* Do we have at least one active particle in the cell ?*/
- if (c->ti_end_min == ti_current) return 1;
+ if (cell_is_active(c, e)) return 1;
/* Loop over the pair tasks that involve this cell */
for (struct link *l = c->density; l != NULL; l = l->next) {
@@ -882,9 +883,9 @@ int cell_is_drift_needed(struct cell *c, int ti_current) {
if (l->t->type != task_type_pair && l->t->type != task_type_sub_pair)
continue;
- /* Does the other cell in the pair have an active particle ? */
- if ((l->t->ci == c && l->t->cj->ti_end_min == ti_current) ||
- (l->t->cj == c && l->t->ci->ti_end_min == ti_current))
+ /* Is the other cell in the pair active ? */
+ if ((l->t->ci == c && cell_is_active(l->t->cj, e)) ||
+ (l->t->cj == c && cell_is_active(l->t->ci, e)))
return 1;
}
diff --git a/src/cell.h b/src/cell.h
index 9e5bed091178b59e1b757c420bc1d5fde0b9ce42..10c67be768352374b47736bed225a136f05dba28 100644
--- a/src/cell.h
+++ b/src/cell.h
@@ -37,6 +37,7 @@
#include "task.h"
/* Avoid cyclic inclusions */
+struct engine;
struct space;
struct scheduler;
@@ -290,7 +291,7 @@ int cell_are_neighbours(const struct cell *restrict ci,
const struct cell *restrict cj);
void cell_check_multipole(struct cell *c, void *data);
void cell_clean(struct cell *c);
-int cell_is_drift_needed(struct cell *c, int ti_current);
+int cell_is_drift_needed(struct cell *c, const struct engine *e);
int cell_unskip_tasks(struct cell *c, struct scheduler *s);
void cell_set_super(struct cell *c, struct cell *super);
diff --git a/src/debug.c b/src/debug.c
index 196c2f7e49afa827f8d80539cdeb3a975cbf31fc..5be0370f64f2c21e6b28c40cc9802520087ae07f 100644
--- a/src/debug.c
+++ b/src/debug.c
@@ -193,6 +193,59 @@ int checkSpacehmax(struct space *s) {
return 0;
}
+/**
+ * @brief Check if the h_max and dx_max values of a cell's hierarchy are
+ * consistent with the particles. Report verbosely if not.
+ *
+ * @param c the top cell of the hierarchy.
+ * @param depth the recursion depth for use in messages. Set to 0 initially.
+ * @result 1 or 0
+ */
+int checkCellhdxmax(const struct cell *c, int *depth) {
+
+ *depth = *depth + 1;
+
+ float h_max = 0.0f;
+ float dx_max = 0.0f;
+ if (!c->split) {
+ const size_t nr_parts = c->count;
+ struct part *parts = c->parts;
+ for (size_t k = 0; k < nr_parts; k++) {
+ h_max = (h_max > parts[k].h) ? h_max : parts[k].h;
+ }
+ } else {
+ for (int k = 0; k < 8; k++)
+ if (c->progeny[k] != NULL) {
+ struct cell *cp = c->progeny[k];
+ checkCellhdxmax(cp, depth);
+ dx_max = max(dx_max, cp->dx_max);
+ h_max = max(h_max, cp->h_max);
+ }
+ }
+
+ /* Check. */
+ int result = 1;
+ if (c->h_max != h_max) {
+ message("%d Inconsistent h_max: cell %f != parts %f", *depth, c->h_max, h_max);
+ message("location: %f %f %f", c->loc[0], c->loc[1], c->loc[2]);
+ result = 0;
+ }
+ if (c->dx_max != dx_max) {
+ message("%d Inconsistent dx_max: %f != %f", *depth, c->dx_max, dx_max);
+ message("location: %f %f %f", c->loc[0], c->loc[1], c->loc[2]);
+ result = 0;
+ }
+
+ /* Check rebuild criterion. */
+ if (h_max > c->dmin) {
+ message("%d Inconsistent c->dmin: %f > %f", *depth, h_max, c->dmin);
+ message("location: %f %f %f", c->loc[0], c->loc[1], c->loc[2]);
+ result = 0;
+ }
+
+ return result;
+}
+
#ifdef HAVE_METIS
/**
diff --git a/src/debug.h b/src/debug.h
index d75ce91f1ea23221626f32e11472f713e9731789..7422a6f7f9815490966f08415e0312876ce0123f 100644
--- a/src/debug.h
+++ b/src/debug.h
@@ -32,6 +32,7 @@ void printParticle_single(const struct part *p, const struct xpart *xp);
void printgParticle_single(struct gpart *gp);
int checkSpacehmax(struct space *s);
+int checkCellhdxmax(const struct cell *c, int *depth);
#ifdef HAVE_METIS
#include "metis.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 cd047514ab6ff4c67e13bd67daabc8f00af5f67a..2d6da4e4aedc9c40d1dade243e605e9aeda86dbe 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) {
@@ -182,8 +183,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);
}
}
@@ -484,12 +484,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) {
@@ -503,7 +503,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);
@@ -516,7 +516,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);
@@ -541,12 +541,12 @@ void runner_do_extra_ghost(struct runner *r, struct cell *c, int timer) {
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;
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) {
@@ -560,7 +560,7 @@ void runner_do_extra_ghost(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 force calculation */
hydro_end_gradient(p);
@@ -588,20 +588,20 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
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) {
@@ -630,10 +630,10 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
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;
@@ -754,20 +754,18 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
*/
static void runner_do_drift(struct cell *c, struct engine *e, int drift) {
- const int ti_current = e->ti_current;
-
/* Unskip any active tasks. */
- if (c->ti_end_min == e->ti_current) {
+ if (cell_is_active(c, e)) {
const int forcerebuild = cell_unskip_tasks(c, &e->sched);
if (forcerebuild) atomic_inc(&e->forcerebuild);
}
/* Do we really need to drift? */
if (drift) {
- if (!e->drift_all && !cell_is_drift_needed(c, ti_current)) return;
+ if (!e->drift_all && !cell_is_drift_needed(c, e)) return;
} else {
- /* Not drifting, but may still need to recurse for task skipping. */
+ /* Not drifting, but may still need to recurse for task un-skipping. */
if (c->split) {
for (int k = 0; k < 8; k++) {
if (c->progeny[k] != NULL) {
@@ -779,8 +777,12 @@ static void runner_do_drift(struct cell *c, struct engine *e, int drift) {
return;
}
+ /* Now, we can drift */
+
+ /* Get some information first */
const double timeBase = e->timeBase;
const int ti_old = c->ti_old;
+ const int ti_current = e->ti_current;
struct part *const parts = c->parts;
struct xpart *const xparts = c->xparts;
struct gpart *const gparts = c->gparts;
@@ -793,7 +795,7 @@ static void runner_do_drift(struct cell *c, struct engine *e, int drift) {
if (!c->split) {
/* Check that we are actually going to move forward. */
- if (ti_current >= ti_old) {
+ if (ti_current > ti_old) {
/* Loop over all the g-particles in the cell */
const size_t nr_gparts = c->gcount;
@@ -837,6 +839,12 @@ static void runner_do_drift(struct cell *c, struct engine *e, int drift) {
dx_max = sqrtf(dx2_max);
} /* Check that we are actually going to move forward. */
+
+ else {
+ /* ti_old == ti_current, just keep the current cell values. */
+ h_max = c->h_max;
+ dx_max = c->dx_max;
+ }
}
/* Otherwise, aggregate data from children. */
@@ -898,18 +906,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;
@@ -930,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);
@@ -961,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);
@@ -1119,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",
@@ -1127,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 "
@@ -1136,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..c067d3bc9a576ee9c7dfb6e910eb1aa01012f755 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/space.c b/src/space.c
index 18f0e21e87b22715e6318be5dfbc3c0cc9a2cc12..5f7e3522eef8c398faba1aee4b03c49ab94232b1 100644
--- a/src/space.c
+++ b/src/space.c
@@ -346,7 +346,6 @@ void space_regrid(struct space *s, int verbose) {
if (verbose)
message("set cell dimensions to [ %i %i %i ].", cdim[0], cdim[1],
cdim[2]);
- fflush(stdout);
#ifdef WITH_MPI
if (oldnodeIDs != NULL) {
@@ -1563,6 +1562,15 @@ void space_split_mapper(void *map_data, int num_cells, void *extra_data) {
else
c->owner = 0; /* Ok, there is really nothing on this rank... */
}
+
+#ifdef SWIFT_DEBUG_CHECKS
+ /* All cells and particles should have consistent h_max values. */
+ for (int ind = 0; ind < num_cells; ind++) {
+ int depth = 0;
+ if (!checkCellhdxmax(&cells_top[ind], &depth))
+ message(" at cell depth %d", depth);
+ }
+#endif
}
/**
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);
diff --git a/tests/test27cells.c b/tests/test27cells.c
index 22619af53c39218da34d771fab6ed2d10993689c..f58b4dc410637f3d91369dab1b442de0b7044c08 100644
--- a/tests/test27cells.c
+++ b/tests/test27cells.c
@@ -167,7 +167,7 @@ void zero_particle_fields(struct cell *c) {
*/
void end_calculation(struct cell *c) {
for (int pid = 0; pid < c->count; pid++) {
- hydro_end_density(&c->parts[pid], 1);
+ hydro_end_density(&c->parts[pid]);
}
}