From d1dfe28b52940dc70e47d374d8979e2991ab1b4b Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <matthieu.schaller@durham.ac.uk>
Date: Tue, 9 Feb 2016 23:28:19 +0000
Subject: [PATCH] Correct smoothing length time derivative in all cases

---
 src/hydro/Gadget2/hydro.h      | 4 +++-
 src/hydro/Gadget2/hydro_iact.h | 6 +++---
 src/runner.c                   | 7 ++++---
 3 files changed, 10 insertions(+), 7 deletions(-)

diff --git a/src/hydro/Gadget2/hydro.h b/src/hydro/Gadget2/hydro.h
index 30d714e253..40b0a93715 100644
--- a/src/hydro/Gadget2/hydro.h
+++ b/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));
 }
diff --git a/src/hydro/Gadget2/hydro_iact.h b/src/hydro/Gadget2/hydro_iact.h
index 260f78ba0e..d8e5559400 100644
--- a/src/hydro/Gadget2/hydro_iact.h
+++ b/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);
diff --git a/src/runner.c b/src/runner.c
index eeee039344..88a94f0dcd 100644
--- a/src/runner.c
+++ b/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 =
-- 
GitLab