From 1d205f82c442ce94cc7cb89bcb1b904e71456fb0 Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <matthieu.schaller@durham.ac.uk>
Date: Mon, 22 Aug 2016 16:52:15 +0100
Subject: [PATCH] Added \gamma=1.4 to the list of supported adiabatic indices.
---
src/adiabatic_index.h | 81 +++++++++++++++++++++++++++++++++++++------
src/const.h | 1 +
src/potentials.h | 2 +-
tests/test125cells.c | 2 +-
tests/test27cells.c | 5 +--
tests/testPair.c | 2 +-
6 files changed, 78 insertions(+), 15 deletions(-)
diff --git a/src/adiabatic_index.h b/src/adiabatic_index.h
index da0ce943e9..a0c9ce09e3 100644
--- a/src/adiabatic_index.h
+++ b/src/adiabatic_index.h
@@ -35,6 +35,7 @@
/* Local headers. */
#include "const.h"
#include "debug.h"
+#include "error.h"
#include "inline.h"
/* First define some constants */
@@ -47,11 +48,25 @@
#define hydro_gamma_minus_one_over_two_gamma 0.2f
#define hydro_gamma_minus_one_over_gamma_plus_one 0.25f
#define hydro_two_over_gamma_plus_one 0.75f
-#define hydro_two_over_gamma_minus_one 3.0f
+#define hydro_two_over_gamma_minus_one 3.f
#define hydro_gamma_minus_one_over_two 0.33333333333333333f
-#define hydro_two_gamma_over_gamma_minus_one 5.0f
+#define hydro_two_gamma_over_gamma_minus_one 5.f
#define hydro_one_over_gamma 0.6f
+#elif defined(HYDRO_GAMMA_7_5)
+
+#define hydro_gamma 1.4f
+#define hydro_gamma_minus_one 0.4f
+#define hydro_one_over_gamma_minus_one 2.5f
+#define hydro_gamma_plus_one_over_two_gamma 0.857142857f
+#define hydro_gamma_minus_one_over_two_gamma 0.142857143f
+#define hydro_gamma_minus_one_over_gamma_plus_one 0.166666667f
+#define hydro_two_over_gamma_plus_one 0.833333333
+#define hydro_two_over_gamma_minus_one 5.f
+#define hydro_gamma_minus_one_over_two 0.2f
+#define hydro_two_gamma_over_gamma_minus_one 7.f
+#define hydro_one_over_gamma 0.714285714f
+
#elif defined(HYDRO_GAMMA_4_3)
#define hydro_gamma 1.33333333333333333f
@@ -61,9 +76,9 @@
#define hydro_gamma_minus_one_over_two_gamma 0.125f
#define hydro_gamma_minus_one_over_gamma_plus_one 0.142857143f
#define hydro_two_over_gamma_plus_one 0.857142857f
-#define hydro_two_over_gamma_minus_one 6.0f
+#define hydro_two_over_gamma_minus_one 6.f
#define hydro_gamma_minus_one_over_two 0.166666666666666666f
-#define hydro_two_gamma_over_gamma_minus_one 8.0f
+#define hydro_two_gamma_over_gamma_minus_one 8.f
#define hydro_one_over_gamma 0.75f
#elif defined(HYDRO_GAMMA_2_1)
@@ -75,9 +90,9 @@
#define hydro_gamma_minus_one_over_two_gamma 0.25f
#define hydro_gamma_minus_one_over_gamma_plus_one 0.33333333333333333f
#define hydro_two_over_gamma_plus_one 0.66666666666666666f
-#define hydro_two_over_gamma_minus_one 2.0f
+#define hydro_two_over_gamma_minus_one 2.f
#define hydro_gamma_minus_one_over_two 0.5f
-#define hydro_two_gamma_over_gamma_minus_one 4.0f
+#define hydro_two_gamma_over_gamma_minus_one 4.f
#define hydro_one_over_gamma 0.5f
#else
@@ -98,6 +113,10 @@ __attribute__((always_inline)) INLINE static float pow_gamma(float x) {
const float cbrt = cbrtf(x); /* x^(1/3) */
return cbrt * cbrt * x; /* x^(5/3) */
+#elif defined(HYDRO_GAMMA_7_5)
+
+ return powf(x, 1.4f); /* x^(7/5) */
+
#elif defined(HYDRO_GAMMA_4_3)
return cbrtf(x) * x; /* x^(4/3) */
@@ -128,6 +147,10 @@ __attribute__((always_inline)) INLINE static float pow_gamma_minus_one(
const float cbrt = cbrtf(x); /* x^(1/3) */
return cbrt * cbrt; /* x^(2/3) */
+#elif defined(HYDRO_GAMMA_7_5)
+
+ return powf(x, 0.4f); /* x^(2/5) */
+
#elif defined(HYDRO_GAMMA_4_3)
return cbrtf(x); /* x^(1/3) */
@@ -158,6 +181,10 @@ __attribute__((always_inline)) INLINE static float pow_minus_gamma_minus_one(
const float cbrt_inv = 1.f / cbrtf(x); /* x^(-1/3) */
return cbrt_inv * cbrt_inv; /* x^(-2/3) */
+#elif defined(HYDRO_GAMMA_7_5)
+
+ return powf(x, -0.4f); /* x^(-2/5) */
+
#elif defined(HYDRO_GAMMA_4_3)
return 1.f / cbrtf(x); /* x^(-1/3) */
@@ -191,6 +218,10 @@ __attribute__((always_inline)) INLINE static float pow_minus_gamma(float x) {
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_7_5)
+
+ return powf(x, -1.4f); /* x^(-7/5) */
+
#elif defined(HYDRO_GAMMA_4_3)
const float cbrt_inv = 1.f / cbrtf(x); /* x^(-1/3) */
@@ -226,9 +257,16 @@ __attribute__((always_inline)) INLINE static float pow_two_over_gamma_minus_one(
return x * x * x; /* x^3 */
+#elif defined(HYDRO_GAMMA_7_5)
+
+ const float x2 = x * x;
+ const float x3 = x2 * x;
+ return x2 * x3;
+
#elif defined(HYDRO_GAMMA_4_3)
- return x * x * x * x * x * x; /* x^6 */
+ const float x3 = x * x * x; /* x^3 */
+ return x3 * x3; /* x^6 */
#elif defined(HYDRO_GAMMA_2_1)
@@ -257,15 +295,26 @@ pow_two_gamma_over_gamma_minus_one(float x) {
#if defined(HYDRO_GAMMA_5_3)
- return x * x * x * x * x; /* x^5 */
+ const float x2 = x * x;
+ const float x3 = x2 * x;
+ return x2 * x3;
+
+#elif defined(HYDRO_GAMMA_7_5)
+
+ const float x2 = x * x;
+ const float x4 = x2 * x2;
+ return x4 * x2 * x;
#elif defined(HYDRO_GAMMA_4_3)
- return x * x * x * x * x * x * x * x; /* x^8 */
+ const float x2 = x * x;
+ const float x4 = x2 * x2;
+ return x4 * x4; /* x^8 */
#elif defined(HYDRO_GAMMA_2_1)
- return x * x * x * x; /* x^4 */
+ const float x2 = x * x;
+ return x2 * x2; /* x^4 */
#else
@@ -292,6 +341,10 @@ pow_gamma_minus_one_over_two_gamma(float x) {
return powf(x, 0.2f); /* x^0.2 */
+#elif defined(HYDRO_GAMMA_7_5)
+
+ return powf(x, hydro_gamma_minus_one_over_two_gamma);
+
#elif defined(HYDRO_GAMMA_4_3)
return powf(x, 0.125f); /* x^0.125 */
@@ -325,6 +378,10 @@ pow_minus_gamma_plus_one_over_two_gamma(float x) {
return powf(x, -0.8f); /* x^-0.8 */
+#elif defined(HYDRO_GAMMA_7_5)
+
+ return powf(x, -hydro_gamma_plus_one_over_two_gamma);
+
#elif defined(HYDRO_GAMMA_4_3)
return powf(x, -0.875f); /* x^-0.875 */
@@ -355,6 +412,10 @@ __attribute__((always_inline)) INLINE static float pow_one_over_gamma(float x) {
return powf(x, 0.6f); /* x^(3/5) */
+#elif defined(HYDRO_GAMMA_7_5)
+
+ return powf(x, hydro_one_over_gamma);
+
#elif defined(HYDRO_GAMMA_4_3)
return powf(x, 0.75f); /* x^(3/4) */
diff --git a/src/const.h b/src/const.h
index 64d07b3f3b..5f7bf0d264 100644
--- a/src/const.h
+++ b/src/const.h
@@ -46,6 +46,7 @@
/* Hydrodynamical adiabatic index. */
#define HYDRO_GAMMA_5_3
+//#define HYDRO_GAMMA_7_5
//#define HYDRO_GAMMA_4_3
//#define HYDRO_GAMMA_2_1
diff --git a/src/potentials.h b/src/potentials.h
index 010eb62784..74f0fd2856 100644
--- a/src/potentials.h
+++ b/src/potentials.h
@@ -146,7 +146,7 @@ external_gravity_disk_patch_timestep(
__attribute__((always_inline)) INLINE static void
external_gravity_disk_patch_potential(
double time, const struct external_potential* restrict potential,
- const struct phys_const* restrict phys_const, struct gpart *restrict g) {
+ const struct phys_const* restrict phys_const, struct gpart* restrict g) {
const float G_newton = phys_const->const_newton_G;
const float dz = g->x[2] - potential->disk_patch_potential.z_disk;
diff --git a/tests/test125cells.c b/tests/test125cells.c
index f3d02a68b4..c7e01693b4 100644
--- a/tests/test125cells.c
+++ b/tests/test125cells.c
@@ -259,7 +259,7 @@ struct cell *make_cell(size_t n, const double offset[3], double size, double h,
part->h = size * h / (float)n;
#ifdef GIZMO_SPH
- part->conserved.mass = density * volume / count;
+ part->conserved.mass = density * volume / count;
#else
part->mass = density * volume / count;
#endif
diff --git a/tests/test27cells.c b/tests/test27cells.c
index 739d316231..1a1ab88748 100644
--- a/tests/test27cells.c
+++ b/tests/test27cells.c
@@ -105,7 +105,7 @@ struct cell *make_cell(size_t n, double *offset, double size, double h,
part->h = size * h / (float)n;
part->id = ++(*partId);
#ifdef GIZMO_SPH
- part->conserved.mass = density * volume / count;
+ part->conserved.mass = density * volume / count;
#else
part->mass = density * volume / count;
#endif
@@ -188,7 +188,8 @@ void dump_particle_fields(char *fileName, struct cell *main_cell,
main_cell->parts[pid].id, main_cell->parts[pid].x[0],
main_cell->parts[pid].x[1], main_cell->parts[pid].x[2],
main_cell->parts[pid].v[0], main_cell->parts[pid].v[1],
- main_cell->parts[pid].v[2], hydro_get_density(&main_cell->parts[pid]),
+ main_cell->parts[pid].v[2],
+ hydro_get_density(&main_cell->parts[pid]),
#if defined(GIZMO_SPH)
0.f,
#else
diff --git a/tests/testPair.c b/tests/testPair.c
index a77409eeeb..efa1e628c2 100644
--- a/tests/testPair.c
+++ b/tests/testPair.c
@@ -64,7 +64,7 @@ struct cell *make_cell(size_t n, double *offset, double size, double h,
part->h = size * h / (float)n;
part->id = ++(*partId);
#ifdef GIZMO_SPH
- part->conserved.mass = density * volume / count;
+ part->conserved.mass = density * volume / count;
#else
part->mass = density * volume / count;
#endif
--
GitLab