Moved h_dt into the force sub-structure

src/hydro/Default/hydro.h

@@ -186,7 +186,7 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration(
 
   /* Reset the time derivatives. */
   p->force.u_dt = 0.0f;
-  p->h_dt = 0.0f;
+  p->force.h_dt = 0.0f;
   p->force.v_sig = 0.0f;
 }
 
@@ -225,7 +225,9 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
  * @param p The particle to act upon
  */
 __attribute__((always_inline)) INLINE static void hydro_end_force(
-    struct part *restrict p) {}
+    struct part *restrict p) {
+  p->force.h_dt *= p->h * 0.333333333f;
+}
 
 /**
  * @brief Kick the additional variables

src/hydro/Default/hydro_iact.h

@@ -437,8 +437,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
   pj->force.u_dt += mi * tc * (pj->u - pi->u);
 
   /* Get the time derivative for h. */
-  pi->h_dt -= mj * dvdr / rhoj * wi_dr;
-  pj->h_dt -= mi * dvdr / rhoi * wj_dr;
+  pi->force.h_dt -= mj * dvdr / rhoj * wi_dr;
+  pj->force.h_dt -= mi * dvdr / rhoi * wj_dr;
 
   /* Update the signal velocity. */
   pi->force.v_sig = fmaxf(pi->force.v_sig, v_sig);
@@ -645,8 +645,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_vec_force(
   for (k = 0; k < VEC_SIZE; k++) {
     pi[k]->force.u_dt += piu_dt.f[k];
     pj[k]->force.u_dt += pju_dt.f[k];
-    pi[k]->h_dt -= pih_dt.f[k];
-    pj[k]->h_dt -= pjh_dt.f[k];
+    pi[k]->force.h_dt -= pih_dt.f[k];
+    pj[k]->force.h_dt -= pjh_dt.f[k];
     pi[k]->force.v_sig = vi_sig.f[k];
     pj[k]->force.v_sig = vj_sig.f[k];
     for (j = 0; j < 3; j++) {
@@ -746,7 +746,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
   pi->force.u_dt += mj * tc * (pi->u - pj->u);
 
   /* Get the time derivative for h. */
-  pi->h_dt -= mj * dvdr / rhoj * wi_dr;
+  pi->force.h_dt -= mj * dvdr / rhoj * wi_dr;
 
   /* Update the signal velocity. */
   pi->force.v_sig = fmaxf(pi->force.v_sig, v_sig);
@@ -944,7 +944,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_vec_force(
   /* Store the forces back on the particles. */
   for (k = 0; k < VEC_SIZE; k++) {
     pi[k]->force.u_dt += piu_dt.f[k];
-    pi[k]->h_dt -= pih_dt.f[k];
+    pi[k]->force.h_dt -= pih_dt.f[k];
     pi[k]->force.v_sig = vi_sig.f[k];
     pj[k]->force.v_sig = vj_sig.f[k];
     for (j = 0; j < 3; j++) pi[k]->a_hydro[j] -= pia[j].f[k];

src/hydro/Default/hydro_part.h

@@ -46,9 +46,6 @@ struct part {
   /* Particle cutoff radius. */
   float h;
 
-  /* Change in smoothing length over time. */
-  float h_dt;
-
   /* Particle time of beginning of time-step. */
   int ti_begin;
 
@@ -104,6 +101,9 @@ struct part {
       /* Sound speed */
       float c;
 
+      /* Change in smoothing length over time. */
+      float h_dt;
+
     } force;
   };
 

src/hydro/Gadget2/hydro.h

@@ -204,7 +204,7 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
       (p->entropy + p->force.entropy_dt * dt_entr) * pow_gamma(p->rho);
 
   const float irho = 1.f / p->rho;
-  
+
   /* Divide the pressure by the density and density gradient */
   const float P_over_rho2 = pressure * irho * irho * p->rho_dh;
 

src/hydro/Gadget2/hydro_iact.h

@@ -445,7 +445,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
 
   /* Now, convolve with the kernel */
   const float visc_term = 0.5f * visc * (wi_dr + wj_dr) * r_inv;
-  const float sph_term = (P_over_rho2_i * wi_dr + P_over_rho2_j * wj_dr) * r_inv;
+  const float sph_term =
+      (P_over_rho2_i * wi_dr + P_over_rho2_j * wj_dr) * r_inv;
 
   /* Eventually got the acceleration */
   const float acc = visc_term + sph_term;
@@ -547,7 +548,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
 
   /* Now, convolve with the kernel */
   const float visc_term = 0.5f * visc * (wi_dr + wj_dr) * r_inv;
-  const float sph_term = (P_over_rho2_i * wi_dr + P_over_rho2_j * wj_dr) * r_inv;
+  const float sph_term =
+      (P_over_rho2_i * wi_dr + P_over_rho2_j * wj_dr) * r_inv;
 
   /* Eventually got the acceleration */
   const float acc = visc_term + sph_term;

src/hydro/Gadget2/hydro_part.h

@@ -45,7 +45,7 @@ struct part {
 
   /* Particle mass. */
   float mass;
-  
+
   /* Particle time of beginning of time-step. */
   int ti_begin;
 
@@ -92,7 +92,7 @@ struct part {
 
       /* Signal velocity. */
       float v_sig;
-      
+
       /* Entropy time derivative */
       float entropy_dt;
 

src/hydro/Minimal/hydro.h

@@ -155,7 +155,7 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration(
 
   /* Reset the time derivatives. */
   p->u_dt = 0.0f;
-  p->h_dt = 0.0f;
+  p->force.h_dt = 0.0f;
   p->force.v_sig = 0.0f;
 }
 
@@ -191,7 +191,10 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
  * @param p The particle to act upon
  */
 __attribute__((always_inline)) INLINE static void hydro_end_force(
-    struct part *restrict p) {}
+    struct part *restrict p) {
+
+  p->force.h_dt *= p->h * 0.333333333f;
+}
 
 /**
  * @brief Kick the additional variables

src/hydro/Minimal/hydro_debug.h

@@ -28,6 +28,6 @@ __attribute__((always_inline)) INLINE static void hydro_debug_particle(
       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],
       xp->u_full, p->u, p->u_dt, p->force.v_sig, p->force.pressure, p->h,
-      p->h_dt, (int)p->density.wcount, p->mass, p->rho_dh, p->rho, p->ti_begin,
-      p->ti_end);
+      p->force.h_dt, (int)p->density.wcount, p->mass, p->rho_dh, p->rho,
+      p->ti_begin, p->ti_end);
 }

src/hydro/Minimal/hydro_iact.h

@@ -173,8 +173,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
   pj->u_dt += P_over_rho_j * mi * dvdr * r_inv * wj_dr;
 
   /* Get the time derivative for h. */
-  pi->h_dt -= mj * dvdr * r_inv / rhoj * wi_dr;
-  pj->h_dt -= mi * dvdr * r_inv / rhoi * wj_dr;
+  pi->force.h_dt -= mj * dvdr * r_inv / rhoj * wi_dr;
+  pj->force.h_dt -= mi * dvdr * r_inv / rhoi * wj_dr;
 
   /* Update the signal velocity. */
   pi->force.v_sig = fmaxf(pi->force.v_sig, v_sig);
@@ -251,7 +251,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
   pi->u_dt += P_over_rho_i * mj * dvdr * r_inv * wi_dr;
 
   /* Get the time derivative for h. */
-  pi->h_dt -= mj * dvdr * r_inv / rhoj * wi_dr;
+  pi->force.h_dt -= mj * dvdr * r_inv / rhoj * wi_dr;
 
   /* Update the signal velocity. */
   pi->force.v_sig = fmaxf(pi->force.v_sig, v_sig);

src/hydro/Minimal/hydro_part.h

@@ -54,8 +54,6 @@ struct part {
 
   float h; /*!< Particle smoothing length. */
 
-  float h_dt; /*!< Time derivative of smoothing length  */
-
   int ti_begin; /*!< Time at the beginning of time-step. */
 
   int ti_end; /*!< Time at the end of time-step. */
@@ -91,7 +89,7 @@ struct part {
      * neighbours.
      *
      * Quantities in this sub-structure should only be accessed in the force
-     * loop over neighbours and the ghost and kick tasks.
+     * loop over neighbours and the ghost, drift and kick tasks.
      */
     struct {
 
@@ -99,6 +97,8 @@ struct part {
 
       float v_sig; /*!< Particle signal velocity */
 
+      float h_dt; /*!< Time derivative of smoothing length  */
+
     } force;
   };