Added extra adiabatic index power functions for use in Riemann solvers.

src/adiabatic_index.h

 /*******************************************************************************
  * This file is part of SWIFT.
  * Copyright (c) 2016   Matthieu Schaller (matthieu.schaller@durham.ac.uk).
+ *                      Bert Vandenbroucke (bert.vandenbroucke@gmail.com).
  *
  * 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
@@ -167,4 +168,105 @@ __attribute__((always_inline)) INLINE static float pow_minus_gamma_minus_one(
 #endif
 }
 
+/**
+ * @brief Returns one over the argument to the power given by the adiabatic
+ * index
+ *
+ * Computes \f$x^{-\gamma}\f$.
+ *
+ * @param x Argument
+ * @return One over the argument to the power given by the adiabatic index
+ */
+__attribute__((always_inline)) INLINE static float pow_minus_gamma(float x) {
+
+#if defined(HYDRO_GAMMA_5_3)
+
+  const float cbrt_inv = 1.f / cbrtf(x);       /* x^(-1/3) */
+  const float cbrt_inv2 = cbrt_inv * cbrt_inv; /* x^(-2/3) */
+  return cbrt_inv * cbrt_inv2 * cbrt_inv2;     /* x^(-5/3) */
+
+#elif defined(HYDRO_GAMMA_4_3)
+
+  const float cbrt_inv = 1.f / cbrtf(x);       /* x^(-1/3) */
+  const float cbrt_inv2 = cbrt_inv * cbrt_inv; /* x^(-2/3) */
+  return cbrt_inv2 * cbrt_inv2;                /* x^(-4/3) */
+
+#elif defined(HYDRO_GAMMA_2_1)
+
+  const float inv = 1.f / x;
+  return inv * inv;
+
+#else
+
+  error("The adiabatic index is not defined !");
+  return 0.f;
+
+#endif
+}
+
+/**
+ * @brief Return the argument to the power given by two divided by the adiabatic
+ * index minus one
+ *
+ * Computes \f$x^{\frac{2}{\gamma - 1}}\f$.
+ *
+ * @param x Argument
+ * @return Argument to the power two divided by the adiabatic index minus one
+ */
+__attribute__((always_inline)) INLINE static float pow_two_over_gamma_minus_one(
+    float x) {
+
+#if defined(HYDRO_GAMMA_5_3)
+
+  return x * x * x; /* x^3 */
+
+#elif defined(HYDRO_GAMMA_4_3)
+
+  return x * x * x * x * x * x; /* x^6 */
+
+#elif defined(HYDRO_GAMMA_2_1)
+
+  return x * x; /* x^2 */
+
+#else
+
+  error("The adiabatic index is not defined !");
+  return 0.f;
+
+#endif
+}
+
+/**
+ * @brief Return the argument to the power given by two times the adiabatic
+ * index divided by the adiabatic index minus one
+ *
+ * Computes \f$x^{\frac{2\gamma}{\gamma - 1}}\f$.
+ *
+ * @param x Argument
+ * @return Argument to the power two times the adiabatic index divided by the
+ * adiabatic index minus one
+ */
+__attribute__((always_inline)) INLINE static float
+pow_two_gamma_over_gamma_minus_one(float x) {
+
+#if defined(HYDRO_GAMMA_5_3)
+
+  return x * x * x * x * x; /* x^5 */
+
+#elif defined(HYDRO_GAMMA_4_3)
+
+  return x * x * x * x * x * x * x * x; /* x^8 */
+
+#elif defined(HYDRO_GAMMA_2_1)
+
+  return x * x * x * x; /* x^4 */
+
+#else
+
+  error("The adiabatic index is not defined !");
+  return 0.f;
+
+#endif
+}
+
 #endif /* SWIFT_ADIABATIC_INDEX_H */

tests/testRiemann.c

@@ -21,7 +21,7 @@
 #include "error.h"
 
 void check_value(float a, float b, const char* s) {
-  if (fabsf(a - b) > 1.e-6f) {
+  if (fabsf(a - b) > 1.e-5f) {
     error("Values are inconsistent: %g %g (%s)!", a, b, s);
   } else {
     message("Values are consistent: %g %g (%s).", a, b, s);
@@ -57,9 +57,28 @@ void check_constants() {
   check_value(val, hydro_one_over_gamma, "1/gamma");
 }
 
+void check_functions() {
+  float val_a, val_b;
+  const float x = 0.4;
+
+  val_a = pow(x, -hydro_gamma);
+  val_b = pow_minus_gamma(x);
+  check_value(val_a, val_b, "x^(-gamma)");
+
+  val_a = pow(x, 2.0f / (hydro_gamma - 1.0f));
+  val_b = pow_two_over_gamma_minus_one(x);
+  check_value(val_a, val_b, "x^(2/(gamma-1))");
+
+  val_a = pow(x, 2.0f * hydro_gamma / (hydro_gamma - 1.0f));
+  val_b = pow_two_gamma_over_gamma_minus_one(x);
+  check_value(val_a, val_b, "x^((2 gamma)/(gamma-1))");
+}
+
 int main() {
 
   check_constants();
 
+  check_functions();
+
   return 0;
 }