Moved h_dt into the force sub-structure

src/drift.h

@@ -77,7 +77,7 @@ __attribute__((always_inline)) INLINE static void drift_part(
   p->v[2] += p->a_hydro[2] * dt;
 
   /* Predict smoothing length */
-  const float w1 = p->h_dt * h_inv * dt;
+  const float w1 = p->force.h_dt * h_inv * dt;
   if (fabsf(w1) < 0.2f)
     p->h *= approx_expf(w1); /* 4th order expansion of exp(w) */
   else

src/hydro/Gadget2/hydro.h

@@ -177,10 +177,9 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration(
   p->a_hydro[1] = 0.0f;
   p->a_hydro[2] = 0.0f;
 
-  p->h_dt = 0.0f;
-
   /* Reset the time derivatives. */
   p->force.entropy_dt = 0.0f;
+  p->force.h_dt = 0.0f;
 
   /* Reset maximal signal velocity */
   p->force.v_sig = 0.0f;
@@ -227,6 +226,8 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
 __attribute__((always_inline)) INLINE static void hydro_end_force(
     struct part *restrict p) {
 
+  p->force.h_dt *= p->h * 0.333333333f;
+
   p->force.entropy_dt *=
       hydro_gamma_minus_one * pow_minus_gamma_minus_one(p->rho);
 }

src/hydro/Gadget2/hydro_debug.h

@@ -35,5 +35,5 @@ __attribute__((always_inline)) INLINE static void hydro_debug_particle(
       p->rho, p->force.P_over_rho2, p->entropy, p->force.entropy_dt,
       p->force.soundspeed, p->density.div_v, p->density.rot_v[0],
       p->density.rot_v[1], p->density.rot_v[2], p->force.balsara,
-      p->force.v_sig, p->h_dt, p->ti_begin, p->ti_end);
+      p->force.v_sig, p->force.h_dt, p->ti_begin, p->ti_end);
 }

src/hydro/Gadget2/hydro_iact.h

@@ -460,8 +460,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
   pj->a_hydro[2] += mi * acc * dx[2];
 
   /* 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);
@@ -558,7 +558,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
   pi->a_hydro[2] -= mj * acc * dx[2];
 
   /* 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/Gadget2/hydro_part.h

@@ -43,9 +43,9 @@ struct part {
   /* Particle cutoff radius. */
   float h;
 
-  /* Time derivative of the smoothing length */
-  float h_dt;
-
+  /* Particle mass. */
+  float mass;
+  
   /* Particle time of beginning of time-step. */
   int ti_begin;
 
@@ -58,13 +58,9 @@ struct part {
   /* Particle entropy. */
   float entropy;
 
-  /* Derivative of the density with respect to this particle's smoothing length.
-   */
+  /* Derivative of the density with respect to  smoothing length. */
   float rho_dh;
 
-  /* Particle mass. */
-  float mass;
-
   union {
 
     struct {
@@ -88,7 +84,7 @@ struct part {
       /* Balsara switch */
       float balsara;
 
-      /* Pressure over density*/
+      /* Pressure over density squared (including drho/dh term) */
       float P_over_rho2;
 
       /* Particle sound speed. */
@@ -100,6 +96,9 @@ struct part {
       /* Entropy time derivative */
       float entropy_dt;
 
+      /* Time derivative of the smoothing length */
+      float h_dt;
+
     } force;
   };
 

src/runner.c

@@ -787,9 +787,6 @@ void runner_do_kick_fixdt(struct runner *r, struct cell *c, int timer) {
       struct xpart *restrict xp = &xparts[k];
 
       /* First, finish the force loop */
-      p->h_dt *= p->h * 0.333333333f;
-
-      /* And do the same of the extra variable */
       hydro_end_force(p);
       if (p->gpart != NULL) gravity_end_force(p->gpart);
 
@@ -908,9 +905,6 @@ void runner_do_kick(struct runner *r, struct cell *c, int timer) {
       if (p->ti_end <= ti_current) {
 
         /* First, finish the force loop */
-        p->h_dt *= p->h * 0.333333333f;
-
-        /* And do the same of the extra variable */
         hydro_end_force(p);
         if (p->gpart != NULL) gravity_end_force(p->gpart);
 

src/timestep.h

@@ -117,8 +117,8 @@ __attribute__((always_inline)) INLINE static int get_part_timestep(
 
   /* Limit change in h */
   const float dt_h_change =
-      (p->h_dt != 0.0f)
-          ? fabsf(e->hydro_properties->log_max_h_change * p->h / p->h_dt)
+      (p->force.h_dt != 0.0f)
+          ? fabsf(e->hydro_properties->log_max_h_change * p->h / p->force.h_dt)
           : FLT_MAX;
 
   new_dt = fminf(new_dt, dt_h_change);