From e3ae510d12d0af94f0738845f31dbbce9c593f91 Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <matthieu.schaller@durham.ac.uk>
Date: Thu, 9 Nov 2017 10:07:13 +0000
Subject: [PATCH] Use the equation of state functions in a few more places in
 the gizmo code to be more robust to changes.

---
 src/hydro/Gizmo/hydro.h    | 36 +++++++++++++++---------------------
 src/hydro/Gizmo/hydro_io.h | 18 +++---------------
 2 files changed, 18 insertions(+), 36 deletions(-)

diff --git a/src/hydro/Gizmo/hydro.h b/src/hydro/Gizmo/hydro.h
index ccbeff6743..6273f5009a 100644
--- a/src/hydro/Gizmo/hydro.h
+++ b/src/hydro/Gizmo/hydro.h
@@ -24,6 +24,7 @@
 #include "approx_math.h"
 #include "equation_of_state.h"
 #include "hydro_gradients.h"
+#include "hydro_properties.h"
 #include "hydro_space.h"
 #include "hydro_unphysical.h"
 #include "hydro_velocities.h"
@@ -130,8 +131,9 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
 /* remember that we store the total thermal energy, not the specific thermal
    energy (as in Gadget) */
 #if defined(EOS_ISOTHERMAL_GAS)
-  /* this overwrites the internal energy from the initial condition file */
-  p->conserved.energy = mass * eos.isothermal_internal_energy;
+  /* this overwrites the internal energy from the initial condition file
+   * Note that we call the EoS function just to get the constant u here. */
+  p->conserved.energy = mass * gas_internal_energy_from_entropy(0.f, 0.f);
 #else
   p->conserved.energy *= mass;
 #endif
@@ -608,7 +610,9 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
   p->conserved.momentum[1] += p->conserved.flux.momentum[1];
   p->conserved.momentum[2] += p->conserved.flux.momentum[2];
 #if defined(EOS_ISOTHERMAL_GAS)
-  p->conserved.energy = p->conserved.mass * eos.isothermal_internal_energy;
+  /* We use the EoS equation in a sneaky way here just to get the constant u */
+  p->conserved.energy =
+      p->conserved.mass * gas_internal_energy_from_entropy(0.f, 0.f);
 #else
   p->conserved.energy += p->conserved.flux.energy;
 #endif
@@ -707,15 +711,11 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
 __attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
     const struct part* restrict p) {
 
-  if (p->primitives.rho > 0.) {
-#ifdef EOS_ISOTHERMAL_GAS
-    return eos.isothermal_internal_energy;
-#else
-    return p->primitives.P / hydro_gamma_minus_one / p->primitives.rho;
-#endif
-  } else {
+  if (p->primitives.rho > 0.)
+    return gas_internal_energy_from_pressure(p->primitives.rho,
+                                             p->primitives.P);
+  else
     return 0.;
-  }
 }
 
 /**
@@ -727,7 +727,7 @@ __attribute__((always_inline)) INLINE static float hydro_get_entropy(
     const struct part* restrict p) {
 
   if (p->primitives.rho > 0.) {
-    return p->primitives.P / pow_gamma(p->primitives.rho);
+    return gas_entropy_from_pressure(p->primitives.rho, p->primitives.P);
   } else {
     return 0.;
   }
@@ -741,16 +741,10 @@ __attribute__((always_inline)) INLINE static float hydro_get_entropy(
 __attribute__((always_inline)) INLINE static float hydro_get_soundspeed(
     const struct part* restrict p) {
 
-  if (p->primitives.rho > 0.) {
-#ifdef EOS_ISOTHERMAL_GAS
-    return sqrtf(eos.isothermal_internal_energy * hydro_gamma *
-                 hydro_gamma_minus_one);
-#else
-    return sqrtf(hydro_gamma * p->primitives.P / p->primitives.rho);
-#endif
-  } else {
+  if (p->primitives.rho > 0.)
+    return gas_soundspeed_from_pressure(p->primitives.rho, p->primitives.P);
+  else
     return 0.;
-  }
 }
 
 /**
diff --git a/src/hydro/Gizmo/hydro_io.h b/src/hydro/Gizmo/hydro_io.h
index 09d667d8ff..ea0a083f37 100644
--- a/src/hydro/Gizmo/hydro_io.h
+++ b/src/hydro/Gizmo/hydro_io.h
@@ -18,7 +18,7 @@
  ******************************************************************************/
 
 #include "adiabatic_index.h"
-#include "equation_of_state.h"
+#include "hydro.h"
 #include "hydro_flux_limiters.h"
 #include "hydro_gradients.h"
 #include "hydro_slope_limiters.h"
@@ -72,15 +72,7 @@ void hydro_read_particles(struct part* parts, struct io_props* list,
  * @return Internal energy of the particle
  */
 float convert_u(struct engine* e, struct part* p) {
-  if (p->primitives.rho > 0.) {
-#ifdef EOS_ISOTHERMAL_GAS
-    return eos.isothermal_internal_energy;
-#else
-    return p->primitives.P / hydro_gamma_minus_one / p->primitives.rho;
-#endif
-  } else {
-    return 0.;
-  }
+  return hydro_get_internal_energy(p);
 }
 
 /**
@@ -91,11 +83,7 @@ float convert_u(struct engine* e, struct part* p) {
  * @return Entropic function of the particle
  */
 float convert_A(struct engine* e, struct part* p) {
-  if (p->primitives.rho > 0.) {
-    return p->primitives.P / pow_gamma(p->primitives.rho);
-  } else {
-    return 0.;
-  }
+  return hydro_get_entropy(p);
 }
 
 /**
-- 
GitLab