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

Merge branch 'master' into gravity_fixes

parents b5e7b69d d408db13
......@@ -100,6 +100,7 @@ tests/testVoronoi3D
tests/testDump
tests/testLogger
tests/benchmarkInteractions
tests/testGravityDerivatives
theory/latex/swift.pdf
theory/SPH/Kernels/kernels.pdf
......
# This file is part of SWIFT.
# tHIS FIle is part of SWIFT.
# Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
# Matthieu Schaller (matthieu.schaller@durham.ac.uk).
#
......@@ -95,9 +95,7 @@ EXTRA_DIST = BigCosmoVolume/makeIC.py \
EXTRA_DIST += parameter_example.yml
# Scripts to plot task graphs
EXTRA_DIST += plot_tasks_MPI.py plot_tasks.py \
analyse_tasks_MPI.py analyse_tasks.py \
process_plot_tasks_MPI process_plot_tasks
EXTRA_DIST += plot_tasks.py analyse_tasks.py process_plot_tasks_MPI process_plot_tasks
# Scripts to plot threadpool 'task' graphs
EXTRA_DIST += analyse_threadpool_tasks.py \
......
......@@ -370,11 +370,11 @@ __attribute__((always_inline)) INLINE void cache_read_two_partial_cells_sorted(
#ifdef SWIFT_DEBUG_CHECKS
const float shift_threshold_x =
4. * ci->width[0] * (1. + 2. * space_maxreldx);
2. * ci->width[0] + 2. * max(ci->dx_max_part, cj->dx_max_part);
const float shift_threshold_y =
4. * ci->width[1] * (1. + 2. * space_maxreldx);
2. * ci->width[1] + 2. * max(ci->dx_max_part, cj->dx_max_part);
const float shift_threshold_z =
4. * ci->width[2] * (1. + 2. * space_maxreldx);
2. * ci->width[2] + 2. * max(ci->dx_max_part, cj->dx_max_part);
/* Make sure that particle positions have been shifted correctly. */
for (int i = 0; i < ci_cache_count; i++) {
......
......@@ -82,73 +82,6 @@ int cell_getsize(struct cell *c) {
return count;
}
/**
* @brief Unpack the data of a given cell and its sub-cells.
*
* @param pc An array of packed #pcell.
* @param c The #cell in which to unpack the #pcell.
* @param s The #space in which the cells are created.
*
* @return The number of cells created.
*/
int cell_unpack(struct pcell *pc, struct cell *c, struct space *s) {
#ifdef WITH_MPI
/* Unpack the current pcell. */
c->h_max = pc->h_max;
c->ti_end_min = pc->ti_end_min;
c->ti_end_max = pc->ti_end_max;
c->ti_old_part = pc->ti_old_part;
c->ti_old_gpart = pc->ti_old_gpart;
c->count = pc->count;
c->gcount = pc->gcount;
c->scount = pc->scount;
c->tag = pc->tag;
/* Number of new cells created. */
int count = 1;
/* Fill the progeny recursively, depth-first. */
for (int k = 0; k < 8; k++)
if (pc->progeny[k] >= 0) {
struct cell *temp;
space_getcells(s, 1, &temp);
temp->count = 0;
temp->gcount = 0;
temp->scount = 0;
temp->loc[0] = c->loc[0];
temp->loc[1] = c->loc[1];
temp->loc[2] = c->loc[2];
temp->width[0] = c->width[0] / 2;
temp->width[1] = c->width[1] / 2;
temp->width[2] = c->width[2] / 2;
temp->dmin = c->dmin / 2;
if (k & 4) temp->loc[0] += temp->width[0];
if (k & 2) temp->loc[1] += temp->width[1];
if (k & 1) temp->loc[2] += temp->width[2];
temp->depth = c->depth + 1;
temp->split = 0;
temp->dx_max_part = 0.f;
temp->dx_max_gpart = 0.f;
temp->dx_max_sort = 0.f;
temp->nodeID = c->nodeID;
temp->parent = c;
c->progeny[k] = temp;
c->split = 1;
count += cell_unpack(&pc[pc->progeny[k]], temp, s);
}
/* Return the total number of unpacked cells. */
c->pcell_size = count;
return count;
#else
error("SWIFT was not compiled with MPI support.");
return 0;
#endif
}
/**
* @brief Link the cells recursively to the given #part array.
*
......@@ -233,7 +166,7 @@ int cell_link_sparts(struct cell *c, struct spart *sparts) {
*
* @return The number of packed cells.
*/
int cell_pack(struct cell *c, struct pcell *pc) {
int cell_pack(struct cell *restrict c, struct pcell *restrict pc) {
#ifdef WITH_MPI
......@@ -267,26 +200,97 @@ int cell_pack(struct cell *c, struct pcell *pc) {
#endif
}
/**
* @brief Unpack the data of a given cell and its sub-cells.
*
* @param pc An array of packed #pcell.
* @param c The #cell in which to unpack the #pcell.
* @param s The #space in which the cells are created.
*
* @return The number of cells created.
*/
int cell_unpack(struct pcell *restrict pc, struct cell *restrict c,
struct space *restrict s) {
#ifdef WITH_MPI
/* Unpack the current pcell. */
c->h_max = pc->h_max;
c->ti_end_min = pc->ti_end_min;
c->ti_end_max = pc->ti_end_max;
c->ti_old_part = pc->ti_old_part;
c->ti_old_gpart = pc->ti_old_gpart;
c->count = pc->count;
c->gcount = pc->gcount;
c->scount = pc->scount;
c->tag = pc->tag;
/* Number of new cells created. */
int count = 1;
/* Fill the progeny recursively, depth-first. */
for (int k = 0; k < 8; k++)
if (pc->progeny[k] >= 0) {
struct cell *temp;
space_getcells(s, 1, &temp);
temp->count = 0;
temp->gcount = 0;
temp->scount = 0;
temp->loc[0] = c->loc[0];
temp->loc[1] = c->loc[1];
temp->loc[2] = c->loc[2];
temp->width[0] = c->width[0] / 2;
temp->width[1] = c->width[1] / 2;
temp->width[2] = c->width[2] / 2;
temp->dmin = c->dmin / 2;
if (k & 4) temp->loc[0] += temp->width[0];
if (k & 2) temp->loc[1] += temp->width[1];
if (k & 1) temp->loc[2] += temp->width[2];
temp->depth = c->depth + 1;
temp->split = 0;
temp->dx_max_part = 0.f;
temp->dx_max_gpart = 0.f;
temp->dx_max_sort = 0.f;
temp->nodeID = c->nodeID;
temp->parent = c;
c->progeny[k] = temp;
c->split = 1;
count += cell_unpack(&pc[pc->progeny[k]], temp, s);
}
/* Return the total number of unpacked cells. */
c->pcell_size = count;
return count;
#else
error("SWIFT was not compiled with MPI support.");
return 0;
#endif
}
/**
* @brief Pack the time information of the given cell and all it's sub-cells.
*
* @param c The #cell.
* @param ti_ends (output) The time information we pack into
* @param pcells (output) The end-of-timestep information we pack into
*
* @return The number of packed cells.
*/
int cell_pack_ti_ends(struct cell *c, integertime_t *ti_ends) {
int cell_pack_end_step(struct cell *restrict c,
struct pcell_step *restrict pcells) {
#ifdef WITH_MPI
/* Pack this cell's data. */
ti_ends[0] = c->ti_end_min;
pcells[0].ti_end_min = c->ti_end_min;
pcells[0].dx_max_part = c->dx_max_part;
pcells[0].dx_max_gpart = c->dx_max_gpart;
/* Fill in the progeny, depth-first recursion. */
int count = 1;
for (int k = 0; k < 8; k++)
if (c->progeny[k] != NULL) {
count += cell_pack_ti_ends(c->progeny[k], &ti_ends[count]);
count += cell_pack_end_step(c->progeny[k], &pcells[count]);
}
/* Return the number of packed values. */
......@@ -302,22 +306,25 @@ int cell_pack_ti_ends(struct cell *c, integertime_t *ti_ends) {
* @brief Unpack the time information of a given cell and its sub-cells.
*
* @param c The #cell
* @param ti_ends The time information to unpack
* @param pcells The end-of-timestep information to unpack
*
* @return The number of cells created.
*/
int cell_unpack_ti_ends(struct cell *c, integertime_t *ti_ends) {
int cell_unpack_end_step(struct cell *restrict c,
struct pcell_step *restrict pcells) {
#ifdef WITH_MPI
/* Unpack this cell's data. */
c->ti_end_min = ti_ends[0];
c->ti_end_min = pcells[0].ti_end_min;
c->dx_max_part = pcells[0].dx_max_part;
c->dx_max_gpart = pcells[0].dx_max_gpart;
/* Fill in the progeny, depth-first recursion. */
int count = 1;
for (int k = 0; k < 8; k++)
if (c->progeny[k] != NULL) {
count += cell_unpack_ti_ends(c->progeny[k], &ti_ends[count]);
count += cell_unpack_end_step(c->progeny[k], &pcells[count]);
}
/* Return the number of packed values. */
......
......@@ -70,24 +70,63 @@ struct link {
struct link *next;
};
/* Packed cell. */
/**
* @brief Packed cell for information correct at rebuild time.
*
* Contains all the information for a tree walk in a non-local cell.
*/
struct pcell {
/* Stats on this cell's particles. */
/*! Maximal smoothing length. */
double h_max;
integertime_t ti_end_min, ti_end_max, ti_beg_max, ti_old_part, ti_old_gpart;
/* Number of particles in this cell. */
int count, gcount, scount;
/*! Minimal integer end-of-timestep in this cell */
integertime_t ti_end_min;
/*! Maximal integer end-of-timestep in this cell */
integertime_t ti_end_max;
/*! Maximal integer beginning-of-timestep in this cell */
integertime_t ti_beg_max;
/*! Integer time of the last drift of the #part in this cell */
integertime_t ti_old_part;
/*! Integer time of the last drift of the #gpart in this cell */
integertime_t ti_old_gpart;
/*! Number of #part in this cell. */
int count;
/*! Number of #gpart in this cell. */
int gcount;
/*! Number of #spart in this cell. */
int scount;
/* tag used for MPI communication. */
/*! tag used for MPI communication. */
int tag;
/* Relative indices of the cell's progeny. */
/*! Relative indices of the cell's progeny. */
int progeny[8];
} SWIFT_STRUCT_ALIGN;
/**
* @brief Cell information at the end of a time-step.
*/
struct pcell_step {
/*! Minimal integer end-of-timestep in this cell */
integertime_t ti_end_min;
/*! Maximal distance any #part has travelled since last rebuild */
float dx_max_part;
/*! Maximal distance any #gpart has travelled since last rebuild */
float dx_max_gpart;
};
/**
* @brief Cell within the tree structure.
*
......@@ -389,8 +428,8 @@ int cell_slocktree(struct cell *c);
void cell_sunlocktree(struct cell *c);
int cell_pack(struct cell *c, struct pcell *pc);
int cell_unpack(struct pcell *pc, struct cell *c, struct space *s);
int cell_pack_ti_ends(struct cell *c, integertime_t *ti_ends);
int cell_unpack_ti_ends(struct cell *c, integertime_t *ti_ends);
int cell_pack_end_step(struct cell *c, struct pcell_step *pcell);
int cell_unpack_end_step(struct cell *c, struct pcell_step *pcell);
int cell_getsize(struct cell *c);
int cell_link_parts(struct cell *c, struct part *parts);
int cell_link_gparts(struct cell *c, struct gpart *gparts);
......
......@@ -1961,7 +1961,7 @@ void *runner_main(void *data) {
break;
case task_type_recv:
if (t->subtype == task_subtype_tend) {
cell_unpack_ti_ends(ci, t->buff);
cell_unpack_end_step(ci, t->buff);
free(t->buff);
} else if (t->subtype == task_subtype_xv) {
runner_do_recv_part(r, ci, 1, 1);
......
......@@ -808,6 +808,15 @@ void DOPAIR1(struct runner *r, struct cell *ci, struct cell *cj, const int sid,
cj->nodeID, ci->nodeID, d, sort_j[pjd].d, cj->dx_max_sort,
cj->dx_max_sort_old);
}
/* Some constants used to checks that the parts are in the right frame */
const float shift_threshold_x =
2. * ci->width[0] + 2. * max(ci->dx_max_part, cj->dx_max_part);
const float shift_threshold_y =
2. * ci->width[1] + 2. * max(ci->dx_max_part, cj->dx_max_part);
const float shift_threshold_z =
2. * ci->width[2] + 2. * max(ci->dx_max_part, cj->dx_max_part);
#endif /* SWIFT_DEBUG_CHECKS */
/* Get some other useful values. */
......@@ -859,6 +868,32 @@ void DOPAIR1(struct runner *r, struct cell *ci, struct cell *cj, const int sid,
const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
#ifdef SWIFT_DEBUG_CHECKS
/* Check that particles are in the correct frame after the shifts */
if (pix > shift_threshold_x || pix < -shift_threshold_x)
error(
"Invalid particle position in X for pi (pix=%e ci->width[0]=%e)",
pix, ci->width[0]);
if (piy > shift_threshold_y || piy < -shift_threshold_y)
error(
"Invalid particle position in Y for pi (piy=%e ci->width[1]=%e)",
piy, ci->width[1]);
if (piz > shift_threshold_z || piz < -shift_threshold_z)
error(
"Invalid particle position in Z for pi (piz=%e ci->width[2]=%e)",
piz, ci->width[2]);
if (pjx > shift_threshold_x || pjx < -shift_threshold_x)
error(
"Invalid particle position in X for pj (pjx=%e ci->width[0]=%e)",
pjx, ci->width[0]);
if (pjy > shift_threshold_y || pjy < -shift_threshold_y)
error(
"Invalid particle position in Y for pj (pjy=%e ci->width[1]=%e)",
pjy, ci->width[1]);
if (pjz > shift_threshold_z || pjz < -shift_threshold_z)
error(
"Invalid particle position in Z for pj (pjz=%e ci->width[2]=%e)",
pjz, ci->width[2]);
/* Check that particles have been drifted to the current time */
if (pi->ti_drift != e->ti_current)
error("Particle pi not drifted to current time");
......@@ -913,6 +948,32 @@ void DOPAIR1(struct runner *r, struct cell *ci, struct cell *cj, const int sid,
const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
#ifdef SWIFT_DEBUG_CHECKS
/* Check that particles are in the correct frame after the shifts */
if (pix > shift_threshold_x || pix < -shift_threshold_x)
error(
"Invalid particle position in X for pi (pix=%e ci->width[0]=%e)",
pix, ci->width[0]);
if (piy > shift_threshold_y || piy < -shift_threshold_y)
error(
"Invalid particle position in Y for pi (piy=%e ci->width[1]=%e)",
piy, ci->width[1]);
if (piz > shift_threshold_z || piz < -shift_threshold_z)
error(
"Invalid particle position in Z for pi (piz=%e ci->width[2]=%e)",
piz, ci->width[2]);
if (pjx > shift_threshold_x || pjx < -shift_threshold_x)
error(
"Invalid particle position in X for pj (pjx=%e ci->width[0]=%e)",
pjx, ci->width[0]);
if (pjy > shift_threshold_y || pjy < -shift_threshold_y)
error(
"Invalid particle position in Y for pj (pjy=%e ci->width[1]=%e)",
pjy, ci->width[1]);
if (pjz > shift_threshold_z || pjz < -shift_threshold_z)
error(
"Invalid particle position in Z for pj (pjz=%e ci->width[2]=%e)",
pjz, ci->width[2]);
/* Check that particles have been drifted to the current time */
if (pi->ti_drift != e->ti_current)
error("Particle pi not drifted to current time");
......@@ -1050,6 +1111,15 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj) {
cj->nodeID, ci->nodeID, d, sort_j[pjd].d, cj->dx_max_sort,
cj->dx_max_sort_old);
}
/* Some constants used to checks that the parts are in the right frame */
const float shift_threshold_x =
2. * ci->width[0] + 2. * max(ci->dx_max_part, cj->dx_max_part);
const float shift_threshold_y =
2. * ci->width[1] + 2. * max(ci->dx_max_part, cj->dx_max_part);
const float shift_threshold_z =
2. * ci->width[2] + 2. * max(ci->dx_max_part, cj->dx_max_part);
#endif /* SWIFT_DEBUG_CHECKS */
/* Get some other useful values. */
......@@ -1154,6 +1224,32 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj) {
const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
#ifdef SWIFT_DEBUG_CHECKS
/* Check that particles are in the correct frame after the shifts */
if (pix > shift_threshold_x || pix < -shift_threshold_x)
error(
"Invalid particle position in X for pi (pix=%e ci->width[0]=%e)",
pix, ci->width[0]);
if (piy > shift_threshold_y || piy < -shift_threshold_y)
error(
"Invalid particle position in Y for pi (piy=%e ci->width[1]=%e)",
piy, ci->width[1]);
if (piz > shift_threshold_z || piz < -shift_threshold_z)
error(
"Invalid particle position in Z for pi (piz=%e ci->width[2]=%e)",
piz, ci->width[2]);
if (pjx > shift_threshold_x || pjx < -shift_threshold_x)
error(
"Invalid particle position in X for pj (pjx=%e ci->width[0]=%e)",
pjx, ci->width[0]);
if (pjy > shift_threshold_y || pjy < -shift_threshold_y)
error(
"Invalid particle position in Y for pj (pjy=%e ci->width[1]=%e)",
pjy, ci->width[1]);
if (pjz > shift_threshold_z || pjz < -shift_threshold_z)
error(
"Invalid particle position in Z for pj (pjz=%e ci->width[2]=%e)",
pjz, ci->width[2]);
/* Check that particles have been drifted to the current time */
if (pi->ti_drift != e->ti_current)
error("Particle pi not drifted to current time");
......@@ -1188,6 +1284,32 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj) {
const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
#ifdef SWIFT_DEBUG_CHECKS
/* Check that particles are in the correct frame after the shifts */
if (pix > shift_threshold_x || pix < -shift_threshold_x)
error(
"Invalid particle position in X for pi (pix=%e ci->width[0]=%e)",
pix, ci->width[0]);
if (piy > shift_threshold_y || piy < -shift_threshold_y)
error(
"Invalid particle position in Y for pi (piy=%e ci->width[1]=%e)",
piy, ci->width[1]);
if (piz > shift_threshold_z || piz < -shift_threshold_z)
error(
"Invalid particle position in Z for pi (piz=%e ci->width[2]=%e)",
piz, ci->width[2]);
if (pjx > shift_threshold_x || pjx < -shift_threshold_x)
error(
"Invalid particle position in X for pj (pjx=%e ci->width[0]=%e)",
pjx, ci->width[0]);
if (pjy > shift_threshold_y || pjy < -shift_threshold_y)
error(
"Invalid particle position in Y for pj (pjy=%e ci->width[1]=%e)",
pjy, ci->width[1]);
if (pjz > shift_threshold_z || pjz < -shift_threshold_z)
error(
"Invalid particle position in Z for pj (pjz=%e ci->width[2]=%e)",
pjz, ci->width[2]);
/* Check that particles have been drifted to the current time */
if (pi->ti_drift != e->ti_current)
error("Particle pi not drifted to current time");
......@@ -1252,6 +1374,32 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj) {
const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
#ifdef SWIFT_DEBUG_CHECKS
/* Check that particles are in the correct frame after the shifts */
if (pix > shift_threshold_x || pix < -shift_threshold_x)
error(
"Invalid particle position in X for pi (pix=%e ci->width[0]=%e)",
pix, ci->width[0]);
if (piy > shift_threshold_y || piy < -shift_threshold_y)
error(
"Invalid particle position in Y for pi (piy=%e ci->width[1]=%e)",
piy, ci->width[1]);
if (piz > shift_threshold_z || piz < -shift_threshold_z)
error(
"Invalid particle position in Z for pi (piz=%e ci->width[2]=%e)",
piz, ci->width[2]);
if (pjx > shift_threshold_x || pjx < -shift_threshold_x)
error(
"Invalid particle position in X for pj (pjx=%e ci->width[0]=%e)",
pjx, ci->width[0]);
if (pjy > shift_threshold_y || pjy < -shift_threshold_y)
error(
"Invalid particle position in Y for pj (pjy=%e ci->width[1]=%e)",
pjy, ci->width[1]);
if (pjz > shift_threshold_z || pjz < -shift_threshold_z)
error(
"Invalid particle position in Z for pj (pjz=%e ci->width[2]=%e)",
pjz, ci->width[2]);
/* Check that particles have been drifted to the current time */
if (pi->ti_drift != e->ti_current)
error("Particle pi not drifted to current time");
......@@ -1288,6 +1436,32 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj) {
const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
#ifdef SWIFT_DEBUG_CHECKS
/* Check that particles are in the correct frame after the shifts */
if (pix > shift_threshold_x || pix < -shift_threshold_x)
error(
"Invalid particle position in X for pi (pix=%e ci->width[0]=%e)",
pix, ci->width[0]);
if (piy > shift_threshold_y || piy < -shift_threshold_y)
error(
"Invalid particle position in Y for pi (piy=%e ci->width[1]=%e)",
piy, ci->width[1]);
if (piz > shift_threshold_z || piz < -shift_threshold_z)
error(
"Invalid particle position in Z for pi (piz=%e ci->width[2]=%e)",
piz, ci->width[2]);
if (pjx > shift_threshold_x || pjx < -shift_threshold_x)
error(
"Invalid particle position in X for pj (pjx=%e ci->width[0]=%e)",
pjx, ci->width[0]);
if (pjy > shift_threshold_y || pjy < -shift_threshold_y)
error(
"Invalid particle position in Y for pj (pjy=%e ci->width[1]=%e)",
pjy, ci->width[1]);
if (pjz > shift_threshold_z || pjz < -shift_threshold_z)
error(
"Invalid particle position in Z for pj (pjz=%e ci->width[2]=%e)",
pjz, ci->width[2]);
/* Check that particles have been drifted to the current time */
if (pi->ti_drift != e->ti_current)
error("Particle pi not drifted to current time");
......
......@@ -325,7 +325,7 @@ __attribute__((always_inline)) INLINE static void populate_max_index_no_cache(
last_pj = active_id;
/* Find the maximum index into cell i for each particle in range in cell j. */
if (last_pj > 0) {
if (last_pj >= 0) {
/* Start from the last particle in cell i. */
temp = ci->count - 1;
......
......@@ -1274,10 +1274,10 @@ void scheduler_enqueue(struct scheduler *s, struct task *t) {
case task_type_recv:
#ifdef WITH_MPI
if (t->subtype == task_subtype_tend) {
t->buff = malloc(sizeof(integertime_t) * t->ci->pcell_size);
err = MPI_Irecv(t->buff, t->ci->pcell_size * sizeof(integertime_t),
MPI_BYTE, t->ci->nodeID, t->flags, MPI_COMM_WORLD,
&t->req);
t->buff = malloc(sizeof(struct pcell_step) * t->ci->pcell_size);
err = MPI_Irecv(
t->buff, t->ci->pcell_size * sizeof(struct pcell_step), MPI_BYTE,
t->ci->nodeID, t->flags, MPI_COMM_WORLD, &t->req);
} else if (t->subtype == task_subtype_xv ||
t->subtype == task_subtype_rho ||
t->subtype == task_subtype_gradient) {
......@@ -1309,11 +1309,11 @@ void scheduler_enqueue(struct scheduler *s, struct task *t) {
case task_type_send:
#ifdef WITH_MPI
if (t->subtype == task_subtype_tend) {
t->buff = malloc(sizeof(integertime_t) * t->ci->pcell_size);
cell_pack_ti_ends(t->ci, t->buff);
err = MPI_Isend(t->buff, t->ci->pcell_size * sizeof(integertime_t),