diff --git a/src/hydro/Minimal/hydro.h b/src/hydro/Minimal/hydro.h
index 935b5d96843c5ed13b167ad1152c96b5ced647bd..571fce786bc14abfdf625b77890605e83745ec3b 100644
--- a/src/hydro/Minimal/hydro.h
+++ b/src/hydro/Minimal/hydro.h
@@ -35,6 +35,7 @@
#include "adiabatic_index.h"
#include "approx_math.h"
+#include "cosmology.h"
#include "dimension.h"
#include "equation_of_state.h"
#include "hydro_properties.h"
@@ -43,7 +44,7 @@
#include "minmax.h"
/**
- * @brief Returns the internal energy of a particle
+ * @brief Returns the comoving internal energy of a particle
*
* For implementations where the main thermodynamic variable
* is not internal energy, this function computes the internal
@@ -51,25 +52,61 @@
*
* @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 p->u;
}
/**
- * @brief Returns the pressure of a particle
+ * @brief Returns the physical internal energy of a particle
+ *
+ * For implementations where the main thermodynamic variable
+ * is not internal energy, this function computes the internal
+ * energy from the thermodynamic variable and converts it to
+ * physical coordinates.
+ *
+ * @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 p->u * cosmo->a_factor_internal_energy;
+}
+
+/**
+ * @brief Returns the comoving pressure of a particle
+ *
+ * Computes the pressure based on the particle's properties.
*
* @param p The particle of interest
*/
-__attribute__((always_inline)) INLINE static float hydro_get_pressure(
+__attribute__((always_inline)) INLINE static float hydro_get_comoving_pressure(
const struct part *restrict p) {
return gas_pressure_from_internal_energy(p->rho, p->u);
}
/**
- * @brief Returns the entropy of a particle
+ * @brief Returns the physical pressure of a particle
+ *
+ * Computes the pressure based on the particle's properties and
+ * convert it to physical coordinates.
+ *
+ * @param p The particle of interest
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_physical_pressure(
+ const struct part *restrict p, const struct cosmology *cosmo) {
+
+ return cosmo->a_factor_pressure *
+ gas_pressure_from_internal_energy(p->rho, p->u);
+}
+
+/**
+ * @brief Returns the comoving entropy of a particle
*
* For implementations where the main thermodynamic variable
* is not entropy, this function computes the entropy from
@@ -77,34 +114,78 @@ __attribute__((always_inline)) INLINE static float hydro_get_pressure(
*
* @param p The particle of interest
*/
-__attribute__((always_inline)) INLINE static float hydro_get_entropy(
+__attribute__((always_inline)) INLINE static float hydro_get_comoving_entropy(
const struct part *restrict p) {
return gas_entropy_from_internal_energy(p->rho, p->u);
}
/**
- * @brief Returns the sound speed of a particle
+ * @brief Returns the physical entropy of a particle
+ *
+ * For implementations where the main thermodynamic variable
+ * is not entropy, this function computes the entropy from
+ * the thermodynamic variable and converts it to
+ * physical coordinates.
*
* @param p The particle of interest
+ * @param cosmo The cosmological model.
*/
-__attribute__((always_inline)) INLINE static float hydro_get_soundspeed(
- const struct part *restrict p) {
+__attribute__((always_inline)) INLINE static float hydro_get_physical_entropy(
+ const struct part *restrict p, const struct cosmology *cosmo) {
+
+ /* Note: no cosmological conversion required here with our choice of
+ * coordinates. */
+ return gas_entropy_from_internal_energy(p->rho, p->u);
+}
+
+/**
+ * @brief Returns the comoving sound speed of a particle
+ *
+ * @param p The particle of interest
+ */
+__attribute__((always_inline)) INLINE static float
+hydro_get_comoving_soundspeed(const struct part *restrict p) {
return p->force.soundspeed;
}
/**
- * @brief Returns the density of a particle
+ * @brief Returns the physical sound speed of a particle
*
* @param p The particle of interest
+ * @param cosmo The cosmological model.
*/
-__attribute__((always_inline)) INLINE static float hydro_get_density(
+__attribute__((always_inline)) INLINE static float
+hydro_get_physical_soundspeed(const struct part *restrict p,
+ const struct cosmology *cosmo) {
+
+ return cosmo->a_factor_sound_speed * p->force.soundspeed;
+}
+
+/**
+ * @brief Returns the comoving density of a particle
+ *
+ * @param p The particle of interest
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_comoving_density(
const struct part *restrict p) {
return p->rho;
}
+/**
+ * @brief Returns the comoving 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 cosmo->a3_inv * p->rho;
+}
+
/**
* @brief Returns the mass of a particle
*
@@ -121,16 +202,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;
}
/**
@@ -168,17 +253,18 @@ __attribute__((always_inline)) INLINE static void hydro_set_internal_energy_dt(
* @param p Pointer to the particle data
* @param xp Pointer to the extended particle data
* @param hydro_properties The SPH parameters
- *
+ * @param cosmo The cosmological model.
*/
__attribute__((always_inline)) INLINE static float hydro_compute_timestep(
const struct part *restrict p, const struct xpart *restrict xp,
- const struct hydro_props *restrict hydro_properties) {
+ const struct hydro_props *restrict hydro_properties,
+ const struct cosmology *restrict cosmo) {
const float CFL_condition = hydro_properties->CFL_condition;
/* CFL condition */
- const float dt_cfl =
- 2.f * kernel_gamma * CFL_condition * p->h / p->force.v_sig;
+ const float dt_cfl = 2.f * kernel_gamma * CFL_condition * cosmo->a * p->h /
+ (cosmo->a_factor_sound_speed * p->force.v_sig);
return dt_cfl;
}
@@ -218,13 +304,15 @@ __attribute__((always_inline)) INLINE static void hydro_init_part(
* Multiplies the density and number of neighbours by the appropiate constants
* and add the self-contribution term.
* Additional quantities such as velocity gradients will also get the final
- *terms
- * added to them here.
+ * terms added to them here.
+ *
+ * Also adds/multiplies the cosmological terms if need be.
*
* @param p The particle to act upon
+ * @param cosmo The cosmological model.
*/
__attribute__((always_inline)) INLINE static void hydro_end_density(
- struct part *restrict p) {
+ struct part *restrict p, const struct cosmology *cosmo) {
/* Some smoothing length multiples. */
const float h = p->h;
@@ -248,11 +336,17 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
/**
* @brief Sets all particle fields to sensible values when the #part has 0 ngbs.
*
+ * In the desperate case where a particle has no neighbours (likely because
+ * of the h_max ceiling), set the particle fields to something sensible to avoid
+ * NaNs in the next calculations.
+ *
* @param p The particle to act upon
* @param xp The extended particle data to act upon
+ * @param cosmo The cosmological model.
*/
__attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours(
- struct part *restrict p, struct xpart *restrict xp) {
+ struct part *restrict p, struct xpart *restrict xp,
+ const struct cosmology *cosmo) {
/* Some smoothing length multiples. */
const float h = p->h;
@@ -278,9 +372,11 @@ __attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours(
*
* @param p The particle to act upon
* @param xp The extended particle data to act upon
+ * @param cosmo The current cosmological model.
*/
__attribute__((always_inline)) INLINE static void hydro_prepare_force(
- struct part *restrict p, struct xpart *restrict xp) {
+ struct part *restrict p, struct xpart *restrict xp,
+ const struct cosmology *cosmo) {
/* Compute the pressure */
const float pressure = gas_pressure_from_internal_energy(p->rho, p->u);
@@ -346,17 +442,22 @@ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values(
* Additional hydrodynamic quantites are drifted forward in time here. These
* include thermal quantities (thermal energy or total energy or entropy, ...).
*
+ * Note the different time-step sizes used for the different quantities as they
+ * include cosmological factors.
+ *
* @param p The particle.
* @param xp The extended data of the particle.
- * @param dt The drift time-step.
+ * @param dt_drift The drift time-step for positions.
+ * @param dt_therm The drift time-step for thermal quantities.
*/
__attribute__((always_inline)) INLINE static void hydro_predict_extra(
- struct part *restrict p, const struct xpart *restrict xp, float dt) {
+ struct part *restrict p, const struct xpart *restrict xp, float dt_drift,
+ float dt_therm) {
const float h_inv = 1.f / p->h;
/* Predict smoothing length */
- const float w1 = p->force.h_dt * h_inv * dt;
+ const float w1 = p->force.h_dt * h_inv * dt_drift;
if (fabsf(w1) < 0.2f)
p->h *= approx_expf(w1); /* 4th order expansion of exp(w) */
else
@@ -370,7 +471,7 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
p->rho *= expf(w2);
/* Predict the internal energy */
- p->u += p->u_dt * dt;
+ p->u += p->u_dt * dt_therm;
/* Compute the new pressure */
const float pressure = gas_pressure_from_internal_energy(p->rho, p->u);
@@ -386,13 +487,16 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
* @brief Finishes the force calculation.
*
* Multiplies the force and accelerations by the appropiate constants
- * and add the self-contribution term. In most cases, there is nothing
+ * and add the self-contribution term. In most cases, there is little
* to do here.
*
+ * Cosmological terms are also added/multiplied here.
+ *
* @param p The particle to act upon
+ * @param cosmo The current cosmological model.
*/
__attribute__((always_inline)) INLINE static void hydro_end_force(
- struct part *restrict p) {
+ struct part *restrict p, const struct cosmology *cosmo) {
p->force.h_dt *= p->h * hydro_dimension_inv;
}
@@ -403,18 +507,18 @@ __attribute__((always_inline)) INLINE static void hydro_end_force(
* Additional hydrodynamic quantites are kicked forward in time here. These
* include thermal quantities (thermal energy or total energy or entropy, ...).
*
- * @param p The particle to act upon
- * @param xp The particle extended data to act upon
- * @param dt The time-step for this kick
+ * @param p The particle to act upon.
+ * @param xp The particle extended data to act upon.
+ * @param dt_therm The time-step for this kick (for thermodynamic quantities).
*/
__attribute__((always_inline)) INLINE static void hydro_kick_extra(
- struct part *restrict p, struct xpart *restrict xp, float dt) {
+ struct part *restrict p, struct xpart *restrict xp, float dt_therm) {
/* Do not decrease the energy by more than a factor of 2*/
- if (dt > 0. && p->u_dt * dt < -0.5f * xp->u_full) {
- p->u_dt = -0.5f * xp->u_full / dt;
+ if (dt_therm > 0. && p->u_dt * dt_therm < -0.5f * xp->u_full) {
+ p->u_dt = -0.5f * xp->u_full / dt_therm;
}
- xp->u_full += p->u_dt * dt;
+ xp->u_full += p->u_dt * dt_therm;
/* Compute the pressure */
const float pressure = gas_pressure_from_internal_energy(p->rho, xp->u_full);
diff --git a/src/hydro/Minimal/hydro_debug.h b/src/hydro/Minimal/hydro_debug.h
index 876ce148824489d4c43358c2c519aa3b90dcf002..541029ee06dd2799443fc89b688d7baca3fae0f8 100644
--- a/src/hydro/Minimal/hydro_debug.h
+++ b/src/hydro/Minimal/hydro_debug.h
@@ -43,8 +43,9 @@ __attribute__((always_inline)) INLINE static void hydro_debug_particle(
"time_bin=%d\n",
p->x[0], p->x[1], p->x[2], p->v[0], p->v[1], p->v[2], xp->v_full[0],
xp->v_full[1], xp->v_full[2], p->a_hydro[0], p->a_hydro[1], p->a_hydro[2],
- p->u, p->u_dt, p->force.v_sig, p->force.pressure, p->h, p->force.h_dt,
- (int)p->density.wcount, p->mass, p->density.rho_dh, p->rho, p->time_bin);
+ p->u, p->u_dt, p->force.v_sig, hydro_get_comoving_pressure(p), p->h,
+ p->force.h_dt, (int)p->density.wcount, p->mass, p->density.rho_dh, p->rho,
+ p->time_bin);
}
#endif /* SWIFT_MINIMAL_HYDRO_DEBUG_H */
diff --git a/src/hydro/Minimal/hydro_iact.h b/src/hydro/Minimal/hydro_iact.h
index 0d4db1123d23e59146a2b026700c7b44a9e99f0c..bfd8f1dff5febbd14b8a1e3a1cb8c01d49bc85ba 100644
--- a/src/hydro/Minimal/hydro_iact.h
+++ b/src/hydro/Minimal/hydro_iact.h
@@ -36,10 +36,20 @@
#include "minmax.h"
/**
- * @brief Density loop
+ * @brief Density interaction between two particles.
+ *
+ * @param r2 Comoving square distance between the two particles.
+ * @param dx Comoving vector separating both particles (pi - pj).
+ * @param hi Comoving smoothing-length of particle i.
+ * @param hj Comoving smoothing-length of particle j.
+ * @param pi First particle.
+ * @param pj Second particle.
+ * @param a Current scale factor.
+ * @param H Current Hubble parameter.
*/
__attribute__((always_inline)) INLINE static void runner_iact_density(
- float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {
+ float r2, const float *dx, float hi, float hj, struct part *restrict pi,
+ struct part *restrict pj, float a, float H) {
float wi, wj, wi_dx, wj_dx;
@@ -71,10 +81,20 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
}
/**
- * @brief Density loop (non-symmetric version)
+ * @brief Density interaction between two particles (non-symmetric).
+ *
+ * @param r2 Comoving square distance between the two particles.
+ * @param dx Comoving vector separating both particles (pi - pj).
+ * @param hi Comoving smoothing-length of particle i.
+ * @param hj Comoving smoothing-length of particle j.
+ * @param pi First particle.
+ * @param pj Second particle (not updated).
+ * @param a Current scale factor.
+ * @param H Current Hubble parameter.
*/
__attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
- float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {
+ float r2, const float *dx, float hi, float hj, struct part *restrict pi,
+ const struct part *restrict pj, float a, float H) {
float wi, wi_dx;
@@ -95,12 +115,24 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
}
/**
- * @brief Force loop
+ * @brief Force interaction between two particles.
+ *
+ * @param r2 Comoving square distance between the two particles.
+ * @param dx Comoving vector separating both particles (pi - pj).
+ * @param hi Comoving smoothing-length of particle i.
+ * @param hj Comoving smoothing-length of particle j.
+ * @param pi First particle.
+ * @param pj Second particle.
+ * @param a Current scale factor.
+ * @param H Current Hubble parameter.
*/
__attribute__((always_inline)) INLINE static void runner_iact_force(
- float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {
+ float r2, const float *dx, float hi, float hj, struct part *restrict pi,
+ struct part *restrict pj, float a, float H) {
- const float fac_mu = 1.f; /* Will change with cosmological integration */
+ /* Cosmological factors entering the EoMs */
+ const float fac_mu = pow_three_gamma_minus_five_over_two(a);
+ const float a2_Hubble = a * a * H;
const float r = sqrtf(r2);
const float r_inv = 1.0f / r;
@@ -136,7 +168,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
/* Compute dv dot r. */
const float dvdr = (pi->v[0] - pj->v[0]) * dx[0] +
(pi->v[1] - pj->v[1]) * dx[1] +
- (pi->v[2] - pj->v[2]) * dx[2];
+ (pi->v[2] - pj->v[2]) * dx[2] + a2_Hubble * r2;
/* Are the particles moving towards each others ? */
const float omega_ij = min(dvdr, 0.f);
@@ -195,12 +227,24 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
}
/**
- * @brief Force loop (non-symmetric version)
+ * @brief Force interaction between two particles (non-symmetric).
+ *
+ * @param r2 Comoving square distance between the two particles.
+ * @param dx Comoving vector separating both particles (pi - pj).
+ * @param hi Comoving smoothing-length of particle i.
+ * @param hj Comoving smoothing-length of particle j.
+ * @param pi First particle.
+ * @param pj Second particle (not updated).
+ * @param a Current scale factor.
+ * @param H Current Hubble parameter.
*/
__attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
- float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {
+ float r2, const float *dx, float hi, float hj, struct part *restrict pi,
+ const struct part *restrict pj, float a, float H) {
- const float fac_mu = 1.f; /* Will change with cosmological integration */
+ /* Cosmological factors entering the EoMs */
+ const float fac_mu = pow_three_gamma_minus_five_over_two(a);
+ const float a2_Hubble = a * a * H;
const float r = sqrtf(r2);
const float r_inv = 1.0f / r;
@@ -236,7 +280,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
/* Compute dv dot r. */
const float dvdr = (pi->v[0] - pj->v[0]) * dx[0] +
(pi->v[1] - pj->v[1]) * dx[1] +
- (pi->v[2] - pj->v[2]) * dx[2];
+ (pi->v[2] - pj->v[2]) * dx[2] + a2_Hubble * r2;
/* Are the particles moving towards each others ? */
const float omega_ij = min(dvdr, 0.f);
diff --git a/src/hydro/Minimal/hydro_io.h b/src/hydro/Minimal/hydro_io.h
index d30d1e6ce6ecc00f5c370239e966ab1391055c3b..5922834e12f6dd74e951b46c2c93f975d7e8b82e 100644
--- a/src/hydro/Minimal/hydro_io.h
+++ b/src/hydro/Minimal/hydro_io.h
@@ -71,12 +71,12 @@ void hydro_read_particles(struct part* parts, struct io_props* list,
void convert_S(const struct engine* e, const struct part* p, float* ret) {
- ret[0] = hydro_get_entropy(p);
+ ret[0] = hydro_get_comoving_entropy(p);
}
void convert_P(const struct engine* e, const struct part* p, float* ret) {
- ret[0] = hydro_get_pressure(p);
+ ret[0] = hydro_get_comoving_pressure(p);
}
void convert_part_pos(const struct engine* e, const struct part* p,
diff --git a/src/hydro/Minimal/hydro_part.h b/src/hydro/Minimal/hydro_part.h
index 7556f415b9a3fb03923c98352322a9922f078f9a..c33f1b9a214cf9839f1acb965b686d4a4962865c 100644
--- a/src/hydro/Minimal/hydro_part.h
+++ b/src/hydro/Minimal/hydro_part.h
@@ -53,6 +53,9 @@ struct xpart {
/*! Velocity at the last full step. */
float v_full[3];
+ /*! Gravitational acceleration at the last full step. */
+ float a_grav[3];
+
/*! Internal energy at the last full step. */
float u_full;