diff --git a/src/cell.c b/src/cell.c
index a8281a292e2b77034998d8f8192e92eb12ff2779..50d3a441568acf8f20da905c0734965ed0a38993 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -686,9 +686,10 @@ void cell_sanitize(struct cell *c) {
void cell_convert_hydro(struct cell *c, void *data) {
struct part *p = c->parts;
+ struct xpart *xp = c->xparts;
for (int i = 0; i < c->count; ++i) {
- hydro_convert_quantities(&p[i]);
+ hydro_convert_quantities(&p[i], &xp[i]);
}
}
diff --git a/src/hydro/Gadget2/hydro.h b/src/hydro/Gadget2/hydro.h
index eba1a03663b007ca31ca16cac1b7f649cf3bb04b..581cb63d829dfec18b2371a128a29cd4f793b53b 100644
--- a/src/hydro/Gadget2/hydro.h
+++ b/src/hydro/Gadget2/hydro.h
@@ -48,7 +48,7 @@
__attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
const struct part *restrict p, float dt) {
- const float entropy = p->entropy + p->entropy_dt * dt;
+ const float entropy = p->entropy + dt * p->entropy_dt;
return gas_internal_energy_from_entropy(p->rho, entropy);
}
@@ -62,7 +62,7 @@ __attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
__attribute__((always_inline)) INLINE static float hydro_get_pressure(
const struct part *restrict p, float dt) {
- const float entropy = p->entropy + p->entropy_dt * dt;
+ const float entropy = p->entropy + dt * p->entropy_dt;
return gas_pressure_from_entropy(p->rho, entropy);
}
@@ -76,7 +76,7 @@ __attribute__((always_inline)) INLINE static float hydro_get_pressure(
__attribute__((always_inline)) INLINE static float hydro_get_entropy(
const struct part *restrict p, float dt) {
- return p->entropy + p->entropy_dt * dt;
+ return p->entropy + dt * p->entropy_dt;
}
/**
@@ -216,6 +216,7 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
xp->v_full[0] = p->v[0];
xp->v_full[1] = p->v[1];
xp->v_full[2] = p->v[2];
+ xp->entropy_full = p->entropy;
}
/**
@@ -228,9 +229,10 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
*/
__attribute__((always_inline)) INLINE static void hydro_init_part(
struct part *restrict p) {
+
+ p->rho = 0.f;
p->density.wcount = 0.f;
p->density.wcount_dh = 0.f;
- p->rho = 0.f;
p->density.rho_dh = 0.f;
p->density.div_v = 0.f;
p->density.rot_v[0] = 0.f;
@@ -302,11 +304,7 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
const float abs_div_v = fabsf(p->density.div_v);
/* Compute the pressure */
- const integertime_t ti_begin =
- get_integer_time_begin(ti_current, p->time_bin);
- const integertime_t ti_end = get_integer_time_end(ti_current, p->time_bin);
- const float half_dt = (ti_current - (ti_begin + ti_end) / 2) * timeBase;
- const float pressure = hydro_get_pressure(p, half_dt);
+ const float pressure = gas_pressure_from_entropy(p->rho, p->entropy);
/* Compute the sound speed */
const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure);
@@ -368,6 +366,9 @@ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values(
p->v[0] = xp->v_full[0];
p->v[1] = xp->v_full[1];
p->v[2] = xp->v_full[2];
+
+ /* Re-set the entropy */
+ p->entropy = xp->entropy_full;
}
/**
@@ -400,11 +401,11 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
else
p->rho *= expf(w2);
- /* Drift the pressure */
- const integertime_t ti_begin = get_integer_time_begin(t0, p->time_bin);
- const integertime_t ti_end = get_integer_time_end(t0, p->time_bin);
- const float dt_entr = (t1 - (ti_begin + ti_end) / 2) * timeBase;
- const float pressure = hydro_get_pressure(p, dt_entr);
+ /* Predict the entropy */
+ p->entropy += p->entropy_dt * dt;
+
+ /* Re-compute the pressure */
+ const float pressure = gas_pressure_from_entropy(p->rho, p->entropy);
/* Compute the new sound speed */
const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure);
@@ -446,13 +447,13 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
/* Do not decrease the entropy (temperature) by more than a factor of 2*/
const float entropy_change = p->entropy_dt * dt;
- if (entropy_change > -0.5f * p->entropy)
- p->entropy += entropy_change;
+ if (entropy_change > -0.5f * xp->entropy_full)
+ xp->entropy_full += entropy_change;
else
- p->entropy *= 0.5f;
+ xp->entropy_full *= 0.5f;
/* Compute the pressure */
- const float pressure = gas_pressure_from_entropy(p->rho, p->entropy);
+ const float pressure = gas_pressure_from_entropy(p->rho, xp->entropy_full);
/* Compute the new sound speed */
const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure);
@@ -473,10 +474,11 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
* @param p The particle to act upon
*/
__attribute__((always_inline)) INLINE static void hydro_convert_quantities(
- struct part *restrict p) {
+ struct part *restrict p, struct xpart *restrict xp) {
/* We read u in the entropy field. We now get S from u */
- p->entropy = gas_entropy_from_internal_energy(p->rho, p->entropy);
+ xp->entropy_full = gas_entropy_from_internal_energy(p->rho, p->entropy);
+ p->entropy = xp->entropy_full;
/* Compute the pressure */
const float pressure = gas_pressure_from_entropy(p->rho, p->entropy);
diff --git a/src/hydro/Gadget2/hydro_part.h b/src/hydro/Gadget2/hydro_part.h
index 56162055f4f200aff376306d121e7cd24a5cd79b..8fb4b1d7b7e60ee5fe80b41b150ff872ae433c60 100644
--- a/src/hydro/Gadget2/hydro_part.h
+++ b/src/hydro/Gadget2/hydro_part.h
@@ -42,6 +42,9 @@ struct xpart {
/* Velocity at the last full step. */
float v_full[3];
+ /* Entropy at the last full step. */
+ float entropy_full;
+
/* Additional data used to record cooling information */
struct cooling_xpart_data cooling_data;
diff --git a/src/statistics.c b/src/statistics.c
index 55b76eba974df19cdcd1e8ef24e5b6f8d0679c31..10ec276f567c62dda78f6fbd1b8b2f43360446a8 100644
--- a/src/statistics.c
+++ b/src/statistics.c
@@ -104,8 +104,8 @@ void stats_collect_part_mapper(void *map_data, int nr_parts, void *extra_data) {
const struct part *restrict parts = (struct part *)map_data;
const struct xpart *restrict xparts =
s->xparts + (ptrdiff_t)(parts - s->parts);
- // const integertime_t ti_current = s->e->ti_current;
- // const double timeBase = s->e->timeBase;
+ const integertime_t ti_current = s->e->ti_current;
+ const double timeBase = s->e->timeBase;
const double time = s->e->time;
struct statistics *const global_stats = data->stats;
@@ -125,32 +125,27 @@ void stats_collect_part_mapper(void *map_data, int nr_parts, void *extra_data) {
const struct xpart *xp = &xparts[k];
const struct gpart *gp = (p->gpart != NULL) ? gp = p->gpart : NULL;
- /* Get useful variables */
- // const integertime_t ti_begin =
- // get_integer_time_begin(ti_current, p->time_bin);
- // const integertime_t ti_end = get_integer_time_end(ti_current,
- // p->time_bin);
- // const integertime_t dti = get_integer_timestep(p->time_bin);
- const float dt =
- 0.f; //(ti_current - (ti_begin + ti_end) / 2) * timeBase; //MATTHIEU
- const double x[3] = {p->x[0], p->x[1], p->x[2]};
+ /* Get useful time variables */
+ const integertime_t ti_begin =
+ get_integer_time_begin(ti_current, p->time_bin);
+ const integertime_t ti_step = get_integer_timestep(p->time_bin);
+ const float dt = (ti_current - (ti_begin + ti_step / 2)) * timeBase;
+
+ /* Get the total acceleration */
float a_tot[3] = {p->a_hydro[0], p->a_hydro[1], p->a_hydro[2]};
if (gp != NULL) {
a_tot[0] += gp->a_grav[0];
a_tot[1] += gp->a_grav[1];
a_tot[2] += gp->a_grav[2];
}
+
+ /* Extrapolate velocities to current time */
const float v[3] = {xp->v_full[0] + a_tot[0] * dt,
xp->v_full[1] + a_tot[1] * dt,
xp->v_full[2] + a_tot[2] * dt};
const float m = hydro_get_mass(p);
-
- /* if (p->id == ICHECK) */
- /* message("bin=%d dti=%lld ti_begin=%lld ti_end=%lld dt=%e v=[%e %e %e]",
- */
- /* p->time_bin, dti, ti_begin, ti_end, */
- /* dt, v[0], v[1], v[2]); */
+ const double x[3] = {p->x[0], p->x[1], p->x[2]};
/* Collect mass */
stats.mass += m;
@@ -167,13 +162,12 @@ void stats_collect_part_mapper(void *map_data, int nr_parts, void *extra_data) {
/* Collect energies. */
stats.E_kin += 0.5f * m * (v[0] * v[0] + v[1] * v[1] + v[2] * v[2]);
+ stats.E_int += m * hydro_get_internal_energy(p, dt);
+ stats.E_rad += cooling_get_radiated_energy(xp);
stats.E_pot_self += 0.f;
if (gp != NULL)
stats.E_pot_ext += m * external_gravity_get_potential_energy(
time, potential, phys_const, gp);
- stats.E_int += m * hydro_get_internal_energy(p, dt);
- stats.E_rad += cooling_get_radiated_energy(xp);
-
/* Collect entropy */
stats.entropy += m * hydro_get_entropy(p, dt);
}