From ef9ba85f020e7c7ac05ca5fba2c846a2ac36b688 Mon Sep 17 00:00:00 2001
From: Jacob Kegerreis <jacob.kegerreis@durham.ac.uk>
Date: Tue, 29 May 2018 16:05:07 +0100
Subject: [PATCH 1/6] Flag explicitly when a particle has no neighbours. Add
 no-neighbour function for chemistry.

---
 src/chemistry/EAGLE/chemistry.h | 16 ++++++++++++++++
 src/chemistry/GEAR/chemistry.h  | 16 ++++++++++++++++
 src/chemistry/none/chemistry.h  | 14 ++++++++++++++
 src/runner.c                    |  8 +++++++-
 4 files changed, 53 insertions(+), 1 deletion(-)

diff --git a/src/chemistry/EAGLE/chemistry.h b/src/chemistry/EAGLE/chemistry.h
index 459de24ef..f438425b5 100644
--- a/src/chemistry/EAGLE/chemistry.h
+++ b/src/chemistry/EAGLE/chemistry.h
@@ -78,6 +78,22 @@ __attribute__((always_inline)) INLINE static void chemistry_end_density(
     struct part* restrict p, const struct chemistry_global_data* cd,
     const struct cosmology* cosmo) {}
 
+/**
+ * @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
+ * @param cd #chemistry_global_data containing chemistry informations.
+ * @param cosmo The current cosmological model.
+ */
+__attribute__((always_inline)) INLINE static void
+chemistry_part_has_no_neighbours(struct part* restrict p,
+                                 struct xpart* restrict xp,
+                                 const struct chemistry_global_data* cd,
+                                 const struct cosmology* cosmo) {
+  error("Needs implementing!");
+}
+
 /**
  * @brief Sets the chemistry properties of the (x-)particles to a valid start
  * state.
diff --git a/src/chemistry/GEAR/chemistry.h b/src/chemistry/GEAR/chemistry.h
index a51051ca3..b138e66a8 100644
--- a/src/chemistry/GEAR/chemistry.h
+++ b/src/chemistry/GEAR/chemistry.h
@@ -134,6 +134,22 @@ __attribute__((always_inline)) INLINE static void chemistry_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
+ * @param cd #chemistry_global_data containing chemistry informations.
+ * @param cosmo The current cosmological model.
+ */
+__attribute__((always_inline)) INLINE static void
+chemistry_part_has_no_neighbours(struct part* restrict p,
+                                 struct xpart* restrict xp,
+                                 const struct chemistry_global_data* cd,
+                                 const struct cosmology* cosmo) {
+  error("Needs implementing!");
+}
+
 /**
  * @brief Sets the chemistry properties of the (x-)particles to a valid start
  * state.
diff --git a/src/chemistry/none/chemistry.h b/src/chemistry/none/chemistry.h
index 3ca51660d..6c2e7a21a 100644
--- a/src/chemistry/none/chemistry.h
+++ b/src/chemistry/none/chemistry.h
@@ -86,6 +86,20 @@ __attribute__((always_inline)) INLINE static void chemistry_end_density(
     struct part* restrict p, const struct chemistry_global_data* cd,
     const struct cosmology* cosmo) {}
 
+/**
+ * @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
+ * @param cd #chemistry_global_data containing chemistry informations.
+ * @param cosmo The current cosmological model.
+ */
+__attribute__((always_inline)) INLINE static void
+chemistry_part_has_no_neighbours(struct part* restrict p,
+                                 struct xpart* restrict xp,
+                                 const struct chemistry_global_data* cd,
+                                 const struct cosmology* cosmo) {}
+
 /**
  * @brief Sets the chemistry properties of the (x-)particles to a valid start
  * state.
diff --git a/src/runner.c b/src/runner.c
index 5a382a916..b86f84b26 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -710,9 +710,13 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
         const float h_old_dim = pow_dimension(h_old);
         const float h_old_dim_minus_one = pow_dimension_minus_one(h_old);
         float h_new;
+        int has_no_neighbours = 0;
 
         if (p->density.wcount == 0.f) { /* No neighbours case */
 
+          /* Flag that there were no neighbours */
+          has_no_neighbours = 1;
+
           /* Double h and try again */
           h_new = 2.f * h_old;
         } else {
@@ -768,8 +772,10 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
             p->h = hydro_h_max;
 
             /* Do some damage control if no neighbours at all were found */
-            if (p->density.wcount == kernel_root * kernel_norm)
+            if (has_no_neighbours) {
               hydro_part_has_no_neighbours(p, xp, cosmo);
+              chemistry_part_has_no_neighbours(p, xp, chemistry, cosmo);
+            }
           }
         }
 
-- 
GitLab


From 4bfc9a0c2fcb8ef153cfba5c4474d927b2e2697b Mon Sep 17 00:00:00 2001
From: Jacob Kegerreis <jacob.kegerreis@durham.ac.uk>
Date: Tue, 29 May 2018 16:07:21 +0100
Subject: [PATCH 2/6] Avoid floating point exception from f_prime = 0

---
 src/runner.c | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/runner.c b/src/runner.c
index b86f84b26..3d595fec9 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -733,7 +733,8 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
               p->density.wcount_dh * h_old_dim +
               hydro_dimension * p->density.wcount * h_old_dim_minus_one;
 
-          h_new = h_old - f / f_prime;
+          /* Avoid floating point exception from f_prime = 0 */
+          h_new = h_old - f / (f_prime + FLT_MIN);
 
 #ifdef SWIFT_DEBUG_CHECKS
           if ((f > 0.f && h_new > h_old) || (f < 0.f && h_new < h_old))
-- 
GitLab


From bb129f3d980da338b44e693e21413f0da1b5dd12 Mon Sep 17 00:00:00 2001
From: Jacob Kegerreis <jacob.kegerreis@durham.ac.uk>
Date: Tue, 29 May 2018 16:11:28 +0100
Subject: [PATCH 3/6] Update hydro no-neighbour and reset functions

---
 src/hydro/Minimal/hydro.h         | 5 +++--
 src/hydro/MinimalMultiMat/hydro.h | 5 +++--
 2 files changed, 6 insertions(+), 4 deletions(-)

diff --git a/src/hydro/Minimal/hydro.h b/src/hydro/Minimal/hydro.h
index 3f9d99683..8bbaacfbe 100644
--- a/src/hydro/Minimal/hydro.h
+++ b/src/hydro/Minimal/hydro.h
@@ -367,9 +367,10 @@ __attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours(
 
   /* Re-set problematic values */
   p->rho = p->mass * kernel_root * h_inv_dim;
-  p->density.wcount = kernel_root * kernel_norm * h_inv_dim;
+  p->density.wcount = kernel_root * h_inv_dim;
   p->density.rho_dh = 0.f;
   p->density.wcount_dh = 0.f;
+  p->force.v_sig = 1.f;
 }
 
 /**
@@ -426,7 +427,7 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration(
   /* Reset the time derivatives. */
   p->u_dt = 0.0f;
   p->force.h_dt = 0.0f;
-  p->force.v_sig = 0.0f;
+  p->force.v_sig = p->force.soundspeed;
 }
 
 /**
diff --git a/src/hydro/MinimalMultiMat/hydro.h b/src/hydro/MinimalMultiMat/hydro.h
index 5383ffda8..405b4b5b0 100644
--- a/src/hydro/MinimalMultiMat/hydro.h
+++ b/src/hydro/MinimalMultiMat/hydro.h
@@ -368,9 +368,10 @@ __attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours(
 
   /* Re-set problematic values */
   p->rho = p->mass * kernel_root * h_inv_dim;
-  p->density.wcount = kernel_root * kernel_norm * h_inv_dim;
+  p->density.wcount = kernel_root * h_inv_dim;
   p->density.rho_dh = 0.f;
   p->density.wcount_dh = 0.f;
+  p->force.v_sig = 1.f;
 }
 
 /**
@@ -429,7 +430,7 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration(
   /* Reset the time derivatives. */
   p->u_dt = 0.0f;
   p->force.h_dt = 0.0f;
-  p->force.v_sig = 0.0f;
+  p->force.v_sig = p->force.soundspeed;
 }
 
 /**
-- 
GitLab


From 8134b7834a7e49e3a9412ac4cf7b381b7e50ee5e Mon Sep 17 00:00:00 2001
From: Jacob Kegerreis <jacob.kegerreis@durham.ac.uk>
Date: Tue, 29 May 2018 16:13:48 +0100
Subject: [PATCH 4/6] Tidy hydro debug

---
 src/hydro/Minimal/hydro_debug.h         | 21 ++++++++++++---------
 src/hydro/MinimalMultiMat/hydro_debug.h | 21 ++++++++++++---------
 2 files changed, 24 insertions(+), 18 deletions(-)

diff --git a/src/hydro/Minimal/hydro_debug.h b/src/hydro/Minimal/hydro_debug.h
index 541029ee0..eb0c99dd4 100644
--- a/src/hydro/Minimal/hydro_debug.h
+++ b/src/hydro/Minimal/hydro_debug.h
@@ -36,15 +36,18 @@
 __attribute__((always_inline)) INLINE static void hydro_debug_particle(
     const struct part* p, const struct xpart* xp) {
   printf(
-      "x=[%.3e,%.3e,%.3e], "
-      "v=[%.3e,%.3e,%.3e],v_full=[%.3e,%.3e,%.3e] \n a=[%.3e,%.3e,%.3e], "
-      "u=%.3e, du/dt=%.3e v_sig=%.3e, P=%.3e\n"
-      "h=%.3e, dh/dt=%.3e wcount=%d, m=%.3e, dh_drho=%.3e, rho=%.3e, "
-      "time_bin=%d\n",
-      p->x[0], p->x[1], p->x[2], p->v[0], p->v[1], p->v[2], xp->v_full[0],
-      xp->v_full[1], xp->v_full[2], p->a_hydro[0], p->a_hydro[1], p->a_hydro[2],
-      p->u, p->u_dt, p->force.v_sig, hydro_get_comoving_pressure(p), p->h,
-      p->force.h_dt, (int)p->density.wcount, p->mass, p->density.rho_dh, p->rho,
+      "\n "
+      "x=[%.6g, %.6g, %.6g], v=[%.3g, %.3g, %.3g], \n "
+      "v_full=[%.3g, %.3g, %.3g], a=[%.3g, %.3g, %.3g], \n "
+      "m=%.3g, u=%.3g, du/dt=%.3g, P=%.3g, c_s=%.3g, \n "
+      "v_sig=%.3g, h=%.3g, dh/dt=%.3g, wcount=%.3g, rho=%.3g, \n "
+      "dh_drho=%.3g, time_bin=%d \n",
+      p->x[0], p->x[1], p->x[2], p->v[0], p->v[1], p->v[2],
+      xp->v_full[0], xp->v_full[1], xp->v_full[2],
+      p->a_hydro[0], p->a_hydro[1], p->a_hydro[2],
+      p->mass, p->u, p->u_dt, hydro_get_comoving_pressure(p),
+      p->force.soundspeed, p->force.v_sig, p->h, p->force.h_dt,
+      p->density.wcount, p->rho, p->density.rho_dh,
       p->time_bin);
 }
 
diff --git a/src/hydro/MinimalMultiMat/hydro_debug.h b/src/hydro/MinimalMultiMat/hydro_debug.h
index d8fe73313..01d4d24b1 100644
--- a/src/hydro/MinimalMultiMat/hydro_debug.h
+++ b/src/hydro/MinimalMultiMat/hydro_debug.h
@@ -38,15 +38,18 @@
 __attribute__((always_inline)) INLINE static void hydro_debug_particle(
     const struct part* p, const struct xpart* xp) {
   printf(
-      "x=[%.3e,%.3e,%.3e], "
-      "v=[%.3e,%.3e,%.3e],v_full=[%.3e,%.3e,%.3e] \n a=[%.3e,%.3e,%.3e], "
-      "u=%.3e, du/dt=%.3e v_sig=%.3e, P=%.3e\n"
-      "h=%.3e, dh/dt=%.3e wcount=%d, m=%.3e, dh_drho=%.3e, rho=%.3e, "
-      "time_bin=%d, mat_id=%d\n",
-      p->x[0], p->x[1], p->x[2], p->v[0], p->v[1], p->v[2], xp->v_full[0],
-      xp->v_full[1], xp->v_full[2], p->a_hydro[0], p->a_hydro[1], p->a_hydro[2],
-      p->u, p->u_dt, p->force.v_sig, hydro_get_comoving_pressure(p), p->h,
-      p->force.h_dt, (int)p->density.wcount, p->mass, p->density.rho_dh, p->rho,
+      "\n "
+      "x=[%.6g, %.6g, %.6g], v=[%.3g, %.3g, %.3g], \n "
+      "v_full=[%.3g, %.3g, %.3g], a=[%.3g, %.3g, %.3g], \n "
+      "m=%.3g, u=%.3g, du/dt=%.3g, P=%.3g, c_s=%.3g, \n "
+      "v_sig=%.3g, h=%.3g, dh/dt=%.3g, wcount=%.3g, rho=%.3g, \n "
+      "dh_drho=%.3g, time_bin=%d, mat_id=%d \n",
+      p->x[0], p->x[1], p->x[2], p->v[0], p->v[1], p->v[2],
+      xp->v_full[0], xp->v_full[1], xp->v_full[2],
+      p->a_hydro[0], p->a_hydro[1], p->a_hydro[2],
+      p->mass, p->u, p->u_dt, hydro_get_comoving_pressure(p),
+      p->force.soundspeed, p->force.v_sig, p->h, p->force.h_dt,
+      p->density.wcount, p->rho, p->density.rho_dh,
       p->time_bin, p->mat_id);
 }
 
-- 
GitLab


From fc43e3851ad03f027c6d66b79cfd376bc466a1e6 Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <matthieu.schaller@durham.ac.uk>
Date: Wed, 30 May 2018 09:32:19 +0200
Subject: [PATCH 5/6] Do not set the signal velocity to 1 in
 hydro_has_no_neighbours(). Code formatting.

---
 src/hydro/Minimal/hydro.h               | 1 -
 src/hydro/Minimal/hydro_debug.h         | 8 +++-----
 src/hydro/MinimalMultiMat/hydro.h       | 1 -
 src/hydro/MinimalMultiMat/hydro_debug.h | 8 +++-----
 4 files changed, 6 insertions(+), 12 deletions(-)

diff --git a/src/hydro/Minimal/hydro.h b/src/hydro/Minimal/hydro.h
index 8bbaacfbe..812f8ad72 100644
--- a/src/hydro/Minimal/hydro.h
+++ b/src/hydro/Minimal/hydro.h
@@ -370,7 +370,6 @@ __attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours(
   p->density.wcount = kernel_root * h_inv_dim;
   p->density.rho_dh = 0.f;
   p->density.wcount_dh = 0.f;
-  p->force.v_sig = 1.f;
 }
 
 /**
diff --git a/src/hydro/Minimal/hydro_debug.h b/src/hydro/Minimal/hydro_debug.h
index eb0c99dd4..73ffc26b8 100644
--- a/src/hydro/Minimal/hydro_debug.h
+++ b/src/hydro/Minimal/hydro_debug.h
@@ -42,13 +42,11 @@ __attribute__((always_inline)) INLINE static void hydro_debug_particle(
       "m=%.3g, u=%.3g, du/dt=%.3g, P=%.3g, c_s=%.3g, \n "
       "v_sig=%.3g, h=%.3g, dh/dt=%.3g, wcount=%.3g, rho=%.3g, \n "
       "dh_drho=%.3g, time_bin=%d \n",
-      p->x[0], p->x[1], p->x[2], p->v[0], p->v[1], p->v[2],
-      xp->v_full[0], xp->v_full[1], xp->v_full[2],
-      p->a_hydro[0], p->a_hydro[1], p->a_hydro[2],
+      p->x[0], p->x[1], p->x[2], p->v[0], p->v[1], p->v[2], xp->v_full[0],
+      xp->v_full[1], xp->v_full[2], p->a_hydro[0], p->a_hydro[1], p->a_hydro[2],
       p->mass, p->u, p->u_dt, hydro_get_comoving_pressure(p),
       p->force.soundspeed, p->force.v_sig, p->h, p->force.h_dt,
-      p->density.wcount, p->rho, p->density.rho_dh,
-      p->time_bin);
+      p->density.wcount, p->rho, p->density.rho_dh, p->time_bin);
 }
 
 #endif /* SWIFT_MINIMAL_HYDRO_DEBUG_H */
diff --git a/src/hydro/MinimalMultiMat/hydro.h b/src/hydro/MinimalMultiMat/hydro.h
index 405b4b5b0..cfad6b2b2 100644
--- a/src/hydro/MinimalMultiMat/hydro.h
+++ b/src/hydro/MinimalMultiMat/hydro.h
@@ -371,7 +371,6 @@ __attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours(
   p->density.wcount = kernel_root * h_inv_dim;
   p->density.rho_dh = 0.f;
   p->density.wcount_dh = 0.f;
-  p->force.v_sig = 1.f;
 }
 
 /**
diff --git a/src/hydro/MinimalMultiMat/hydro_debug.h b/src/hydro/MinimalMultiMat/hydro_debug.h
index 01d4d24b1..17b624ad0 100644
--- a/src/hydro/MinimalMultiMat/hydro_debug.h
+++ b/src/hydro/MinimalMultiMat/hydro_debug.h
@@ -44,13 +44,11 @@ __attribute__((always_inline)) INLINE static void hydro_debug_particle(
       "m=%.3g, u=%.3g, du/dt=%.3g, P=%.3g, c_s=%.3g, \n "
       "v_sig=%.3g, h=%.3g, dh/dt=%.3g, wcount=%.3g, rho=%.3g, \n "
       "dh_drho=%.3g, time_bin=%d, mat_id=%d \n",
-      p->x[0], p->x[1], p->x[2], p->v[0], p->v[1], p->v[2],
-      xp->v_full[0], xp->v_full[1], xp->v_full[2],
-      p->a_hydro[0], p->a_hydro[1], p->a_hydro[2],
+      p->x[0], p->x[1], p->x[2], p->v[0], p->v[1], p->v[2], xp->v_full[0],
+      xp->v_full[1], xp->v_full[2], p->a_hydro[0], p->a_hydro[1], p->a_hydro[2],
       p->mass, p->u, p->u_dt, hydro_get_comoving_pressure(p),
       p->force.soundspeed, p->force.v_sig, p->h, p->force.h_dt,
-      p->density.wcount, p->rho, p->density.rho_dh,
-      p->time_bin, p->mat_id);
+      p->density.wcount, p->rho, p->density.rho_dh, p->time_bin, p->mat_id);
 }
 
 #endif /* SWIFT_MINIMAL_MULTI_MAT_HYDRO_DEBUG_H */
-- 
GitLab


From f0329b64dca75e76e936f1df60d8edf221e46646 Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <matthieu.schaller@durham.ac.uk>
Date: Wed, 30 May 2018 09:40:37 +0200
Subject: [PATCH 6/6] Apply the same bug fix to the other hydro schemes.

---
 src/hydro/Default/hydro.h         | 4 ++--
 src/hydro/Gadget2/hydro.h         | 2 +-
 src/hydro/GizmoMFM/hydro.h        | 2 +-
 src/hydro/GizmoMFV/hydro.h        | 2 +-
 src/hydro/PressureEnergy/hydro.h  | 4 ++--
 src/hydro/PressureEntropy/hydro.h | 4 ++--
 src/hydro/Shadowswift/hydro.h     | 2 +-
 7 files changed, 10 insertions(+), 10 deletions(-)

diff --git a/src/hydro/Default/hydro.h b/src/hydro/Default/hydro.h
index b1a999b63..2c3a9c46f 100644
--- a/src/hydro/Default/hydro.h
+++ b/src/hydro/Default/hydro.h
@@ -339,7 +339,7 @@ __attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours(
 
   /* Re-set problematic values */
   p->rho = p->mass * kernel_root * h_inv_dim;
-  p->density.wcount = kernel_root * kernel_norm * h_inv_dim;
+  p->density.wcount = kernel_root * h_inv_dim;
   p->rho_dh = 0.f;
   p->density.wcount_dh = 0.f;
   p->density.div_v = 0.f;
@@ -423,7 +423,7 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration(
   /* Reset the time derivatives. */
   p->force.u_dt = 0.0f;
   p->force.h_dt = 0.0f;
-  p->force.v_sig = 0.0f;
+  p->force.v_sig = p->force.soundspeed;
 }
 
 /**
diff --git a/src/hydro/Gadget2/hydro.h b/src/hydro/Gadget2/hydro.h
index bc06a24e2..26e3bf97d 100644
--- a/src/hydro/Gadget2/hydro.h
+++ b/src/hydro/Gadget2/hydro.h
@@ -349,7 +349,7 @@ __attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours(
 
   /* Re-set problematic values */
   p->rho = p->mass * kernel_root * h_inv_dim;
-  p->density.wcount = kernel_root * kernel_norm * h_inv_dim;
+  p->density.wcount = kernel_root * h_inv_dim;
   p->density.rho_dh = 0.f;
   p->density.wcount_dh = 0.f;
   p->density.div_v = 0.f;
diff --git a/src/hydro/GizmoMFM/hydro.h b/src/hydro/GizmoMFM/hydro.h
index 9c4be6af3..ae0fe13bc 100644
--- a/src/hydro/GizmoMFM/hydro.h
+++ b/src/hydro/GizmoMFM/hydro.h
@@ -398,7 +398,7 @@ __attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours(
   const float h_inv_dim = pow_dimension(h_inv); /* 1/h^d */
 
   /* Re-set problematic values */
-  p->density.wcount = kernel_root * kernel_norm * h_inv_dim;
+  p->density.wcount = kernel_root * h_inv_dim;
   p->density.wcount_dh = 0.f;
   p->geometry.volume = 1.0f;
   p->geometry.matrix_E[0][0] = 1.0f;
diff --git a/src/hydro/GizmoMFV/hydro.h b/src/hydro/GizmoMFV/hydro.h
index 1d5abeaaf..e94a87408 100644
--- a/src/hydro/GizmoMFV/hydro.h
+++ b/src/hydro/GizmoMFV/hydro.h
@@ -398,7 +398,7 @@ __attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours(
   const float h_inv_dim = pow_dimension(h_inv); /* 1/h^d */
 
   /* Re-set problematic values */
-  p->density.wcount = kernel_root * kernel_norm * h_inv_dim;
+  p->density.wcount = kernel_root * h_inv_dim;
   p->density.wcount_dh = 0.f;
   p->geometry.volume = 1.0f;
   p->geometry.matrix_E[0][0] = 1.0f;
diff --git a/src/hydro/PressureEnergy/hydro.h b/src/hydro/PressureEnergy/hydro.h
index 1e0d8208e..ea086daee 100644
--- a/src/hydro/PressureEnergy/hydro.h
+++ b/src/hydro/PressureEnergy/hydro.h
@@ -403,7 +403,7 @@ __attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours(
   p->rho = p->mass * kernel_root * h_inv_dim;
   p->pressure_bar =
       p->mass * p->u * hydro_gamma_minus_one * kernel_root * h_inv_dim;
-  p->density.wcount = kernel_root * kernel_norm * h_inv_dim;
+  p->density.wcount = kernel_root * h_inv_dim;
   p->density.rho_dh = 0.f;
   p->density.wcount_dh = 0.f;
   p->density.pressure_bar_dh = 0.f;
@@ -480,7 +480,7 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration(
   /* Reset the time derivatives. */
   p->u_dt = 0.0f;
   p->force.h_dt = 0.0f;
-  p->force.v_sig = 0.0f;
+  p->force.v_sig = p->force.soundspeed;
 }
 
 /**
diff --git a/src/hydro/PressureEntropy/hydro.h b/src/hydro/PressureEntropy/hydro.h
index 87d46c6d4..e4b7cf06e 100644
--- a/src/hydro/PressureEntropy/hydro.h
+++ b/src/hydro/PressureEntropy/hydro.h
@@ -357,7 +357,7 @@ __attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours(
   /* Re-set problematic values */
   p->rho = p->mass * kernel_root * h_inv_dim;
   p->rho_bar = p->mass * kernel_root * h_inv_dim;
-  p->density.wcount = kernel_root * kernel_norm * h_inv_dim;
+  p->density.wcount = kernel_root * h_inv_dim;
   p->density.rho_dh = 0.f;
   p->density.wcount_dh = 0.f;
   p->density.pressure_dh = 0.f;
@@ -441,7 +441,7 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration(
   p->force.h_dt = 0.0f;
 
   /* Reset maximal signal velocity */
-  p->force.v_sig = 0.0f;
+  p->force.v_sig = p->force.soundspeed;
 }
 
 /**
diff --git a/src/hydro/Shadowswift/hydro.h b/src/hydro/Shadowswift/hydro.h
index 36078798c..025779f17 100644
--- a/src/hydro/Shadowswift/hydro.h
+++ b/src/hydro/Shadowswift/hydro.h
@@ -257,7 +257,7 @@ __attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours(
   const float h_inv_dim = pow_dimension(h_inv); /* 1/h^d */
 
   /* Re-set problematic values */
-  p->density.wcount = kernel_root * kernel_norm * h_inv_dim;
+  p->density.wcount = kernel_root * h_inv_dim;
   p->density.wcount_dh = 0.f;
 }
 
-- 
GitLab