Commit 1d205f82 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Added \gamma=1.4 to the list of supported adiabatic indices.

parent fa51a0e6
......@@ -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) */
......
......@@ -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
......
......@@ -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;
......
......@@ -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
......
......@@ -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
......
......@@ -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
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment