Done !

src/hydro/Gadget2/hydro.h

@@ -31,7 +31,7 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
   float ac = sqrtf(p->a[0]*p->a[0] + p->a[1]*p->a[1] + p->a[2]*p->a[2]);
   ac = fmaxf(ac, 1e-30);
 
-  const float dt_accel = sqrtf(2.f);
+  const float dt_accel = sqrtf(2.f); //MATTHIEU
   
   /* CFL condition */
   const float dt_cfl = 2.f * const_cfl * p->h / p->v_sig;
@@ -214,6 +214,27 @@ __attribute__((always_inline)) INLINE static void hydro_end_force(struct part* p
   
 }
 
+/**
+ * @brief Kick the additional variables
+ *
+ * @param p The particle to act upon
+ */
+__attribute__((always_inline)) INLINE static void hydro_kick_extra(struct part* p, float dt) {
+
+  /* Do not decrease the entropy (temperature) by more than a factor of 2*/
+  const float entropy_change = p->entropy_dt * dt;
+  if( entropy_change > -0.5f * p->entropy)
+    p->entropy += entropy_change;
+  else
+    p->entropy *= 0.5f;
+
+  /* Do not 'overcool' when timestep increases */
+  if(p->entropy + 0.5f * p->entropy_dt * dt < 0.5f * p->entropy)
+    p->entropy_dt = -0.5f * p->entropy / dt;
+  
+}
+
+
 /**
  * @brief Converts hydro quantity of a particle 
  *

src/hydro/Gadget2/hydro_iact.h

@@ -158,8 +158,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
   const float dvdr = dv[0] * dx[0] + dv[1] * dx[1] + dv[2] * dx[2];
   pi->div_v -= fac * dvdr;
 
-  if(pi->id == 1000 && pj->id == 1203)
-    message("Interacting with %lld. r=%e hi=%e u=%e W=%e dW/dx=%e dh_drho1=%e dh_drho2=%e\n fac=%e dvdr=%e",
+  if(pi->id == 515050 && pj->id == 504849)
+    message("Interacting with %lld. r=%e hi=%e u=%e W=%e dW/dx=%e dh_drho1=%e dh_drho2=%e\n fac=%e dvdr=%e pj->v=[%.3e,%.3e,%.3e]",
 	    pj->id,
 	    r,
 	    hi,
@@ -169,7 +169,10 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
 	    -mj * (3.f * kernel_igamma * wi) * ih4,
 	    -mj * u * wi_dx * kernel_igamma * ih4,
 	    fac * ih4,
-	    dvdr
+	    dvdr,
+	    pj->v[0],
+	    pj->v[1],
+	    pj->v[2]
 	    );
 
   
@@ -252,7 +255,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
   const float acc = visc_term + sph_term;
 
 
-  if(pi->id == 1000 && pj->id == 1203)
+  //if(pi->id == 1000 && pj->id == 1100)
+  if(pi->id == 515050 && pj->id == 504849)
     message("Interacting with %lld. r=%e hi=%e hj=%e dWi/dx=%e dWj/dx=%3e dvdr=%e visc=%e sph=%e",
 	    pj->id,
 	    r,
@@ -264,7 +268,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
 	    visc_term,
 	    sph_term
 	    );
-  if(pi->id == 1203 && pj->id == 1000)
+  if(pi->id == 1100 && pj->id == 1000)
     message("oO");
 
   
@@ -282,8 +286,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
   pj->v_sig = fmaxf(pj->v_sig, v_sig) ;
   
   /* Change in entropy */
-  pi->entropy_dt -= 0.5f * visc_term * dvdr;
-  pj->entropy_dt += 0.5f * visc_term * dvdr;
+  pi->entropy_dt += 0.5f * visc_term * dvdr;
+  pj->entropy_dt -= 0.5f * visc_term * dvdr;
 }
 
 
@@ -363,7 +367,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
   pi->v_sig = fmaxf(pi->v_sig, v_sig) ;
   
   /* Change in entropy */
-  pi->entropy_dt -= 0.5f * visc_term * dvdr;
+  pi->entropy_dt += 0.5f * visc_term * dvdr;
 }
 
 

src/runner.c

@@ -758,6 +758,17 @@ void runner_dodrift(struct runner *r, struct cell *c, int timer) {
       /* Predict the values of the extra fields */
       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",
+		p->id,
+		r->e->time,
+		r->e->timeOld,
+		r->e->time,
+		p->v[0],
+		p->v[1],
+		p->v[2]);
+
+      
       /* Compute motion since last cell construction */
       const float dx = sqrtf((p->x[0] - xp->x_old[0]) * (p->x[0] - xp->x_old[0]) +
 			     (p->x[1] - xp->x_old[1]) * (p->x[1] - xp->x_old[1]) +