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

Fix to #314 (closed).

src/hydro/Default/hydro.h

@@ -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;
@@ -223,6 +220,31 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
 }
 
 /**
+ * @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) {
+
+  /* Some smoothing length multiples. */
+  const float h = p->h;
+  const float h_inv = 1.0f / h;                 /* 1/h */
+  const float h_inv_dim = pow_dimension(h_inv); /* 1/h^d */
+
+  /* Re-set problematic values */
+  p->rho = p->mass * kernel_root * h_inv_dim;
+  p->density.wcount = kernel_root * kernel_norm * h_inv_dim;
+  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.
  *
  * Computes viscosity term, conduction term and smoothing length gradient terms.
@@ -249,6 +271,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);