Commit 4a6150af authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Time-step re-ordered

parent fa9029e0
......@@ -6,9 +6,6 @@ InternalUnitSystem:
UnitCurrent_in_cgs: 1 # Amperes
UnitTemp_in_cgs: 1 # Kelvin
#Scheduler:
# max_top_level_cells: 20
# Parameters governing the time integration
TimeIntegration:
time_begin: 0. # The starting time of the simulation (in internal units).
......
......@@ -6,6 +6,9 @@ InternalUnitSystem:
UnitCurrent_in_cgs: 1 # Amperes
UnitTemp_in_cgs: 1 # Kelvin
#Scheduler:
# max_top_level_cells: 3
# Parameters governing the time integration
TimeIntegration:
time_begin: 0. # The starting time of the simulation (in internal units).
......
......@@ -143,7 +143,6 @@ __attribute__((always_inline)) INLINE static int gpart_is_active(
return (ti_end == ti_current);
}
/* Are cells / particles active for kick1 tasks ? */
/**
......@@ -179,7 +178,8 @@ __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);
const integertime_t ti_beg =
get_integer_time_begin(ti_current + 1, p->time_bin);
#ifdef SWIFT_DEBUG_CHECKS
if (ti_beg > ti_current)
......@@ -203,7 +203,8 @@ __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);
const integertime_t ti_beg =
get_integer_time_begin(ti_current + 1, gp->time_bin);
#ifdef SWIFT_DEBUG_CHECKS
if (ti_beg > ti_current)
......@@ -216,5 +217,4 @@ __attribute__((always_inline)) INLINE static int gpart_is_starting(
return (ti_beg == ti_current);
}
#endif /* SWIFT_ACTIVE_H */
......@@ -912,6 +912,11 @@ int cell_unskip_tasks(struct cell *c, struct scheduler *s) {
struct engine *e = s->space->e;
#endif
//message("unskip! c=%p c->kick1=%p c->drift=%p", c, c->kick1, c->drift);
//task_print(c->drift);
int ret = 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;
......@@ -937,7 +942,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;
ret = 1;
#ifdef WITH_MPI
/* Activate the send/recv flags. */
......@@ -1034,7 +1039,9 @@ 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;
//task_print(c->drift);
return ret;
}
/**
......@@ -1077,7 +1084,8 @@ void cell_drift(struct cell *c, const struct engine *e) {
float dx_max = 0.f, dx2_max = 0.f, h_max = 0.f;
/* if (c->gcount > 0) */
/* message("dt=%e, ti_old=%lld ti_current=%lld", dt, c->ti_old, e->ti_current); */
/* message("dt=%e, ti_old=%lld ti_current=%lld", dt, c->ti_old,
* e->ti_current); */
/* Check that we are actually going to move forward. */
if (ti_current < ti_old) error("Attempt to drift to the past");
......
......@@ -226,7 +226,7 @@ struct cell {
/*! 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;
......
......@@ -53,8 +53,8 @@ __attribute__((always_inline)) INLINE static void drift_gpart(
gp->ti_drift = ti_current;
#endif
//message("dt= %e", dt);
//fprintf(files_timestep[gp->id_or_neg_offset], "drift: dt=%e\n", dt);
// message("dt= %e", dt);
// fprintf(files_timestep[gp->id_or_neg_offset], "drift: dt=%e\n", dt);
/* Drift... */
gp->x[0] += gp->v_full[0] * dt;
......
......@@ -71,7 +71,7 @@
#include "units.h"
#include "version.h"
FILE* files_timestep[num_files];
FILE *files_timestep[num_files];
/* Particle cache size. */
#define CACHE_SIZE 512
......@@ -2162,10 +2162,11 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
}
/* Kick/Drift/Init? */
else if (t->type == task_type_kick1 || t->type == task_type_kick2 || t->type == task_type_drift || t->type == task_type_init) {
else if (t->type == task_type_kick1 || t->type == task_type_kick2 ||
t->type == task_type_drift || t->type == task_type_init) {
if (cell_is_active(t->ci, e)) scheduler_activate(s, t);
}
/* Time-step? */
else if (t->type == task_type_timestep) {
t->ci->updated = 0;
......@@ -2192,6 +2193,8 @@ int engine_marktasks(struct engine *e) {
const ticks tic = getticks();
int rebuild_space = 0;
message("marktasks");
/* Run through the tasks and mark as skip or not. */
size_t extra_data[3] = {(size_t)e, rebuild_space, (size_t)&e->sched};
threadpool_map(&e->threadpool, engine_marktasks_mapper, s->tasks, s->nr_tasks,
......@@ -2280,7 +2283,7 @@ void engine_rebuild(struct engine *e) {
error("engine_marktasks failed after space_rebuild.");
/* Print the status of the system */
//if (e->verbose) engine_print_task_counts(e);
// if (e->verbose) engine_print_task_counts(e);
if (e->verbose)
message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
......@@ -2318,60 +2321,6 @@ void engine_prepare(struct engine *e) {
clocks_from_ticks(getticks() - tic), clocks_getunit());
}
/* void engine_prepare(struct engine *e, int drift_all, int postrepart) { */
/* TIMER_TIC; */
/* // message("e->forcerebuild: %d", e->forcerebuild); */
/* /\* 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; */
/* // message("e->forcerebuild: %d", 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 */
/* * previously been active on this rank. *\/ */
/* if (!postrepart) space_check_drift_point(e->s, e->ti_current); */
/* #endif */
/* engine_rebuild(e); */
/* } */
/* if (postrepart) engine_unskip(e); */
/* /\* Re-rank the tasks every now and then. *\/ */
/* if (e->tasks_age % engine_tasksreweight == 1) { */
/* scheduler_reweight(&e->sched, e->verbose); */
/* } */
/* e->tasks_age += 1; */
/* TIMER_TOC(timer_prepare); */
/* if (e->verbose) */
/* message("took %.3f %s (including drift all, rebuild and reweight).", */
/* clocks_from_ticks(getticks() - tic), clocks_getunit()); */
/* } */
/**
* @brief Implements a barrier for the #runner threads.
*
......@@ -2468,7 +2417,7 @@ void engine_collect_timestep(struct engine *e) {
const ticks tic = getticks();
int updates = 0, g_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. */
......@@ -2481,6 +2430,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;
......@@ -2514,6 +2465,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;
......@@ -2668,7 +2621,7 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs) {
/* 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);
......@@ -2701,7 +2654,7 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs) {
/* Print the number of active tasks ? */
if (e->verbose) engine_print_task_counts(e);
engine_launch(e, e->nr_threads);
/* Recover the (integer) end of the next time-step */
......@@ -2715,7 +2668,7 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs) {
/* Ready to go */
e->step = 0;
e->forcerebuild = 0;
e->forcerebuild = 1;
e->wallclock_time = (float)clocks_diff(&time1, &time2);
if (e->verbose) message("took %.3f %s.", e->wallclock_time, clocks_getunit());
......@@ -2746,8 +2699,8 @@ void engine_step(struct engine *e) {
if (e->nodeID == 0) {
/* Print some information to the screen */
printf(" %6d %14e %14e %10zu %10zu %21.3f\n", e->step, e->time,
e->timeStep, e->updates, e->g_updates, e->wallclock_time);
printf(" %6d %lld %14e %14e %10zu %10zu %21.3f\n", e->step, e->ti_current,
e->time, e->timeStep, e->updates, e->g_updates, e->wallclock_time);
fflush(stdout);
fprintf(e->file_timesteps, " %6d %14e %14e %10zu %10zu %21.3f\n", e->step,
......@@ -2755,24 +2708,26 @@ void engine_step(struct engine *e) {
fflush(e->file_timesteps);
/* for(int i=0; i<num_files; ++i) { */
/* fprintf(files_timestep[i], " %6d %14e %14e %10zu %10zu %21.3f\n", e->step, */
/* e->time, e->timeStep, e->updates, e->g_updates, e->wallclock_time); */
/* fprintf(files_timestep[i], " %6d %14e %14e %10zu %10zu %21.3f\n",
* e->step, */
/* e->time, e->timeStep, e->updates, e->g_updates,
* e->wallclock_time); */
/* fflush(files_timestep[i]); */
/* } */
}
/* Do we need repartitioning ? */
if(e->forcerepart) engine_repartition(e);
if (e->forcerepart != REPART_NONE) engine_repartition(e);
/* Do we need rebuilding ? */
if(e->forcerebuild) engine_rebuild(e);
if (e->forcerebuild) engine_rebuild(e);
/* Prepare the tasks to be launched. */
engine_prepare(e);
if (e->step % 99 == 0) e->forcerebuild = 1;
if (e->step == 1) e->forcerebuild = 1;
/* Print the number of active tasks ? */
if (e->verbose) engine_print_task_counts(e);
......@@ -2792,8 +2747,7 @@ void engine_step(struct engine *e) {
/* Drift everybody (i.e. what has not yet been drifted) */
/* to the current time */
if(e->dump_snapshot || e->forcerebuild || e->forcerepart != REPART_NONE) {
message("snap=%d, rebuild=%d repart=%d", e->dump_snapshot, e->forcerebuild, e->forcerepart);
if (e->dump_snapshot || e->forcerebuild || e->forcerepart != REPART_NONE) {
engine_drift_all(e);
......@@ -2803,11 +2757,14 @@ void engine_step(struct engine *e) {
* previously been active on this rank. */
space_check_drift_point(e->s, e->ti_current);
#endif
}
message("snap=%d, rebuild=%d repart=%d", e->dump_snapshot, e->forcerebuild,
e->forcerepart);
/* Write a snapshot ? */
if(e->dump_snapshot) {
if (e->dump_snapshot) {
/* Dump... */
engine_dump_snapshot(e);
......@@ -2815,49 +2772,11 @@ void engine_step(struct engine *e) {
/* ... and find the next output time */
engine_compute_next_snapshot_time(e);
}
/* Recover the (integer) end of the next time-step */
engine_collect_timestep(e);
/* /\* Check for snapshot dump *\/ */
/* while { */
/* /\* Drift everybody (i.e. what has not yet been drifted) *\/ */
/* /\* to the current time *\/ */
/* engine_drift_all(e); */
/* /\* Dump... *\/ */
/* engine_dump_snapshot(e); */
/* /\* ... and find the next output time *\/ */
/* engine_compute_next_snapshot_time(e); */
/* } */
/* /\* Repartition if necessary *\/ */
/* if (e->forcerepart != REPART_NONE) { */
/* /\* Drift all particles to the current time if needed. *\/ */
/* engine_drift_all(e); */
/* engine_repartition(e); */
/* } */
/* /\* Rebuild if necessary *\/ */
/* if (e->forcerebuild) { */
/* /\* Drift all particles to the current time if needed. *\/ */
/* engine_drift_all(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_current); */
/* #endif */
/* engine_rebuild(e); */
/* } */
TIMER_TOC2(timer_step);
......@@ -2867,91 +2786,6 @@ void engine_step(struct engine *e) {
e->toc_step = getticks();
}
/* void engine_step(struct engine *e) { */
/* TIMER_TIC2; */
/* struct clocks_time time1, time2; */
/* clocks_gettime(&time1); */
/* e->tic_step = getticks(); */
/* /\* Recover the (integer) end of the next time-step *\/ */
/* engine_collect_timestep(e); */
/* /\* Move forward in time *\/ */
/* e->ti_old = e->ti_current; */
/* e->ti_current = e->ti_end_min; */
/* e->step += 1; */
/* 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; */
/* // if(e->step % 100 == 0) e->forcerebuild = 1; */
/* if (e->nodeID == 0) { */
/* /\* Print some information to the screen *\/ */
/* printf(" %6d %14e %14e %10zu %10zu %21.3f\n", e->step, e->time, */
/* e->timeStep, e->updates, e->g_updates, e->wallclock_time); */
/* fflush(stdout); */
/* fprintf(e->file_timesteps, " %6d %14e %14e %10zu %10zu %21.3f\n",
* e->step, */
/* e->time, e->timeStep, e->updates, e->g_updates,
* e->wallclock_time); */
/* fflush(e->file_timesteps); */
/* } */
/* /\* Drift only the necessary particles, that means all particles */
/* * if we are about to repartition. *\/ */
/* const int repart = (e->forcerepart != REPART_NONE); */
/* // message("repart: %d", repart); */
/* const int drift_all = (e->policy & engine_policy_drift_all); */
/* // message("drift_all: %d", drift_all); */
/* if (repart || drift_all) engine_drift_all(e); */
/* /\* Re-distribute the particles amongst the nodes? *\/ */
/* if (repart) engine_repartition(e); */
/* /\* Prepare the space. *\/ */
/* engine_prepare(e, !(drift_all || repart), repart); */
/* if (e->verbose) engine_print_task_counts(e); */
/* /\* Save some statistics *\/ */
/* if (e->time - e->timeLastStatistics >= e->deltaTimeStatistics) { */
/* engine_print_stats(e); */
/* e->timeLastStatistics += e->deltaTimeStatistics; */
/* } */
/* /\* Send off the runners. *\/ */
/* TIMER_TIC; */
/* engine_launch(e, e->nr_threads); */
/* TIMER_TOC(timer_runners); */
/* /\* Check for output *\/ */
/* while (e->ti_end_min >= e->ti_nextSnapshot && e->ti_nextSnapshot > 0) { */
/* /\* Drift everybody (i.e. what has not yet been drifted) 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); */
/* } */
/* TIMER_TOC2(timer_step); */
/* clocks_gettime(&time2); */
/* e->wallclock_time = (float)clocks_diff(&time1, &time2); */
/* e->toc_step = getticks(); */
/* } */
/**
* @brief Returns 1 if the simulation has reached its end point, 0 otherwise
*/
......@@ -2966,6 +2800,8 @@ int engine_is_done(struct engine *e) {
*/
void engine_unskip(struct engine *e) {
message("unskip");
const ticks tic = getticks();
threadpool_map(&e->threadpool, runner_do_unskip_mapper, e->s->cells_top,
e->s->nr_cells, sizeof(struct cell), 1, e);
......@@ -2988,7 +2824,6 @@ void engine_drift_all(struct engine *e) {
/* fprintf(files_timestep[i], "drift all\n"); */
/* } */
const ticks tic = getticks();
threadpool_map(&e->threadpool, runner_do_drift_mapper, e->s->cells_top,
e->s->nr_cells, sizeof(struct cell), 1, e);
......@@ -3210,7 +3045,6 @@ void engine_dump_snapshot(struct engine *e) {
/* fprintf(files_timestep[i], "dump\n"); */
/* } */
message("dump");
/* Dump... */
......@@ -3230,7 +3064,7 @@ void engine_dump_snapshot(struct engine *e) {
#endif
e->dump_snapshot = 0;
clocks_gettime(&time2);
if (e->verbose)
message("writing particle properties took %.3f %s.",
......@@ -3714,12 +3548,12 @@ void engine_init(struct engine *e, struct space *s,
free(buf);
#endif
for(int i=0; i<num_files; ++i) {
for (int i = 0; i < num_files; ++i) {
char name[10];
sprintf(name, "dt_%d.txt", i);
files_timestep[i] = fopen(name, "w");
}
/* Wait for the runner threads to be in place. */
while (e->barrier_running || e->barrier_launch)
if (pthread_cond_wait(&e->barrier_cond, &e->barrier_mutex) != 0)
......@@ -3787,9 +3621,8 @@ void engine_compute_next_snapshot_time(struct engine *e) {
*/
void engine_clean(struct engine *e) {
for(int i=0; i<num_files; ++i)
fclose(files_timestep[i]);
for (int i = 0; i < num_files; ++i) fclose(files_timestep[i]);
free(e->snapshotUnits);
free(e->links);
scheduler_clean(&e->sched);
......
......@@ -130,6 +130,12 @@ struct engine {
/* Minimal ti_end for the next time-step */
integertime_t ti_end_min;
/* Maximal ti_end for the next time-step */
integertime_t ti_end_max;
/* Maximal ti_beg for the next time-step */
integertime_t ti_beg_max;
/* Number of particles updated */
size_t updates, g_updates;
......@@ -182,7 +188,7 @@ struct engine {
/* Need to dump a snapshot ? */
int dump_snapshot;
/* How many steps have we done with the same set of tasks? */
int tasks_age;
......
......@@ -48,13 +48,12 @@ __attribute__((always_inline)) INLINE static void kick_gpart(
"Particle has not been kicked to the current time gp->ti_kick=%lld, "
"ti_start=%lld, ti_end=%lld",
gp->ti_kick, ti_start, ti_end);
gp->ti_kick = ti_end;
#endif
//message("dt= %e gp->ti_kick=%lld", dt, gp->ti_kick);
//fprintf(files_timestep[gp->id_or_neg_offset], "kick: dt=%e\n", dt);
// message("dt= %e gp->ti_kick=%lld", dt, gp->ti_kick);
// fprintf(files_timestep[gp->id_or_neg_offset], "kick: dt=%e\n", dt);
/* Kick particles in momentum space */
gp->v_full[0] += gp->a_grav[0] * dt;
......
......@@ -80,14 +80,17 @@ __attribute__((always_inline)) INLINE static float external_gravity_timestep(
/* (-g->v_full[1] + 3.f * rinv2 * drdv * dy); */
/* const float dota_z = G_newton * potential->mass * rinv3 * */
/* (-g->v_full[2] + 3.f * rinv2 * drdv * dz); */
/* const float dota_2 = dota_x * dota_x + dota_y * dota_y + dota_z * dota_z; */
/* const float a_2 = g->a_grav[0] * g->a_grav[0] + g->a_grav[1] * g->a_grav[1] */
/* const float dota_2 = dota_x * dota_x + dota_y * dota_y + dota_z * dota_z;
*/
/* const float a_2 = g->a_grav[0] * g->a_grav[0] + g->a_grav[1] * g->a_grav[1]
*/
/* + g->a_grav[2] * g->a_grav[2]; */
/* return potential->timestep_mult * sqrtf(a_2 / dota_2); */
const double r = sqrt(dx * dx + dy * dy + dz * dz);
return potential->timestep_mult * sqrt(r * r * r / (G_newton * potential->mass));
return potential->timestep_mult *
sqrt(r * r * r / (G_newton * potential->mass));
}
/**
......
......@@ -783,6 +783,8 @@ static void runner_do_unskip(struct cell *c, struct engine *e) {
if (forcerebuild) atomic_inc(&e->forcerebuild);
}
// message("c->depth=%d c->split=%d c->count=%.5d c->ti_end_min=%lld c->ti_end_max=%lld c->ti_beg_max=%lld ti_current=%lld", c->depth, c->split, c->count, c->ti_end_min, c->ti_end_max, c->ti_beg_max ,e->ti_current);
/* Recurse */
if (c->split) {
for (int k = 0; k < 8; k++) {
......@@ -869,7 +871,7 @@ void runner_do_kick1(struct runner *r, struct cell *c, int timer) {
/* Anything to do here? */
if (!cell_is_starting(c, e)) return;
/* Recurse? */
if (c->split) {
for (int k = 0; k < 8; k++)
......@@ -890,16 +892,16 @@ void runner_do_kick1(struct runner *r, struct cell *c, int timer) {
const integertime_t ti_begin =
get_integer_time_begin(ti_current + 1, p->time_bin);
/* #ifdef SWIFT_DEBUG_CHECKS */
/* const integertime_t ti_end = */
/* get_integer_time_end(ti_current, p->time_bin); */
#ifdef SWIFT_DEBUG_CHECKS
const integertime_t ti_end =
get_integer_time_end(ti_current + 1, p->time_bin);
/* if (ti_end - ti_begin != ti_step) */
/* error( */
/* "Particle in wrong time-bin, ti_end=%lld, ti_begin=%lld, " */
/* "ti_step=%lld time_bin=%d ti_current=%lld", */
/* ti_end, ti_begin, ti_step, p->time_bin, ti_current); */
/* #endif */
if (ti_begin != ti_current)
error(
"Particle in wrong time-bin, ti_end=%lld, ti_begin=%lld, "
"ti_step=%lld time_bin=%d ti_current=%lld",
ti_end, ti_begin, ti_step, p->time_bin, ti_current);
#endif