diff --git a/src/adiabatic_index.h b/src/adiabatic_index.h
index de7c3871cfb6e42739edf45a7d5b4882547d3cc2..db94a074d624e9107c0474543ccbc3d78f149db9 100644
--- a/src/adiabatic_index.h
+++ b/src/adiabatic_index.h
@@ -545,4 +545,39 @@ pow_three_gamma_minus_five_over_two(float x) {
#endif
}
+
+/**
+ * @brief Return the argument to the power three (adiabatic index - 1).
+ *
+ * Computes \f$x^{3(\gamma - 1)}\f$.
+ *
+ * @param x Argument
+ */
+__attribute__((always_inline, const)) INLINE static float
+pow_three_gamma_minus_one(float x) {
+
+#if defined(HYDRO_GAMMA_5_3)
+
+ return x * x; /* x^(2) */
+
+#elif defined(HYDRO_GAMMA_7_5)
+
+ return powf(x, 1.2f); /* x^(6/5) */
+
+#elif defined(HYDRO_GAMMA_4_3)
+
+ return x; /* x^(1) */
+
+#elif defined(HYDRO_GAMMA_2_1)
+
+ return x * x * x; /* x^(3) */
+
+#else
+
+ error("The adiabatic index is not defined !");
+ return 0.f;
+
+#endif
+}
+
#endif /* SWIFT_ADIABATIC_INDEX_H */
diff --git a/src/hydro/Gadget2/hydro.h b/src/hydro/Gadget2/hydro.h
index 5e13e6c558f6767c12846d94b55e8765b512d667..0beeafa87cc265bf2d686f2b78b03f8dc578fbb3 100644
--- a/src/hydro/Gadget2/hydro.h
+++ b/src/hydro/Gadget2/hydro.h
@@ -390,7 +390,7 @@ hydro_set_drifted_physical_internal_energy(struct part *p,
/* Compute the pressure */
float comoving_pressure = gas_pressure_from_entropy(p->rho, p->entropy);
- comoving_pressure = pressure_floor_get_comoving_pressure(p, p->rho, comoving_pressure);
+ comoving_pressure = pressure_floor_get_comoving_pressure(p, comoving_pressure, cosmo);
/* Compute the sound speed */
const float soundspeed =
@@ -595,7 +595,7 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
/* Compute the pressure */
float comoving_pressure = gas_pressure_from_entropy(p->rho, p->entropy);
- comoving_pressure = pressure_floor_get_comoving_pressure(p, p->rho, comoving_pressure);
+ comoving_pressure = pressure_floor_get_comoving_pressure(p, comoving_pressure, cosmo);
/* Compute the sound speed */
const float soundspeed =
@@ -657,9 +657,11 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration(
*
* @param p The particle.
* @param xp The extended data of this particle.
+ * @param cosmo The cosmological model.
*/
__attribute__((always_inline)) INLINE static void hydro_reset_predicted_values(
- struct part *restrict p, const struct xpart *restrict xp) {
+ struct part *restrict p, const struct xpart *restrict xp,
+ const struct cosmology *cosmo) {
/* Re-set the predicted velocities */
p->v[0] = xp->v_full[0];
@@ -671,7 +673,7 @@ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values(
/* Re-compute the pressure */
float comoving_pressure = gas_pressure_from_entropy(p->rho, p->entropy);
- comoving_pressure = pressure_floor_get_comoving_pressure(p, p->rho, comoving_pressure);
+ comoving_pressure = pressure_floor_get_comoving_pressure(p, comoving_pressure, cosmo);
/* Compute the new sound speed */
const float soundspeed =
@@ -738,7 +740,7 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
/* Re-compute the pressure */
float comoving_pressure = gas_pressure_from_entropy(p->rho, p->entropy);
- comoving_pressure = pressure_floor_get_comoving_pressure(p, p->rho, comoving_pressure);
+ comoving_pressure = pressure_floor_get_comoving_pressure(p, comoving_pressure, cosmo);
/* Compute the new sound speed */
const float soundspeed =
@@ -850,7 +852,7 @@ __attribute__((always_inline)) INLINE static void hydro_convert_quantities(
/* Compute the pressure */
float comoving_pressure = gas_pressure_from_entropy(p->rho, p->entropy);
- comoving_pressure = pressure_floor_get_comoving_pressure(p, p->rho, comoving_pressure);
+ comoving_pressure = pressure_floor_get_comoving_pressure(p, comoving_pressure, cosmo);
/* Compute the sound speed */
const float soundspeed =
diff --git a/src/pressure_floor/GEAR/pressure_floor.h b/src/pressure_floor/GEAR/pressure_floor.h
index 0870182bae80228f0ed27396e6b93f4338fac188..873615cf7fb12cc0b45a1c0bce5a7cded397e858 100644
--- a/src/pressure_floor/GEAR/pressure_floor.h
+++ b/src/pressure_floor/GEAR/pressure_floor.h
@@ -21,9 +21,9 @@
/* Forward declaration */
__attribute__((always_inline)) static INLINE float pressure_floor_get_comoving_pressure(
- const struct part* p, const float rho, const float pressure);
+ const struct part* p, const float pressure, const struct cosmology *cosmo);
__attribute__((always_inline)) static INLINE float pressure_floor_get_physical_pressure(
- const struct part* p, const float rho, const float pressure, const struct cosmology *cosmo);
+ const struct part* p, const float pressure, const struct cosmology *cosmo);
#include "adiabatic_index.h"
#include "cosmology.h"
@@ -58,21 +58,23 @@ struct pressure_floor_properties {
* Note that the particle is not updated!!
*
* @param p The #part.
- * @param rho The physical or comoving density.
- * @param pressure The physical pressure without any pressure floor.
+ * @param pressure_physical The physical pressure without any pressure floor.
* @param cosmo The #cosmology model.
*
* @return The physical pressure with the floor.
*/
__attribute__((always_inline)) static INLINE float pressure_floor_get_physical_pressure(
- const struct part* p, const float rho, const float pressure, const struct cosmology *cosmo) {
+ const struct part* p, const float pressure_physical, const struct cosmology *cosmo) {
- /* Compute pressure floor */
- float floor = p->h * p->h * rho * pressure_floor_props.constants -
- p->pressure_floor_data.sigma2 * cosmo->a * cosmo->a;
+ const float h_phys = p->h * cosmo->a_inv;
+ const float rho = hydro_get_physical_density(p, cosmo);
+
+ /* Compute the pressure floor */
+ float floor = h_phys * h_phys * rho * pressure_floor_props.constants -
+ p->pressure_floor_data.sigma2;
floor *= rho * hydro_one_over_gamma;
- return fmaxf(pressure, floor);
+ return fmaxf(pressure_physical, floor);
}
/**
@@ -81,20 +83,22 @@ __attribute__((always_inline)) static INLINE float pressure_floor_get_physical_p
* Note that the particle is not updated!!
*
* @param p The #part.
- * @param rho The comoving density.
- * @param pressure The comoving pressure without any pressure floor.
+ * @param pressure_comoving The comoving pressure without any pressure floor.
*
* @return The physical or comoving pressure with the floor.
*/
__attribute__((always_inline)) static INLINE float pressure_floor_get_comoving_pressure(
- const struct part* p, const float rho, const float pressure) {
+ const struct part* p, const float pressure_comoving, const struct cosmology *cosmo) {
- /* Compute pressure floor */
+ const float a_coef = pow_three_gamma_minus_one(cosmo->a);
+ const float rho = hydro_get_comoving_density(p);
+
+ /* Compute the pressure floor */
float floor = p->h * p->h * rho * pressure_floor_props.constants -
- p->pressure_floor_data.sigma2;
- floor *= rho * hydro_one_over_gamma;
+ p->pressure_floor_data.sigma2 * cosmo->a * cosmo->a;
+ floor *= a_coef * rho * hydro_one_over_gamma;
- return fmaxf(pressure, floor);
+ return fmaxf(pressure_comoving, floor);
}
/**
@@ -160,6 +164,9 @@ __attribute__((always_inline)) INLINE static void pressure_floor_end_density(
/* To finish the turbulence estimation we devide by the density */
p->pressure_floor_data.sigma2 /=
pow_dimension(p->h) * hydro_get_comoving_density(p);
+
+ /* Add the cosmological term */
+ p->pressure_floor_data.sigma2 *= cosmo->a2_inv;
}
/**
diff --git a/src/pressure_floor/GEAR/pressure_floor_iact.h b/src/pressure_floor/GEAR/pressure_floor_iact.h
index 1ecf8c29ea34089ab0e62e3bf6114d018bfecac3..a1928cc1503062276f32e245f38620013a5cdb3f 100644
--- a/src/pressure_floor/GEAR/pressure_floor_iact.h
+++ b/src/pressure_floor/GEAR/pressure_floor_iact.h
@@ -53,13 +53,17 @@ __attribute__((always_inline)) INLINE static void runner_iact_pressure_floor(
/* Delta v */
float dv[3] = {pi->v[0] - pj->v[0], pi->v[1] - pj->v[1], pi->v[2] - pj->v[2]};
- /* Norms */
- const float norm_v2 = dv[0] * dv[0] + dv[1] * dv[1] + dv[2] * dv[2];
- const float norm_v = sqrtf(norm_v2);
- const float r = sqrtf(r2);
-
/* Compute the velocity dispersion */
- const float sigma2 = norm_v2 + H * H * r2 + 2 * H * r * norm_v;
+ const float a2H = a * a * H;
+ const float sigma[3] = {
+ dv[0] + a2H * dx[0],
+ dv[1] + a2H * dx[1],
+ dv[2] + a2H * dx[2]
+ };
+ const float sigma2 =
+ sigma[0] * sigma[0] +
+ sigma[1] * sigma[1] +
+ sigma[2] * sigma[2];
/* Compute the velocity dispersion */
pi->pressure_floor_data.sigma2 += sigma2 * wi * hydro_get_mass(pj);
@@ -91,13 +95,17 @@ runner_iact_nonsym_pressure_floor(float r2, const float *dx, float hi, float hj,
/* Delta v */
float dv[3] = {pi->v[0] - pj->v[0], pi->v[1] - pj->v[1], pi->v[2] - pj->v[2]};
- /* Norms */
- const float norm_v2 = dv[0] * dv[0] + dv[1] * dv[1] + dv[2] * dv[2];
- const float norm_v = sqrtf(norm_v2);
- const float r = sqrtf(r2);
-
/* Compute the velocity dispersion */
- const float sigma2 = norm_v2 + H * H * r2 + 2 * H * r * norm_v;
+ const float a2H = a * a * H;
+ const float sigma[3] = {
+ dv[0] + a2H * dx[0],
+ dv[1] + a2H * dx[1],
+ dv[2] + a2H * dx[2]
+ };
+ const float sigma2 =
+ sigma[0] * sigma[0] +
+ sigma[1] * sigma[1] +
+ sigma[2] * sigma[2];
/* Compute the velocity dispersion */
pi->pressure_floor_data.sigma2 += sigma2 * wi * hydro_get_mass(pj);
diff --git a/src/pressure_floor/GEAR/pressure_floor_struct.h b/src/pressure_floor/GEAR/pressure_floor_struct.h
index b33540a8ce55186c4c53bf674ccc40f7831b482c..1eb70b86dcfb79ce9818e1d7c598d49d5e654efd 100644
--- a/src/pressure_floor/GEAR/pressure_floor_struct.h
+++ b/src/pressure_floor/GEAR/pressure_floor_struct.h
@@ -25,7 +25,7 @@
*/
struct pressure_floor_part_data {
/*! Estimation of local turbulence (squared)
- * Units: length^2 / time^2 (comoving) */
+ * Units: length^2 / time^2 (physical) */
float sigma2;
};
diff --git a/src/runner.c b/src/runner.c
index 2f274992aba0864355e9fc2fa0c47730e1130c06..599521830aa9d4ceb2f5b580d62e62e14dcfed07 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -2961,7 +2961,7 @@ void runner_do_kick2(struct runner *r, struct cell *c, int timer) {
#endif
/* Prepare the values to be drifted */
- hydro_reset_predicted_values(p, xp);
+ hydro_reset_predicted_values(p, xp, cosmo);
}
}