Merge branch 'equation_of_state' into 'master'

Equation of state

This implements @tt's suggestion (#195) to allow for changes in the equation of state of the gas. The transformation between thermodynamic variables are now all located in one single file. Users can supply a different set of transformations if they want to simulate a different kind of gas. The SPH schemes themselves remain unchanged.

This is "simply" a harmless re-factoring of the code gathering common equations into a single file. It passes all the tests but you may want to cross-check things. Thanks!
Branch can be removed.

See merge request !215

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  \

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
 
 /**

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

src/equation_of_state.h

+/*******************************************************************************
+ * 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 */

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;
-}

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]);

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);
 }

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