Commit 8f560a8a authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Re-arranged the time-step sequence. Still some debugging checks in place.

parent 26129ed0
......@@ -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")
......
......@@ -6,6 +6,9 @@ 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).
......@@ -42,5 +45,5 @@ PointMassPotential:
position_y: 50.
position_z: 50.
mass: 1e10 # mass of external point mass in internal units
timestep_mult: 0.03 # controls time step
timestep_mult: 0.001 # controls time step
......@@ -96,6 +96,10 @@ grp1 = file.create_group("/PartType1")
#generate particle positions
radius = max_radius * (numpy.random.rand(numPart))**(1./3.)
#radius = numpy.zeros(3)
#radius[0] = 8.
#radius[1] = 10.
#radius[2] = 20.
print '---------------------'
print 'Radius: minimum = ',min(radius)
print 'Radius: maximum = ',max(radius)
......
......@@ -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.
......@@ -141,4 +143,78 @@ __attribute__((always_inline)) INLINE static int gpart_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);
}
#endif /* SWIFT_ACTIVE_H */
......@@ -60,6 +60,8 @@
#include "space.h"
#include "timers.h"
int counter = 0;
/* Global variables. */
int cell_next_tag = 0;
......@@ -1074,6 +1076,9 @@ void cell_drift(struct cell *c, const struct engine *e) {
const double dt = (ti_current - ti_old) * timeBase;
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); */
/* Check that we are actually going to move forward. */
if (ti_current < ti_old) error("Attempt to drift to the past");
......@@ -1101,6 +1106,16 @@ void cell_drift(struct cell *c, const struct engine *e) {
/* Drift... */
drift_gpart(gp, dt, timeBase, ti_old, ti_current);
#ifdef ICHECK
if (gp->id_or_neg_offset == ICHECK) {
message("--- ti_current=%lld time=%e dt=%e---", e->ti_current, e->time,
dt);
counter++;
message("Drift counter: %d", counter);
printgParticle_single(gp);
}
#endif
/* 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] +
......
......@@ -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;
......@@ -224,6 +224,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;
......
......@@ -20,6 +20,8 @@
#ifndef SWIFT_CONST_H
#define SWIFT_CONST_H
#include <stdio.h>
/* SPH Viscosity constants. */
#define const_viscosity_alpha 0.8f
#define const_viscosity_alpha_min \
......@@ -59,4 +61,9 @@
#define SOURCETERMS_NONE
//#define SOURCETERMS_SN_FEEDBACK
#define ICHECK 0
#define num_files 4
extern FILE* files_timestep[num_files];
#endif /* SWIFT_CONST_H */
......@@ -39,8 +39,23 @@
* @param ti_current Integer end of time-step
*/
__attribute__((always_inline)) INLINE static void drift_gpart(
struct gpart *restrict gp, float dt, double timeBase, integertime_t ti_old,
struct gpart *restrict gp, double 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
//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;
gp->x[1] += gp->v_full[1] * dt;
......@@ -63,7 +78,7 @@ __attribute__((always_inline)) INLINE static void drift_gpart(
* @param ti_current Integer end of time-step
*/
__attribute__((always_inline)) INLINE static void drift_part(
struct part *restrict p, struct xpart *restrict xp, float dt,
struct part *restrict p, struct xpart *restrict xp, double dt,
double timeBase, integertime_t ti_old, integertime_t ti_current) {
#ifdef SWIFT_DEBUG_CHECKS
......
......@@ -71,6 +71,8 @@
#include "units.h"
#include "version.h"
FILE* files_timestep[num_files];
/* Particle cache size. */
#define CACHE_SIZE 512
......@@ -156,12 +158,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. */
......@@ -2160,8 +2162,7 @@ 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);
}
......@@ -2255,6 +2256,8 @@ void engine_rebuild(struct engine *e) {
const ticks tic = getticks();
message("REBUILD !!!!");
/* Clear the forcerebuild flag, whatever it was. */
e->forcerebuild = 0;
......@@ -2277,7 +2280,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),
......@@ -2295,41 +2298,12 @@ void engine_rebuild(struct engine *e) {
* 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
* 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);
engine_unskip(e);
/* Re-rank the tasks every now and then. */
if (e->tasks_age % engine_tasksreweight == 1) {
......@@ -2340,10 +2314,64 @@ 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());
}
/* 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.
*
......@@ -2399,7 +2427,7 @@ void engine_collect_kick(struct cell *c) {
/* Counters for the different quantities. */
int updated = 0, g_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++) {
......@@ -2411,6 +2439,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;
......@@ -2422,6 +2452,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;
}
......@@ -2554,7 +2586,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;
......@@ -2564,7 +2596,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;
}
}
......@@ -2623,15 +2655,20 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs) {
struct clocks_time time1, time2;
clocks_gettime(&time1);
if (e->nodeID == 0) message("Computing initial gas densities.");
if (e->nodeID == 0) message("Computing initial gas densities.\n\n\n");
engine_prepare(e, 0, 0);
engine_rebuild(e);
// engine_prepare(e);
engine_marktasks(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);
......@@ -2660,10 +2697,16 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs) {
engine_marktasks(e);
engine_skip_drift_and_kick(e);
engine_skip_drift(e);
/* 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 */
engine_collect_timestep(e);
clocks_gettime(&time2);
#ifdef SWIFT_DEBUG_CHECKS
......@@ -2672,7 +2715,7 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs) {
/* Ready to go */
e->step = 0;
e->forcerebuild = 1;
e->forcerebuild = 0;
e->wallclock_time = (float)clocks_diff(&time1, &time2);
if (e->verbose) message("took %.3f %s.", e->wallclock_time, clocks_getunit());
......@@ -2685,8 +2728,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;
......@@ -2694,36 +2735,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);
}
/* 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 + snapshot_drift_time;
e->timeStep = (e->ti_current - e->ti_old) * e->timeBase;
if (e->nodeID == 0) {
......@@ -2735,33 +2753,112 @@ void engine_step(struct engine *e) {
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);
/* 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); */
/* fflush(files_timestep[i]); */
/* } */
}
/* Drift only the necessary particles, that means all particles
* if we are about to repartition. */
const int repart = (e->forcerepart != REPART_NONE);
const int drift_all = (e->policy & engine_policy_drift_all);
if (repart || drift_all) engine_drift_all(e);
/* Do we need repartitioning ? */
if(e->forcerepart) engine_repartition(e);
/* Do we need rebuilding ? */
if(e->forcerebuild) engine_rebuild(e);
/* Re-distribute the particles amongst the nodes? */
if (repart) engine_repartition(e);
/* Prepare the tasks to be launched. */
engine_prepare(e);
/* Prepare the space. */
engine_prepare(e, !(drift_all || repart), repart);
if (e->step % 99 == 0) e->forcerebuild = 1;
/* Print the number of active tasks ? */
if (e->verbose) engine_print_task_counts(e);
/* Save some statistics */
/* Save some statistics ? */
if (e->time - e->timeLastStatistics >= e->deltaTimeStatistics) {
engine_print_stats(e);
e->timeLastStatistics += e->deltaTimeStatistics;
}
/* Send off the runners. */
/* Start all the tasks. */
TIMER_TIC;
engine_launch(e, e->nr_threads);
TIMER_TOC(timer_runners);
if (e->ti_end_min >= e->ti_nextSnapshot && e->ti_nextSnapshot > 0)
e->dump_snapshot = 1;
/* 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);
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
}
/* Write a snapshot ? */
if(e->dump_snapshot) {
/* Dump... */
engine_dump_snapshot(e);
/* ... and find the next output time */