diff --git a/src/Makefile.am b/src/Makefile.am
index 4aa7278ba5162ec106fe25f10529f766f9a703a8..7e2ceeb6067bb831620472b50ccd51fdfd50b494 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -56,6 +56,7 @@ AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
nobase_noinst_HEADERS = approx_math.h atomic.h cycle.h error.h inline.h kernel_hydro.h kernel_gravity.h \
kernel_long_gravity.h vector.h runner_doiact.h runner_doiact_grav.h runner_doiact_fft.h \
units.h intrinsics.h minmax.h kick.h timestep.h drift.h adiabatic_index.h io_properties.h \
+ equation_of_state.h \
gravity.h gravity_io.h \
gravity/Default/gravity.h gravity/Default/gravity_iact.h gravity/Default/gravity_io.h \
gravity/Default/gravity_debug.h gravity/Default/gravity_part.h \
diff --git a/src/adiabatic_index.h b/src/adiabatic_index.h
index ba7a1a64d5c97f1f6d276e5969782351762be4b1..50d24e881585abfb9c482b936780b5e1935df314 100644
--- a/src/adiabatic_index.h
+++ b/src/adiabatic_index.h
@@ -49,6 +49,10 @@
#define hydro_gamma_minus_one 1.f
#define hydro_one_over_gamma_minus_one 1.f
+#else
+
+#error "An adiabatic index needs to be chosen in const.h !"
+
#endif
/**
diff --git a/src/const.h b/src/const.h
index 736ba0155e03dc0f595ac06edfd0763b38302e9c..498b60b711e88e8054b5e0e73a18475d98277448 100644
--- a/src/const.h
+++ b/src/const.h
@@ -41,6 +41,10 @@
//#define HYDRO_GAMMA_4_3
//#define HYDRO_GAMMA_2_1
+/* Equation of state choice */
+#define EOS_IDEAL_GAS
+//#define EOS_ISOTHERMAL_GAS
+
/* Kernel function to use */
#define CUBIC_SPLINE_KERNEL
//#define QUARTIC_SPLINE_KERNEL
diff --git a/src/equation_of_state.h b/src/equation_of_state.h
new file mode 100644
index 0000000000000000000000000000000000000000..d9f1a46f240a8bdf965749a5513ca3448d08aeb9
--- /dev/null
+++ b/src/equation_of_state.h
@@ -0,0 +1,197 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk).
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_EQUATION_OF_STATE_H
+#define SWIFT_EQUATION_OF_STATE_H
+
+/* Config parameters. */
+#include "../config.h"
+
+/* Some standard headers. */
+#include <math.h>
+
+/* Local headers. */
+#include "adiabatic_index.h"
+#include "const.h"
+#include "debug.h"
+#include "inline.h"
+
+/* ------------------------------------------------------------------------- */
+#if defined(EOS_IDEAL_GAS)
+
+/**
+ * @brief Returns the internal energy given density and entropy
+ *
+ * Computes \f$u = \frac{S\rho^{\gamma-1} }{\gamma - 1}\f$.
+ *
+ * @param density The density \f$\rho\f$.
+ * @param entropy The entropy \f$S\f$.
+ */
+__attribute__((always_inline)) INLINE static float
+gas_internal_energy_from_entropy(float density, float entropy) {
+
+ return entropy * pow_gamma_minus_one(density) *
+ hydro_one_over_gamma_minus_one;
+}
+
+/**
+ * @brief Returns the pressure given density and entropy
+ *
+ * Computes \f$P = S\rho^\gamma\f$.
+ *
+ * @param density The density \f$\rho\f$.
+ * @param entropy The entropy \f$S\f$.
+ */
+__attribute__((always_inline)) INLINE static float gas_pressure_from_entropy(
+ float density, float entropy) {
+
+ return entropy * pow_gamma(density);
+}
+
+/**
+ * @brief Returns the sound speed given density and entropy
+ *
+ * Computes \f$c = \sqrt{\gamma S \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(hydro_gamma * pow_gamma_minus_one(density) * entropy);
+}
+
+/**
+ * @brief Returns the entropy given density and internal energy
+ *
+ * Computes \f$S = \frac{(\gamma - 1)u}{\rho^{\gamma-1}}\f$.
+ *
+ * @param density The density \f$\rho\f$
+ * @param u The internal energy \f$u\f$
+ */
+__attribute__((always_inline)) INLINE static float
+gas_entropy_from_internal_energy(float density, float u) {
+
+ return hydro_gamma_minus_one * u * pow_minus_gamma_minus_one(density);
+}
+
+/**
+ * @brief Returns the pressure given density and internal energy
+ *
+ * Computes \f$P = (\gamma - 1)u\rho\f$.
+ *
+ * @param density The density \f$\rho\f$
+ * @param u The internal energy \f$u\f$
+ */
+__attribute__((always_inline)) INLINE static float
+gas_pressure_from_internal_energy(float density, float u) {
+
+ return hydro_gamma_minus_one * u * density;
+}
+
+/**
+ * @brief Returns the sound speed given density and internal energy
+ *
+ * Computes \f$c = \sqrt{\gamma (\gamma - 1) u }\f$.
+ *
+ * @param density The density \f$\rho\f$
+ * @param u The internal energy \f$u\f$
+ */
+__attribute__((always_inline)) INLINE static float
+gas_soundspeed_from_internal_energy(float density, float u) {
+
+ return sqrtf(u * hydro_gamma * hydro_gamma_minus_one);
+}
+
+/* ------------------------------------------------------------------------- */
+#elif defined(EOS_ISOTHERMAL_GAS)
+
+/**
+ * @brief Returns the internal energy given density and entropy
+ *
+ * @param density The density
+ * @param entropy The entropy
+ */
+__attribute__((always_inline)) INLINE static float
+gas_internal_energy_from_entropy(float density, float entropy) {
+
+ error("Missing definition !");
+ return 0.f;
+}
+
+/**
+ * @brief Returns the pressure given density and entropy
+ *
+ * @param density The density
+ * @param entropy The entropy
+ */
+__attribute__((always_inline)) INLINE static float gas_pressure_from_entropy(
+ float density, float entropy) {
+
+ error("Missing definition !");
+ return 0.f;
+}
+
+/**
+ * @brief Returns the entropy given density and internal energy
+ *
+ * @param density The density
+ * @param u The internal energy
+ */
+__attribute__((always_inline)) INLINE static float
+gas_entropy_from_internal_energy(float density, float u) {
+
+ error("Missing definition !");
+ return 0.f;
+}
+
+/**
+ * @brief Returns the pressure given density and internal energy
+ *
+ * @param density The density
+ * @param u The internal energy
+ */
+__attribute__((always_inline)) INLINE static float
+gas_pressure_from_internal_energy(float density, float u) {
+
+ error("Missing definition !");
+ return 0.f;
+}
+
+/**
+ * @brief Returns the sound speed given density and internal energy
+ *
+ * @param density The density
+ * @param u The internal energy
+ */
+__attribute__((always_inline)) INLINE static float
+gas_soundspeed_from_internal_energy(float density, float u) {
+
+ error("Missing definition !");
+ return 0.f;
+}
+
+/* ------------------------------------------------------------------------- */
+#else
+
+#error "An Equation of state needs to be chosen in const.h !"
+
+#endif
+
+#endif /* SWIFT_EQUATION_OF_STATE_H */
diff --git a/src/hydro/Default/hydro.h b/src/hydro/Default/hydro.h
index 45b21531e610391fed6693daa07b8907bcdd8a5b..7fe78f75fc00569ee8d34ea30a822bac394531e2 100644
--- a/src/hydro/Default/hydro.h
+++ b/src/hydro/Default/hydro.h
@@ -19,9 +19,89 @@
#include "adiabatic_index.h"
#include "approx_math.h"
+#include "equation_of_state.h"
#include <float.h>
+/**
+ * @brief Returns the internal energy of a particle
+ *
+ * @param p The particle of interest
+ * @param dt Time since the last kick
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
+ const struct part *restrict p, float dt) {
+
+ return p->u;
+}
+
+/**
+ * @brief Returns the pressure of a particle
+ *
+ * @param p The particle of interest
+ * @param dt Time since the last kick
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_pressure(
+ const struct part *restrict p, float dt) {
+
+ return gas_pressure_from_internal_energy(p->rho, p->u);
+}
+
+/**
+ * @brief Returns the entropy of a particle
+ *
+ * @param p The particle of interest
+ * @param dt Time since the last kick
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_entropy(
+ const struct part *restrict p, float dt) {
+
+ return gas_entropy_from_internal_energy(p->rho, p->u);
+}
+
+/**
+ * @brief Returns the sound speed of a particle
+ *
+ * @param p The particle of interest
+ * @param dt Time since the last kick
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_soundspeed(
+ const struct part *restrict p, float dt) {
+
+ return p->force.soundspeed;
+}
+
+/**
+ * @brief Modifies the thermal state of a particle to the imposed internal
+ * energy
+ *
+ * This overrides the current state of the particle but does *not* change its
+ * time-derivatives
+ *
+ * @param p The particle
+ * @param u The new internal energy
+ */
+__attribute__((always_inline)) INLINE static void hydro_set_internal_energy(
+ struct part *restrict p, float u) {
+
+ p->u = u;
+}
+
+/**
+ * @brief Modifies the thermal state of a particle to the imposed entropy
+ *
+ * This overrides the current state of the particle but does *not* change its
+ * time-derivatives
+ *
+ * @param p The particle
+ * @param S The new entropy
+ */
+__attribute__((always_inline)) INLINE static void hydro_set_entropy(
+ struct part *restrict p, float S) {
+
+ p->u = gas_internal_energy_from_entropy(p->rho, S);
+}
+
/**
* @brief Computes the hydro time-step of a given particle
*
@@ -253,51 +333,3 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
*/
__attribute__((always_inline)) INLINE static void hydro_convert_quantities(
struct part *restrict p) {}
-
-/**
- * @brief Returns the internal energy of a particle
- *
- * @param p The particle of interest
- * @param dt Time since the last kick
- */
-__attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
- const struct part *restrict p, float dt) {
-
- return p->u;
-}
-
-/**
- * @brief Returns the pressure of a particle
- *
- * @param p The particle of interest
- * @param dt Time since the last kick
- */
-__attribute__((always_inline)) INLINE static float hydro_get_pressure(
- const struct part *restrict p, float dt) {
-
- return p->force.P_over_rho2 * p->rho * p->rho / p->rho_dh;
-}
-
-/**
- * @brief Returns the entropy of a particle
- *
- * @param p The particle of interest
- * @param dt Time since the last kick
- */
-__attribute__((always_inline)) INLINE static float hydro_get_entropy(
- const struct part *restrict p, float dt) {
-
- return hydro_gamma_minus_one * p->u * pow_minus_gamma_minus_one(p->rho);
-}
-
-/**
- * @brief Returns the sound speed of a particle
- *
- * @param p The particle of interest
- * @param dt Time since the last kick
- */
-__attribute__((always_inline)) INLINE static float hydro_get_soundspeed(
- const struct part *restrict p, float dt) {
-
- return p->force.soundspeed;
-}
diff --git a/src/hydro/Default/hydro_iact.h b/src/hydro/Default/hydro_iact.h
index 5db3604e1717ba2c14d50c531c0537409501226d..141ab8f884da7f120f2b636933904f0d795a5d8b 100644
--- a/src/hydro/Default/hydro_iact.h
+++ b/src/hydro/Default/hydro_iact.h
@@ -468,7 +468,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_vec_force(
vector pia[3], pja[3];
vector piu_dt, pju_dt;
vector pih_dt, pjh_dt;
- vector ci, cj, v_sig, vi_sig, vj_sig;
+ vector ci, cj, v_sig;
vector omega_ij, Pi_ij, balsara;
vector pialpha, pjalpha, alpha_ij, v_sig_u, tc;
int j, k;
@@ -503,12 +503,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_vec_force(
pj[2]->force.soundspeed, pj[3]->force.soundspeed,
pj[4]->force.soundspeed, pj[5]->force.soundspeed,
pj[6]->force.soundspeed, pj[7]->force.soundspeed);
- vi_sig.v = vec_set(pi[0]->force.v_sig, pi[1]->force.v_sig, pi[2]->force.v_sig,
- pi[3]->force.v_sig, pi[4]->force.v_sig, pi[5]->force.v_sig,
- pi[6]->force.v_sig, pi[7]->force.v_sig);
- vj_sig.v = vec_set(pj[0]->force.v_sig, pj[1]->force.v_sig, pj[2]->force.v_sig,
- pj[3]->force.v_sig, pj[4]->force.v_sig, pj[5]->force.v_sig,
- pj[6]->force.v_sig, pj[7]->force.v_sig);
for (k = 0; k < 3; k++) {
vi[k].v = vec_set(pi[0]->v[k], pi[1]->v[k], pi[2]->v[k], pi[3]->v[k],
pi[4]->v[k], pi[5]->v[k], pi[6]->v[k], pi[7]->v[k]);
@@ -544,10 +538,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_vec_force(
pi[2]->force.soundspeed, pi[3]->force.soundspeed);
cj.v = vec_set(pj[0]->force.soundspeed, pj[1]->force.soundspeed,
pj[2]->force.soundspeed, pj[3]->force.soundspeed);
- vi_sig.v = vec_set(pi[0]->force.v_sig, pi[1]->force.v_sig, pi[2]->force.v_sig,
- pi[3]->force.v_sig);
- vj_sig.v = vec_set(pj[0]->force.v_sig, pj[1]->force.v_sig, pj[2]->force.v_sig,
- pj[3]->force.v_sig);
for (k = 0; k < 3; k++) {
vi[k].v = vec_set(pi[0]->v[k], pi[1]->v[k], pi[2]->v[k], pi[3]->v[k]);
vj[k].v = vec_set(pj[0]->v[k], pj[1]->v[k], pj[2]->v[k], pj[3]->v[k]);
@@ -773,7 +763,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_vec_force(
vector pia[3];
vector piu_dt;
vector pih_dt;
- vector ci, cj, v_sig, vi_sig, vj_sig;
+ vector ci, cj, v_sig;
vector omega_ij, Pi_ij, balsara;
vector pialpha, pjalpha, alpha_ij, v_sig_u, tc;
int j, k;
@@ -806,12 +796,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_vec_force(
pj[2]->force.soundspeed, pj[3]->force.soundspeed,
pj[4]->force.soundspeed, pj[5]->force.soundspeed,
pj[6]->force.soundspeed, pj[7]->force.soundspeed);
- vi_sig.v = vec_set(pi[0]->force.v_sig, pi[1]->force.v_sig, pi[2]->force.v_sig,
- pi[3]->force.v_sig, pi[4]->force.v_sig, pi[5]->force.v_sig,
- pi[6]->force.v_sig, pi[7]->force.v_sig);
- vj_sig.v = vec_set(pj[0]->force.v_sig, pj[1]->force.v_sig, pj[2]->force.v_sig,
- pj[3]->force.v_sig, pj[4]->force.v_sig, pj[5]->force.v_sig,
- pj[6]->force.v_sig, pj[7]->force.v_sig);
for (k = 0; k < 3; k++) {
vi[k].v = vec_set(pi[0]->v[k], pi[1]->v[k], pi[2]->v[k], pi[3]->v[k],
pi[4]->v[k], pi[5]->v[k], pi[6]->v[k], pi[7]->v[k]);
@@ -846,10 +830,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_vec_force(
pi[2]->force.soundspeed, pi[3]->force.soundspeed);
cj.v = vec_set(pj[0]->force.soundspeed, pj[1]->force.soundspeed,
pj[2]->force.soundspeed, pj[3]->force.soundspeed);
- vi_sig.v = vec_set(pi[0]->force.v_sig, pi[1]->force.v_sig, pi[2]->force.v_sig,
- pi[3]->force.v_sig);
- vj_sig.v = vec_set(pj[0]->force.v_sig, pj[1]->force.v_sig, pj[2]->force.v_sig,
- pj[3]->force.v_sig);
for (k = 0; k < 3; k++) {
vi[k].v = vec_set(pi[0]->v[k], pi[1]->v[k], pi[2]->v[k], pi[3]->v[k]);
vj[k].v = vec_set(pj[0]->v[k], pj[1]->v[k], pj[2]->v[k], pj[3]->v[k]);
diff --git a/src/hydro/Gadget2/hydro.h b/src/hydro/Gadget2/hydro.h
index 5efce17d78ddc22911ca297a1ef9024d31971ba9..a1448d8686a30de88e24629558f01661daca223a 100644
--- a/src/hydro/Gadget2/hydro.h
+++ b/src/hydro/Gadget2/hydro.h
@@ -18,6 +18,90 @@
******************************************************************************/
#include "adiabatic_index.h"
+#include "equation_of_state.h"
+
+/**
+ * @brief Returns the internal energy of a particle
+ *
+ * @param p The particle of interest
+ * @param dt Time since the last kick
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
+ const struct part *restrict p, float dt) {
+
+ const float entropy = p->entropy + p->entropy_dt * dt;
+
+ return gas_internal_energy_from_entropy(p->rho, entropy);
+}
+
+/**
+ * @brief Returns the pressure of a particle
+ *
+ * @param p The particle of interest
+ * @param dt Time since the last kick
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_pressure(
+ const struct part *restrict p, float dt) {
+
+ const float entropy = p->entropy + p->entropy_dt * dt;
+
+ return gas_pressure_from_entropy(p->rho, entropy);
+}
+
+/**
+ * @brief Returns the entropy of a particle
+ *
+ * @param p The particle of interest
+ * @param dt Time since the last kick
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_entropy(
+ const struct part *restrict p, float dt) {
+
+ return p->entropy + p->entropy_dt * dt;
+}
+
+/**
+ * @brief Returns the sound speed of a particle
+ *
+ * @param p The particle of interest
+ * @param dt Time since the last kick
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_soundspeed(
+ const struct part *restrict p, float dt) {
+
+ return p->force.soundspeed;
+}
+
+/**
+ * @brief Modifies the thermal state of a particle to the imposed internal
+ * energy
+ *
+ * This overrides the current state of the particle but does *not* change its
+ * time-derivatives
+ *
+ * @param p The particle
+ * @param u The new internal energy
+ */
+__attribute__((always_inline)) INLINE static void hydro_set_internal_energy(
+ struct part *restrict p, float u) {
+
+ p->entropy = gas_entropy_from_internal_energy(p->rho, u);
+}
+
+/**
+ * @brief Modifies the thermal state of a particle to the imposed entropy
+ *
+ * This overrides the current state of the particle but does *not* change its
+ * time-derivatives
+ *
+ * @param p The particle
+ * @param S The new entropy
+ */
+__attribute__((always_inline)) INLINE static void hydro_set_entropy(
+ struct part *restrict p, float S) {
+
+ p->entropy = S;
+}
/**
* @brief Computes the hydro time-step of a given particle
@@ -140,8 +224,7 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
/* Compute the pressure */
const float half_dt = (ti_current - (p->ti_begin + p->ti_end) / 2) * timeBase;
- const float pressure =
- (p->entropy + p->entropy_dt * half_dt) * pow_gamma(p->rho);
+ const float pressure = hydro_get_pressure(p, half_dt);
const float irho = 1.f / p->rho;
@@ -200,8 +283,7 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
/* Drift the pressure */
const float dt_entr = (t1 - (p->ti_begin + p->ti_end) / 2) * timeBase;
- const float pressure =
- (p->entropy + p->entropy_dt * dt_entr) * pow_gamma(p->rho);
+ const float pressure = hydro_get_pressure(p, dt_entr);
const float irho = 1.f / p->rho;
@@ -265,56 +347,6 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
__attribute__((always_inline)) INLINE static void hydro_convert_quantities(
struct part *restrict p) {
- p->entropy =
- hydro_gamma_minus_one * p->entropy * pow_minus_gamma_minus_one(p->rho);
-}
-
-/**
- * @brief Returns the internal energy of a particle
- *
- * @param p The particle of interest
- * @param dt Time since the last kick
- */
-__attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
- const struct part *restrict p, float dt) {
-
- const float entropy = p->entropy + p->entropy_dt * dt;
-
- return entropy * pow_gamma_minus_one(p->rho) * hydro_one_over_gamma_minus_one;
-}
-
-/**
- * @brief Returns the pressure of a particle
- *
- * @param p The particle of interest
- * @param dt Time since the last kick
- */
-__attribute__((always_inline)) INLINE static float hydro_get_pressure(
- const struct part *restrict p, float dt) {
-
- return p->force.P_over_rho2 * p->rho * p->rho / p->rho_dh;
-}
-
-/**
- * @brief Returns the entropy of a particle
- *
- * @param p The particle of interest
- * @param dt Time since the last kick
- */
-__attribute__((always_inline)) INLINE static float hydro_get_entropy(
- const struct part *restrict p, float dt) {
-
- return p->entropy + p->entropy_dt * dt;
-}
-
-/**
- * @brief Returns the sound speed of a particle
- *
- * @param p The particle of interest
- * @param dt Time since the last kick
- */
-__attribute__((always_inline)) INLINE static float hydro_get_soundspeed(
- const struct part *restrict p, float dt) {
-
- return p->force.soundspeed;
+ /* We read u in the entropy field. We now get S from u */
+ p->entropy = gas_entropy_from_internal_energy(p->rho, p->entropy);
}
diff --git a/src/hydro/Minimal/hydro.h b/src/hydro/Minimal/hydro.h
index f1f284d5d2b56af9d20b24c91eaf4a481f94851f..444920ab2bec940af9c10ea0a6ba471b2ab79762 100644
--- a/src/hydro/Minimal/hydro.h
+++ b/src/hydro/Minimal/hydro.h
@@ -19,6 +19,94 @@
#include "adiabatic_index.h"
#include "approx_math.h"
+#include "equation_of_state.h"
+
+/**
+ * @brief Returns the 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.
+ *
+ * @param p The particle of interest
+ * @param dt Time since the last kick
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
+ const struct part *restrict p, float dt) {
+
+ return p->u;
+}
+
+/**
+ * @brief Returns the pressure of a particle
+ *
+ * @param p The particle of interest
+ * @param dt Time since the last kick
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_pressure(
+ const struct part *restrict p, float dt) {
+
+ return gas_pressure_from_internal_energy(p->rho, p->u);
+}
+
+/**
+ * @brief Returns the entropy of a particle
+ *
+ * For implementations where the main thermodynamic variable
+ * is not entropy, this function computes the entropy from
+ * the thermodynamic variable.
+ *
+ * @param p The particle of interest
+ * @param dt Time since the last kick
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_entropy(
+ const struct part *restrict p, float dt) {
+
+ return gas_entropy_from_internal_energy(p->rho, p->u);
+}
+
+/**
+ * @brief Returns the sound speed of a particle
+ *
+ * @param p The particle of interest
+ * @param dt Time since the last kick
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_soundspeed(
+ const struct part *restrict p, float dt) {
+
+ return gas_soundspeed_from_internal_energy(p->rho, p->u);
+}
+
+/**
+ * @brief Modifies the thermal state of a particle to the imposed internal
+ * energy
+ *
+ * This overrides the current state of the particle but does *not* change its
+ * time-derivatives
+ *
+ * @param p The particle
+ * @param u The new internal energy
+ */
+__attribute__((always_inline)) INLINE static void hydro_set_internal_energy(
+ struct part *restrict p, float u) {
+
+ p->u = u;
+}
+
+/**
+ * @brief Modifies the thermal state of a particle to the imposed entropy
+ *
+ * This overrides the current state of the particle but does *not* change its
+ * time-derivatives
+ *
+ * @param p The particle
+ * @param S The new entropy
+ */
+__attribute__((always_inline)) INLINE static void hydro_set_entropy(
+ struct part *restrict p, float S) {
+
+ p->u = gas_internal_energy_from_entropy(p->rho, S);
+}
/**
* @brief Computes the hydro time-step of a given particle
@@ -230,55 +318,3 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
*/
__attribute__((always_inline)) INLINE static void hydro_convert_quantities(
struct part *restrict p) {}
-
-/**
- * @brief Returns the 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.
- *
- * @param p The particle of interest
- * @param dt Time since the last kick
- */
-__attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
- const struct part *restrict p, float dt) {
-
- return p->u;
-}
-
-/**
- * @brief Returns the pressure of a particle
- *
- * @param p The particle of interest
- * @param dt Time since the last kick
- */
-__attribute__((always_inline)) INLINE static float hydro_get_pressure(
- const struct part *restrict p, float dt) {
-
- return p->u * p->rho * hydro_gamma_minus_one;
-}
-
-/**
- * @brief Returns the entropy of a particle
- *
- * @param p The particle of interest
- * @param dt Time since the last kick
- */
-__attribute__((always_inline)) INLINE static float hydro_get_entropy(
- const struct part *restrict p, float dt) {
-
- return hydro_gamma_minus_one * p->u * pow_minus_gamma_minus_one(p->rho);
-}
-
-/**
- * @brief Returns the sound speed of a particle
- *
- * @param p The particle of interest
- * @param dt Time since the last kick
- */
-__attribute__((always_inline)) INLINE static float hydro_get_soundspeed(
- const struct part *restrict p, float dt) {
-
- return sqrt(p->u * hydro_gamma * hydro_gamma_minus_one);
-}
diff --git a/src/timestep.h b/src/timestep.h
index 99747484b2fad2d1dfadad749232f77848e56f35..569120cf9cf989b633da35529f7693c8f62a1910 100644
--- a/src/timestep.h
+++ b/src/timestep.h
@@ -108,7 +108,7 @@ __attribute__((always_inline)) INLINE static int get_part_timestep(
e->external_potential, e->physical_constants, p->gpart);
/* const float new_dt_self = */
/* gravity_compute_timestep_self(e->physical_constants, p->gpart); */
- const float new_dt_self = FLT_MAX; // MATTHIEU
+ const float new_dt_self = FLT_MAX; // MATTHIEU
new_dt_grav = fminf(new_dt_external, new_dt_self);
}