Commit 907b0ce7 authored by Peter W. Draper's avatar Peter W. Draper
Browse files

Merge branch 'new_time_line' into 'master'

Improve readability of active/inactive condition

This, as is, is rather minor but it's on the way towards improving the time-stepping. 

I have replaced all explicit checks that a particle or cell is active by calls to inlined functions. No change in functionality thus far.

The plan is to slowly introduce a better time-stepping mechanism where the particles carry less information (lower footprint) and make the whole time-step limiter problem easier. With this merge request in then I only need to make the changes in one single location (file active.h) and not everywhere. 

See merge request !279
parents fb49f1ee fa3be81e
......@@ -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 "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 */
......@@ -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;
}
......
......@@ -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);
......
......@@ -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
/**
......
......@@ -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"
......
......@@ -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) {
......@@ -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",
......
......@@ -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;