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

Cleaner way to handle the iosthermal equation of state constant

parent fb1c657f
......@@ -329,7 +329,6 @@ int main(int argc, char *argv[]) {
/* Initialise the hydro properties */
struct hydro_props hydro_properties;
hydro_props_init(&hydro_properties, params);
if (with_hydro && myrank == 0) eos_print();
/* Initialise the external potential properties */
struct external_potential potential;
......
......@@ -36,6 +36,10 @@
/* Time integration constants. */
#define const_max_u_change 0.1f
/* Thermal energy per unit mass to use as a constant when using an isothermal
* EoS */
#define const_isothermal_internal_energy 20.2615290634f
/* Dimensionality of the problem */
#define HYDRO_DIMENSION_3D
//#define HYDRO_DIMENSION_2D
......@@ -47,10 +51,8 @@
//#define HYDRO_GAMMA_2_1
/* Equation of state choice */
#define EOS_IDEAL_GAS
/* if EOS_ISOTHERMAL_GAS is defined, keep thermal energy per unit mass equal to
* EOS_ISOTHERMAL_GAS in programme units */
//#define EOS_ISOTHERMAL_GAS (20.2615290634)
//#define EOS_IDEAL_GAS
#define EOS_ISOTHERMAL_GAS
/* Kernel function to use */
#define CUBIC_SPLINE_KERNEL
......
......@@ -40,9 +40,6 @@
/* ------------------------------------------------------------------------- */
#if defined(EOS_IDEAL_GAS)
#if defined(EOS_ISOTHERMAL_GAS)
#error : cannot have ideal and isothermal gas at the same time!
#endif
/**
* @brief Returns the internal energy given density and entropy
......@@ -128,18 +125,7 @@ gas_soundspeed_from_internal_energy(float density, float u) {
return sqrtf(u * hydro_gamma * hydro_gamma_minus_one);
}
__attribute__((always_inline)) INLINE static void eos_print() {
message(" assuming an equation of state of an ideal gas, with gamma = %e",
hydro_gamma);
}
__attribute__((always_inline)) INLINE static float hydro_update_entropy(
float density, float entropy) {
return entropy;
}
__attribute__((always_inline)) INLINE static float hydro_update_entropy_rate(
float density, float entropy) {
return hydro_gamma_minus_one * pow_minus_gamma_minus_one(density);
}
/* ------------------------------------------------------------------------- */
#elif defined(EOS_ISOTHERMAL_GAS)
......@@ -152,7 +138,7 @@ __attribute__((always_inline)) INLINE static float hydro_update_entropy_rate(
__attribute__((always_inline)) INLINE static float
gas_internal_energy_from_entropy(float density, float entropy) {
return EOS_ISOTHERMAL_GAS;
return const_isothermal_internal_energy;
}
/**
* @brief Returns the pressure given density and entropy
......@@ -163,7 +149,22 @@ gas_internal_energy_from_entropy(float density, float entropy) {
__attribute__((always_inline)) INLINE static float gas_pressure_from_entropy(
float density, float entropy) {
return hydro_gamma_minus_one * EOS_ISOTHERMAL_GAS * density;
return hydro_gamma_minus_one * const_isothermal_internal_energy * density;
}
/**
* @brief Returns the sound speed given density and entropy
*
* Computes \f$c = \sqrt{u \gamma \rho^{\gamma-1}}\f$.
*
* @param density The density \f$\rho\f$.
* @param entropy The entropy \f$S\f$.
*/
__attribute__((always_inline)) INLINE static float gas_soundspeed_from_entropy(
float density, float entropy) {
return sqrtf(const_isothermal_internal_energy * hydro_gamma *
hydro_gamma_minus_one);
}
/**
......@@ -175,7 +176,7 @@ __attribute__((always_inline)) INLINE static float gas_pressure_from_entropy(
__attribute__((always_inline)) INLINE static float
gas_entropy_from_internal_energy(float density, float u) {
return hydro_gamma_minus_one * EOS_ISOTHERMAL_GAS /
return hydro_gamma_minus_one * const_isothermal_internal_energy /
pow(density, hydro_gamma - 1);
}
......@@ -188,8 +189,7 @@ gas_entropy_from_internal_energy(float density, float u) {
__attribute__((always_inline)) INLINE static float
gas_pressure_from_internal_energy(float density, float u) {
error("Missing definition !");
return 0.f;
return hydro_gamma_minus_one * const_isothermal_internal_energy * density;
}
/**
......@@ -201,24 +201,10 @@ gas_pressure_from_internal_energy(float density, float u) {
__attribute__((always_inline)) INLINE static float
gas_soundspeed_from_internal_energy(float density, float u) {
error("Missing definition !");
return 0.f;
}
__attribute__((always_inline)) INLINE static void eos_print() {
message(
" assuming an equation of state for an isothermal gas with utherm= %e "
"and gamma = %e",
EOS_ISOTHERMAL_GAS, hydro_gamma);
}
__attribute__((always_inline)) INLINE static float hydro_update_entropy(
float density, float entropy) {
float thermal_energy = EOS_ISOTHERMAL_GAS;
return gas_entropy_from_internal_energy(density, thermal_energy);
}
__attribute__((always_inline)) INLINE static float hydro_update_entropy_rate(
float density, float entropy) {
return 0;
return sqrtf(const_isothermal_internal_energy * hydro_gamma *
hydro_gamma_minus_one);
}
/* ------------------------------------------------------------------------- */
#else
......
......@@ -56,6 +56,15 @@ void hydro_props_init(struct hydro_props *p,
void hydro_props_print(const struct hydro_props *p) {
#if defined(EOS_IDEAL_GAS)
message("Equation of state: Ideal gas.");
#elif defined(EOS_ISOTHERMAL_GAS)
message(
"Equation of state: Isothermal with internal energy "
"per unit mass set to %f.",
const_isothermal_internal_energy);
#endif
message("Adiabatic index gamma: %f.", hydro_gamma);
message("Hydrodynamic scheme: %s in %dD.", SPH_IMPLEMENTATION,
......
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