Re-introduced rho_dh in the Pressure-Energy SPH hydro case to allow for the...

Re-introduced rho_dh in the Pressure-Energy SPH hydro case to allow for the alternative neighbour number definition.

src/hydro/PressureEnergy/hydro.h

@@ -517,6 +517,7 @@ __attribute__((always_inline)) INLINE static void hydro_init_part(
   p->density.wcount = 0.f;
   p->density.wcount_dh = 0.f;
   p->rho = 0.f;
+  p->density.rho_dh = 0.f;
   p->pressure_bar = 0.f;
   p->density.pressure_bar_dh = 0.f;
 
@@ -550,6 +551,7 @@ __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->pressure_bar += p->mass * p->u * kernel_root;
   p->density.pressure_bar_dh -= hydro_dimension * p->mass * p->u * kernel_root;
   p->density.wcount += kernel_root;
@@ -557,6 +559,7 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
 
   /* Finish the calculation by inserting the missing h-factors */
   p->rho *= h_inv_dim;
+  p->density.rho_dh *= h_inv_dim_plus_one;
   p->pressure_bar *= (h_inv_dim * hydro_gamma_minus_one);
   p->density.pressure_bar_dh *= (h_inv_dim_plus_one * hydro_gamma_minus_one);
   p->density.wcount *= h_inv_dim;
@@ -601,6 +604,7 @@ __attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours(
   p->pressure_bar =
       p->mass * p->u * hydro_gamma_minus_one * kernel_root * 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;
 

src/hydro/PressureEnergy/hydro_iact.h

@@ -69,6 +69,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
   kernel_deval(ui, &wi, &wi_dx);
 
   pi->rho += mj * wi;
+  pi->density.rho_dh -= mj * (hydro_dimension * wi + ui * wi_dx);
 
   pi->pressure_bar += mj * wi * pj->u;
   pi->density.pressure_bar_dh -=
@@ -82,6 +83,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
   kernel_deval(uj, &wj, &wj_dx);
 
   pj->rho += mi * wj;
+  pj->density.rho_dh -= mi * (hydro_dimension * wj + uj * wj_dx);
+
   pj->pressure_bar += mi * wj * pi->u;
   pj->density.pressure_bar_dh -=
       mi * pi->u * (hydro_dimension * wj + uj * wj_dx);
@@ -147,6 +150,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
   kernel_deval(ui, &wi, &wi_dx);
 
   pi->rho += mj * wi;
+  pi->density.rho_dh -= mj * (hydro_dimension * wi + ui * wi_dx);
 
   pi->pressure_bar += mj * wi * pj->u;
 

src/hydro/PressureEnergy/hydro_part.h

@@ -133,6 +133,9 @@ struct part {
       /*! Derivative of the neighbour number with respect to h. */
       float wcount_dh;
 
+      /*! Derivative of density with respect to h */
+      float rho_dh;
+
       /*! Derivative of the weighted pressure with respect to h */
       float pressure_bar_dh;