Commit 8def1372 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Merge branch 'fix_santa_barbara' into 'master'

Fixes to the cosmo hydro

See merge request !641
parents fd9a5416 ae6f5ce2
......@@ -41,7 +41,7 @@
*
* @param T The temperature of the gas [K].
* @param H_mass_fraction The hydrogen mass fraction of the gas.
* @param T_trans The temperature of the transition from HII to HI [K].
* @param T_transition The temperature of the transition from HII to HI [K].
*/
__attribute__((always_inline, const)) INLINE static double
mean_molecular_weight(const double T, const double H_mass_fraction,
......@@ -59,7 +59,7 @@ mean_molecular_weight(const double T, const double H_mass_fraction,
*
* @param u_cgs The internal energy per unit mass of the gas [erg * g^-1].
* @param H_mass_fraction The hydrogen mass fraction of the gas.
* @param T_trans The temperature of the transition from HII to HI [K].
* @param T_transition The temperature of the transition from HII to HI [K].
* @param m_H_cgs The mass of the Hydorgen atom [g].
* @param k_B_cgs The Boltzmann constant in cgs units [erg * K^-1]
* @return The temperature of the gas [K]
......
......@@ -157,6 +157,8 @@ void cosmology_update(struct cosmology *c, const struct phys_const *phys_const,
pow(a, -3. * hydro_gamma + 2.); /* 1 / a^(3*gamma - 2) */
c->a_factor_mu =
pow(a, 0.5 * (3. * hydro_gamma - 5.)); /* a^{(3*gamma - 5) / 2} */
c->a_factor_Balsara_eps =
pow(a, 0.5 * (1. - 3. * hydro_gamma)); /* a^{(1 - 3*gamma) / 2} */
/* Redshift */
c->z = a_inv - 1.;
......@@ -507,6 +509,10 @@ void cosmology_init(struct swift_params *params, const struct unit_system *us,
c->H0 = H0_cgs * units_cgs_conversion_factor(us, UNIT_CONV_TIME);
c->Hubble_time = 1. / c->H0;
/* Critical density at present day */
c->critical_density_0 =
3. * c->H0 * c->H0 / (8. * M_PI * phys_const->const_newton_G);
/* Initialise the interpolation tables */
c->drift_fac_interp_table = NULL;
c->grav_kick_fac_interp_table = NULL;
......@@ -548,6 +554,7 @@ void cosmology_init_no_cosmo(struct cosmology *c) {
c->log_a_end = 0.;
c->H = 0.;
c->H0 = 0.;
c->a = 1.;
c->z = 0.;
c->a_inv = 1.;
......@@ -557,10 +564,12 @@ void cosmology_init_no_cosmo(struct cosmology *c) {
c->a_factor_pressure = 1.;
c->a_factor_sound_speed = 1.;
c->a_factor_mu = 1.;
c->a_factor_Balsara_eps = 1.;
c->a_factor_hydro_accel = 1.;
c->a_factor_grav_accel = 1.;
c->critical_density = 0.;
c->critical_density_0 = 0.;
c->time_step_factor = 1.;
......
......@@ -54,10 +54,12 @@ struct cosmology {
/*! Power of the scale-factor used for sound-speed conversion to physical */
double a_factor_sound_speed;
/*! Power of the scale-factor used for relative velocities in viscosity term
*/
/*! Power of the scale-factor used for relative velocities in visc. terms */
double a_factor_mu;
/*! {ower of the scale-factor used for epsilon term in the Balsara switch */
double a_factor_Balsara_eps;
/*! Power of the scale-factor used for gravity accelerations */
double a_factor_grav_accel;
......@@ -73,6 +75,9 @@ struct cosmology {
/*! The critical density at the current redshift (in internal units) */
double critical_density;
/*! The critical density at redshift 0 (in internal units) */
double critical_density_0;
/*! Conversion factor from internal time-step size to cosmological step */
double time_step_factor;
......
......@@ -481,18 +481,22 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
const struct cosmology *cosmo, const struct hydro_props *hydro_props,
const float dt_alpha) {
const float fac_mu = cosmo->a_factor_mu;
const float fac_Balsara_eps = cosmo->a_factor_Balsara_eps;
/* Inverse of the physical density */
/* Inverse of the co-moving density */
const float rho_inv = 1.f / p->rho;
/* Inverse of the smoothing length */
const float h_inv = 1.f / p->h;
/* Compute the norm of the curl */
const float curl_v = sqrtf(p->density.rot_v[0] * p->density.rot_v[0] +
p->density.rot_v[1] * p->density.rot_v[1] +
p->density.rot_v[2] * p->density.rot_v[2]);
/* Compute the norm of div v */
const float abs_div_v = fabsf(p->density.div_v);
/* Compute the norm of div v including the Hubble flow term */
const float div_physical_v = p->density.div_v + 3.f * cosmo->H;
const float abs_div_physical_v = fabsf(div_physical_v);
/* Compute the pressure */
const float pressure = gas_pressure_from_entropy(p->rho, p->entropy);
......@@ -506,9 +510,9 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
/* Compute the Balsara switch */
/* Pre-multiply in the AV factor; hydro_props are not passed to the iact
* functions */
const float balsara =
hydro_props->viscosity.alpha * abs_div_v /
(abs_div_v + curl_v + 0.0001f * fac_mu * soundspeed / p->h);
const float balsara = hydro_props->viscosity.alpha * abs_div_physical_v /
(abs_div_physical_v + curl_v +
0.0001f * fac_Balsara_eps * soundspeed * h_inv);
/* Compute the "grad h" term */
const float omega_inv =
......@@ -656,12 +660,13 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
xp->entropy_full += p->entropy_dt * dt_therm;
/* Apply the minimal energy limit */
const float density = p->rho * cosmo->a3_inv;
const float min_energy = hydro_props->minimal_internal_energy;
const float min_entropy =
gas_entropy_from_internal_energy(density, min_energy);
if (xp->entropy_full < min_entropy) {
xp->entropy_full = min_entropy;
const float physical_density = p->rho * cosmo->a3_inv;
const float min_physical_energy = hydro_props->minimal_internal_energy;
const float min_physical_entropy =
gas_entropy_from_internal_energy(physical_density, min_physical_energy);
const float min_comoving_entropy = min_physical_entropy; /* A' = A */
if (xp->entropy_full < min_comoving_entropy) {
xp->entropy_full = min_comoving_entropy;
p->entropy_dt = 0.f;
}
......@@ -693,20 +698,21 @@ __attribute__((always_inline)) INLINE static void hydro_convert_quantities(
struct part *restrict p, struct xpart *restrict xp,
const struct cosmology *cosmo, const struct hydro_props *hydro_props) {
/* We read u in the entropy field. We now get (comoving) S from (physical) u
* and (physical) rho. Note that comoving S == physical S */
/* We read u in the entropy field. We now get (comoving) A from (physical) u
* and (physical) rho. Note that comoving A (A') == physical A */
xp->entropy_full =
gas_entropy_from_internal_energy(p->rho * cosmo->a3_inv, p->entropy);
p->entropy = xp->entropy_full;
/* Apply the minimal energy limit */
const float density = p->rho * cosmo->a3_inv;
const float min_energy = hydro_props->minimal_internal_energy;
const float min_entropy =
gas_entropy_from_internal_energy(density, min_energy);
if (xp->entropy_full < min_entropy) {
xp->entropy_full = min_entropy;
p->entropy = min_entropy;
const float physical_density = p->rho * cosmo->a3_inv;
const float min_physical_energy = hydro_props->minimal_internal_energy;
const float min_physical_entropy =
gas_entropy_from_internal_energy(physical_density, min_physical_energy);
const float min_comoving_entropy = min_physical_entropy; /* A' = A */
if (xp->entropy_full < min_comoving_entropy) {
xp->entropy_full = min_comoving_entropy;
p->entropy = min_comoving_entropy;
p->entropy_dt = 0.f;
}
......
......@@ -525,8 +525,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
pj->force.v_sig = max(pj->force.v_sig, v_sig);
/* Change in entropy */
pi->entropy_dt += mj * visc_term * dvdr;
pj->entropy_dt += mi * visc_term * dvdr;
pi->entropy_dt += mj * visc_term * dvdr_Hubble;
pj->entropy_dt += mi * visc_term * dvdr_Hubble;
#ifdef DEBUG_INTERACTIONS_SPH
/* Update ngb counters */
......@@ -641,7 +641,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
pi->force.v_sig = max(pi->force.v_sig, v_sig);
/* Change in entropy */
pi->entropy_dt += mj * visc_term * dvdr;
pi->entropy_dt += mj * visc_term * dvdr_Hubble;
#ifdef DEBUG_INTERACTIONS_SPH
/* Update ngb counters */
......@@ -769,7 +769,7 @@ runner_iact_nonsym_1_vec_force(
vec_div(vec_mul(mj.v, vec_mul(dvdr.v, vec_mul(ri.v, wi_dr.v))), pjrho.v);
/* Change in entropy */
entropy_dt.v = vec_mul(mj.v, vec_mul(visc_term.v, dvdr.v));
entropy_dt.v = vec_mul(mj.v, vec_mul(visc_term.v, dvdr_Hubble.v));
/* Store the forces back on the particles. */
a_hydro_xSum->v = vec_mask_sub(a_hydro_xSum->v, piax.v, mask);
......@@ -981,8 +981,8 @@ runner_iact_nonsym_2_vec_force(
pjrho_2.v);
/* Change in entropy */
entropy_dt.v = vec_mul(mj.v, vec_mul(visc_term.v, dvdr.v));
entropy_dt_2.v = vec_mul(mj_2.v, vec_mul(visc_term_2.v, dvdr_2.v));
entropy_dt.v = vec_mul(mj.v, vec_mul(visc_term.v, dvdr_Hubble.v));
entropy_dt_2.v = vec_mul(mj_2.v, vec_mul(visc_term_2.v, dvdr_Hubble_2.v));
/* Store the forces back on the particles. */
if (mask_cond) {
......
......@@ -494,15 +494,19 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
const struct cosmology *cosmo, const struct hydro_props *hydro_props,
const float dt_alpha) {
const float fac_mu = cosmo->a_factor_mu;
const float fac_Balsara_eps = cosmo->a_factor_Balsara_eps;
/* Inverse of the smoothing length */
const float h_inv = 1.f / p->h;
/* Compute the norm of the curl */
const float curl_v = sqrtf(p->density.rot_v[0] * p->density.rot_v[0] +
p->density.rot_v[1] * p->density.rot_v[1] +
p->density.rot_v[2] * p->density.rot_v[2]);
/* Compute the norm of div v */
const float abs_div_v = fabsf(p->density.div_v);
/* Compute the norm of div v including the Hubble flow term */
const float div_physical_v = p->density.div_v + 3.f * cosmo->H;
const float abs_div_physical_v = fabsf(div_physical_v);
/* Compute the pressure */
const float pressure = gas_pressure_from_internal_energy(p->rho, p->u);
......@@ -518,9 +522,9 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
/* Compute the Balsara switch */
/* Pre-multiply in the AV factor; hydro_props are not passed to the iact
* functions */
const float balsara =
hydro_props->viscosity.alpha * abs_div_v /
(abs_div_v + curl_v + 0.0001f * fac_mu * soundspeed / p->h);
const float balsara = hydro_props->viscosity.alpha * abs_div_physical_v /
(abs_div_physical_v + curl_v +
0.0001f * fac_Balsara_eps * soundspeed * h_inv);
/* Update variables. */
p->force.f = grad_h_term;
......@@ -662,10 +666,10 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
xp->u_full += p->u_dt * dt_therm;
/* Apply the minimal energy limit */
const float min_energy =
const float min_comoving_energy =
hydro_props->minimal_internal_energy / cosmo->a_factor_internal_energy;
if (xp->u_full < min_energy) {
xp->u_full = min_energy;
if (xp->u_full < min_comoving_energy) {
xp->u_full = min_comoving_energy;
p->u_dt = 0.f;
}
......@@ -704,11 +708,11 @@ __attribute__((always_inline)) INLINE static void hydro_convert_quantities(
xp->u_full = p->u;
/* Apply the minimal energy limit */
const float min_energy =
const float min_comoving_energy =
hydro_props->minimal_internal_energy / cosmo->a_factor_internal_energy;
if (xp->u_full < min_energy) {
xp->u_full = min_energy;
p->u = min_energy;
if (xp->u_full < min_comoving_energy) {
xp->u_full = min_comoving_energy;
p->u = min_comoving_energy;
p->u_dt = 0.f;
}
......
......@@ -266,7 +266,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
const float sph_du_term_j = P_over_rho2_j * dvdr * r_inv * wj_dr;
/* Viscosity term */
const float visc_du_term = 0.5f * visc_acc_term * dvdr;
const float visc_du_term = 0.5f * visc_acc_term * dvdr_Hubble;
/* Assemble the energy equation term */
const float du_dt_i = sph_du_term_i + visc_du_term;
......@@ -381,7 +381,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
const float sph_du_term_i = P_over_rho2_i * dvdr * r_inv * wi_dr;
/* Viscosity term */
const float visc_du_term = 0.5f * visc_acc_term * dvdr;
const float visc_du_term = 0.5f * visc_acc_term * dvdr_Hubble;
/* Assemble the energy equation term */
const float du_dt_i = sph_du_term_i + visc_du_term;
......
......@@ -69,13 +69,13 @@ struct hydro_props {
/*! Minimal temperature allowed */
float minimal_temperature;
/*! Minimal internal energy per unit mass */
/*! Minimal physical internal energy per unit mass */
float minimal_internal_energy;
/*! Initial temperature */
float initial_temperature;
/*! Initial internal energy per unit mass */
/*! Initial physical internal energy per unit mass */
float initial_internal_energy;
/*! Primordial hydrogen mass fraction for initial energy conversion */
......
......@@ -95,7 +95,7 @@ const double const_earth_mass_cgs = 5.9724e27;
/*! Temperature of the CMB at present day [K] */
const double const_T_CMB_0_cgs = 2.7255;
/*! Primordial Helium fraction */
/*! Primordial Helium fraction [-] */
const double const_primordial_He_fraction_cgs = 0.245;
#endif /* SWIFT_PHYSICAL_CONSTANTS_CGS_H */
......@@ -2867,7 +2867,7 @@ void space_init(struct space *s, struct swift_params *params,
/* Are we generating gas from the DM-only ICs? */
if (generate_gas_in_ics) {
space_generate_gas(s, cosmo, verbose);
space_generate_gas(s, cosmo, periodic, dim, verbose);
parts = s->parts;
gparts = s->gparts;
Npart = s->nr_parts;
......@@ -3180,10 +3180,12 @@ void space_replicate(struct space *s, int replicate, int verbose) {
*
* @param s The #space to create the particles in.
* @param cosmo The current #cosmology model.
* @param periodic Are we using periodic boundary conditions?
* @param dim The size of the box (for periodic wrapping).
* @param verbose Are we talkative?
*/
void space_generate_gas(struct space *s, const struct cosmology *cosmo,
int verbose) {
int periodic, const double dim[3], int verbose) {
/* Check that this is a sensible ting to do */
if (!s->hydro)
......@@ -3225,7 +3227,7 @@ void space_generate_gas(struct space *s, const struct cosmology *cosmo,
/* Compute some constants */
const double mass_ratio = cosmo->Omega_b / cosmo->Omega_m;
const double bg_density = cosmo->Omega_m * cosmo->critical_density;
const double bg_density = cosmo->Omega_m * cosmo->critical_density_0;
const double bg_density_inv = 1. / bg_density;
/* Update the particle properties */
......@@ -3239,9 +3241,11 @@ void space_generate_gas(struct space *s, const struct cosmology *cosmo,
p->id = gp_gas->id_or_neg_offset * 2 + 1;
gp_dm->id_or_neg_offset *= 2;
if (gp_dm->id_or_neg_offset <= 0) error("DM particle ID overflowd");
if (gp_dm->id_or_neg_offset < 0)
error("DM particle ID overflowd (DM id=%lld gas id=%lld)",
gp_dm->id_or_neg_offset, p->id);
if (p->id <= 0) error("gas particle ID overflowd");
if (p->id < 0) error("gas particle ID overflowd (id=%lld)", p->id);
/* Set the links correctly */
p->gpart = gp_gas;
......@@ -3250,8 +3254,8 @@ void space_generate_gas(struct space *s, const struct cosmology *cosmo,
/* Compute positions shift */
const double d = cbrt(gp_dm->mass * bg_density_inv);
const double shift_dm = d * mass_ratio;
const double shift_gas = d * (1. - mass_ratio);
const double shift_dm = 0.5 * d * mass_ratio;
const double shift_gas = 0.5 * d * (1. - mass_ratio);
/* Set the masses */
gp_dm->mass *= (1. - mass_ratio);
......@@ -3262,20 +3266,35 @@ void space_generate_gas(struct space *s, const struct cosmology *cosmo,
gp_dm->x[0] += shift_dm;
gp_dm->x[1] += shift_dm;
gp_dm->x[2] += shift_dm;
gp_gas->x[0] += shift_gas;
gp_gas->x[1] += shift_gas;
gp_gas->x[2] += shift_gas;
gp_gas->x[0] -= shift_gas;
gp_gas->x[1] -= shift_gas;
gp_gas->x[2] -= shift_gas;
/* Make sure the positions are identical between linked particles */
p->x[0] = gp_gas->x[0];
p->x[1] = gp_gas->x[1];
p->x[2] = gp_gas->x[2];
/* Box-wrap the whole thing to be safe */
if (periodic) {
gp_dm->x[0] = box_wrap(gp_dm->x[0], 0., dim[0]);
gp_dm->x[1] = box_wrap(gp_dm->x[1], 0., dim[1]);
gp_dm->x[2] = box_wrap(gp_dm->x[2], 0., dim[2]);
gp_gas->x[0] = box_wrap(gp_gas->x[0], 0., dim[0]);
gp_gas->x[1] = box_wrap(gp_gas->x[1], 0., dim[1]);
gp_gas->x[2] = box_wrap(gp_gas->x[2], 0., dim[2]);
p->x[0] = box_wrap(p->x[0], 0., dim[0]);
p->x[1] = box_wrap(p->x[1], 0., dim[1]);
p->x[2] = box_wrap(p->x[2], 0., dim[2]);
}
/* Also copy the velocities */
p->v[0] = gp_gas->v_full[0];
p->v[1] = gp_gas->v_full[1];
p->v[2] = gp_gas->v_full[2];
/* Set the smoothing length to the mean inter-particle separation */
p->h = 30. * d;
p->h = d;
/* Note that the thermodynamic properties (u, S, ...) will be set later */
}
......
......@@ -268,7 +268,7 @@ void space_check_top_multipoles_drift_point(struct space *s,
void space_check_timesteps(struct space *s);
void space_replicate(struct space *s, int replicate, int verbose);
void space_generate_gas(struct space *s, const struct cosmology *cosmo,
int verbose);
int periodic, const double dim[3], int verbose);
void space_check_cosmology(struct space *s, const struct cosmology *cosmo,
int rank);
void space_reset_task_counters(struct space *s);
......
......@@ -27,7 +27,10 @@
void print_bytes(void *p, size_t len) {
printf("(");
for (size_t i = 0; i < len; ++i) printf("%02x", ((unsigned char *)p)[i]);
for (size_t i = 0; i < len; ++i) {
printf("%02x", ((unsigned char *)p)[i]);
if (i % 4 == 3) printf("|");
}
printf(")\n");
}
......@@ -162,8 +165,8 @@ void test(void) {
if (i_not_ok) {
printParticle_single(&pi, &xpi);
printParticle_single(&pi2, &xpi);
print_bytes(&pj, sizeof(struct part));
print_bytes(&pj2, sizeof(struct part));
print_bytes(&pi, sizeof(struct part));
print_bytes(&pi2, sizeof(struct part));
error("Particles 'pi' do not match after density (byte = %d)", i_not_ok);
}
if (j_not_ok) {
......@@ -220,17 +223,15 @@ void test(void) {
j_not_ok |= c_is_d;
}
#else
i_not_ok =
strncmp((const char *)&pi, (const char *)&pi2, sizeof(struct part));
j_not_ok =
strncmp((const char *)&pj, (const char *)&pj2, sizeof(struct part));
i_not_ok = memcmp((char *)&pi, (char *)&pi2, sizeof(struct part));
j_not_ok = memcmp((char *)&pj, (char *)&pj2, sizeof(struct part));
#endif
if (i_not_ok) {
printParticle_single(&pi, &xpi);
printParticle_single(&pi2, &xpi);
print_bytes(&pj, sizeof(struct part));
print_bytes(&pj2, sizeof(struct part));
print_bytes(&pi, sizeof(struct part));
print_bytes(&pi2, sizeof(struct part));
error("Particles 'pi' do not match after force (byte = %d)", i_not_ok);
}
if (j_not_ok) {
......
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