Do the same for pow(x, gamma-1)

src/adiabatic_index.h

@@ -35,16 +35,19 @@
 
 #define hydro_gamma 1.66666666666666667f
 #define hydro_gamma_minus_one 0.66666666666666667f
+#define hydro_one_over_gamma_minus_one 1.5f
 
 #elif defined(HYDRO_GAMMA_4_3)
 
 #define hydro_gamma 1.33333333333333333f
 #define hydro_gamma_minus_one 0.33333333333333333f
+#define hydro_one_over_gamma_minus_one 3.f
 
 #elif defined(HYDRO_GAMMA_2_1)
 
 #define hydro_gamma 2.f
 #define hydro_gamma_minus_one 1.f
+#define hydro_one_over_gamma_minus_one 1.f
 
 #endif
 
@@ -72,8 +75,39 @@ __attribute__((always_inline)) INLINE static float pow_gamma(float x) {
   return x * x;
 
 #else
+
+  error("The adiabatic index is not defined !");
+  return 0.f;
+
+#endif
+}
+
+/**
+ * @brief Returns the argument to the power given by the adiabatic index minus
+ * one
+ *
+ * Computes $x^(\gamma - 1)$.
+ */
+__attribute__((always_inline)) INLINE static float pow_gamma_minus_one(
+    float x) {
+
+#if defined(HYDRO_GAMMA_5_3)
+
+  return cbrtf(x * x);
+
+#elif defined(HYDRO_GAMMA_4_3)
+
+  return cbrtf(x);
+
+#elif defined(HYDRO_GAMMA_2_1)
+
+  return x;
+
+#else
+
   error("The adiabatic index is not defined !");
   return 0.f;
+
 #endif
 }
 
@@ -99,8 +133,10 @@ __attribute__((always_inline)) INLINE static float pow_minus_gamma_minus_one(
   return 1.f / x;
 
 #else
+
   error("The adiabatic index is not defined !");
   return 0.f;
+
 #endif
 }
 

src/hydro/Default/hydro_iact.h

@@ -914,7 +914,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_vec_force(
   Pi_ij.v *= (wi_dr.v + wj_dr.v);
 
   /* Thermal conductivity */
-  v_sig_u.v = vec_sqrt(vec_set1(2.f * (const_hydro_gamma - 1.f)) *
+  v_sig_u.v = vec_sqrt(vec_set1(2.f * hydro_gamma_minus_one) *
                        vec_fabs(pirho.v * piu.v - pjrho.v * pju.v) /
                        (pirho.v + pjrho.v));
   tc.v = vec_set1(const_conductivity_alpha) * v_sig_u.v / (pirho.v + pjrho.v);

src/hydro/Gadget2/hydro.h

@@ -247,7 +247,5 @@ __attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
 
   const float entropy = p->entropy + p->entropy_dt * dt;
 
-  return entropy * powf(p->rho, hydro_gamma - 1.f) *
-         (1.f / (hydro_gamma - 1.f));
-  return 1;
+  return entropy * pow_gamma_minus_one(p->rho) * hydro_one_over_gamma_minus_one;
 }

src/hydro/Minimal/hydro.h

@@ -17,8 +17,8 @@
  *
  ******************************************************************************/
 
-#include "approx_math.h"
 #include "adiabatic_index.h"
+#include "approx_math.h"
 
 /**
  * @brief Computes the hydro time-step of a given particle