From 3922dc4ff43b2b9ab38fcdec81e8ee9430064ccf Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <matthieu.schaller@durham.ac.uk>
Date: Wed, 13 Jul 2016 15:19:37 +0100
Subject: [PATCH] Faster implementation of pow(x, 5/3) and more generic version
for various implementations
---
src/Makefile.am | 2 +-
src/const.h | 8 +--
src/gamma.h | 107 ++++++++++++++++++++++++++++++++++++++
src/hydro/Gadget2/hydro.h | 23 ++++----
src/units.c | 9 ++--
5 files changed, 130 insertions(+), 19 deletions(-)
create mode 100644 src/gamma.h
diff --git a/src/Makefile.am b/src/Makefile.am
index bc263d6fb0..3b78a428f6 100644
--- a/src/Makefile.am
+++ b/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 \
diff --git a/src/const.h b/src/const.h
index 673aa3bf00..65404d6a3e 100644
--- a/src/const.h
+++ b/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
diff --git a/src/gamma.h b/src/gamma.h
new file mode 100644
index 0000000000..f51f16b5c0
--- /dev/null
+++ b/src/gamma.h
@@ -0,0 +1,107 @@
+/*******************************************************************************
+ * 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 */
diff --git a/src/hydro/Gadget2/hydro.h b/src/hydro/Gadget2/hydro.h
index 0973acb0fb..4081679d49 100644
--- a/src/hydro/Gadget2/hydro.h
+++ b/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;
}
diff --git a/src/units.c b/src/units.c
index 8e0ff08f1d..ca75fc7ee6 100644
--- a/src/units.c
+++ b/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;
--
GitLab