diff --git a/src/hydro/Default/hydro.h b/src/hydro/Default/hydro.h index 4a6b1900374eb6b422e6fa7422901293b36fd5eb..a670515381a5fce6bcd524bc12e72f28d597687d 100644 --- a/src/hydro/Default/hydro.h +++ b/src/hydro/Default/hydro.h @@ -18,6 +18,7 @@ ******************************************************************************/ #include "approx_math.h" +#include "gamma.h" /** * @brief Computes the hydro time-step of a given particle @@ -138,12 +139,11 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force( /* Compute this particle's sound speed. */ const float u = p->u; - const float fc = p->force.c = - sqrtf(const_hydro_gamma * (const_hydro_gamma - 1.0f) * u); + const float fc = p->force.c = sqrtf(hydro_gamma * hydro_gamma_minus_one * u); /* Compute the P/Omega/rho2. */ xp->omega = 1.0f + 0.3333333333f * h * p->rho_dh / p->rho; - p->force.POrho2 = u * (const_hydro_gamma - 1.0f) / (p->rho * xp->omega); + p->force.POrho2 = u * hydro_gamma_minus_one / (p->rho * xp->omega); /* Balsara switch */ p->force.balsara = normDiv_v / (normDiv_v + normCurl_v + 0.0001f * fc * ih); @@ -207,7 +207,7 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra( u = p->u *= expf(w); /* Predict gradient term */ - p->force.POrho2 = u * (const_hydro_gamma - 1.0f) / (p->rho * xp->omega); + p->force.POrho2 = u * hydro_gamma_minus_one / (p->rho * xp->omega); } /** diff --git a/src/hydro/Default/hydro_iact.h b/src/hydro/Default/hydro_iact.h index 0a577931b5e0ca67ab07dfe414a548da66e82cdd..cf5bce5c69c625f0536837df3ce9bea7405a5444 100644 --- a/src/hydro/Default/hydro_iact.h +++ b/src/hydro/Default/hydro_iact.h @@ -20,6 +20,8 @@ #ifndef SWIFT_RUNNER_IACT_H #define SWIFT_RUNNER_IACT_H +#include "gamma.h" + /** * @brief SPH interaction functions following the Gadget-2 version of SPH. * @@ -408,7 +410,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_force( Pi_ij *= (pi->force.balsara + pj->force.balsara); /* Thermal conductivity */ - v_sig_u = sqrtf(2.f * (const_hydro_gamma - 1.f) * + v_sig_u = sqrtf(2.f * hydro_gamma_minus_one * fabs(rhoi * pi->u - rhoj * pj->u) / (rhoi + rhoj)); tc = const_conductivity_alpha * v_sig_u / (rhoi + rhoj); tc *= (wi_dr + wj_dr); @@ -608,7 +610,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_vec_force( Pi_ij.v *= (wi_dr.v + wj_dr.v); /* Thermal conductivity */ - v_sig_u.v = vec_sqrt(vec_set1(2.f * (const_hydro_gamma - 1.f)) * + v_sig_u.v = vec_sqrt(vec_set1(2.f * hydro_gamma_minus_one) * vec_fabs(pirho.v * piu.v - pjrho.v * pju.v) / (pirho.v + pjrho.v)); tc.v = vec_set1(const_conductivity_alpha) * v_sig_u.v / (pirho.v + pjrho.v); @@ -721,7 +723,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force( Pi_ij *= (pi->force.balsara + pj->force.balsara); /* Thermal conductivity */ - v_sig_u = sqrtf(2.f * (const_hydro_gamma - 1.f) * + v_sig_u = sqrtf(2.f * hydro_gamma_minus_one * fabs(rhoi * pi->u - rhoj * pj->u) / (rhoi + rhoj)); tc = const_conductivity_alpha * v_sig_u / (rhoi + rhoj); tc *= (wi_dr + wj_dr); diff --git a/src/hydro/Minimal/hydro.h b/src/hydro/Minimal/hydro.h index 4222daafe82e7dc977cd87f57a5b9a235d505f00..da4699d2ebced5650edea8608d6ec97db112090d 100644 --- a/src/hydro/Minimal/hydro.h +++ b/src/hydro/Minimal/hydro.h @@ -18,6 +18,7 @@ ******************************************************************************/ #include "approx_math.h" +#include "gamma.h" /** * @brief Computes the hydro time-step of a given particle @@ -132,7 +133,7 @@ __attribute__((always_inline)) INLINE static void hydro_end_density( __attribute__((always_inline)) INLINE static void hydro_prepare_force( struct part* p, struct xpart* xp, int ti_current, double timeBase) { - p->force.pressure = p->rho * p->u * (const_hydro_gamma - 1.f); + p->force.pressure = p->rho * p->u * hydro_gamma_minus_one; } /** @@ -175,7 +176,7 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra( p->u = xp->u_full; /* Need to recompute the pressure as well */ - p->force.pressure = p->rho * p->u * (const_hydro_gamma - 1.f); + p->force.pressure = p->rho * p->u * hydro_gamma_minus_one; } /** diff --git a/src/hydro/Minimal/hydro_iact.h b/src/hydro/Minimal/hydro_iact.h index c9da185b8a29eafe2a58420ae5de3a05ff043225..dd22aebd7d7f03080c50b079c4ffa87a5ae0ab7d 100644 --- a/src/hydro/Minimal/hydro_iact.h +++ b/src/hydro/Minimal/hydro_iact.h @@ -19,6 +19,8 @@ #ifndef SWIFT_RUNNER_IACT_MINIMAL_H #define SWIFT_RUNNER_IACT_MINIMAL_H +#include "gamma.h" + /** * @brief Minimal conservative implementation of SPH * @@ -128,8 +130,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_force( (pi->v[2] - pj->v[2]) * dx[2]; /* Compute sound speeds */ - const float ci = sqrtf(const_hydro_gamma * pressurei / rhoi); - const float cj = sqrtf(const_hydro_gamma * pressurej / rhoj); + const float ci = sqrtf(hydro_gamma * pressurei / rhoi); + const float cj = sqrtf(hydro_gamma * pressurej / rhoj); const float v_sig = ci + cj; /* SPH acceleration term */ @@ -200,8 +202,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force( (pi->v[2] - pj->v[2]) * dx[2]; /* Compute sound speeds */ - const float ci = sqrtf(const_hydro_gamma * pressurei / rhoi); - const float cj = sqrtf(const_hydro_gamma * pressurej / rhoj); + const float ci = sqrtf(hydro_gamma * pressurei / rhoi); + const float cj = sqrtf(hydro_gamma * pressurej / rhoj); const float v_sig = ci + cj; /* SPH acceleration term */ diff --git a/src/hydro/Minimal/hydro_io.h b/src/hydro/Minimal/hydro_io.h index bed6b8f8221dc3a6e94c65d101c68aa9b3bdea91..d5e52e18ef9fd933904565fe5bbd9c7dd328a14e 100644 --- a/src/hydro/Minimal/hydro_io.h +++ b/src/hydro/Minimal/hydro_io.h @@ -104,7 +104,6 @@ void writeSPHflavour(hid_t h_grpsph) { /* Kernel function description */ writeAttribute_s(h_grpsph, "Kernel", kernel_name); - writeAttribute_f(h_grpsph, "Hydro gamma", const_hydro_gamma); /* Viscosity and thermal conduction */ /* Nothing in this minimal model... */ diff --git a/src/hydro_properties.c b/src/hydro_properties.c index 16216f81a5b505fc3a887e86ca4898bc4179e4d5..45106d1a86651c150e74b352306f6e4e2f2edfb2 100644 --- a/src/hydro_properties.c +++ b/src/hydro_properties.c @@ -26,6 +26,7 @@ /* Local headers. */ #include "error.h" +#include "gamma.h" #include "hydro.h" #include "kernel_hydro.h" @@ -54,6 +55,7 @@ void hydro_props_init(struct hydro_props *p, void hydro_props_print(const struct hydro_props *p) { + message("Adiabatic index gamma: %f.", hydro_gamma); message("Hydrodynamic scheme: %s.", SPH_IMPLEMENTATION); message("Hydrodynamic kernel: %s with %.2f +/- %.2f neighbours (eta=%f).", kernel_name, p->target_neighbours, p->delta_neighbours,