diff --git a/src/adiabatic_index.h b/src/adiabatic_index.h index da0ce943e966c8061bd8732c0335f63766a8621c..a0c9ce09e3e004af07e8b208ef9f1af5f46c9e81 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 64d07b3f3b883dd9f574ee12f4bd98bdb53675a1..5f7bf0d26431887c8e37d96b20b0eba460de55a6 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 010eb627848a9bc2b35dc50e22e5dbd94a0ec3fb..74f0fd28566f355962c83e5d743aeae9afe09c59 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 f3d02a68b4ae1e89ebd88ddd988b9be96c8c36ac..c7e01693b45f76fa21ffd397289fc06ad36af03d 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 739d31623168b0933dd9aaa7d6ed44e34a9ace07..1a1ab88748d922b3e7fbb30a73a10809dca10863 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 a77409eeeb763c3ca0fd002783784c4796600e6e..efa1e628c2d57bf7922be8affdd5436ebca2f9cf 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