diff --git a/src/hydro/AnarchyDU/hydro.h b/src/hydro/AnarchyDU/hydro.h
index 2cce0b0336f43acedf20b12867bfc31afa1cec6c..ff8eac12510cce6f5e5a578fbf91479aa67b3f5c 100644
--- a/src/hydro/AnarchyDU/hydro.h
+++ b/src/hydro/AnarchyDU/hydro.h
@@ -501,7 +501,7 @@ __attribute__((always_inline)) INLINE static void hydro_remove_part(
__attribute__((always_inline)) INLINE static void hydro_init_part(
struct part *restrict p, const struct hydro_space *hs) {
- p->density.wcount = 1.f;
+ p->density.wcount = 0.f;
p->density.wcount_dh = 0.f;
p->rho = 0.f;
p->density.rho_dh = 0.f;
@@ -539,17 +539,14 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
/* Final operation on the density (add self-contribution). */
p->rho += p->mass * kernel_root;
p->density.rho_dh -= hydro_dimension * p->mass * kernel_root;
- // p->density.wcount += kernel_root;
- // p->density.wcount_dh -= hydro_dimension * kernel_root;
+ p->density.wcount += kernel_root;
+ p->density.wcount_dh -= hydro_dimension * kernel_root;
/* Finish the calculation by inserting the missing h-factors */
p->rho *= h_inv_dim;
p->density.rho_dh *= h_inv_dim_plus_one;
- // p->density.wcount *= h_inv_dim;
- // p->density.wcount_dh *= h_inv_dim_plus_one;
-
- p->density.wcount = p->rho / p->mass;
- p->density.wcount_dh = p->density.rho_dh / p->mass;
+ p->density.wcount *= h_inv_dim;
+ p->density.wcount_dh *= h_inv_dim_plus_one;
const float rho_inv = 1.f / p->rho;
const float a_inv2 = cosmo->a2_inv;
diff --git a/src/hydro/AnarchyDU/hydro_iact.h b/src/hydro/AnarchyDU/hydro_iact.h
index 49106e4ee89e9c31e8fdb06bad406a264d99fcbb..09ecde01d44cfef1c0edbfe0d3f1aa909d553a95 100644
--- a/src/hydro/AnarchyDU/hydro_iact.h
+++ b/src/hydro/AnarchyDU/hydro_iact.h
@@ -65,8 +65,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
pi->rho += mj * wi;
pi->density.rho_dh -= mj * (hydro_dimension * wi + ui * wi_dx);
- // pi->density.wcount += wi;
- // pi->density.wcount_dh -= (hydro_dimension * wi + ui * wi_dx);
+
+ pi->density.wcount += wi;
+ pi->density.wcount_dh -= (hydro_dimension * wi + ui * wi_dx);
/* Compute density of pj. */
const float hj_inv = 1.f / hj;
@@ -75,8 +76,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
pj->rho += mi * wj;
pj->density.rho_dh -= mi * (hydro_dimension * wj + uj * wj_dx);
- // pj->density.wcount += wj;
- // pj->density.wcount_dh -= (hydro_dimension * wj + uj * wj_dx);
+ pj->density.wcount += wj;
+ pj->density.wcount_dh -= (hydro_dimension * wj + uj * wj_dx);
/* Now we need to compute the div terms */
const float r_inv = 1.f / r;
@@ -138,8 +139,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
pi->rho += mj * wi;
pi->density.rho_dh -= mj * (hydro_dimension * wi + ui * wi_dx);
- // pi->density.wcount += wi;
- // pi->density.wcount_dh -= (hydro_dimension * wi + ui * wi_dx);
+
+ pi->density.wcount += wi;
+ pi->density.wcount_dh -= (hydro_dimension * wi + ui * wi_dx);
const float r_inv = 1.f / r;
const float faci = mj * wi_dx * r_inv;
diff --git a/src/hydro/Minimal/hydro.h b/src/hydro/Minimal/hydro.h
index bcdb5306c58859d991cf8cb7fdbe01c31b0e9461..cdfbdbf3bcf3790b8cba84ce74084d3885cf75bf 100644
--- a/src/hydro/Minimal/hydro.h
+++ b/src/hydro/Minimal/hydro.h
@@ -470,7 +470,7 @@ __attribute__((always_inline)) INLINE static void hydro_timestep_extra(
__attribute__((always_inline)) INLINE static void hydro_init_part(
struct part *restrict p, const struct hydro_space *hs) {
- p->density.wcount = 1.f;
+ p->density.wcount = 0.f;
p->density.wcount_dh = 0.f;
p->rho = 0.f;
p->density.rho_dh = 0.f;
@@ -505,17 +505,14 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
/* Final operation on the density (add self-contribution). */
p->rho += p->mass * kernel_root;
p->density.rho_dh -= hydro_dimension * p->mass * kernel_root;
- // p->density.wcount += kernel_root;
- // p->density.wcount_dh -= hydro_dimension * kernel_root;
+ p->density.wcount += kernel_root;
+ p->density.wcount_dh -= hydro_dimension * kernel_root;
/* Finish the calculation by inserting the missing h-factors */
p->rho *= h_inv_dim;
p->density.rho_dh *= h_inv_dim_plus_one;
- // p->density.wcount *= h_inv_dim;
- // p->density.wcount_dh *= h_inv_dim_plus_one;
-
- p->density.wcount = p->rho / p->mass;
- p->density.wcount_dh = p->density.rho_dh / p->mass;
+ p->density.wcount *= h_inv_dim;
+ p->density.wcount_dh *= h_inv_dim_plus_one;
const float rho_inv = 1.f / p->rho;
const float a_inv2 = cosmo->a2_inv;
diff --git a/src/hydro/Minimal/hydro_iact.h b/src/hydro/Minimal/hydro_iact.h
index e2b41852e71053066b903c09f3e91572d8707e2f..8929ede0ac8ef736776e745323d33bb0b2265021 100644
--- a/src/hydro/Minimal/hydro_iact.h
+++ b/src/hydro/Minimal/hydro_iact.h
@@ -77,8 +77,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
pi->rho += mj * wi;
pi->density.rho_dh -= mj * (hydro_dimension * wi + ui * wi_dx);
- // pi->density.wcount += wi;
- // pi->density.wcount_dh -= (hydro_dimension * wi + ui * wi_dx);
+ pi->density.wcount += wi;
+ pi->density.wcount_dh -= (hydro_dimension * wi + ui * wi_dx);
/* Compute density of pj. */
const float hj_inv = 1.f / hj;
@@ -87,8 +87,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
pj->rho += mi * wj;
pj->density.rho_dh -= mi * (hydro_dimension * wj + uj * wj_dx);
- // pj->density.wcount += wj;
- // pj->density.wcount_dh -= (hydro_dimension * wj + uj * wj_dx);
+ pj->density.wcount += wj;
+ pj->density.wcount_dh -= (hydro_dimension * wj + uj * wj_dx);
/* Compute dv dot r */
float dv[3], curlvr[3];
@@ -156,8 +156,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
pi->rho += mj * wi;
pi->density.rho_dh -= mj * (hydro_dimension * wi + ui * wi_dx);
- // pi->density.wcount += wi;
- // pi->density.wcount_dh -= (hydro_dimension * wi + ui * wi_dx);
+ pi->density.wcount += wi;
+ pi->density.wcount_dh -= (hydro_dimension * wi + ui * wi_dx);
/* Compute dv dot r */
float dv[3], curlvr[3];