Drift is now Gadget-2 identical

src/engine.c

@@ -1737,8 +1737,6 @@ void engine_init_particles(struct engine *e) {
   printParticle(e->s->parts, 1000, e->s->nr_parts);
   printParticle(e->s->parts, 515050, e->s->nr_parts);
   
-  exit(0);
-  
   /* Ready to go */
   e->step = -1;
 }
@@ -1812,9 +1810,6 @@ if ( e->nodeID == 0 )
   printf("%d %f %f %d\n", e->step, e->time, e->timeStep, updates);
   fflush(stdout);
 
-  printParticle(e->s->parts, 1000, e->s->nr_parts);
-  printParticle(e->s->parts, 515050, e->s->nr_parts);
-
   message("\nDRIFT\n");
   
   /* Drift everybody */
@@ -1823,6 +1818,8 @@ if ( e->nodeID == 0 )
   printParticle(e->s->parts, 1000, e->s->nr_parts);
   printParticle(e->s->parts, 515050, e->s->nr_parts);
 
+  exit(0);
+  
   /* Move forward in time */
   e->timeOld = e->time;
   e->time = t_end_min;

src/hydro/Gadget2/hydro.h

@@ -183,12 +183,22 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration(struc
  *
  * @param p The particle
  * @param xp The extended data of the particle
- * @param dt The time-step over which to drift
+ * @param t0 The time at the start of the drift
+ * @param t1 The time at the end of the drift
  */
 __attribute__((always_inline)) INLINE static void hydro_predict_extra(struct part* p,
 								      struct xpart* xp,
-								      float dt) {
- 
+								      float t0,
+								      float t1) {
+
+  const float dt = t1 - t0;
+  
+  p->rho *= expf(-p->div_v * dt);
+  p->h *= expf( 0.33333333f * p->div_v * dt);
+
+  const float dt_entr = t1 - 0.5f*(p->t_begin + p->t_end);
+  p->pressure = (p->entropy + p->entropy_dt * dt_entr) * powf(p->rho, const_hydro_gamma);
+  
 }
 
 

src/runner.c

@@ -707,7 +707,6 @@ void runner_dodrift(struct runner *r, struct cell *c, int timer) {
   const float dt = r->e->time - r->e->timeOld;
   struct part *restrict p, *restrict parts = c->parts;
   struct xpart *restrict xp, *restrict xparts = c->xparts;
-  float w;
   float dx_max = 0.f, h_max = 0.f;
   
   TIMER_TIC
@@ -722,10 +721,6 @@ void runner_dodrift(struct runner *r, struct cell *c, int timer) {
       p = &parts[k];
       xp = &xparts[k];
 
-      /* Get local copies of particle data. */
-      const float h = p->h;
-      const float ih = 1.0f / h;
-
       /* Drift... */
       p->x[0] += xp->v_full[0] * dt;
       p->x[1] += xp->v_full[1] * dt;
@@ -736,26 +731,30 @@ void runner_dodrift(struct runner *r, struct cell *c, int timer) {
       p->v[1] += p->a[1] * dt;
       p->v[2] += p->a[2] * dt;
 
-      /* Predict smoothing length */
-      w = p->force.h_dt * ih * dt;
-      if (fabsf(w) < 0.01f) /* 1st order expansion of exp(w) */
-	p->h *=
-	  1.0f +
-	  w * (1.0f + w * (0.5f + w * (1.0f / 6.0f + 1.0f / 24.0f * w)));
-      else
-	p->h *= expf(w);
+      /* /\* Predict smoothing length *\/ */
+      /* w = p->force.h_dt * ih * dt; */
+      /* if (fabsf(w) < 0.01f) /\* 1st order expansion of exp(w) *\/ */
+      /* 	p->h *= */
+      /* 	  1.0f + */
+      /* 	  w * (1.0f + w * (0.5f + w * (1.0f / 6.0f + 1.0f / 24.0f * w))); */
+      /* else */
+      /* 	p->h *= expf(w); */
+
+      //MATTHIEU
       
-      /* Predict density */
-      w = -3.0f * p->force.h_dt * ih * dt;
-      if (fabsf(w) < 0.1f)
-	p->rho *=
-	  1.0f +
-	  w * (1.0f + w * (0.5f + w * (1.0f / 6.0f + 1.0f / 24.0f * w)));
-      else
-	p->rho *= expf(w);
+      /* /\* Predict density *\/ */
+      /* w = -3.0f * p->force.h_dt * ih * dt; */
+      /* if (fabsf(w) < 0.1f) */
+      /* 	p->rho *= */
+      /* 	  1.0f + */
+      /* 	  w * (1.0f + w * (0.5f + w * (1.0f / 6.0f + 1.0f / 24.0f * w))); */
+      /* else */
+      /* 	p->rho *= expf(w); */
+
+
       
       /* Predict the values of the extra fields */
-      hydro_predict_extra(p, xp, dt);
+      hydro_predict_extra(p, xp, r->e->timeOld, r->e->time);
 
       /* Compute motion since last cell construction */
       const float dx = sqrtf((p->x[0] - xp->x_old[0]) * (p->x[0] - xp->x_old[0]) +