Skip to content
Snippets Groups Projects

If h is h_max then set rho_dh to 0, prevent insanely high grad_h_term

Merged Jacob Kegerreis requested to merge planetary into master
@@ -518,17 +518,15 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
@@ -518,17 +518,15 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
/* Compute the "grad h" term */
/* Compute the "grad h" term */
const float rho_inv = 1.f / p->rho;
const float rho_inv = 1.f / p->rho;
float grad_h_term;
float rho_dh = p->density.rho_dh;
const float grad_h_term_inv =
/* Ignore changing-kernel effects when h is h_max */
1.f + hydro_dimension_inv * p->h * p->density.rho_dh * rho_inv;
if (p->h == hydro_props->h_max) {
/* Avoid 1/0 from only having one neighbour right at the edge of the kernel */
rho_dh = 0.f;
if (grad_h_term_inv != 0.f) {
grad_h_term = 1.f / grad_h_term_inv;
} else {
grad_h_term = 0.f;
}
}
 
const float grad_h_term =
 
1.f / (1.f + hydro_dimension_inv * p->h * rho_dh * rho_inv);
/* Compute the Balsara switch */
/* Compute the Balsara switch */
#ifdef PLANETARY_SPH_NO_BALSARA
#ifdef PLANETARY_SPH_NO_BALSARA
const float balsara = hydro_props->viscosity.alpha;
const float balsara = hydro_props->viscosity.alpha;
#else
#else
Loading