Skip to content
Snippets Groups Projects
Commit a08bc6b1 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Port changes to other hydro schemes. Adevertize what we do at start-up

parent 3922dc4f
No related branches found
No related tags found
1 merge request!195Gamma power speed up
...@@ -18,6 +18,7 @@ ...@@ -18,6 +18,7 @@
******************************************************************************/ ******************************************************************************/
#include "approx_math.h" #include "approx_math.h"
#include "gamma.h"
/** /**
* @brief Computes the hydro time-step of a given particle * @brief Computes the hydro time-step of a given particle
...@@ -138,12 +139,11 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force( ...@@ -138,12 +139,11 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
/* Compute this particle's sound speed. */ /* Compute this particle's sound speed. */
const float u = p->u; const float u = p->u;
const float fc = p->force.c = const float fc = p->force.c = sqrtf(hydro_gamma * hydro_gamma_minus_one * u);
sqrtf(const_hydro_gamma * (const_hydro_gamma - 1.0f) * u);
/* Compute the P/Omega/rho2. */ /* Compute the P/Omega/rho2. */
xp->omega = 1.0f + 0.3333333333f * h * p->rho_dh / p->rho; 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 */ /* Balsara switch */
p->force.balsara = normDiv_v / (normDiv_v + normCurl_v + 0.0001f * fc * ih); 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( ...@@ -207,7 +207,7 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
u = p->u *= expf(w); u = p->u *= expf(w);
/* Predict gradient term */ /* 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);
} }
/** /**
......
...@@ -20,6 +20,8 @@ ...@@ -20,6 +20,8 @@
#ifndef SWIFT_RUNNER_IACT_H #ifndef SWIFT_RUNNER_IACT_H
#define SWIFT_RUNNER_IACT_H #define SWIFT_RUNNER_IACT_H
#include "gamma.h"
/** /**
* @brief SPH interaction functions following the Gadget-2 version of SPH. * @brief SPH interaction functions following the Gadget-2 version of SPH.
* *
...@@ -408,7 +410,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_force( ...@@ -408,7 +410,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
Pi_ij *= (pi->force.balsara + pj->force.balsara); Pi_ij *= (pi->force.balsara + pj->force.balsara);
/* Thermal conductivity */ /* 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)); fabs(rhoi * pi->u - rhoj * pj->u) / (rhoi + rhoj));
tc = const_conductivity_alpha * v_sig_u / (rhoi + rhoj); tc = const_conductivity_alpha * v_sig_u / (rhoi + rhoj);
tc *= (wi_dr + wj_dr); tc *= (wi_dr + wj_dr);
...@@ -608,7 +610,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_vec_force( ...@@ -608,7 +610,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_vec_force(
Pi_ij.v *= (wi_dr.v + wj_dr.v); Pi_ij.v *= (wi_dr.v + wj_dr.v);
/* Thermal conductivity */ /* 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) / vec_fabs(pirho.v * piu.v - pjrho.v * pju.v) /
(pirho.v + pjrho.v)); (pirho.v + pjrho.v));
tc.v = vec_set1(const_conductivity_alpha) * v_sig_u.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( ...@@ -721,7 +723,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
Pi_ij *= (pi->force.balsara + pj->force.balsara); Pi_ij *= (pi->force.balsara + pj->force.balsara);
/* Thermal conductivity */ /* 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)); fabs(rhoi * pi->u - rhoj * pj->u) / (rhoi + rhoj));
tc = const_conductivity_alpha * v_sig_u / (rhoi + rhoj); tc = const_conductivity_alpha * v_sig_u / (rhoi + rhoj);
tc *= (wi_dr + wj_dr); tc *= (wi_dr + wj_dr);
......
...@@ -18,6 +18,7 @@ ...@@ -18,6 +18,7 @@
******************************************************************************/ ******************************************************************************/
#include "approx_math.h" #include "approx_math.h"
#include "gamma.h"
/** /**
* @brief Computes the hydro time-step of a given particle * @brief Computes the hydro time-step of a given particle
...@@ -132,7 +133,7 @@ __attribute__((always_inline)) INLINE static void hydro_end_density( ...@@ -132,7 +133,7 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
__attribute__((always_inline)) INLINE static void hydro_prepare_force( __attribute__((always_inline)) INLINE static void hydro_prepare_force(
struct part* p, struct xpart* xp, int ti_current, double timeBase) { 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( ...@@ -175,7 +176,7 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
p->u = xp->u_full; p->u = xp->u_full;
/* Need to recompute the pressure as well */ /* 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;
} }
/** /**
......
...@@ -19,6 +19,8 @@ ...@@ -19,6 +19,8 @@
#ifndef SWIFT_RUNNER_IACT_MINIMAL_H #ifndef SWIFT_RUNNER_IACT_MINIMAL_H
#define SWIFT_RUNNER_IACT_MINIMAL_H #define SWIFT_RUNNER_IACT_MINIMAL_H
#include "gamma.h"
/** /**
* @brief Minimal conservative implementation of SPH * @brief Minimal conservative implementation of SPH
* *
...@@ -128,8 +130,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_force( ...@@ -128,8 +130,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
(pi->v[2] - pj->v[2]) * dx[2]; (pi->v[2] - pj->v[2]) * dx[2];
/* Compute sound speeds */ /* Compute sound speeds */
const float ci = sqrtf(const_hydro_gamma * pressurei / rhoi); const float ci = sqrtf(hydro_gamma * pressurei / rhoi);
const float cj = sqrtf(const_hydro_gamma * pressurej / rhoj); const float cj = sqrtf(hydro_gamma * pressurej / rhoj);
const float v_sig = ci + cj; const float v_sig = ci + cj;
/* SPH acceleration term */ /* SPH acceleration term */
...@@ -200,8 +202,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force( ...@@ -200,8 +202,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
(pi->v[2] - pj->v[2]) * dx[2]; (pi->v[2] - pj->v[2]) * dx[2];
/* Compute sound speeds */ /* Compute sound speeds */
const float ci = sqrtf(const_hydro_gamma * pressurei / rhoi); const float ci = sqrtf(hydro_gamma * pressurei / rhoi);
const float cj = sqrtf(const_hydro_gamma * pressurej / rhoj); const float cj = sqrtf(hydro_gamma * pressurej / rhoj);
const float v_sig = ci + cj; const float v_sig = ci + cj;
/* SPH acceleration term */ /* SPH acceleration term */
......
...@@ -104,7 +104,6 @@ void writeSPHflavour(hid_t h_grpsph) { ...@@ -104,7 +104,6 @@ void writeSPHflavour(hid_t h_grpsph) {
/* Kernel function description */ /* Kernel function description */
writeAttribute_s(h_grpsph, "Kernel", kernel_name); writeAttribute_s(h_grpsph, "Kernel", kernel_name);
writeAttribute_f(h_grpsph, "Hydro gamma", const_hydro_gamma);
/* Viscosity and thermal conduction */ /* Viscosity and thermal conduction */
/* Nothing in this minimal model... */ /* Nothing in this minimal model... */
......
...@@ -26,6 +26,7 @@ ...@@ -26,6 +26,7 @@
/* Local headers. */ /* Local headers. */
#include "error.h" #include "error.h"
#include "gamma.h"
#include "hydro.h" #include "hydro.h"
#include "kernel_hydro.h" #include "kernel_hydro.h"
...@@ -54,6 +55,7 @@ void hydro_props_init(struct hydro_props *p, ...@@ -54,6 +55,7 @@ void hydro_props_init(struct hydro_props *p,
void hydro_props_print(const 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 scheme: %s.", SPH_IMPLEMENTATION);
message("Hydrodynamic kernel: %s with %.2f +/- %.2f neighbours (eta=%f).", message("Hydrodynamic kernel: %s with %.2f +/- %.2f neighbours (eta=%f).",
kernel_name, p->target_neighbours, p->delta_neighbours, kernel_name, p->target_neighbours, p->delta_neighbours,
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment