Correct smoothing length time derivative in all cases

src/hydro/Gadget2/hydro.h

@@ -168,6 +168,8 @@ __attribute__((always_inline))
   p->a[1] = 0.0f;
   p->a[2] = 0.0f;
 
+  p->force.h_dt = 0.0f;
+
   /* Reset the time derivatives. */
   p->entropy_dt = 0.0f;
 
@@ -203,7 +205,7 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
  */
 __attribute__((always_inline))
     INLINE static void hydro_end_force(struct part* p) {
-  
+
   p->entropy_dt *=
       (const_hydro_gamma - 1.f) * powf(p->rho, -(const_hydro_gamma - 1.f));
 }

src/hydro/Gadget2/hydro_iact.h

@@ -281,8 +281,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
   pj->a[2] += acc * dx[2];
 
   /* Get the time derivative for h. */
-  pi->force.h_dt -= mj * dvdr / rhoj * wi_dr;
-  pj->force.h_dt -= mi * dvdr / 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->v_sig = fmaxf(pi->v_sig, v_sig);
@@ -371,7 +371,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
   pi->a[2] -= acc * dx[2];
 
   /* Get the time derivative for h. */
-  pi->force.h_dt -= mj * dvdr / rhoj * wi_dr;
+  pi->force.h_dt -= mj * dvdr * r_inv / rhoj * wi_dr;
   
   /* Update the signal velocity. */
   pi->v_sig = fmaxf(pi->v_sig, v_sig);

src/runner.c

@@ -751,15 +751,16 @@ void runner_dodrift(struct runner *r, struct cell *c, int timer) {
       hydro_predict_extra(p, xp, r->e->timeOld, r->e->time);
 
       /* if(p->id == 1000 || p->id == 515050 || p->id == 504849) */
-      /* 	message("%lld: current_t=%f t0=%f t1=%f  v=[%.3e %.3e %.3e]\n",
-       */
+      /* 	message("%lld: current_t=%f t0=%f t1=%f  v=[%.3e %.3e %.3e] dh/dt=%.3e div_v=%.3e\n", */
       /* 		p->id, */
       /* 		r->e->time, */
       /* 		r->e->timeOld, */
       /* 		r->e->time, */
       /* 		p->v[0], */
       /* 		p->v[1], */
-      /* 		p->v[2]); */
+      /* 		p->v[2], */
+      /* 		p->force.h_dt * h_inv * dt, */
+      /* 		0.333333f * p->div_v * dt); */
 
       /* Compute motion since last cell construction */
       const float dx =