diff --git a/src/cell.h b/src/cell.h
index a69a1a74648d76aeb13126d37a3c1265397dc2bf..cde4d75b7000a8d61b5391a7e7bf3f33a62ffad7 100644
--- a/src/cell.h
+++ b/src/cell.h
@@ -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;
diff --git a/src/engine.c b/src/engine.c
index 8441bdd596483e9eb4678a9827da32267291948b..d8121f0a9373b9d7e55986f562f0221809051ba0 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -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 */
diff --git a/src/engine.h b/src/engine.h
index 6c9ee32249f56246a917a7b347f3f5983cabecc5..92c7527e3c704f762e89bd65f59fe43286ba6239 100644
--- a/src/engine.h
+++ b/src/engine.h
@@ -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;
diff --git a/src/runner.c b/src/runner.c
index 9eaa942f5d1fa6c916d33070ebcd1e6b60b2e02a..148308cc66ccfd71899f766e7399e020273f1ba8 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -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);
}