diff --git a/src/equation_of_state.h b/src/equation_of_state.h index af59d8a2cad1632c67b6d377b5ed9dfe9484b4aa..5b19f99ab85d8a2b4e5e6f0b4eeb542925b4ee50 100644 --- a/src/equation_of_state.h +++ b/src/equation_of_state.h @@ -126,6 +126,21 @@ gas_soundspeed_from_internal_energy(float density, float u) { return sqrtf(u * hydro_gamma * hydro_gamma_minus_one); } +/** + * @brief Returns the sound speed given density and pressure + * + * Computes \f$c = \sqrt{\frac{\gamma P}{\rho} }\f$. + * + * @param density The density \f$\rho\f$ + * @param P The pressure \f$P\f$ + */ +__attribute__((always_inline)) INLINE static float gas_soundspeed_from_pressure( + float density, float P) { + + const float density_inv = 1.f / density; + return sqrtf(hydro_gamma * P * density_inv); +} + /* ------------------------------------------------------------------------- */ #elif defined(EOS_ISOTHERMAL_GAS) @@ -221,6 +236,22 @@ gas_soundspeed_from_internal_energy(float density, float u) { hydro_gamma_minus_one); } +/** + * @brief Returns the sound speed given density and pressure + * + * Since we are using an isothermal EoS, the pressure value is ignored + * Computes \f$c = \sqrt{u_{cst} \gamma \rho^{\gamma-1}}\f$. + * + * @param density The density \f$\rho\f$ + * @param P The pressure \f$P\f$ + */ +__attribute__((always_inline)) INLINE static float gas_soundspeed_from_pressure( + float density, float P) { + + return sqrtf(const_isothermal_internal_energy * hydro_gamma * + hydro_gamma_minus_one); +} + /* ------------------------------------------------------------------------- */ #else diff --git a/src/hydro/Minimal/hydro.h b/src/hydro/Minimal/hydro.h index d4dfc68fd35f235806c32cd0a39603371ddd2233..42dc2aa3d8899ce01d1b7f0bc55f5faa43e45a68 100644 --- a/src/hydro/Minimal/hydro.h +++ b/src/hydro/Minimal/hydro.h @@ -65,9 +65,7 @@ __attribute__((always_inline)) INLINE static float hydro_get_internal_energy( __attribute__((always_inline)) INLINE static float hydro_get_pressure( const struct part *restrict p, float dt) { - const float u = p->u + p->u_dt * dt; - - return gas_pressure_from_internal_energy(p->rho, u); + return p->force.pressure; } /** @@ -97,9 +95,7 @@ __attribute__((always_inline)) INLINE static float hydro_get_entropy( __attribute__((always_inline)) INLINE static float hydro_get_soundspeed( const struct part *restrict p, float dt) { - const float u = p->u + p->u_dt * dt; - - return gas_soundspeed_from_internal_energy(p->rho, u); + return p->force.soundspeed; } /** @@ -139,12 +135,11 @@ __attribute__((always_inline)) INLINE static void hydro_set_internal_energy( p->u = u; - /* Compute the pressure */ + /* Compute the new pressure */ const float pressure = gas_pressure_from_internal_energy(p->rho, p->u); - /* Compute the sound speed from the pressure*/ - const float rho_inv = 1.f / p->rho; - const float soundspeed = sqrtf(hydro_gamma * pressure * rho_inv); + /* Compute the new sound speed */ + const float soundspeed = gas_soundspeed_from_internal_energy(p->rho, p->u); p->force.soundspeed = soundspeed; p->force.pressure = pressure; @@ -167,9 +162,8 @@ __attribute__((always_inline)) INLINE static void hydro_set_entropy( /* Compute the pressure */ const float pressure = gas_pressure_from_internal_energy(p->rho, p->u); - /* Compute the sound speed from the pressure*/ - const float rho_inv = 1.f / p->rho; - const float soundspeed = sqrtf(hydro_gamma * pressure * rho_inv); + /* Compute the new sound speed */ + const float soundspeed = gas_soundspeed_from_internal_energy(p->rho, p->u); p->force.soundspeed = soundspeed; p->force.pressure = pressure; @@ -293,12 +287,11 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force( const float half_dt = (ti_current - (p->ti_begin + p->ti_end) / 2) * timeBase; const float pressure = hydro_get_pressure(p, half_dt); - const float rho_inv = 1.f / p->rho; - /* Compute the sound speed */ - const float soundspeed = sqrtf(hydro_gamma * pressure * rho_inv); + const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure); /* Compute the "grad h" term */ + const float rho_inv = 1.f / p->rho; const float grad_h_term = 1.f / (1.f + hydro_dimension_inv * p->h * p->density.rho_dh * rho_inv); @@ -367,10 +360,8 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra( const float dt_entr = (t1 - (p->ti_begin + p->ti_end) / 2) * timeBase; const float pressure = hydro_get_pressure(p, dt_entr); - const float rho_inv = 1.f / p->rho; - /* Compute the new sound speed */ - const float soundspeed = sqrtf(hydro_gamma * pressure * rho_inv); + const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure); p->force.pressure = pressure; p->force.soundspeed = soundspeed; @@ -419,12 +410,11 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra( /* Compute the pressure */ const float pressure = gas_pressure_from_internal_energy(p->rho, p->u); - /* Compute the sound speed from the pressure*/ - const float rho_inv = 1.f / p->rho; - const float soundspeed = sqrtf(hydro_gamma * pressure * rho_inv); + /* Compute the sound speed */ + const float soundspeed = gas_soundspeed_from_internal_energy(p->rho, p->u); - p->force.soundspeed = soundspeed; p->force.pressure = pressure; + p->force.soundspeed = soundspeed; } /** @@ -442,7 +432,12 @@ __attribute__((always_inline)) INLINE static void hydro_convert_quantities( /* Compute the pressure */ const float pressure = gas_pressure_from_internal_energy(p->rho, p->u); + + /* Compute the sound speed */ + const float soundspeed = gas_soundspeed_from_internal_energy(p->rho, p->u); + p->force.pressure = pressure; + p->force.soundspeed = soundspeed; } #endif /* SWIFT_MINIMAL_HYDRO_H */