diff --git a/src/cell.h b/src/cell.h index 9f35f4386cb161517b1b0c57468f1c0317c0c6a7..b9e060050827e6e0bbd62c87024bbf066d3807f8 100644 --- a/src/cell.h +++ b/src/cell.h @@ -152,7 +152,7 @@ struct cell { 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; + double mass, e_pot, e_int, e_kin, entropy; /* Number of particles updated in this cell. */ int updated, g_updated; diff --git a/src/engine.c b/src/engine.c index 8e471466fad5bac7c17603a42a59fba50d8c024e..6609f7f0d53453c6649a8089d2b656f555e00f75 100644 --- a/src/engine.c +++ b/src/engine.c @@ -2098,7 +2098,7 @@ void engine_collect_drift(struct cell *c) { 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 e_kin = 0.0, e_int = 0.0, e_pot = 0.0, entropy = 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 */ @@ -2120,6 +2120,7 @@ void engine_collect_drift(struct cell *c) { e_kin += cp->e_kin; e_int += cp->e_int; e_pot += cp->e_pot; + entropy += cp->entropy; mom[0] += cp->mom[0]; mom[1] += cp->mom[1]; mom[2] += cp->mom[2]; @@ -2135,6 +2136,7 @@ void engine_collect_drift(struct cell *c) { c->e_kin = e_kin; c->e_int = e_int; c->e_pot = e_pot; + c->entropy = entropy; c->mom[0] = mom[0]; c->mom[1] = mom[1]; c->mom[2] = mom[2]; @@ -2151,7 +2153,7 @@ 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 e_kin = 0.0, e_int = 0.0, e_pot = 0.0, entropy = 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. */ @@ -2167,6 +2169,7 @@ void engine_print_stats(struct engine *e) { e_kin += c->e_kin; e_int += c->e_int; e_pot += c->e_pot; + entropy += c->entropy; mom[0] += c->mom[0]; mom[1] += c->mom[1]; mom[2] += c->mom[2]; @@ -2178,8 +2181,8 @@ void engine_print_stats(struct engine *e) { /* 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]; + double in[11] = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}; + double out[11]; out[0] = e_kin; out[1] = e_int; out[2] = e_pot; @@ -2190,7 +2193,8 @@ void engine_print_stats(struct engine *e) { 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) != + out[10] = entropy; + if (MPI_Reduce(out, in, 11, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD) != MPI_SUCCESS) error("Failed to aggregate stats."); e_kin = out[0]; @@ -2203,6 +2207,7 @@ void engine_print_stats(struct engine *e) { ang_mom[1] = out[7]; ang_mom[2] = out[8]; mass = out[9]; + entropy = out[10]; } #endif @@ -2210,10 +2215,11 @@ void engine_print_stats(struct engine *e) { /* Print info */ if (e->nodeID == 0) { - fprintf(e->file_stats, - " %14e %14e %14e %14e %14e %14e %14e %14e %14e %14e %14e %14e\n", - 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]); + fprintf( + e->file_stats, + " %14e %14e %14e %14e %14e %14e %14e %14e %14e %14e %14e %14e %14e\n", + e->time, mass, e_tot, e_kin, e_int, e_pot, entropy, mom[0], mom[1], + mom[2], ang_mom[0], ang_mom[1], ang_mom[2]); fflush(e->file_stats); } } @@ -2973,10 +2979,11 @@ void engine_init(struct engine *e, struct space *s, engine_default_energy_file_name); sprintf(energyfileName + strlen(energyfileName), ".txt"); e->file_stats = fopen(energyfileName, "w"); - fprintf(e->file_stats, - "# %14s %14s %14s %14s %14s %14s %14s %14s %14s %14s %14s %14s\n", - "Time", "Mass", "E_tot", "E_kin", "E_int", "E_pot", "p_x", "p_y", - "p_z", "ang_x", "ang_y", "ang_z"); + fprintf( + e->file_stats, + "#%14s %14s %14s %14s %14s %14s %14s %14s %14s %14s %14s %14s %14s\n", + "Time", "Mass", "E_tot", "E_kin", "E_int", "E_pot", "Entropy", "p_x", + "p_y", "p_z", "ang_x", "ang_y", "ang_z"); fflush(e->file_stats); char timestepsfileName[200] = ""; diff --git a/src/hydro/Gadget2/hydro.h b/src/hydro/Gadget2/hydro.h index e998fae3e6078d4fa31141f1b31a3dc037313ada..7ff0ca03997b1526a4d741f31ff02f56bac7f7f1 100644 --- a/src/hydro/Gadget2/hydro.h +++ b/src/hydro/Gadget2/hydro.h @@ -282,3 +282,15 @@ __attribute__((always_inline)) INLINE static float hydro_get_internal_energy( return entropy * pow_gamma_minus_one(p->rho) * hydro_one_over_gamma_minus_one; } + +/** + * @brief Returns the entropy of a particle + * + * @param p The particle of interest + * @param dt Time since the last kick + */ +__attribute__((always_inline)) INLINE static float hydro_get_entropy( + const struct part *restrict p, float dt) { + + return p->entropy + p->entropy_dt * dt; +} diff --git a/src/runner.c b/src/runner.c index e48d2ccaba69315f4295a0dbacd36e9c7d26a0c4..e845d29288ecbd2b0d76587b42555c3777639864 100644 --- a/src/runner.c +++ b/src/runner.c @@ -593,7 +593,7 @@ void runner_do_drift(struct runner *r, struct cell *c, int timer) { struct gpart *restrict gparts = c->gparts; float dx_max = 0.f, dx2_max = 0.f, h_max = 0.f; - double e_kin = 0.0, e_int = 0.0, e_pot = 0.0, mass = 0.0; + double e_kin = 0.0, e_int = 0.0, e_pot = 0.0, entropy = 0.0, mass = 0.0; double mom[3] = {0.0, 0.0, 0.0}; double ang_mom[3] = {0.0, 0.0, 0.0}; @@ -670,6 +670,9 @@ void runner_do_drift(struct runner *r, struct cell *c, int timer) { e_kin += 0.5 * m * (v[0] * v[0] + v[1] * v[1] + v[2] * v[2]); e_pot += 0.; e_int += m * hydro_get_internal_energy(p, half_dt); + + /* Collect entropy */ + entropy += hydro_get_entropy(p, half_dt); } /* Now, get the maximal particle motion from its square */ @@ -694,6 +697,7 @@ void runner_do_drift(struct runner *r, struct cell *c, int timer) { e_kin += cp->e_kin; e_int += cp->e_int; e_pot += cp->e_pot; + entropy += cp->entropy; mom[0] += cp->mom[0]; mom[1] += cp->mom[1]; mom[2] += cp->mom[2]; @@ -710,6 +714,7 @@ void runner_do_drift(struct runner *r, struct cell *c, int timer) { c->e_kin = e_kin; c->e_int = e_int; c->e_pot = e_pot; + c->entropy = entropy; c->mom[0] = mom[0]; c->mom[1] = mom[1]; c->mom[2] = mom[2];