diff --git a/src/adiabatic_index.h b/src/adiabatic_index.h index e741a6d1fb33d01003e7916bd0fdc32dc581625f..e7798f7be6e28536191475f3c940c6f84a92bdad 100644 --- a/src/adiabatic_index.h +++ b/src/adiabatic_index.h @@ -437,7 +437,6 @@ __attribute__((always_inline)) INLINE static float pow_one_over_gamma(float x) { * Computes \f$x^{3\gamma - 2}\f$. * * @param x Argument - * @return Argument to the power one over the adiabatic index */ __attribute__((always_inline)) INLINE static float pow_three_gamma_minus_two( float x) { @@ -466,4 +465,39 @@ __attribute__((always_inline)) INLINE static float pow_three_gamma_minus_two( #endif } +/** + * @brief Return the argument to the power three adiabatic index minus five over + * two. + * + * Computes \f$x^{(3\gamma - 5)/2}\f$. + * + * @param x Argument + */ +__attribute__((always_inline)) INLINE static float +pow_three_gamma_minus_five_over_two(float x) { + +#if defined(HYDRO_GAMMA_5_3) + + return 1.f; /* x^(0) */ + +#elif defined(HYDRO_GAMMA_7_5) + + return powf(x, -0.4f); /* x^(-2/5) */ + +#elif defined(HYDRO_GAMMA_4_3) + + return 1.f / sqrtf(x); /* x^(-1/2) */ + +#elif defined(HYDRO_GAMMA_2_1) + + return sqrtf(x); /* x^(1/2) */ + +#else + + error("The adiabatic index is not defined !"); + return 0.f; + +#endif +} + #endif /* SWIFT_ADIABATIC_INDEX_H */ diff --git a/src/hydro/Gadget2/hydro_iact.h b/src/hydro/Gadget2/hydro_iact.h index e5396988333d77801122e2508f9c2ae6c1ab69fc..56efbe945d7ebe8c0d21364153ecc443cbf51c55 100644 --- a/src/hydro/Gadget2/hydro_iact.h +++ b/src/hydro/Gadget2/hydro_iact.h @@ -436,9 +436,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_force( float wi, wj, wi_dx, wj_dx; - /* Will change with cosmological integration */ - const float fac_mu = 1.f; - const float a2_Hubble = 0.f; + /* Cosmological factors entering the EoMs */ + const float fac_mu = pow_three_gamma_minus_five_over_two(a); + const float a2_Hubble = a * a * H; const float r = sqrtf(r2); const float r_inv = 1.0f / r; @@ -555,9 +555,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force( float wi, wj, wi_dx, wj_dx; - /* Will change with cosmological integration */ - const float fac_mu = 1.f; - const float a2_Hubble = 0.f; + /* Cosmological factors entering the EoMs */ + const float fac_mu = pow_three_gamma_minus_five_over_two(a); + const float a2_Hubble = a * a * H; const float r = sqrtf(r2); const float r_inv = 1.0f / r;