Commit 76469999 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Reset the value of rho_dh to something sensible if no neighbours are found when imposing h=h_max

parent ec42caed
......@@ -210,9 +210,6 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
const float irho = 1.f / p->rho;
/* Compute the derivative term */
p->rho_dh = 1.f / (1.f + hydro_dimension_inv * p->h * p->rho_dh * irho);
/* Finish calculation of the velocity curl components */
p->density.rot_v[0] *= h_inv_dim_plus_one * irho;
p->density.rot_v[1] *= h_inv_dim_plus_one * irho;
......@@ -222,6 +219,24 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
p->density.div_v *= h_inv_dim_plus_one * irho;
}
/**
* @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) {
/* Re-set problematic values */
p->rho_dh = 0.f;
p->density.wcount_dh = 0.f;
p->density.div_v = 0.f;
p->density.rot_v[0] = 0.f;
p->density.rot_v[1] = 0.f;
p->density.rot_v[2] = 0.f;
}
/**
* @brief Prepare a particle for the force calculation.
*
......@@ -249,6 +264,9 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
const float fc = p->force.soundspeed =
sqrtf(hydro_gamma * hydro_gamma_minus_one * u);
/* Compute the derivative term */
p->rho_dh = 1.f / (1.f + hydro_dimension_inv * p->h * p->rho_dh / p->rho);
/* Compute the P/Omega/rho2. */
xp->omega = 1.0f + hydro_dimension_inv * h * p->rho_dh / p->rho;
p->force.P_over_rho2 = u * hydro_gamma_minus_one / (p->rho * xp->omega);
......
......@@ -224,6 +224,24 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
p->density.div_v *= h_inv_dim_plus_one * rho_inv;
}
/**
* @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) {
/* Re-set problematic values */
p->density.rho_dh = 0.f;
p->density.wcount_dh = 0.f;
p->density.div_v = 0.f;
p->density.rot_v[0] = 0.f;
p->density.rot_v[1] = 0.f;
p->density.rot_v[2] = 0.f;
}
/**
* @brief Prepare a particle for the force calculation.
*
......
......@@ -366,6 +366,19 @@ __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) {
/* Re-set problematic values */
p->density.wcount_dh = 0.f;
}
/**
* @brief Prepare a particle for the gradient calculation.
*
......
......@@ -227,6 +227,20 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
p->density.wcount_dh *= h_inv * kernel_gamma * kernel_norm;
}
/**
* @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) {
/* Re-set problematic values */
p->density.rho_dh = 0.f;
p->density.wcount_dh = 0.f;
}
/**
* @brief Prepare a particle for the force calculation.
*
......
......@@ -236,6 +236,25 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
p->density.div_v *= h_inv_dim_plus_one * rho_inv;
}
/**
* @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) {
/* Re-set problematic values */
p->density.rho_dh = 0.f;
p->density.wcount_dh = 0.f;
p->density.pressure_dh = 0.f; // MATTHIEU: to be checked
p->density.div_v = 0.f;
p->density.rot_v[0] = 0.f;
p->density.rot_v[1] = 0.f;
p->density.rot_v[2] = 0.f;
}
/**
* @brief Prepare a particle for the force calculation.
*
......
......@@ -238,6 +238,19 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
#endif
}
/**
* @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) {
/* Re-set problematic values */
p->density.wcount_dh = 0.f;
}
/**
* @brief Prepare a particle for the gradient calculation.
*
......
......@@ -708,6 +708,10 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
/* Ok, this particle is a lost cause... */
p->h = hydro_h_max;
/* Do some damage control if no neighbours at all were found */
if (p->density.wcount == kernel_root * kernel_norm)
hydro_part_has_no_neighbours(p, xp);
}
}
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment