diff --git a/src/cosmology.c b/src/cosmology.c
index a818f0e5ea34fc1cfae27001d846b0fa54f39e81..c85a23788cd4e7300ccaeb7a2787974f43db767a 100644
--- a/src/cosmology.c
+++ b/src/cosmology.c
@@ -410,7 +410,7 @@ void cosmology_init(const struct swift_params *params,
/* Set remaining variables to alid values */
cosmology_update(c, 0);
-
+
/* Initialise the interpolation tables */
c->drift_fac_interp_table = NULL;
c->grav_kick_fac_interp_table = NULL;
@@ -461,7 +461,7 @@ void cosmology_init_no_cosmo(const struct swift_params *params,
c->a_factor_sig_vel = 1.;
c->a_factor_hydro_accel = 1.;
c->a_factor_grav_accel = 1.;
-
+
c->a_dot = 0.;
c->time = 0.;
c->universe_age_at_present_day = 0.;
diff --git a/src/engine.c b/src/engine.c
index 856c7dd0edd238fb587c82817ad6a02a1fe8745e..79c00681df15460ee2e69b634b5d8661b04b3c0f 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -4386,11 +4386,11 @@ void engine_step(struct engine *e) {
e->step += 1;
e->step_props = engine_step_prop_none;
- if(e->policy & engine_policy_cosmology) {
+ if (e->policy & engine_policy_cosmology) {
e->time_old = e->time;
cosmology_update(e->cosmology, e->ti_current);
e->time = e->cosmology->time;
- e->time_step = e->time - e->time_old;
+ e->time_step = e->time - e->time_old;
} else {
e->time = e->ti_current * e->time_base + e->time_begin;
e->time_old = e->ti_old * e->time_base + e->time_begin;
@@ -5262,7 +5262,7 @@ void engine_init(
e->time_end = e->cosmology->time_end;
e->time_old = e->time_begin;
e->time = e->time_begin;
-
+
/* Copy the relevent information from the cosmology model */
e->time_base = e->cosmology->time_base;
e->time_base_inv = e->cosmology->time_base_inv;
@@ -5768,25 +5768,25 @@ void engine_compute_next_snapshot_time(struct engine *e) {
/* Find upper-bound on last output */
double time_end;
- if(e->policy & engine_policy_cosmology)
+ if (e->policy & engine_policy_cosmology)
time_end = e->cosmology->a_end * e->deltaTimeSnapshot;
else
time_end = e->time_end + e->deltaTimeSnapshot;
/* Find next snasphot above current time */
- double time = e->timeFirstSnapshot;
- while(time < time_end) {
+ double time = e->timeFirstSnapshot;
+ while (time < time_end) {
/* Output time on the integer timeline */
- if(e->policy & engine_policy_cosmology)
+ if (e->policy & engine_policy_cosmology)
e->ti_nextSnapshot = log(time / e->cosmology->a_begin) / e->time_base;
else
e->ti_nextSnapshot = (time - e->time_begin) / e->time_base;
/* Found it? */
if (e->ti_nextSnapshot > e->ti_current) break;
-
- if(e->policy & engine_policy_cosmology)
+
+ if (e->policy & engine_policy_cosmology)
time *= e->deltaTimeSnapshot;
else
time += e->deltaTimeSnapshot;
@@ -5799,13 +5799,13 @@ void engine_compute_next_snapshot_time(struct engine *e) {
} else {
/* Be nice, talk... */
- if(e->policy & engine_policy_cosmology) {
+ if (e->policy & engine_policy_cosmology) {
const float next_snapshot_time =
- exp(e->ti_nextSnapshot * e->time_base) * e->cosmology->a_begin;
+ exp(e->ti_nextSnapshot * e->time_base) * e->cosmology->a_begin;
message("Next output time set to a=%e.", next_snapshot_time);
} else {
const float next_snapshot_time =
- e->ti_nextSnapshot * e->time_base + e->time_begin;
+ e->ti_nextSnapshot * e->time_base + e->time_begin;
message("Next output time set to t=%e.", next_snapshot_time);
}
}
diff --git a/src/hydro/Gadget2/hydro.h b/src/hydro/Gadget2/hydro.h
index 7af47b48924c4c013a3f4dd571838755a39861e3..7598f7b454bcda45c51c6090e97dd6843e9bb147 100644
--- a/src/hydro/Gadget2/hydro.h
+++ b/src/hydro/Gadget2/hydro.h
@@ -46,12 +46,25 @@
*
* @param p The particle of interest
*/
-__attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
- const struct part *restrict p) {
+__attribute__((always_inline)) INLINE static float
+hydro_get_comoving_internal_energy(const struct part *restrict p) {
return gas_internal_energy_from_entropy(p->rho, p->entropy);
}
+/**
+ * @brief Returns the internal energy of a particle
+ *
+ * @param p The particle of interest.
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static float
+hydro_get_physical_internal_energy(const struct part *restrict p,
+ const struct cosmology *cosmo) {
+
+ return gas_internal_energy_from_entropy(p->rho * cosmo->a3_inv, p->entropy);
+}
+
/**
* @brief Returns the pressure of a particle
*
@@ -66,10 +79,11 @@ __attribute__((always_inline)) INLINE static float hydro_get_pressure(
/**
* @brief Returns the entropy of a particle
*
- * @param p The particle of interest
+ * @param p The particle of interest.
+ * @param cosmo The cosmological model.
*/
__attribute__((always_inline)) INLINE static float hydro_get_entropy(
- const struct part *restrict p) {
+ const struct part *restrict p, const struct cosmology *cosmo) {
return p->entropy;
}
@@ -86,16 +100,28 @@ __attribute__((always_inline)) INLINE static float hydro_get_soundspeed(
}
/**
- * @brief Returns the density of a particle
+ * @brief Returns the comoving density of a particle
*
* @param p The particle of interest
*/
-__attribute__((always_inline)) INLINE static float hydro_get_density(
+__attribute__((always_inline)) INLINE static float hydro_get_comoving_density(
const struct part *restrict p) {
return p->rho;
}
+/**
+ * @brief Returns the physical density of a particle
+ *
+ * @param p The particle of interest.
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_physical_density(
+ const struct part *restrict p, const struct cosmology *cosmo) {
+
+ return p->rho * cosmo->a3_inv;
+}
+
/**
* @brief Returns the mass of a particle
*
@@ -112,16 +138,20 @@ __attribute__((always_inline)) INLINE static float hydro_get_mass(
*
* @param p The particle of interest
* @param xp The extended data of the particle.
- * @param dt The time since the last kick.
+ * @param dt_kick_hydro The time (for hydro accelerations) since the last kick.
+ * @param dt_kick_grav The time (for gravity accelerations) since the last kick.
* @param v (return) The velocities at the current time.
*/
__attribute__((always_inline)) INLINE static void hydro_get_drifted_velocities(
- const struct part *restrict p, const struct xpart *xp, float dt,
- float v[3]) {
-
- v[0] = xp->v_full[0] + p->a_hydro[0] * dt;
- v[1] = xp->v_full[1] + p->a_hydro[1] * dt;
- v[2] = xp->v_full[2] + p->a_hydro[2] * dt;
+ const struct part *restrict p, const struct xpart *xp, float dt_kick_hydro,
+ float dt_kick_grav, float v[3]) {
+
+ v[0] = xp->v_full[0] + p->a_hydro[0] * dt_kick_hydro +
+ xp->a_grav[0] * dt_kick_grav;
+ v[1] = xp->v_full[1] + p->a_hydro[1] * dt_kick_hydro +
+ xp->a_grav[1] * dt_kick_grav;
+ v[2] = xp->v_full[2] + p->a_hydro[2] * dt_kick_hydro +
+ xp->a_grav[2] * dt_kick_grav;
}
/**
diff --git a/src/hydro/Gadget2/hydro_io.h b/src/hydro/Gadget2/hydro_io.h
index 917e260fbd16526c7f2f6983b674df597555a41f..e4766a468ec66f3cbc2c5e83be694c0a1c298754 100644
--- a/src/hydro/Gadget2/hydro_io.h
+++ b/src/hydro/Gadget2/hydro_io.h
@@ -57,7 +57,7 @@ void hydro_read_particles(struct part* parts, struct io_props* list,
void convert_u(const struct engine* e, const struct part* p, float* ret) {
- ret[0] = hydro_get_internal_energy(p);
+ ret[0] = hydro_get_comoving_internal_energy(p);
}
void convert_P(const struct engine* e, const struct part* p, float* ret) {
diff --git a/src/statistics.c b/src/statistics.c
index 95217022a718829802395a807682c181ad69030e..66b0308f69deab84411d7647168537837d55305e 100644
--- a/src/statistics.c
+++ b/src/statistics.c
@@ -105,17 +105,24 @@ void stats_collect_part_mapper(void *map_data, int nr_parts, void *extra_data) {
/* Unpack the data */
const struct index_data *data = (struct index_data *)extra_data;
const struct space *s = data->s;
+ const struct engine *e = s->e;
+ const int with_cosmology = (e->policy & engine_policy_cosmology);
+ const integertime_t ti_current = e->ti_current;
+ const double time_base = e->time_base;
+ const double time = e->time;
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 time_base = s->e->time_base;
- const double time = s->e->time;
struct statistics *const global_stats = data->stats;
- /* Required for external potential energy */
- const struct external_potential *potential = s->e->external_potential;
- const struct phys_const *phys_const = s->e->physical_constants;
+ /* Some information about the physical model */
+ const struct external_potential *potential = e->external_potential;
+ const struct phys_const *phys_const = e->physical_constants;
+ const struct cosmology *cosmo = e->cosmology;
+
+ /* Some constants from cosmology */
+ const float a_inv = cosmo->a_inv;
+ const float a_inv2 = a_inv * a_inv;
/* Local accumulator */
struct statistics stats;
@@ -129,20 +136,36 @@ 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;
- // MATTHIEU Add cosmology here!
-
/* Get useful time variables */
- const integertime_t ti_begin =
+ const integertime_t ti_beg =
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 dt = (ti_current - ((ti_begin + ti_end) / 2)) * time_base;
+
+ /* Get time-step since the last kick */
+ float dt_kick_grav, dt_kick_hydro, dt_therm;
+ if (with_cosmology) {
+ dt_kick_grav = cosmology_get_grav_kick_factor(cosmo, ti_beg, ti_current);
+ dt_kick_grav -=
+ cosmology_get_grav_kick_factor(cosmo, ti_beg, (ti_beg + ti_end) / 2);
+ dt_kick_hydro =
+ cosmology_get_hydro_kick_factor(cosmo, ti_beg, ti_current);
+ dt_kick_hydro -=
+ cosmology_get_hydro_kick_factor(cosmo, ti_beg, (ti_beg + ti_end) / 2);
+ dt_therm = cosmology_get_therm_kick_factor(cosmo, ti_beg, ti_current);
+ dt_therm -=
+ cosmology_get_therm_kick_factor(cosmo, ti_beg, (ti_beg + ti_end) / 2);
+ } else {
+ dt_kick_grav = (ti_current - ((ti_beg + ti_end) / 2)) * time_base;
+ dt_kick_hydro = (ti_current - ((ti_beg + ti_end) / 2)) * time_base;
+ dt_therm = (ti_current - ((ti_beg + ti_end) / 2)) * time_base;
+ }
float v[3];
- hydro_get_drifted_velocities(p, xp, dt, v);
+ hydro_get_drifted_velocities(p, xp, dt_kick_hydro, dt_kick_grav, v);
const double x[3] = {p->x[0], p->x[1], p->x[2]};
const float m = hydro_get_mass(p);
- const float entropy = hydro_get_entropy(p);
- const float u_int = hydro_get_internal_energy(p);
+ const float entropy = hydro_get_entropy(p, cosmo);
+ const float u_int = hydro_get_physical_internal_energy(p, cosmo);
/* Collect mass */
stats.mass += m;
@@ -163,11 +186,12 @@ void stats_collect_part_mapper(void *map_data, int nr_parts, void *extra_data) {
stats.ang_mom[2] += m * (x[0] * v[1] - x[1] * v[0]);
/* Collect energies. */
- stats.E_kin += 0.5f * m * (v[0] * v[0] + v[1] * v[1] + v[2] * v[2]);
+ stats.E_kin += 0.5f * m * (v[0] * v[0] + v[1] * v[1] + v[2] * v[2]) *
+ a_inv2; /* 1/2 m a^2 \dot{r}^2 */
stats.E_int += m * u_int;
stats.E_rad += cooling_get_radiated_energy(xp);
if (gp != NULL) {
- stats.E_pot_self += m * gravity_get_potential(gp);
+ stats.E_pot_self += m * gravity_get_potential(gp) * a_inv;
stats.E_pot_ext += m * external_gravity_get_potential_energy(
time, potential, phys_const, gp);
}
@@ -194,15 +218,22 @@ void stats_collect_gpart_mapper(void *map_data, int nr_gparts,
/* Unpack the data */
const struct index_data *data = (struct index_data *)extra_data;
const struct space *s = data->s;
+ const struct engine *e = s->e;
+ const int with_cosmology = (e->policy & engine_policy_cosmology);
+ const integertime_t ti_current = e->ti_current;
+ const double time_base = e->time_base;
+ const double time = e->time;
const struct gpart *restrict gparts = (struct gpart *)map_data;
- const integertime_t ti_current = s->e->ti_current;
- const double time_base = s->e->time_base;
- const double time = s->e->time;
struct statistics *const global_stats = data->stats;
- /* Required for external potential energy */
- const struct external_potential *potential = s->e->external_potential;
- const struct phys_const *phys_const = s->e->physical_constants;
+ /* Some information about the physical model */
+ const struct external_potential *potential = e->external_potential;
+ const struct phys_const *phys_const = e->physical_constants;
+ const struct cosmology *cosmo = e->cosmology;
+
+ /* Some constants from cosmology */
+ const float a_inv = cosmo->a_inv;
+ const float a_inv2 = a_inv * a_inv;
/* Local accumulator */
struct statistics stats;
@@ -217,18 +248,25 @@ void stats_collect_gpart_mapper(void *map_data, int nr_gparts,
/* If the g-particle has a counterpart, ignore it */
if (gp->id_or_neg_offset < 0) continue;
- // MATTHIEU Add cosmology here!
-
/* Get useful variables */
- const integertime_t ti_begin =
+ const integertime_t ti_beg =
get_integer_time_begin(ti_current, gp->time_bin);
const integertime_t ti_end = get_integer_time_end(ti_current, gp->time_bin);
- const float dt = (ti_current - ((ti_begin + ti_end) / 2)) * time_base;
+
+ /* Get time-step since the last kick */
+ float dt_kick_grav;
+ if (with_cosmology) {
+ dt_kick_grav = cosmology_get_grav_kick_factor(cosmo, ti_beg, ti_current);
+ dt_kick_grav -=
+ cosmology_get_grav_kick_factor(cosmo, ti_beg, (ti_beg + ti_end) / 2);
+ } else {
+ dt_kick_grav = (ti_current - ((ti_beg + ti_end) / 2)) * time_base;
+ }
/* Extrapolate velocities */
- const float v[3] = {gp->v_full[0] + gp->a_grav[0] * dt,
- gp->v_full[1] + gp->a_grav[1] * dt,
- gp->v_full[2] + gp->a_grav[2] * dt};
+ const float v[3] = {gp->v_full[0] + gp->a_grav[0] * dt_kick_grav,
+ gp->v_full[1] + gp->a_grav[1] * dt_kick_grav,
+ gp->v_full[2] + gp->a_grav[2] * dt_kick_grav};
const float m = gravity_get_mass(gp);
const double x[3] = {gp->x[0], gp->x[1], gp->x[2]};
@@ -252,8 +290,9 @@ void stats_collect_gpart_mapper(void *map_data, int nr_gparts,
stats.ang_mom[2] += m * (x[0] * v[1] - x[1] * v[0]);
/* Collect energies. */
- stats.E_kin += 0.5f * m * (v[0] * v[0] + v[1] * v[1] + v[2] * v[2]);
- stats.E_pot_self += m * gravity_get_potential(gp);
+ stats.E_kin += 0.5f * m * (v[0] * v[0] + v[1] * v[1] + v[2] * v[2]) *
+ a_inv2; /* 1/2 m a^2 \dot{r}^2 */
+ stats.E_pot_self += m * gravity_get_potential(gp) * a_inv;
stats.E_pot_ext += m * external_gravity_get_potential_energy(
time, potential, phys_const, gp);
}
diff --git a/src/tools.c b/src/tools.c
index 89d89e62ea092ba6c6ec661e423d3e0ee44eb7fe..5c89b435313ec4287a8427661876172c8b165003 100644
--- a/src/tools.c
+++ b/src/tools.c
@@ -448,7 +448,7 @@ void engine_single_density(double *dim, long long int pid,
/* Dump the result. */
hydro_end_density(&p);
message("part %lli (h=%e) has wcount=%e, rho=%e.", p.id, p.h,
- p.density.wcount, hydro_get_density(&p));
+ p.density.wcount, hydro_get_comoving_density(&p));
fflush(stdout);
}