Define thermodynamical relations for the current EoS and use them for the Gadget scheme

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 IDEAL_GAS
+//#define 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(IDEAL_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) {
+
+  return entropy * pow_gamma_minus_one(density) *
+         hydro_one_over_gamma_minus_one;
+}
+
+/**
+ * @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) {
+
+  return entropy * pow_gamma(density);
+}
+
+/**
+ * @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) {
+
+  return hydro_gamma_minus_one * u * pow_minus_gamma_minus_one(density);
+}
+
+/**
+ * @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) {
+
+  return hydro_gamma_minus_one * u * density;
+}
+
+
+#elif defined(ISOTHERMAL_GAS)
+
+#error "Missing definitions !"
+
+#else
+
+#error "An Equation of state needs to be chosen in const.h !"
+
+#endif
+
+#endif /* SWIFT_EQUATION_OF_STATE_H */

src/hydro/Gadget2/hydro.h

@@ -18,6 +18,59 @@
  ******************************************************************************/
 
 #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 Computes the hydro time-step of a given particle
@@ -140,8 +193,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 +252,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 +316,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/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);
   }