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

Merge branch 'timestep_order' into 'master'

Timestep order

This is my last batch of changes regarding time integration.
 
 - Additional checks here and there that the particles are synchronized. Especially gpart.
 - Re-arranged the time-step sequence and simplified the logic in engine_step().
 - Moved the decision of repartitioning from the main() to engine_step().

The re-arrangement of the time-step sequence allows for perfect energy conservation even when we rebuild, repartition or dump snapshots. Previously that was not the case as a spurious drift had to be introduced. 

See merge request !314
parents c707ab6b f181430f
......@@ -91,6 +91,9 @@ for i in range(402):
E_tot_snap[i] = E_kin_snap[i] + E_pot_snap[i]
Lz_snap[i] = np.sum(Lz)
print "Starting energy:", E_kin_stats[0], E_pot_stats[0], E_tot_stats[0]
print "Ending energy:", E_kin_stats[-1], E_pot_stats[-1], E_tot_stats[-1]
# Plot energy evolution
figure()
plot(time_stats, E_kin_stats, "r-", lw=0.5, label="Kinetic energy")
......
......@@ -333,10 +333,10 @@ int main(int argc, char *argv[]) {
MPI_Bcast(params, sizeof(struct swift_params), MPI_BYTE, 0, MPI_COMM_WORLD);
#endif
/* Prepare the domain decomposition scheme */
/* Prepare the domain decomposition scheme */
enum repartition_type reparttype = REPART_NONE;
#ifdef WITH_MPI
struct partition initial_partition;
enum repartition_type reparttype;
partition_init(&initial_partition, &reparttype, params, nr_nodes);
/* Let's report what we did */
......@@ -514,8 +514,9 @@ int main(int argc, char *argv[]) {
if (myrank == 0) clocks_gettime(&tic);
struct engine e;
engine_init(&e, &s, params, nr_nodes, myrank, nr_threads, with_aff,
engine_policies, talking, &us, &prog_const, &hydro_properties,
&gravity_properties, &potential, &cooling_func, &sourceterms);
engine_policies, talking, reparttype, &us, &prog_const,
&hydro_properties, &gravity_properties, &potential, &cooling_func,
&sourceterms);
if (myrank == 0) {
clocks_gettime(&toc);
message("engine_init took %.3f %s.", clocks_diff(&tic, &toc),
......@@ -589,11 +590,6 @@ int main(int argc, char *argv[]) {
/* Main simulation loop */
for (int j = 0; !engine_is_done(&e) && e.step != nsteps; j++) {
/* Repartition the space amongst the nodes? */
#ifdef WITH_MPI
if (j % 100 == 2) e.forcerepart = reparttype;
#endif
/* Reset timers */
timers_reset(timers_mask_all);
......@@ -695,6 +691,7 @@ int main(int argc, char *argv[]) {
#endif
/* Write final output. */
engine_drift_all(&e);
engine_dump_snapshot(&e);
#ifdef WITH_MPI
......
......@@ -50,8 +50,10 @@ __attribute__((always_inline)) INLINE static int cell_is_drifted(
return (c->ti_old == e->ti_current);
}
/* Are cells / particles active for regular tasks ? */
/**
* @brief Does a cell contain any active particle ?
* @brief Does a cell contain any particle finishing their time-step now ?
*
* @param c The #cell.
* @param e The #engine containing information about the current time.
......@@ -73,7 +75,7 @@ __attribute__((always_inline)) INLINE static int cell_is_active(
}
/**
* @brief Are *all* particles in a cell active ?
* @brief Are *all* particles in a cell finishing their time-step now ?
*
* @param c The #cell.
* @param e The #engine containing information about the current time.
......@@ -94,7 +96,7 @@ __attribute__((always_inline)) INLINE static int cell_is_all_active(
}
/**
* @brief Is this particle active ?
* @brief Is this particle finishing its time-step now ?
*
* @param p The #part.
* @param e The #engine containing information about the current time.
......@@ -118,7 +120,7 @@ __attribute__((always_inline)) INLINE static int part_is_active(
}
/**
* @brief Is this g-particle active ?
* @brief Is this g-particle finishing its time-step now ?
*
* @param gp The #gpart.
* @param e The #engine containing information about the current time.
......@@ -142,7 +144,7 @@ __attribute__((always_inline)) INLINE static int gpart_is_active(
}
/**
* @brief Is this s-particle active ?
* @brief Is this s-particle finishing its time-step now ?
*
* @param sp The #spart.
* @param e The #engine containing information about the current time.
......@@ -157,7 +159,7 @@ __attribute__((always_inline)) INLINE static int spart_is_active(
#ifdef SWIFT_DEBUG_CHECKS
if (ti_end < ti_current)
error(
"s-particle in an impossible time-zone! gp->ti_end=%lld "
"s-particle in an impossible time-zone! sp->ti_end=%lld "
"e->ti_current=%lld",
ti_end, ti_current);
#endif
......@@ -165,4 +167,102 @@ __attribute__((always_inline)) INLINE static int spart_is_active(
return (ti_end == ti_current);
}
/* Are cells / particles active for kick1 tasks ? */
/**
* @brief Does a cell contain any particle starting their time-step now ?
*
* @param c The #cell.
* @param e The #engine containing information about the current time.
* @return 1 if the #cell contains at least an active particle, 0 otherwise.
*/
__attribute__((always_inline)) INLINE static int cell_is_starting(
const struct cell *c, const struct engine *e) {
#ifdef SWIFT_DEBUG_CHECKS
if (c->ti_beg_max > e->ti_current)
error(
"cell in an impossible time-zone! c->ti_beg_max=%lld (t=%e) and "
"e->ti_current=%lld (t=%e)",
c->ti_beg_max, c->ti_beg_max * e->timeBase, e->ti_current,
e->ti_current * e->timeBase);
#endif
return (c->ti_beg_max == e->ti_current);
}
/**
* @brief Is this particle starting its time-step now ?
*
* @param p The #part.
* @param e The #engine containing information about the current time.
* @return 1 if the #part is active, 0 otherwise.
*/
__attribute__((always_inline)) INLINE static int part_is_starting(
const struct part *p, const struct engine *e) {
const integertime_t ti_current = e->ti_current;
const integertime_t ti_beg =
get_integer_time_begin(ti_current + 1, p->time_bin);
#ifdef SWIFT_DEBUG_CHECKS
if (ti_beg > ti_current)
error(
"particle in an impossible time-zone! p->ti_beg=%lld "
"e->ti_current=%lld",
ti_beg, ti_current);
#endif
return (ti_beg == ti_current);
}
/**
* @brief Is this g-particle starting its time-step now ?
*
* @param gp The #gpart.
* @param e The #engine containing information about the current time.
* @return 1 if the #gpart is active, 0 otherwise.
*/
__attribute__((always_inline)) INLINE static int gpart_is_starting(
const struct gpart *gp, const struct engine *e) {
const integertime_t ti_current = e->ti_current;
const integertime_t ti_beg =
get_integer_time_begin(ti_current + 1, gp->time_bin);
#ifdef SWIFT_DEBUG_CHECKS
if (ti_beg > ti_current)
error(
"g-particle in an impossible time-zone! gp->ti_beg=%lld "
"e->ti_current=%lld",
ti_beg, ti_current);
#endif
return (ti_beg == ti_current);
}
/**
* @brief Is this s-particle starting its time-step now ?
*
* @param sp The #spart.
* @param e The #engine containing information about the current time.
* @return 1 if the #spart is active, 0 otherwise.
*/
__attribute__((always_inline)) INLINE static int spart_is_starting(
const struct spart *sp, const struct engine *e) {
const integertime_t ti_current = e->ti_current;
const integertime_t ti_beg =
get_integer_time_begin(ti_current + 1, sp->time_bin);
#ifdef SWIFT_DEBUG_CHECKS
if (ti_beg > ti_current)
error(
"s-particle in an impossible time-zone! sp->ti_beg=%lld "
"e->ti_current=%lld",
ti_beg, ti_current);
#endif
return (ti_beg == ti_current);
}
#endif /* SWIFT_ACTIVE_H */
......@@ -943,11 +943,34 @@ void cell_clean_links(struct cell *c, void *data) {
*/
void cell_check_drift_point(struct cell *c, void *data) {
integertime_t ti_current = *(integertime_t *)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 != ti_drift)
error("Cell in an incorrect time-zone! c->ti_old=%lld ti_drift=%lld",
c->ti_old, 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);
for (int i = 0; i < c->gcount; ++i)
if (c->gparts[i].ti_drift != ti_drift)
error("g-part in an incorrect time-zone! gp->ti_drift=%lld ti_drift=%lld",
c->gparts[i].ti_drift, ti_drift);
if (c->ti_old != ti_current && c->nodeID == engine_rank)
error("Cell in an incorrect time-zone! c->ti_old=%lld ti_current=%lld",
c->ti_old, ti_current);
for (int i = 0; i < c->scount; ++i)
if (c->sparts[i].ti_drift != ti_drift)
error("s-part in an incorrect time-zone! sp->ti_drift=%lld ti_drift=%lld",
c->sparts[i].ti_drift, ti_drift);
#else
error("Calling debugging code without debugging flag activated.");
#endif
}
/**
......@@ -1095,6 +1118,8 @@ int cell_unskip_tasks(struct cell *c, struct scheduler *s) {
struct engine *e = s->space->e;
#endif
int rebuild = 0;
/* Un-skip the density tasks involved with this cell. */
for (struct link *l = c->density; l != NULL; l = l->next) {
struct task *t = l->t;
......@@ -1120,7 +1145,7 @@ int cell_unskip_tasks(struct cell *c, struct scheduler *s) {
(max(ci->h_max, cj->h_max) + ci->dx_max + cj->dx_max > cj->dmin ||
ci->dx_max > space_maxreldx * ci->h_max ||
cj->dx_max > space_maxreldx * cj->h_max))
return 1;
rebuild = 1;
#ifdef WITH_MPI
/* Activate the send/recv flags. */
......@@ -1217,7 +1242,7 @@ int cell_unskip_tasks(struct cell *c, struct scheduler *s) {
if (c->cooling != NULL) scheduler_activate(s, c->cooling);
if (c->sourceterms != NULL) scheduler_activate(s, c->sourceterms);
return 0;
return rebuild;
}
/**
......@@ -1364,5 +1389,7 @@ void cell_check_timesteps(struct cell *c) {
if (c->parts[i].time_bin == 0)
error("Particle without assigned time-bin");
}
#else
error("Calling debugging code without debugging flag activated.");
#endif
}
......@@ -74,7 +74,7 @@ struct pcell {
/* Stats on this cell's particles. */
double h_max;
integertime_t ti_end_min, ti_end_max, ti_old;
integertime_t ti_end_min, ti_end_max, ti_beg_max, ti_old;
/* Number of particles in this cell. */
int count, gcount, scount;
......@@ -227,6 +227,9 @@ struct cell {
/*! Maximum end of (integer) time step in this cell. */
integertime_t ti_end_max;
/*! Maximum beginning of (integer) time step in this cell. */
integertime_t ti_beg_max;
/*! Last (integer) time the cell's content was drifted forward in time. */
integertime_t ti_old;
......
......@@ -41,6 +41,18 @@
__attribute__((always_inline)) INLINE static void drift_gpart(
struct gpart *restrict gp, float dt, double timeBase, integertime_t ti_old,
integertime_t ti_current) {
#ifdef SWIFT_DEBUG_CHECKS
if (gp->ti_drift != ti_old)
error(
"g-particle has not been drifted to the current time "
"gp->ti_drift=%lld, "
"c->ti_old=%lld, ti_current=%lld",
gp->ti_drift, ti_old, ti_current);
gp->ti_drift = ti_current;
#endif
/* Drift... */
gp->x[0] += gp->v_full[0] * dt;
gp->x[1] += gp->v_full[1] * dt;
......@@ -69,7 +81,7 @@ __attribute__((always_inline)) INLINE static void drift_part(
#ifdef SWIFT_DEBUG_CHECKS
if (p->ti_drift != ti_old)
error(
"Particle has not been drifted to the current time p->ti_drift=%lld, "
"particle has not been drifted to the current time p->ti_drift=%lld, "
"c->ti_old=%lld, ti_current=%lld",
p->ti_drift, ti_old, ti_current);
......@@ -108,6 +120,17 @@ __attribute__((always_inline)) INLINE static void drift_spart(
struct spart *restrict sp, float dt, double timeBase, integertime_t ti_old,
integertime_t ti_current) {
#ifdef SWIFT_DEBUG_CHECKS
if (sp->ti_drift != ti_old)
error(
"s-particle has not been drifted to the current time "
"sp->ti_drift=%lld, "
"c->ti_old=%lld, ti_current=%lld",
sp->ti_drift, ti_old, ti_current);
sp->ti_drift = ti_current;
#endif
/* Drift... */
sp->x[0] += sp->v[0] * dt;
sp->x[1] += sp->v[1] * dt;
......
......@@ -157,12 +157,12 @@ void engine_make_hierarchical_tasks(struct engine *e, struct cell *c) {
0, 0, c, NULL, 0);
scheduler_addunlock(s, c->kick2, c->timestep);
scheduler_addunlock(s, c->timestep, c->kick1);
/* Add the drift task and its dependencies. */
c->drift = scheduler_addtask(s, task_type_drift, task_subtype_none, 0, 0,
c, NULL, 0);
scheduler_addunlock(s, c->kick1, c->drift);
scheduler_addunlock(s, c->drift, c->init);
/* Generate the ghost task. */
......@@ -822,19 +822,18 @@ void engine_repartition(struct engine *e) {
fflush(stdout);
/* Check that all cells have been drifted to the current time */
space_check_drift_point(e->s, e->ti_current);
space_check_drift_point(e->s, e->ti_old);
#endif
/* Clear the repartition flag. */
enum repartition_type reparttype = e->forcerepart;
e->forcerepart = REPART_NONE;
e->forcerepart = 0;
/* Nothing to do if only using a single node. Also avoids METIS
* bug that doesn't handle this case well. */
if (e->nr_nodes == 1) return;
/* Do the repartitioning. */
partition_repartition(reparttype, e->nodeID, e->nr_nodes, e->s,
partition_repartition(e->reparttype, e->nodeID, e->nr_nodes, e->s,
e->sched.tasks, e->sched.nr_tasks);
/* Now comes the tricky part: Exchange particles between all nodes.
......@@ -2592,8 +2591,15 @@ void engine_rebuild(struct engine *e) {
if (engine_marktasks(e))
error("engine_marktasks failed after space_rebuild.");
/* Print the status of the system */
if (e->verbose) engine_print_task_counts(e);
/* Print the status of the system */
// if (e->verbose) engine_print_task_counts(e);
#ifdef SWIFT_DEBUG_CHECKS
/* Check that all cells have been drifted to the current time.
* That can include cells that have not
* previously been active on this rank. */
space_check_drift_point(e->s, e->ti_old);
#endif
if (e->verbose)
message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
......@@ -2604,48 +2610,28 @@ void engine_rebuild(struct engine *e) {
* @brief Prepare the #engine by re-building the cells and tasks.
*
* @param e The #engine to prepare.
* @param drift_all Whether to drift particles before rebuilding or not. Will
* not be necessary if all particles have already been
* drifted (before repartitioning for instance).
* @param postrepart If we have just repartitioned, if so we need to defer the
* skip until after the rebuild and not check the if all
* cells have been drifted.
*/
void engine_prepare(struct engine *e, int drift_all, int postrepart) {
void engine_prepare(struct engine *e) {
TIMER_TIC;
/* Unskip active tasks and check for rebuild */
if (!postrepart) engine_unskip(e);
/* Run through the tasks and mark as skip or not. */
int rebuild = e->forcerebuild;
/* Collect the values of rebuild from all nodes. */
#ifdef WITH_MPI
int buff = 0;
if (MPI_Allreduce(&rebuild, &buff, 1, MPI_INT, MPI_MAX, MPI_COMM_WORLD) !=
MPI_SUCCESS)
error("Failed to aggregate the rebuild flag across nodes.");
rebuild = buff;
#endif
/* And rebuild if necessary. */
if (rebuild) {
/* Drift all particles to the current time if needed. */
if (drift_all) engine_drift_all(e);
#ifdef SWIFT_DEBUG_CHECKS
/* Check that all cells have been drifted to the current time, unless
* we have just repartitioned, that can include cells that have not
if (e->forcerepart || e->forcerebuild) {
/* Check that all cells have been drifted to the current time.
* That can include cells that have not
* previously been active on this rank. */
if (!postrepart) space_check_drift_point(e->s, e->ti_current);
space_check_drift_point(e->s, e->ti_old);
}
#endif
engine_rebuild(e);
}
if (postrepart) engine_unskip(e);
/* Do we need repartitioning ? */
if (e->forcerepart) engine_repartition(e);
/* Do we need rebuilding ? */
if (e->forcerebuild) engine_rebuild(e);
/* Unskip active tasks and check for rebuild */
engine_unskip(e);
/* Re-rank the tasks every now and then. */
if (e->tasks_age % engine_tasksreweight == 1) {
......@@ -2656,7 +2642,7 @@ void engine_prepare(struct engine *e, int drift_all, int postrepart) {
TIMER_TOC(timer_prepare);
if (e->verbose)
message("took %.3f %s (including drift all, rebuild and reweight).",
message("took %.3f %s (including unskip and reweight).",
clocks_from_ticks(getticks() - tic), clocks_getunit());
}
......@@ -2715,7 +2701,7 @@ void engine_collect_kick(struct cell *c) {
/* Counters for the different quantities. */
int updated = 0, g_updated = 0, s_updated = 0;
integertime_t ti_end_min = max_nr_timesteps;
integertime_t ti_end_min = max_nr_timesteps, ti_end_max = 0, ti_beg_max = 0;
/* Collect the values from the progeny. */
for (int k = 0; k < 8; k++) {
......@@ -2727,6 +2713,8 @@ void engine_collect_kick(struct cell *c) {
/* And update */
ti_end_min = min(ti_end_min, cp->ti_end_min);
ti_end_max = max(ti_end_max, cp->ti_end_max);
ti_beg_max = max(ti_beg_max, cp->ti_beg_max);
updated += cp->updated;
g_updated += cp->g_updated;
s_updated += cp->s_updated;
......@@ -2740,6 +2728,8 @@ void engine_collect_kick(struct cell *c) {
/* Store the collected values in the cell. */
c->ti_end_min = ti_end_min;
c->ti_end_max = ti_end_max;
c->ti_beg_max = ti_beg_max;
c->updated = updated;
c->g_updated = g_updated;
c->s_updated = s_updated;
......@@ -2755,7 +2745,7 @@ void engine_collect_timestep(struct engine *e) {
const ticks tic = getticks();
int updates = 0, g_updates = 0, s_updates = 0;
integertime_t ti_end_min = max_nr_timesteps;
integertime_t ti_end_min = max_nr_timesteps, ti_end_max = 0, ti_beg_max = 0;
const struct space *s = e->s;
/* Collect the cell data. */
......@@ -2768,6 +2758,8 @@ void engine_collect_timestep(struct engine *e) {
/* And aggregate */
ti_end_min = min(ti_end_min, c->ti_end_min);
ti_end_max = max(ti_end_max, c->ti_end_max);
ti_beg_max = max(ti_beg_max, c->ti_beg_max);
updates += c->updated;
g_updates += c->g_updated;
s_updates += c->s_updated;
......@@ -2805,6 +2797,8 @@ void engine_collect_timestep(struct engine *e) {
#endif
e->ti_end_min = ti_end_min;
e->ti_end_max = ti_end_max;
e->ti_beg_max = ti_beg_max;
e->updates = updates;
e->g_updates = g_updates;
e->s_updates = s_updates;
......@@ -2878,7 +2872,7 @@ void engine_skip_force_and_kick(struct engine *e) {
*
* @param e The #engine to act on.
*/
void engine_skip_drift_and_kick(struct engine *e) {
void engine_skip_drift(struct engine *e) {
struct task *tasks = e->sched.tasks;
const int nr_tasks = e->sched.nr_tasks;
......@@ -2888,7 +2882,7 @@ void engine_skip_drift_and_kick(struct engine *e) {
struct task *t = &tasks[i];
/* Skip everything that updates the particles */
if (t->type == task_type_drift || t->type == task_type_kick1) t->skip = 1;
if (t->type == task_type_drift) t->skip = 1;
}
}
......@@ -2949,13 +2943,15 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs) {
if (e->nodeID == 0) message("Computing initial gas densities.");
engine_prepare(e, 0, 0);
engine_marktasks(e);
/* Construct all cells and tasks to start everything */
engine_rebuild(e);
/* No time integration. We just want the density and ghosts */
engine_skip_force_and_kick(e);
/* Print the number of active tasks ? */
if (e->verbose) engine_print_task_counts(e);
/* Now, launch the calculation */
TIMER_TIC;
engine_launch(e, e->nr_threads);
......@@ -2967,7 +2963,6 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs) {
if (e->nodeID == 0) message("Converting internal energy variable.");
/* Apply the conversion */
// space_map_cells_pre(s, 0, cell_convert_hydro, NULL);
for (size_t i = 0; i < s->nr_parts; ++i)
hydro_convert_quantities(&s->parts[i], &s->xparts[i]);
......@@ -2982,12 +2977,21 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs) {
/* Now time to get ready for the first time-step */
if (e->nodeID == 0) message("Running initial fake time-step.");
/* Prepare all the tasks again for a new round */
engine_marktasks(e);
engine_skip_drift_and_kick(e);
/* No drift this time */
engine_skip_drift(e);
/* Print the number of active tasks ? */
if (e->verbose) engine_print_task_counts(e);
/* Run the 0th time-step */
engine_launch(e, e->nr_threads);
/* Recover the (integer) end of the next time-step */
engine_collect_timestep(e);
clocks_gettime(&time2);
#ifdef SWIFT_DEBUG_CHECKS
......@@ -2997,7 +3001,7 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs) {
#endif
/* Ready to go */
e->step = 0;
e->step = -1;
e->forcerebuild = 1;
e->wallclock_time = (float)clocks_diff(&time1, &time2);
......@@ -3011,8 +3015,6 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs) {
*/
void engine_step(struct engine *e) {
double snapshot_drift_time = 0.;
TIMER_TIC2;
struct clocks_time time1, time2;
......@@ -3020,36 +3022,13 @@ void engine_step(struct engine *e) {
e->tic_step = getticks();
/* Recover the (integer) end of the next time-step */
engine_collect_timestep(e);
/* Check for output */
while (e->ti_end_min >= e->ti_nextSnapshot && e->ti_nextSnapshot > 0) {
e->ti_old = e->ti_current;
e->ti_current = e->ti_nextSnapshot;
e->time = e->ti_current * e->timeBase + e->timeBegin;
e->timeOld = e->ti_old * e->timeBase + e->timeBegin;
e->timeStep = (e->ti_current - e->ti_old) * e->timeBase;
snapshot_drift_time = e->timeStep;
/* Drift everybody to the snapshot position */
engine_drift_all(e);
/* Dump... */
engine_dump_snapshot(e);
/* ... and find the next output time */
engine_compute_next_snapshot_time(e);
}