Commit 35337e9f authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Split the drift task between the part and gparts

parent 5d8e19b5
......@@ -29,25 +29,48 @@
#include "timeline.h"
/**
* @brief Check that a cell been drifted to the current time.
* @brief Check that the #part in a #cell have been drifted to the current time.
*
* @param c The #cell.
* @param e The #engine containing information about the current time.
* @return 1 if the #cell has been drifted to the current time, 0 otherwise.
*/
__attribute__((always_inline)) INLINE static int cell_is_drifted(
__attribute__((always_inline)) INLINE static int cell_are_part_drifted(
const struct cell *c, const struct engine *e) {
#ifdef SWIFT_DEBUG_CHECKS
if (c->ti_old > e->ti_current)
if (c->ti_old_part > e->ti_current)
error(
"Cell has been drifted too far forward in time! c->ti_old=%lld (t=%e) "
"and e->ti_current=%lld (t=%e)",
c->ti_old, c->ti_old * e->timeBase, e->ti_current,
c->ti_old_part, c->ti_old_part * e->timeBase, e->ti_current,
e->ti_current * e->timeBase);
#endif
return (c->ti_old == e->ti_current);
return (c->ti_old_part == e->ti_current);
}
/**
* @brief Check that the #gpart in a #cell have been drifted to the current
* time.
*
* @param c The #cell.
* @param e The #engine containing information about the current time.
* @return 1 if the #cell has been drifted to the current time, 0 otherwise.
*/
__attribute__((always_inline)) INLINE static int cell_are_gpart_drifted(
const struct cell *c, const struct engine *e) {
#ifdef SWIFT_DEBUG_CHECKS
if (c->ti_old_gpart > e->ti_current)
error(
"Cell has been drifted too far forward in time! c->ti_old=%lld (t=%e) "
"and e->ti_current=%lld (t=%e)",
c->ti_old_gpart, c->ti_old_gpart * e->timeBase, e->ti_current,
e->ti_current * e->timeBase);
#endif
return (c->ti_old_gpart == e->ti_current);
}
/* Are cells / particles active for regular tasks ? */
......
......@@ -99,7 +99,8 @@ int cell_unpack(struct pcell *pc, struct cell *c, struct space *s) {
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 = pc->ti_old;
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;
......@@ -128,7 +129,8 @@ int cell_unpack(struct pcell *pc, struct cell *c, struct space *s) {
if (k & 1) temp->loc[2] += temp->width[2];
temp->depth = c->depth + 1;
temp->split = 0;
temp->dx_max = 0.f;
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;
......@@ -239,7 +241,8 @@ int cell_pack(struct cell *c, struct pcell *pc) {
pc->h_max = c->h_max;
pc->ti_end_min = c->ti_end_min;
pc->ti_end_max = c->ti_end_max;
pc->ti_old = c->ti_old;
pc->ti_old_part = c->ti_old_part;
pc->ti_old_gpart = c->ti_old_gpart;
pc->count = c->count;
pc->gcount = c->gcount;
pc->scount = c->scount;
......@@ -1018,7 +1021,7 @@ void cell_clean_links(struct cell *c, void *data) {
}
/**
* @brief Checks that the particles in a cell are at the
* @brief Checks that the #part in a cell are at the
* current point in time
*
* Calls error() if the cell is not at the current time.
......@@ -1026,7 +1029,7 @@ void cell_clean_links(struct cell *c, void *data) {
* @param c Cell to act upon
* @param data The current time on the integer time-line
*/
void cell_check_particle_drift_point(struct cell *c, void *data) {
void cell_check_part_drift_point(struct cell *c, void *data) {
#ifdef SWIFT_DEBUG_CHECKS
......@@ -1035,14 +1038,40 @@ void cell_check_particle_drift_point(struct cell *c, void *data) {
/* Only check local cells */
if (c->nodeID != engine_rank) return;
if (c->ti_old != ti_drift)
error("Cell in an incorrect time-zone! c->ti_old=%lld ti_drift=%lld",
c->ti_old, ti_drift);
if (c->ti_old_part != ti_drift)
error("Cell in an incorrect time-zone! c->ti_old_part=%lld ti_drift=%lld",
c->ti_old_part, ti_drift);
for (int i = 0; i < c->count; ++i)
if (c->parts[i].ti_drift != ti_drift)
error("part in an incorrect time-zone! p->ti_drift=%lld ti_drift=%lld",
c->parts[i].ti_drift, ti_drift);
#else
error("Calling debugging code without debugging flag activated.");
#endif
}
/**
* @brief Checks that the #gpart and #spart in a cell are 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_gpart_drift_point(struct cell *c, void *data) {
#ifdef SWIFT_DEBUG_CHECKS
const integertime_t ti_drift = *(integertime_t *)data;
/* Only check local cells */
if (c->nodeID != engine_rank) return;
if (c->ti_old_gpart != ti_drift)
error("Cell in an incorrect time-zone! c->ti_old_gpart=%lld ti_drift=%lld",
c->ti_old_gpart, ti_drift);
for (int i = 0; i < c->gcount; ++i)
if (c->gparts[i].ti_drift != ti_drift)
......@@ -1327,7 +1356,7 @@ int cell_unskip_tasks(struct cell *c, struct scheduler *s) {
error("bad flags in sort task.");
#endif
scheduler_activate(s, ci->sorts);
if (ci->nodeID == engine_rank) scheduler_activate(s, ci->drift);
if (ci->nodeID == engine_rank) scheduler_activate(s, ci->drift_part);
}
if (!(cj->sorted & (1 << t->flags))) {
#ifdef SWIFT_DEBUG_CHECKS
......@@ -1335,7 +1364,7 @@ int cell_unskip_tasks(struct cell *c, struct scheduler *s) {
error("bad flags in sort task.");
#endif
scheduler_activate(s, cj->sorts);
if (cj->nodeID == engine_rank) scheduler_activate(s, cj->drift);
if (cj->nodeID == engine_rank) scheduler_activate(s, cj->drift_part);
}
}
......@@ -1344,7 +1373,8 @@ int cell_unskip_tasks(struct cell *c, struct scheduler *s) {
/* Check whether there was too much particle motion, i.e. the
cell neighbour conditions were violated. */
if (max(ci->h_max, cj->h_max) + ci->dx_max + cj->dx_max > cj->dmin)
if (max(ci->h_max, cj->h_max) + ci->dx_max_part + cj->dx_max_part >
cj->dmin)
rebuild = 1;
#ifdef WITH_MPI
......@@ -1367,11 +1397,11 @@ int cell_unskip_tasks(struct cell *c, struct scheduler *s) {
scheduler_activate(s, l->t);
/* Drift both cells, the foreign one at the level which it is sent. */
if (l->t->ci->drift)
scheduler_activate(s, l->t->ci->drift);
if (l->t->ci->drift_part)
scheduler_activate(s, l->t->ci->drift_part);
else
error("Drift task missing !");
if (t->type == task_type_pair) scheduler_activate(s, cj->drift);
if (t->type == task_type_pair) scheduler_activate(s, cj->drift_part);
if (cell_is_active(cj, e)) {
for (l = cj->send_rho; l != NULL && l->t->cj->nodeID != ci->nodeID;
......@@ -1405,11 +1435,11 @@ int cell_unskip_tasks(struct cell *c, struct scheduler *s) {
scheduler_activate(s, l->t);
/* Drift both cells, the foreign one at the level which it is sent. */
if (l->t->ci->drift)
scheduler_activate(s, l->t->ci->drift);
if (l->t->ci->drift_part)
scheduler_activate(s, l->t->ci->drift_part);
else
error("Drift task missing !");
if (t->type == task_type_pair) scheduler_activate(s, ci->drift);
if (t->type == task_type_pair) scheduler_activate(s, ci->drift_part);
if (cell_is_active(ci, e)) {
for (l = ci->send_rho; l != NULL && l->t->cj->nodeID != cj->nodeID;
......@@ -1425,13 +1455,13 @@ int cell_unskip_tasks(struct cell *c, struct scheduler *s) {
scheduler_activate(s, l->t);
}
} else if (t->type == task_type_pair) {
scheduler_activate(s, ci->drift);
scheduler_activate(s, cj->drift);
scheduler_activate(s, ci->drift_part);
scheduler_activate(s, cj->drift_part);
}
#else
if (t->type == task_type_pair) {
scheduler_activate(s, ci->drift);
scheduler_activate(s, cj->drift);
scheduler_activate(s, ci->drift_part);
scheduler_activate(s, cj->drift_part);
}
#endif
}
......@@ -1447,7 +1477,8 @@ int cell_unskip_tasks(struct cell *c, struct scheduler *s) {
if (c->extra_ghost != NULL) scheduler_activate(s, c->extra_ghost);
if (c->ghost != NULL) scheduler_activate(s, c->ghost);
if (c->init_grav != NULL) scheduler_activate(s, c->init_grav);
if (c->drift != NULL) scheduler_activate(s, c->drift);
if (c->drift_part != NULL) scheduler_activate(s, c->drift_part);
if (c->drift_gpart != NULL) scheduler_activate(s, c->drift_gpart);
if (c->kick1 != NULL) scheduler_activate(s, c->kick1);
if (c->kick2 != NULL) scheduler_activate(s, c->kick2);
if (c->timestep != NULL) scheduler_activate(s, c->timestep);
......@@ -1480,30 +1511,28 @@ void cell_set_super(struct cell *c, struct cell *super) {
}
/**
* @brief Recursively drifts particles of all kinds in a cell hierarchy.
* @brief Recursively drifts the #part in a cell hierarchy.
*
* @param c The #cell.
* @param e The #engine (to get ti_current).
*/
void cell_drift_particles(struct cell *c, const struct engine *e) {
void cell_drift_part(struct cell *c, const struct engine *e) {
const float hydro_h_max = e->hydro_properties->h_max;
const double timeBase = e->timeBase;
const integertime_t ti_old = c->ti_old;
const integertime_t ti_old_part = c->ti_old_part;
const integertime_t ti_current = e->ti_current;
struct part *const parts = c->parts;
struct xpart *const xparts = c->xparts;
struct gpart *const gparts = c->gparts;
struct spart *const sparts = c->sparts;
/* Drift from the last time the cell was drifted to the current time */
const double dt = (ti_current - ti_old) * timeBase;
const double dt = (ti_current - ti_old_part) * timeBase;
float dx_max = 0.f, dx2_max = 0.f;
float dx_max_sort = 0.0f, dx2_max_sort = 0.f;
float cell_h_max = 0.f;
/* Check that we are actually going to move forward. */
if (ti_current < ti_old) error("Attempt to drift to the past");
if (ti_current < ti_old_part) error("Attempt to drift to the past");
/* Are we not in a leaf ? */
if (c->split) {
......@@ -1514,37 +1543,15 @@ void cell_drift_particles(struct cell *c, const struct engine *e) {
struct cell *cp = c->progeny[k];
/* Collect */
cell_drift_particles(cp, e);
cell_drift_part(cp, e);
/* Update */
dx_max = max(dx_max, cp->dx_max);
dx_max = max(dx_max, cp->dx_max_part);
dx_max_sort = max(dx_max_sort, cp->dx_max_sort);
cell_h_max = max(cell_h_max, cp->h_max);
}
} else if (ti_current > ti_old) {
/* Loop over all the g-particles in the cell */
const size_t nr_gparts = c->gcount;
for (size_t k = 0; k < nr_gparts; k++) {
/* Get a handle on the gpart. */
struct gpart *const gp = &gparts[k];
/* Drift... */
drift_gpart(gp, dt, timeBase, ti_old, ti_current);
/* Compute (square of) motion since last cell construction */
const float dx2 = gp->x_diff[0] * gp->x_diff[0] +
gp->x_diff[1] * gp->x_diff[1] +
gp->x_diff[2] * gp->x_diff[2];
dx2_max = max(dx2_max, dx2);
/* Init gravity force fields. */
if (gpart_is_active(gp, e)) {
gravity_init_gpart(gp);
}
}
} else if (ti_current > ti_old_part) {
/* Loop over all the gas particles in the cell */
const size_t nr_parts = c->count;
......@@ -1555,7 +1562,7 @@ void cell_drift_particles(struct cell *c, const struct engine *e) {
struct xpart *const xp = &xparts[k];
/* Drift... */
drift_part(p, xp, dt, timeBase, ti_old, ti_current);
drift_part(p, xp, dt, timeBase, ti_old_part, ti_current);
/* Limit h to within the allowed range */
p->h = min(p->h, hydro_h_max);
......@@ -1579,6 +1586,86 @@ void cell_drift_particles(struct cell *c, const struct engine *e) {
}
}
/* Now, get the maximal particle motion from its square */
dx_max = sqrtf(dx2_max);
dx_max_sort = sqrtf(dx2_max_sort);
} else {
cell_h_max = c->h_max;
dx_max = c->dx_max_part;
dx_max_sort = c->dx_max_sort;
}
/* Store the values */
c->h_max = cell_h_max;
c->dx_max_part = dx_max;
c->dx_max_sort = dx_max_sort;
/* Update the time of the last drift */
c->ti_old_part = ti_current;
}
/**
* @brief Recursively drifts the #gpart in a cell hierarchy.
*
* @param c The #cell.
* @param e The #engine (to get ti_current).
*/
void cell_drift_gpart(struct cell *c, const struct engine *e) {
const double timeBase = e->timeBase;
const integertime_t ti_old_gpart = c->ti_old_gpart;
const integertime_t ti_current = e->ti_current;
struct gpart *const gparts = c->gparts;
struct spart *const sparts = c->sparts;
/* Drift from the last time the cell was drifted to the current time */
const double dt = (ti_current - ti_old_gpart) * timeBase;
float dx_max = 0.f, dx2_max = 0.f;
/* Check that we are actually going to move forward. */
if (ti_current < ti_old_gpart) error("Attempt to drift to the past");
/* Are we not in a leaf ? */
if (c->split) {
/* Loop over the progeny and collect their data. */
for (int k = 0; k < 8; k++)
if (c->progeny[k] != NULL) {
struct cell *cp = c->progeny[k];
/* Recurse */
cell_drift_gpart(cp, e);
/* Update */
dx_max = max(dx_max, cp->dx_max_gpart);
}
} else if (ti_current > ti_old_gpart) {
/* Loop over all the g-particles in the cell */
const size_t nr_gparts = c->gcount;
for (size_t k = 0; k < nr_gparts; k++) {
/* Get a handle on the gpart. */
struct gpart *const gp = &gparts[k];
/* Drift... */
drift_gpart(gp, dt, timeBase, ti_old_gpart, ti_current);
/* Compute (square of) motion since last cell construction */
const float dx2 = gp->x_diff[0] * gp->x_diff[0] +
gp->x_diff[1] * gp->x_diff[1] +
gp->x_diff[2] * gp->x_diff[2];
dx2_max = max(dx2_max, dx2);
/* Init gravity force fields. */
if (gpart_is_active(gp, e)) {
gravity_init_gpart(gp);
}
}
/* Loop over all the star particles in the cell */
const size_t nr_sparts = c->scount;
for (size_t k = 0; k < nr_sparts; k++) {
......@@ -1587,29 +1674,24 @@ void cell_drift_particles(struct cell *c, const struct engine *e) {
struct spart *const sp = &sparts[k];
/* Drift... */
drift_spart(sp, dt, timeBase, ti_old, ti_current);
drift_spart(sp, dt, timeBase, ti_old_gpart, ti_current);
/* Note: no need to compute dx_max as all spart have a gpart */
}
/* Now, get the maximal particle motion from its square */
dx_max = sqrtf(dx2_max);
dx_max_sort = sqrtf(dx2_max_sort);
} else {
cell_h_max = c->h_max;
dx_max = c->dx_max;
dx_max_sort = c->dx_max_sort;
dx_max = c->dx_max_gpart;
}
/* Store the values */
c->h_max = cell_h_max;
c->dx_max = dx_max;
c->dx_max_sort = dx_max_sort;
c->dx_max_gpart = dx_max;
/* Update the time of the last drift */
c->ti_old = ti_current;
c->ti_old_gpart = ti_current;
}
/**
......
......@@ -74,7 +74,7 @@ struct pcell {
/* Stats on this cell's particles. */
double h_max;
integertime_t ti_end_min, ti_end_max, ti_beg_max, ti_old;
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;
......@@ -157,8 +157,11 @@ struct cell {
/*! The extra ghost task for complex hydro schemes */
struct task *extra_ghost;
/*! The drift task */
struct task *drift;
/*! The drift task for parts */
struct task *drift_part;
/*! The drift task for gparts */
struct task *drift_gpart;
/*! The first kick task */
struct task *kick1;
......@@ -230,24 +233,30 @@ struct cell {
/*! Maximum beginning of (integer) time step in this cell. */
integertime_t ti_beg_max;
/*! Last (integer) time the cell's particle was drifted forward in time. */
integertime_t ti_old;
/*! Last (integer) time the cell's sort arrays were updated. */
integertime_t ti_sort;
/*! Last (integer) time the cell's part were drifted forward in time. */
integertime_t ti_old_part;
/*! Last (integer) time the cell's gpart were drifted forward in time. */
integertime_t ti_old_gpart;
/*! Last (integer) time the cell's multipole was drifted forward in time. */
integertime_t ti_old_multipole;
/*! Minimum dimension, i.e. smallest edge of this cell (min(width)). */
float dmin;
/*! Maximum particle movement in this cell since last construction. */
float dx_max;
/*! Maximum particle movement in this cell since the last sort. */
float dx_max_sort;
/*! Maximum part movement in this cell since last construction. */
float dx_max_part;
/*! Maximum gpart movement in this cell since last construction. */
float dx_max_gpart;
/*! Nr of #part in this cell. */
int count;
......@@ -354,13 +363,15 @@ void cell_clean_links(struct cell *c, void *data);
void cell_make_multipoles(struct cell *c, integertime_t ti_current);
void cell_check_multipole(struct cell *c, void *data);
void cell_clean(struct cell *c);
void cell_check_particle_drift_point(struct cell *c, void *data);
void cell_check_part_drift_point(struct cell *c, void *data);
void cell_check_gpart_drift_point(struct cell *c, void *data);
void cell_check_multipole_drift_point(struct cell *c, void *data);
void cell_reset_task_counters(struct cell *c);
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);
void cell_drift_particles(struct cell *c, const struct engine *e);
void cell_drift_part(struct cell *c, const struct engine *e);
void cell_drift_gpart(struct cell *c, const struct engine *e);
void cell_drift_multipole(struct cell *c, const struct engine *e);
void cell_drift_all_multipoles(struct cell *c, const struct engine *e);
void cell_check_timesteps(struct cell *c);
......
......@@ -259,8 +259,8 @@ int checkCellhdxmax(const struct cell *c, int *depth) {
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);
if (c->dx_max_part != dx_max) {
message("%d Inconsistent dx_max: %f != %f", *depth, c->dx_max_part, dx_max);
message("location: %f %f %f", c->loc[0], c->loc[1], c->loc[2]);
result = 0;
}
......
......@@ -1043,10 +1043,10 @@ void engine_addtasks_send(struct engine *e, struct cell *ci, struct cell *cj,
#endif
/* Drift before you send */
if (ci->drift == NULL)
ci->drift = scheduler_addtask(s, task_type_drift, task_subtype_none, 0,
0, ci, NULL);
scheduler_addunlock(s, ci->drift, t_xv);
if (ci->drift_part == NULL)
ci->drift_part = scheduler_addtask(s, task_type_drift_part,
task_subtype_none, 0, 0, ci, NULL);
scheduler_addunlock(s, ci->drift_part, t_xv);
/* The super-cell's timestep task should unlock the send_ti task. */
scheduler_addunlock(s, ci->super->timestep, t_ti);
......@@ -1816,10 +1816,11 @@ void engine_count_and_link_tasks(struct engine *e) {
}
/* Link drift tasks to the next-higher drift task. */
else if (t->type == task_type_drift) {
else if (t->type == task_type_drift_part) {
struct cell *finger = ci->parent;
while (finger != NULL && finger->drift == NULL) finger = finger->parent;
if (finger != NULL) scheduler_addunlock(sched, t, finger->drift);
while (finger != NULL && finger->drift_part == NULL)
finger = finger->parent;
if (finger != NULL) scheduler_addunlock(sched, t, finger->drift_part);
}
/* Link self tasks to cells. */
......@@ -1910,7 +1911,7 @@ static inline void engine_make_external_gravity_dependencies(
struct scheduler *sched, struct task *gravity, struct cell *c) {
/* init --> external gravity --> kick */
scheduler_addunlock(sched, c->drift, gravity);
scheduler_addunlock(sched, c->drift_gpart, gravity);
scheduler_addunlock(sched, gravity, c->super->kick2);
}
......@@ -2076,14 +2077,14 @@ void engine_make_extra_hydroloop_tasks(struct engine *e) {
/* Sort tasks depend on the drift of the cell. */
if (t->type == task_type_sort && t->ci->nodeID == engine_rank) {
scheduler_addunlock(sched, t->ci->drift, t);
scheduler_addunlock(sched, t->ci->drift_part, t);
}
/* Self-interaction? */
else if (t->type == task_type_self && t->subtype == task_subtype_density) {
/* Make all density tasks depend on the drift. */
scheduler_addunlock(sched, t->ci->drift, t);
scheduler_addunlock(sched, t->ci->drift_part, t);
#ifdef EXTRA_HYDRO_LOOP
/* Start by constructing the task for the second and third hydro loop */
......@@ -2119,9 +2120,9 @@ void engine_make_extra_hydroloop_tasks(struct engine *e) {
/* Make all density tasks depend on the drift. */
if (t->ci->nodeID == engine_rank)
scheduler_addunlock(sched, t->ci->drift, t);
scheduler_addunlock(sched, t->ci->drift_part, t);
if (t->cj->nodeID == engine_rank)
scheduler_addunlock(sched, t->cj->drift, t);
scheduler_addunlock(sched, t->cj->drift_part, t);
#ifdef EXTRA_HYDRO_LOOP
/* Start by constructing the task for the second and third hydro loop */
......@@ -2478,7 +2479,8 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
if (t->subtype != task_subtype_density) continue;
/* Too much particle movement? */
if (max(ci->h_max, cj->h_max) + ci->dx_max + cj->dx_max > cj->dmin)
if (max(ci->h_max, cj->h_max) + ci->dx_max_part + cj->dx_max_part >
cj->dmin)
*rebuild_space = 1;
/* Set the correct sorting flags */
......@@ -2499,7 +2501,7 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
error("bad flags in sort task.");
#endif
scheduler_activate(s, ci->sorts);
if (ci->nodeID == engine_rank) scheduler_activate(s, ci->drift);
if (ci->nodeID == engine_rank) scheduler_activate(s, ci->drift_part);
}
if (!(cj->sorted & (1 << t->flags))) {
#ifdef SWIFT_DEBUG_CHECKS
......@@ -2507,7 +2509,7 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
error("bad flags in sort task.");
#endif
scheduler_activate(s, cj->sorts);
if (cj->nodeID == engine_rank) scheduler_activate(s, cj->drift);
if (cj->nodeID == engine_rank) scheduler_activate(s, cj->drift_part);
}
}
......@@ -2531,11 +2533,11 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
scheduler_activate(s, l->t);
/* Drift both cells, the foreign one at the level which it is sent. */
if (l->t->ci->drift)
scheduler