Commit ce252df4 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Added cosmological terms to the global statistics

parent ab39ee28
......@@ -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.;
......
......@@ -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);
}
}
......
......@@ -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;
}
/**
......
......@@ -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) {
......
......@@ -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);
}
......
......@@ -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);
}
......
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