Commit 54523c79 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Energy statistics is now done after drift

parent f025e4c9
......@@ -142,7 +142,7 @@ struct cell {
int owner;
/* Momentum of particles in cell. */
float mom[3], ang[3];
double mom[3], ang_mom[3];
/* Mass, potential, internal and kinetic energy of particles in this cell. */
double mass, e_pot, e_int, e_kin;
......
......@@ -1796,8 +1796,9 @@ void engine_barrier(struct engine *e, int tid) {
/**
* @brief Mapping function to collect the data from the kick.
*
* @param c A super-cell.
*/
void engine_collect_kick(struct cell *c) {
/* Skip super-cells (Their values are already set) */
......@@ -1805,15 +1806,13 @@ void engine_collect_kick(struct cell *c) {
/* Counters for the different quantities. */
int updated = 0, g_updated = 0;
double e_kin = 0.0, e_int = 0.0, e_pot = 0.0;
float mom[3] = {0.0f, 0.0f, 0.0f}, ang[3] = {0.0f, 0.0f, 0.0f};
int ti_end_min = max_nr_timesteps, ti_end_max = 0;
int ti_end_min = max_nr_timesteps;
/* Only do something is the cell is non-empty */
if (c->count != 0 || c->gcount != 0) {
/* If this cell is not split, I'm in trouble. */
if (!c->split) error("Cell has no super-cell.");
if (!c->split) error("Cell is not split.");
/* Collect the values from the progeny. */
for (int k = 0; k < 8; k++) {
......@@ -1825,36 +1824,195 @@ 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);
updated += cp->updated;
g_updated += cp->g_updated;
}
}
}
/* Store the collected values in the cell. */
c->ti_end_min = ti_end_min;
c->updated = updated;
c->g_updated = g_updated;
}
/**
* @brief Collects the next time-step by making each super-cell recurse
* to collect the minimal of ti_end and the number of updated particles.
*
* @param The #engine.
*/
void engine_collect_timestep(struct engine *e) {
int updates = 0, g_updates = 0;
int ti_end_min = max_nr_timesteps;
const struct space *s = e->s;
/* Collect the cell data. */
for (int k = 0; k < s->nr_cells; k++)
if (s->cells[k].nodeID == e->nodeID) {
struct cell *c = &s->cells[k];
/* Make the top-cells recurse */
engine_collect_kick(c);
/* And aggregate */
ti_end_min = min(ti_end_min, c->ti_end_min);
updates += c->updated;
g_updates += c->g_updated;
}
/* Aggregate the data from the different nodes. */
#ifdef WITH_MPI
{
int in_i[1], out_i[1];
in_i[0] = 0;
out_i[0] = ti_end_min;
if (MPI_Allreduce(out_i, in_i, 1, MPI_INT, MPI_MIN, MPI_COMM_WORLD) !=
MPI_SUCCESS)
error("Failed to aggregate t_end_min.");
ti_end_min = in_i[0];
}
{
unsigned long long in_ll[2], out_ll[2];
out_ll[0] = updates;
out_ll[1] = g_updates;
if (MPI_Allreduce(out_ll, in_ll, 2, MPI_LONG_LONG_INT, MPI_SUM,
MPI_COMM_WORLD) != MPI_SUCCESS)
error("Failed to aggregate energies.");
updates = in_ll[0];
g_updates = in_ll[1];
}
#endif
e->ti_end_min = ti_end_min;
e->updates = updates;
e->g_updates = g_updates;
}
/**
* @brief Mapping function to collect the data from the drift.
*
* @param c A super-cell.
*/
void engine_collect_drift(struct cell *c) {
/* Skip super-cells (Their values are already set) */
if (c->drift != NULL) return;
/* Counters for the different quantities. */
double e_kin = 0.0, e_int = 0.0, e_pot = 0.0, mass = 0.0;
double mom[3] = {0.0, 0.0, 0.0}, ang_mom[3] = {0.0, 0.0, 0.0};
/* Only do something is the cell is non-empty */
if (c->count != 0 || c->gcount != 0) {
/* If this cell is not split, I'm in trouble. */
if (!c->split) error("Cell has no super-cell.");
/* Collect the values from the progeny. */
for (int k = 0; k < 8; k++) {
struct cell *cp = c->progeny[k];
if (cp != NULL) {
/* Recurse */
engine_collect_drift(cp);
/* And update */
mass += cp->mass;
e_kin += cp->e_kin;
e_int += cp->e_int;
e_pot += cp->e_pot;
mom[0] += cp->mom[0];
mom[1] += cp->mom[1];
mom[2] += cp->mom[2];
ang[0] += cp->ang[0];
ang[1] += cp->ang[1];
ang[2] += cp->ang[2];
ang_mom[0] += cp->ang_mom[0];
ang_mom[1] += cp->ang_mom[1];
ang_mom[2] += cp->ang_mom[2];
}
}
}
/* Store the collected values in the cell. */
c->ti_end_min = ti_end_min;
c->ti_end_max = ti_end_max;
c->updated = updated;
c->g_updated = g_updated;
c->mass = mass;
c->e_kin = e_kin;
c->e_int = e_int;
c->e_pot = e_pot;
c->mom[0] = mom[0];
c->mom[1] = mom[1];
c->mom[2] = mom[2];
c->ang[0] = ang[0];
c->ang[1] = ang[1];
c->ang[2] = ang[2];
c->ang_mom[0] = ang_mom[0];
c->ang_mom[1] = ang_mom[1];
c->ang_mom[2] = ang_mom[2];
}
void engine_print_stats(struct engine *e) {
const struct space *s = e->s;
double e_kin = 0.0, e_int = 0.0, e_pot = 0.0, mass = 0.0;
double mom[3] = {0.0, 0.0, 0.0}, ang_mom[3] = {0.0, 0.0, 0.0};
/* Collect the cell data. */
for (int k = 0; k < s->nr_cells; k++)
if (s->cells[k].nodeID == e->nodeID) {
struct cell *c = &s->cells[k];
/* Make the top-cells recurse */
engine_collect_drift(c);
/* And aggregate */
mass += c->mass;
e_kin += c->e_kin;
e_int += c->e_int;
e_pot += c->e_pot;
mom[0] += c->mom[0];
mom[1] += c->mom[1];
mom[2] += c->mom[2];
ang_mom[0] += c->ang_mom[0];
ang_mom[1] += c->ang_mom[1];
ang_mom[2] += c->ang_mom[2];
}
/* Aggregate the data from the different nodes. */
#ifdef WITH_MPI
{
double in[10] = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0.};
double out[10];
out[0] = e_kin;
out[1] = e_int;
out[2] = e_pot;
out[3] = mom[0];
out[4] = mom[1];
out[5] = mom[2];
out[6] = ang_mom[0];
out[7] = ang_mom[1];
out[8] = ang_mom[2];
out[9] = mass;
if (MPI_Allreduce(out, in, 10, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD) !=
MPI_SUCCESS)
error("Failed to aggregate stats.");
e_kin = out[0];
e_int = out[1];
e_pot = out[2];
mom[0] = out[3];
mom[1] = out[4];
mom[2] = out[5];
ang_mom[0] = out[6];
ang_mom[1] = out[7];
ang_mom[2] = out[8];
mass = out[9];
}
#endif
const double e_tot = e_kin + e_int + e_pot;
/* Print info */
fprintf(e->file_stats,
" %6d %14e %14e %14e %14e %14e %14e %14e %14e %14e %14e %14e %14e\n",
e->step, e->time, mass, e_tot, e_kin, e_int, e_pot, mom[0], mom[1], mom[2],
ang_mom[0], ang_mom[1], ang_mom[2]);
fflush(e->file_stats);
}
/**
......@@ -1865,7 +2023,6 @@ void engine_collect_kick(struct cell *c) {
* @param mask The task mask to launch.
* @param submask The sub-task mask to launch.
*/
void engine_launch(struct engine *e, int nr_runners, unsigned int mask,
unsigned int submask) {
......@@ -1984,13 +2141,7 @@ void engine_init_particles(struct engine *e) {
*/
void engine_step(struct engine *e) {
int updates = 0, g_updates = 0;
int ti_end_min = max_nr_timesteps, ti_end_max = 0;
double e_pot = 0.0, e_int = 0.0, e_kin = 0.0;
float mom[3] = {0.0, 0.0, 0.0};
float ang[3] = {0.0, 0.0, 0.0};
double snapshot_drift_time = 0.;
struct space *s = e->s;
TIMER_TIC2;
......@@ -1998,67 +2149,12 @@ void engine_step(struct engine *e) {
clocks_gettime(&time1);
e->tic_step = getticks();
/* Collect the cell data. */
for (int k = 0; k < s->nr_cells; k++)
if (s->cells[k].nodeID == e->nodeID) {
struct cell *c = &s->cells[k];
/* Recurse */
engine_collect_kick(c);
/* And aggregate */
ti_end_min = min(ti_end_min, c->ti_end_min);
ti_end_max = max(ti_end_max, c->ti_end_max);
e_kin += c->e_kin;
e_int += c->e_int;
e_pot += c->e_pot;
updates += c->updated;
g_updates += c->g_updated;
mom[0] += c->mom[0];
mom[1] += c->mom[1];
mom[2] += c->mom[2];
ang[0] += c->ang[0];
ang[1] += c->ang[1];
ang[2] += c->ang[2];
}
/* Aggregate the data from the different nodes. */
#ifdef WITH_MPI
{
int in_i[1], out_i[1];
in_i[0] = 0;
out_i[0] = ti_end_min;
if (MPI_Allreduce(out_i, in_i, 1, MPI_INT, MPI_MIN, MPI_COMM_WORLD) !=
MPI_SUCCESS)
error("Failed to aggregate t_end_min.");
ti_end_min = in_i[0];
out_i[0] = ti_end_max;
if (MPI_Allreduce(out_i, in_i, 1, MPI_INT, MPI_MAX, MPI_COMM_WORLD) !=
MPI_SUCCESS)
error("Failed to aggregate t_end_max.");
ti_end_max = in_i[0];
}
{
double in_d[5], out_d[5];
out_d[0] = updates;
out_d[1] = g_updates;
out_d[2] = e_kin;
out_d[3] = e_int;
out_d[4] = e_pot;
if (MPI_Allreduce(out_d, in_d, 5, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD) !=
MPI_SUCCESS)
error("Failed to aggregate energies.");
updates = in_d[0];
g_updates = in_d[1];
e_kin = in_d[2];
e_int = in_d[3];
e_pot = in_d[4];
}
#endif
/* Recover the (integer) end of the next time-step */
engine_collect_timestep(e);
/* Check for output */
while (ti_end_min >= e->ti_nextSnapshot && e->ti_nextSnapshot > 0) {
while (e->ti_end_min >= e->ti_nextSnapshot && e->ti_nextSnapshot > 0) {
e->ti_old = e->ti_current;
e->ti_current = e->ti_nextSnapshot;
......@@ -2079,7 +2175,7 @@ void engine_step(struct engine *e) {
/* Move forward in time */
e->ti_old = e->ti_current;
e->ti_current = ti_end_min;
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;
......@@ -2091,17 +2187,14 @@ void engine_step(struct engine *e) {
if (e->nodeID == 0) {
/* Print some information to the screen */
printf(" %6d %14e %14e %10d %10d %21.3f\n", e->step, e->time, e->timeStep,
updates, g_updates, e->wallclock_time);
printf(" %6d %14e %14e %10zd %10zd %21.3f\n", e->step, e->time,
e->timeStep, e->updates, e->g_updates, e->wallclock_time);
fflush(stdout);
/* Write some energy statistics */
fprintf(e->file_stats, "%d %f %f %f %f %f %f %f %f %f %f %f\n", e->step,
e->time, e_kin, e_int, e_pot, e_kin + e_int + e_pot, mom[0], mom[1],
mom[2], ang[0], ang[1], ang[2]);
fflush(e->file_stats);
}
/* Save some statistics */
engine_print_stats(e);
/* Re-distribute the particles amongst the nodes? */
if (e->forcerepart != REPART_NONE) engine_repartition(e);
......@@ -2574,8 +2667,10 @@ void engine_init(struct engine *e, struct space *s,
if (e->nodeID == 0) {
e->file_stats = fopen("energy.txt", "w");
fprintf(e->file_stats,
"# Step Time E_kin E_int E_pot E_tot "
"p_x p_y p_z ang_x ang_y ang_z\n");
"# %6s %14s %14s %14s %14s %14s %14s %14s %14s %14s %14s %14s %14s\n",
"Step", "Time", "Mass", "E_tot", "E_kin", "E_int", "E_pot", "p_x", "p_y", "p_z",
"ang_x", "ang_y", "ang_z");
fflush(e->file_stats);
}
/* Print policy */
......
......@@ -134,6 +134,12 @@ struct engine {
double timeBase;
double timeBase_inv;
/* Minimal ti_end for the next time-step */
int ti_end_min;
/* Number of particles updated */
size_t updates, g_updates;
/* Snapshot information */
double timeFirstSnapshot;
double deltaTimeSnapshot;
......
......@@ -872,13 +872,33 @@ void runner_do_drift(struct runner *r, struct cell *c, int timer) {
/* Collect */
dx_max = fmaxf(dx_max, cp->dx_max);
h_max = fmaxf(h_max, cp->h_max);
mass += cp->mass;
e_kin += cp->e_kin;
e_int += cp->e_int;
e_pot += cp->e_pot;
mom[0] += cp->mom[0];
mom[1] += cp->mom[1];
mom[2] += cp->mom[2];
ang_mom[0] += cp->ang_mom[0];
ang_mom[1] += cp->ang_mom[1];
ang_mom[2] += cp->ang_mom[2];
}
}
/* Store the values */
c->h_max = h_max;
c->dx_max = dx_max;
c->mass = mass;
c->e_kin = e_kin;
c->e_int = e_int;
c->e_pot = e_pot;
c->mom[0] = mom[0];
c->mom[1] = mom[1];
c->mom[2] = mom[2];
c->ang_mom[0] = ang_mom[0];
c->ang_mom[1] = ang_mom[1];
c->ang_mom[2] = ang_mom[2];
if (timer) TIMER_TOC(timer_drift);
}
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment