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

Correct initialisation of the internal energy in comoving coordinates for the Minimal-SPH scheme.

parent a4738f0e
......@@ -34,7 +34,15 @@ Gravity:
comoving_softening: 0.0889 # 1/25th of the mean inter-particle separation: 88.9 kpc
max_physical_softening: 0.0889 # 1/25th of the mean inter-particle separation: 88.9 kpc
mesh_side_length: 64
# Parameters of the hydro scheme
SPH:
resolution_eta: 1.2348 # "48 Ngb" with the cubic spline kernel
CFL_condition: 0.1
initial_temperature: 7075. # (1 + z_ini)^2 * 2.72K
minimal_temperature: 100.
# Parameters governing the snapshots
Snapshots:
basename: snap
......@@ -55,3 +63,5 @@ InitialConditions:
file_name: small_cosmo_volume.hdf5
cleanup_h_factors: 1
cleanup_velocity_factors: 1
generate_gas_in_ics: 1 # Generate gas particles from the DM-only ICs
cleanup_smoothing_lengths: 1 # Since we generate gas, make use of the (expensive) cleaning-up procedure.
......@@ -520,7 +520,7 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
*/
__attribute__((always_inline)) INLINE static void hydro_convert_quantities(
struct part *restrict p, struct xpart *restrict xp,
const struct cosmology *cosmo) {}
const struct cosmology *cosmo, const struct hydro_props *hydro_props) {}
/**
* @brief Initialises the particles for the first time
......
......@@ -571,10 +571,11 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
* @param p The particle to act upon.
* @param xp The extended data.
* @param cosmo The cosmological model.
* @param hydro_props The constants used in the scheme.
*/
__attribute__((always_inline)) INLINE static void hydro_convert_quantities(
struct part *restrict p, struct xpart *restrict xp,
const struct cosmology *cosmo) {
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 */
......@@ -589,7 +590,7 @@ __attribute__((always_inline)) INLINE static void hydro_convert_quantities(
gas_entropy_from_internal_energy(density, min_energy);
if (xp->entropy_full < min_entropy) {
xp->entropy_full = min_entropy;
p->entopy = min_entropy;
p->entropy = min_entropy;
p->entropy_dt = 0.f;
}
......
......@@ -501,7 +501,8 @@ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values(
* @param p The particle to act upon.
*/
__attribute__((always_inline)) INLINE static void hydro_convert_quantities(
struct part* p, struct xpart* xp, const struct cosmology* cosmo) {
struct part* p, struct xpart* xp, const struct cosmology* cosmo,
const struct hydro_props* hydro_props) {
p->conserved.energy /= cosmo->a_factor_internal_energy;
}
......
......@@ -543,7 +543,8 @@ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values(
* @param p The particle to act upon.
*/
__attribute__((always_inline)) INLINE static void hydro_convert_quantities(
struct part* p, struct xpart* xp, const struct cosmology* cosmo) {
struct part* p, struct xpart* xp, const struct cosmology* cosmo,
const struct hydro_props* hydro_props) {
p->conserved.energy /= cosmo->a_factor_internal_energy;
}
......
......@@ -541,7 +541,7 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
/* Apply the minimal energy limit */
const float min_energy =
hydro_props->minimal_internal_energy * cosmo->a_factor_internal_energy;
hydro_props->minimal_internal_energy / cosmo->a_factor_internal_energy;
if (xp->u_full < min_energy) {
xp->u_full = min_energy;
p->u_dt = 0.f;
......@@ -551,7 +551,8 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
const float pressure = gas_pressure_from_internal_energy(p->rho, xp->u_full);
/* Compute the sound speed */
const float soundspeed = gas_soundspeed_from_internal_energy(p->rho, p->u);
const float soundspeed =
gas_soundspeed_from_internal_energy(p->rho, xp->u_full);
p->force.pressure = pressure;
p->force.soundspeed = soundspeed;
......@@ -568,20 +569,21 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
* @param p The particle to act upon
* @param xp The extended particle to act upon
* @param cosmo The cosmological model.
* @param hydro_props The constants used in the scheme.
*/
__attribute__((always_inline)) INLINE static void hydro_convert_quantities(
struct part *restrict p, struct xpart *restrict xp,
const struct cosmology *cosmo) {
const struct cosmology *cosmo, const struct hydro_props *hydro_props) {
/* Convert the physcial internal energy to the comoving one. */
/* u' = a^(3(g-1)) u */
const float factor = 1.f / cosmo->a_factor_internal_energy;
p->u *= factor;
xp->u_full *= factor;
xp->u_full = p->u;
/* Apply the minimal energy limit */
const float min_energy =
hydro_props->minimal_internal_energy * cosmo->a_factor_internal_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;
......
......@@ -611,7 +611,7 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
*/
__attribute__((always_inline)) INLINE static void hydro_convert_quantities(
struct part *restrict p, struct xpart *restrict xp,
const struct cosmology *cosmo) {
const struct cosmology *cosmo, const struct hydro_props *hydro_props) {
/* Compute the pressure */
const float pressure =
......
......@@ -610,7 +610,7 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
*/
__attribute__((always_inline)) INLINE static void hydro_convert_quantities(
struct part *restrict p, struct xpart *restrict xp,
const struct cosmology *cosmo) {}
const struct cosmology *cosmo, const struct hydro_props *hydro_props) {}
/**
* @brief Initialises the particles for the first time
......
......@@ -587,7 +587,7 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
*/
__attribute__((always_inline)) INLINE static void hydro_convert_quantities(
struct part *restrict p, struct xpart *restrict xp,
const struct cosmology *cosmo) {
const struct cosmology *cosmo, const struct hydro_props *hydro_props) {
/* We read u in the entropy field. We now get S from u */
xp->entropy_full =
......
......@@ -411,7 +411,8 @@ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values(
* @param xp The extended particle data to act upon.
*/
__attribute__((always_inline)) INLINE static void hydro_convert_quantities(
struct part* p, struct xpart* xp, const struct cosmology* cosmo) {}
struct part* p, struct xpart* xp, const struct cosmology* cosmo,
const struct hydro_props* hydro_props) {}
/**
* @brief Extra operations to be done during the drift
......
......@@ -2705,12 +2705,13 @@ void space_convert_quantities_mapper(void *restrict map_data, int count,
void *restrict extra_data) {
struct space *s = (struct space *)extra_data;
const struct cosmology *cosmo = s->e->cosmology;
const struct hydro_props *hydro_props = s->e->hydro_properties;
struct part *restrict parts = (struct part *)map_data;
const ptrdiff_t index = parts - s->parts;
struct xpart *restrict xparts = s->xparts + index;
for (int k = 0; k < count; k++)
hydro_convert_quantities(&parts[k], &xparts[k], cosmo);
hydro_convert_quantities(&parts[k], &xparts[k], cosmo, hydro_props);
}
/**
......
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