Faster implementation of pow(x, 5/3) and more generic version for various implementations

src/Makefile.am

@@ -58,7 +58,7 @@ nobase_noinst_HEADERS = approx_math.h atomic.h cycle.h error.h inline.h kernel_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  \
-	 	 hydro.h hydro_io.h \
+	 	 hydro.h hydro_io.h gamma.h \
 		 hydro/Minimal/hydro.h hydro/Minimal/hydro_iact.h hydro/Minimal/hydro_io.h \
                  hydro/Minimal/hydro_debug.h hydro/Minimal/hydro_part.h \
 		 hydro/Default/hydro.h hydro/Default/hydro_iact.h hydro/Default/hydro_io.h \

src/const.h

@@ -20,9 +20,6 @@
 #ifndef SWIFT_CONST_H
 #define SWIFT_CONST_H
 
-/* Hydrodynamical constants. */
-#define const_hydro_gamma (5.0f / 3.0f)
-
 /* SPH Viscosity constants. */
 #define const_viscosity_alpha 0.8f
 #define const_viscosity_alpha_min \
@@ -44,6 +41,11 @@
 #define const_G 6.672e-8f     /* Gravitational constant. */
 #define const_epsilon 0.0014f /* Gravity blending distance. */
 
+/* Hydrodynamical adiabatic index. */
+#define HYDRO_GAMMA_5_3
+//#define HYDRO_GAMMA_4_3
+//#define HYDRO_GAMMA_2_1
+
 /* Kernel function to use */
 #define CUBIC_SPLINE_KERNEL
 //#define QUARTIC_SPLINE_KERNEL

src/gamma.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_GAMMA_H
+#define SWIFT_GAMMA_H
+
+/* Config parameters. */
+#include "../config.h"
+
+/* Some standard headers. */
+#include <math.h>
+
+/* Local headers. */
+#include "const.h"
+#include "debug.h"
+#include "inline.h"
+
+/* First define some constants */
+#if defined(HYDRO_GAMMA_5_3)
+
+#define hydro_gamma 1.66666666666666667f
+#define hydro_gamma_minus_one 0.66666666666666667f
+
+#elif defined(HYDRO_GAMMA_4_3)
+
+#define hydro_gamma 1.33333333333333333f
+#define hydro_gamma_minus_one 0.33333333333333333f
+
+#elif defined(HYDRO_GAMMA_2_1)
+
+#define hydro_gamma 2.f
+#define hydro_gamma_minus_one 1.f
+
+#endif
+
+/**
+ * @brief Returns the argument to the power given by the adiabatic index
+ *
+ * Computes $x^\gamma$.
+ */
+__attribute__((always_inline)) INLINE static float pow_gamma(float x) {
+
+#if defined(HYDRO_GAMMA_5_3)
+
+  const float x2 = x * x;
+  const float x5 = x2 * x2 * x;
+  return cbrtf(x5);
+
+#elif defined(HYDRO_GAMMA_4_3)
+
+  const float x2 = x * x;
+  const float x4 = x2 * x2;
+  return cbrtf(x4);
+
+#elif defined(HYDRO_GAMMA_2_1)
+
+  return x * x;
+
+#else
+  error("The adiabatic index is not defined !");
+  return 0.f;
+#endif
+}
+
+/**
+ * @brief Returns one over the argument to the power given by the adiabatic
+ * index minus one
+ *
+ * Computes $x^{-(\gamma - 1)}$.
+ */
+__attribute__((always_inline)) INLINE static float pow_minus_gamma_minus_one(
+    float x) {
+
+#if defined(HYDRO_GAMMA_5_3)
+
+  return 1.f / cbrtf(x * x);
+
+#elif defined(HYDRO_GAMMA_4_3)
+
+  return 1.f / cbrtf(x);
+
+#elif defined(HYDRO_GAMMA_2_1)
+
+  return 1.f / x;
+
+#else
+  error("The adiabatic index is not defined !");
+  return 0.f;
+#endif
+}
+
+#endif /* SWIFT_GAMMA_H */

src/hydro/Gadget2/hydro.h

@@ -17,6 +17,8 @@
  *
  ******************************************************************************/
 
+#include "gamma.h"
+
 /**
  * @brief Computes the hydro time-step of a given particle
  *
@@ -132,11 +134,10 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
 
   /* Compute the pressure */
   const float dt = (ti_current - (p->ti_begin + p->ti_end) / 2) * timeBase;
-  p->force.pressure =
-      (p->entropy + p->entropy_dt * dt) * powf(p->rho, const_hydro_gamma);
+  p->force.pressure = (p->entropy + p->entropy_dt * dt) * pow_gamma(p->rho);
 
   /* Compute the sound speed */
-  p->force.soundspeed = sqrtf(const_hydro_gamma * p->force.pressure / p->rho);
+  p->force.soundspeed = sqrtf(hydro_gamma * p->force.pressure / p->rho);
 }
 
 /**
@@ -179,10 +180,10 @@ __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;
   p->force.pressure =
-      (p->entropy + p->entropy_dt * dt_entr) * powf(p->rho, const_hydro_gamma);
+      (p->entropy + p->entropy_dt * dt_entr) * pow_gamma(p->rho);
 
   /* Compute the new sound speed */
-  p->force.soundspeed = sqrtf(const_hydro_gamma * p->force.pressure / p->rho);
+  p->force.soundspeed = sqrtf(hydro_gamma * p->force.pressure / p->rho);
 }
 
 /**
@@ -195,8 +196,7 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
 __attribute__((always_inline)) INLINE static void hydro_end_force(
     struct part* p) {
 
-  p->entropy_dt *=
-      (const_hydro_gamma - 1.f) * powf(p->rho, -(const_hydro_gamma - 1.f));
+  p->entropy_dt *= hydro_gamma_minus_one * pow_minus_gamma_minus_one(p->rho);
 }
 
 /**
@@ -232,8 +232,8 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
 __attribute__((always_inline)) INLINE static void hydro_convert_quantities(
     struct part* p) {
 
-  p->entropy = (const_hydro_gamma - 1.f) * p->entropy *
-               powf(p->rho, -(const_hydro_gamma - 1.f));
+  p->entropy =
+      hydro_gamma_minus_one * p->entropy * pow_minus_gamma_minus_one(p->rho);
 }
 
 /**
@@ -247,6 +247,7 @@ __attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
 
   const float entropy = p->entropy + p->entropy_dt * dt;
 
-  return entropy * powf(p->rho, const_hydro_gamma - 1.f) *
-         (1.f / (const_hydro_gamma - 1.f));
+  return entropy * powf(p->rho, hydro_gamma - 1.f) *
+         (1.f / (hydro_gamma - 1.f));
+  return 1;
 }

src/units.c

@@ -39,6 +39,7 @@
 /* Includes. */
 #include "const.h"
 #include "error.h"
+#include "gamma.h"
 
 /**
  * @brief Initialises the UnitSystem structure with the constants given in
@@ -188,14 +189,14 @@ void units_get_base_unit_exponants_array(float baseUnitsExp[5],
       break;
 
     case UNIT_CONV_ENTROPY:
-      baseUnitsExp[UNIT_MASS] = 1.f - const_hydro_gamma;
-      baseUnitsExp[UNIT_LENGTH] = 3.f * const_hydro_gamma - 1.f;
+      baseUnitsExp[UNIT_MASS] = 1.f - hydro_gamma;
+      baseUnitsExp[UNIT_LENGTH] = 3.f * hydro_gamma - 1.f;
       baseUnitsExp[UNIT_TIME] = -2.f;
       break;
 
     case UNIT_CONV_ENTROPY_PER_UNIT_MASS:
-      baseUnitsExp[UNIT_MASS] = -const_hydro_gamma;
-      baseUnitsExp[UNIT_LENGTH] = 3.f * const_hydro_gamma - 1.f;
+      baseUnitsExp[UNIT_MASS] = -hydro_gamma;
+      baseUnitsExp[UNIT_LENGTH] = 3.f * hydro_gamma - 1.f;
       baseUnitsExp[UNIT_TIME] = -2.f;
       break;