Commit 4f7fb709 authored by James Willis's avatar James Willis
Browse files

Merge branch 'master' into space_rebuild

parents 1fbd0435 2474a3c2
......@@ -285,6 +285,11 @@ int main(int argc, char *argv[]) {
message("Running on: %s", hostname());
#endif
/* Do we have debugging checks ? */
#ifdef SWIFT_DEBUG_CHECKS
message("WARNING: Debugging checks activated. Code will be slower !");
#endif
/* Do we choke on FP-exceptions ? */
if (with_fp_exceptions) {
feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
......
......@@ -160,6 +160,9 @@ if test "$ac_test_CFLAGS" != "set"; then
# note that we enable "unsafe" fp optimization with other compilers, too
AX_CHECK_COMPILE_FLAG(-ffast-math, CFLAGS="$CFLAGS -ffast-math")
# not all codes will benefit from this.
AX_CHECK_COMPILE_FLAG(-funroll-loops, CFLAGS="$CFLAGS -funroll-loops")
AX_GCC_ARCHFLAG($acx_maxopt_portable)
;;
......
......@@ -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 Check that a cell been drifted to the current time.
*
* Only used for debugging. Calls error() if the cell has not
* been drifted. Does nothing if SWIFT_DEBUG_CHECKS is not defined.
*
* @param c The #cell.
* @param e The #engine containing information about the current time.
*/
__attribute__((always_inline)) INLINE static void cell_is_drifted(
const struct cell *c, const struct engine *e) {
#ifdef SWIFT_DEBUG_CHECKS
if (c->ti_old > e->ti_current)
error(
"Cell has been drifted too far forward in time! c->ti_old=%d "
"e->ti_current=%d",
c->ti_old, e->ti_current);
if (c->ti_old != e->ti_current) {
error(
"Cell has not been drifted to the current time c->ti_old=%d, "
"e->ti_current=%d",
c->ti_old, e->ti_current);
}
#endif
}
/**
* @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"
......@@ -700,6 +701,23 @@ void cell_clean_links(struct cell *c, void *data) {
c->grav = NULL;
}
/**
* @brief Checks that a cell is at the current point in time
*
* Calls error() if the cell is not at the current time.
*
* @param c Cell to act upon
* @param data The current time on the integer time-line
*/
void cell_check_drift_point(struct cell *c, void *data) {
const int ti_current = *(int *)data;
if (c->ti_old != ti_current)
error("Cell in an incorrect time-zone! c->ti_old=%d ti_current=%d",
c->ti_old, ti_current);
}
/**
* @brief Checks whether the cells are direct neighbours ot not. Both cells have
* to be of the same size
......@@ -805,14 +823,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) {
......@@ -820,9 +838,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,8 @@ 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);
void cell_check_drift_point(struct cell *c, void *data);
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,60 @@ 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"
......
......@@ -557,6 +557,11 @@ void engine_repartition(struct engine *e) {
ticks tic = getticks();
#ifdef SWIFT_DEBUG_CHECKS
/* Check that all cells have been drifted to the current time */
space_check_drift_point(e->s, e->ti_current);
#endif
/* Clear the repartition flag. */
enum repartition_type reparttype = e->forcerepart;
e->forcerepart = REPART_NONE;
......@@ -2234,6 +2239,11 @@ void engine_prepare(struct engine *e, int nodrift) {
e->drift_all = (e->policy & engine_policy_drift_all);
}
#ifdef SWIFT_DEBUG_CHECKS
/* Check that all cells have been drifted to the current time */
space_check_drift_point(e->s, e->ti_current);
#endif
engine_rebuild(e);
}
......
......@@ -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;
......
This diff is collapsed.
......@@ -48,11 +48,13 @@ struct runner {
};
/* Function prototypes. */
void runner_do_ghost(struct runner *r, struct cell *c);
void runner_do_ghost(struct runner *r, struct cell *c, int timer);
void runner_do_extra_ghost(struct runner *r, struct cell *c, int timer);
void runner_do_sort(struct runner *r, struct cell *c, int flag, int clock);
void runner_do_kick(struct runner *r, struct cell *c, int timer);
void runner_do_init(struct runner *r, struct cell *c, int timer);
void runner_do_cooling(struct runner *r, struct cell *c, int timer);
void runner_do_grav_external(struct runner *r, struct cell *c, int timer);
void *runner_main(void *data);
void runner_do_drift_mapper(void *map_data, int num_elements, void *extra_data);
......
......@@ -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,12 @@ 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;
#ifdef SWIFT_DEBUG_CHECKS
cell_is_drifted(ci, e);
cell_is_drifted(cj, e);
#endif
/* Get the sort ID. */
double shift[3] = {0.0, 0.0, 0.0};
......@@ -756,7 +759,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 +821,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 +897,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 +916,12 @@ 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;
#ifdef SWIFT_DEBUG_CHECKS
cell_is_drifted(ci, e);
cell_is_drifted(cj, e);
#endif
/* Get the shift ID. */
double shift[3] = {0.0, 0.0, 0.0};
......@@ -946,28 +953,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 +995,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 +1069,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 +1077,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 +1139,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 +1212,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 +1220,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 +1277,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