Also modified the GIZMO scheme.

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.
  *

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;

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.
  *