From 364c9eb1b22441c6c134a902c39b08d6ddba4c54 Mon Sep 17 00:00:00 2001 From: Jacob Kegerreis Date: Fri, 7 Dec 2018 13:09:44 +0000 Subject: [PATCH 1/3] If h is h_max then set rho_dh to 0, prevent insanely high grad_h_term --- src/hydro/Planetary/hydro.h | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/src/hydro/Planetary/hydro.h b/src/hydro/Planetary/hydro.h index 667bc9e29..286ccf4de 100644 --- a/src/hydro/Planetary/hydro.h +++ b/src/hydro/Planetary/hydro.h @@ -519,16 +519,21 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force( /* Compute the "grad h" term */ const float rho_inv = 1.f / p->rho; float grad_h_term; + float rho_dh = p->density.rho_dh; + /* Ignore changing-kernel effects when h is h_max */ + if (p->h == hydro_props->h_max) { + rho_dh = 0.f; + } const float grad_h_term_inv = - 1.f + hydro_dimension_inv * p->h * p->density.rho_dh * rho_inv; + 1.f + hydro_dimension_inv * p->h * rho_dh * rho_inv; /* Avoid 1/0 from only having one neighbour right at the edge of the kernel */ - if (grad_h_term_inv != 0.f) { - grad_h_term = 1.f / grad_h_term_inv; - } else { + if (grad_h_term_inv == 0.f) { grad_h_term = 0.f; + } else { + grad_h_term = 1.f / grad_h_term_inv; } - /* Compute the Balsara switch */ + /* Compute the Balsara switch */ #ifdef PLANETARY_SPH_NO_BALSARA const float balsara = hydro_props->viscosity.alpha; #else -- GitLab From 8be370e62a8316a8e11d23453e6b6d0ff1d28477 Mon Sep 17 00:00:00 2001 From: Jacob Kegerreis Date: Wed, 19 Dec 2018 13:27:36 +0000 Subject: [PATCH 2/3] Remove old check that is now unnecessary with the more general fix --- src/hydro/Planetary/hydro.h | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) diff --git a/src/hydro/Planetary/hydro.h b/src/hydro/Planetary/hydro.h index 286ccf4de..67a2d475a 100644 --- a/src/hydro/Planetary/hydro.h +++ b/src/hydro/Planetary/hydro.h @@ -526,12 +526,7 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force( } const float grad_h_term_inv = 1.f + hydro_dimension_inv * p->h * rho_dh * rho_inv; - /* Avoid 1/0 from only having one neighbour right at the edge of the kernel */ - if (grad_h_term_inv == 0.f) { - grad_h_term = 0.f; - } else { - grad_h_term = 1.f / grad_h_term_inv; - } + grad_h_term = 1.f / grad_h_term_inv; /* Compute the Balsara switch */ #ifdef PLANETARY_SPH_NO_BALSARA -- GitLab From 57fba300540dc4245ded13390426317a012aae39 Mon Sep 17 00:00:00 2001 From: Jacob Kegerreis Date: Wed, 19 Dec 2018 13:45:54 +0000 Subject: [PATCH 3/3] Tidy up unnecessary inverse term --- src/hydro/Planetary/hydro.h | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/src/hydro/Planetary/hydro.h b/src/hydro/Planetary/hydro.h index 67a2d475a..887693acc 100644 --- a/src/hydro/Planetary/hydro.h +++ b/src/hydro/Planetary/hydro.h @@ -518,15 +518,13 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force( /* Compute the "grad h" term */ const float rho_inv = 1.f / p->rho; - float grad_h_term; float rho_dh = p->density.rho_dh; /* Ignore changing-kernel effects when h is h_max */ if (p->h == hydro_props->h_max) { rho_dh = 0.f; } - const float grad_h_term_inv = - 1.f + hydro_dimension_inv * p->h * rho_dh * rho_inv; - grad_h_term = 1.f / grad_h_term_inv; + const float grad_h_term = + 1.f / (1.f + hydro_dimension_inv * p->h * rho_dh * rho_inv); /* Compute the Balsara switch */ #ifdef PLANETARY_SPH_NO_BALSARA -- GitLab