diff --git a/src/hydro/Gizmo/hydro.h b/src/hydro/Gizmo/hydro.h index 6d39c54d2ddc3571ac34c54fc9eede6f7dee6ac5..6e494b43488cac29f12a4520c9c0d6b3ac5d767b 100644 --- a/src/hydro/Gizmo/hydro.h +++ b/src/hydro/Gizmo/hydro.h @@ -225,14 +225,15 @@ __attribute__((always_inline)) INLINE static void hydro_end_density( /* Some smoothing length multiples. */ const float h = p->h; const float ih = 1.0f / h; + const float ihdim = pow_dimension(ih); + const float ihdim_plus_one = ihdim * ih; /* Final operation on the density. */ p->density.wcount += kernel_root; - p->density.wcount *= kernel_norm; - - p->density.wcount_dh *= ih * kernel_gamma * kernel_norm; + p->density.wcount *= ihdim; - const float ihdim = pow_dimension(ih); + p->density.wcount_dh -= hydro_dimension * kernel_root; + p->density.wcount_dh *= ihdim_plus_one; /* Final operation on the geometry. */ /* we multiply with the smoothing kernel normalization ih3 and calculate the @@ -366,6 +367,17 @@ __attribute__((always_inline)) INLINE static void hydro_end_density( p->density.wcount_dh *= p->density.wcorr; } +/** + * @brief Sets all particle fields to sensible values when the #part has 0 ngbs. + * + * @param p The particle to act upon + * @param xp The extended particle data to act upon + */ +__attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours( + struct part* restrict p, struct xpart* restrict xp) { + error("To be implemented"); +} + /** * @brief Prepare a particle for the gradient calculation. * diff --git a/src/hydro/Gizmo/hydro_iact.h b/src/hydro/Gizmo/hydro_iact.h index 8798dc859a790a83ab7a3b6f1709b1302f574581..c7e86c22da5a8beeece0d7a1fcc9cff560105e6e 100644 --- a/src/hydro/Gizmo/hydro_iact.h +++ b/src/hydro/Gizmo/hydro_iact.h @@ -57,7 +57,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_density( kernel_deval(xi, &wi, &wi_dx); pi->density.wcount += wi; - pi->density.wcount_dh -= xi * wi_dx; + pi->density.wcount_dh -= (hydro_dimension * wi + xi * wi_dx); /* these are eqns. (1) and (2) in the summary */ pi->geometry.volume += wi; @@ -74,7 +74,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_density( kernel_deval(xj, &wj, &wj_dx); pj->density.wcount += wj; - pj->density.wcount_dh -= xj * wj_dx; + pj->density.wcount_dh -= (hydro_dimension * wj + xj * wj_dx); /* these are eqns. (1) and (2) in the summary */ pj->geometry.volume += wj; @@ -121,7 +121,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density( kernel_deval(xi, &wi, &wi_dx); pi->density.wcount += wi; - pi->density.wcount_dh -= xi * wi_dx; + pi->density.wcount_dh -= (hydro_dimension * wi + xi * wi_dx); /* these are eqns. (1) and (2) in the summary */ pi->geometry.volume += wi; diff --git a/src/hydro/PressureEntropy/hydro.h b/src/hydro/PressureEntropy/hydro.h index f92f21ba53b4ad112f7bec154ca08a0dac42df61..6d8d06abb126a574fd82c71a809235eb833494c6 100644 --- a/src/hydro/PressureEntropy/hydro.h +++ b/src/hydro/PressureEntropy/hydro.h @@ -263,7 +263,6 @@ __attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours( p->density.rot_v[2] = 0.f; } - /** * @brief Prepare a particle for the force calculation. *