Move the drifting of the entropy to be like Gadget-3

src/cell.c

@@ -686,9 +686,10 @@ void cell_sanitize(struct cell *c) {
 void cell_convert_hydro(struct cell *c, void *data) {
 
   struct part *p = c->parts;
+  struct xpart *xp = c->xparts;
 
   for (int i = 0; i < c->count; ++i) {
-    hydro_convert_quantities(&p[i]);
+    hydro_convert_quantities(&p[i], &xp[i]);
   }
 }
 

src/hydro/Gadget2/hydro.h

@@ -48,7 +48,7 @@
 __attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
     const struct part *restrict p, float dt) {
 
-  const float entropy = p->entropy + p->entropy_dt * dt;
+  const float entropy = p->entropy + dt * p->entropy_dt;
 
   return gas_internal_energy_from_entropy(p->rho, entropy);
 }
@@ -62,7 +62,7 @@ __attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
 __attribute__((always_inline)) INLINE static float hydro_get_pressure(
     const struct part *restrict p, float dt) {
 
-  const float entropy = p->entropy + p->entropy_dt * dt;
+  const float entropy = p->entropy + dt * p->entropy_dt;
 
   return gas_pressure_from_entropy(p->rho, entropy);
 }
@@ -76,7 +76,7 @@ __attribute__((always_inline)) INLINE static float hydro_get_pressure(
 __attribute__((always_inline)) INLINE static float hydro_get_entropy(
     const struct part *restrict p, float dt) {
 
-  return p->entropy + p->entropy_dt * dt;
+  return p->entropy + dt * p->entropy_dt;
 }
 
 /**
@@ -216,6 +216,7 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
   xp->v_full[0] = p->v[0];
   xp->v_full[1] = p->v[1];
   xp->v_full[2] = p->v[2];
+  xp->entropy_full = p->entropy;
 }
 
 /**
@@ -228,9 +229,10 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
  */
 __attribute__((always_inline)) INLINE static void hydro_init_part(
     struct part *restrict p) {
+
+  p->rho = 0.f;
   p->density.wcount = 0.f;
   p->density.wcount_dh = 0.f;
-  p->rho = 0.f;
   p->density.rho_dh = 0.f;
   p->density.div_v = 0.f;
   p->density.rot_v[0] = 0.f;
@@ -302,11 +304,7 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
   const float abs_div_v = fabsf(p->density.div_v);
 
   /* Compute the pressure */
-  const integertime_t ti_begin =
-      get_integer_time_begin(ti_current, p->time_bin);
-  const integertime_t ti_end = get_integer_time_end(ti_current, p->time_bin);
-  const float half_dt = (ti_current - (ti_begin + ti_end) / 2) * timeBase;
-  const float pressure = hydro_get_pressure(p, half_dt);
+  const float pressure = gas_pressure_from_entropy(p->rho, p->entropy);
 
   /* Compute the sound speed */
   const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure);
@@ -368,6 +366,9 @@ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values(
   p->v[0] = xp->v_full[0];
   p->v[1] = xp->v_full[1];
   p->v[2] = xp->v_full[2];
+
+  /* Re-set the entropy */
+  p->entropy = xp->entropy_full;
 }
 
 /**
@@ -400,11 +401,11 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
   else
     p->rho *= expf(w2);
 
-  /* Drift the pressure */
-  const integertime_t ti_begin = get_integer_time_begin(t0, p->time_bin);
-  const integertime_t ti_end = get_integer_time_end(t0, p->time_bin);
-  const float dt_entr = (t1 - (ti_begin + ti_end) / 2) * timeBase;
-  const float pressure = hydro_get_pressure(p, dt_entr);
+  /* Predict the entropy */
+  p->entropy += p->entropy_dt * dt;
+
+  /* Re-compute the pressure */
+  const float pressure = gas_pressure_from_entropy(p->rho, p->entropy);
 
   /* Compute the new sound speed */
   const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure);
@@ -446,13 +447,13 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
 
   /* 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;
+  if (entropy_change > -0.5f * xp->entropy_full)
+    xp->entropy_full += entropy_change;
   else
-    p->entropy *= 0.5f;
+    xp->entropy_full *= 0.5f;
 
   /* Compute the pressure */
-  const float pressure = gas_pressure_from_entropy(p->rho, p->entropy);
+  const float pressure = gas_pressure_from_entropy(p->rho, xp->entropy_full);
 
   /* Compute the new sound speed */
   const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure);
@@ -473,10 +474,11 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
  * @param p The particle to act upon
  */
 __attribute__((always_inline)) INLINE static void hydro_convert_quantities(
-    struct part *restrict p) {
+    struct part *restrict p, struct xpart *restrict xp) {
 
   /* We read u in the entropy field. We now get S from u */
-  p->entropy = gas_entropy_from_internal_energy(p->rho, p->entropy);
+  xp->entropy_full = gas_entropy_from_internal_energy(p->rho, p->entropy);
+  p->entropy = xp->entropy_full;
 
   /* Compute the pressure */
   const float pressure = gas_pressure_from_entropy(p->rho, p->entropy);

src/hydro/Gadget2/hydro_part.h

@@ -42,6 +42,9 @@ struct xpart {
   /* Velocity at the last full step. */
   float v_full[3];
 
+  /* Entropy at the last full step. */
+  float entropy_full;
+
   /* Additional data used to record cooling information */
   struct cooling_xpart_data cooling_data;
 

src/statistics.c

@@ -104,8 +104,8 @@ void stats_collect_part_mapper(void *map_data, int nr_parts, void *extra_data) {
   const struct part *restrict parts = (struct part *)map_data;
   const struct xpart *restrict xparts =
       s->xparts + (ptrdiff_t)(parts - s->parts);
-  // const integertime_t ti_current = s->e->ti_current;
-  // const double timeBase = s->e->timeBase;
+  const integertime_t ti_current = s->e->ti_current;
+  const double timeBase = s->e->timeBase;
   const double time = s->e->time;
   struct statistics *const global_stats = data->stats;
 
@@ -125,32 +125,27 @@ void stats_collect_part_mapper(void *map_data, int nr_parts, void *extra_data) {
     const struct xpart *xp = &xparts[k];
     const struct gpart *gp = (p->gpart != NULL) ? gp = p->gpart : NULL;
 
-    /* Get useful variables */
-    // const integertime_t ti_begin =
-    //    get_integer_time_begin(ti_current, p->time_bin);
-    // const integertime_t ti_end = get_integer_time_end(ti_current,
-    // p->time_bin);
-    // const integertime_t dti = get_integer_timestep(p->time_bin);
-    const float dt =
-        0.f;  //(ti_current - (ti_begin + ti_end) / 2) * timeBase; //MATTHIEU
-    const double x[3] = {p->x[0], p->x[1], p->x[2]};
+    /* Get useful time variables */
+    const integertime_t ti_begin =
+        get_integer_time_begin(ti_current, p->time_bin);
+    const integertime_t ti_step = get_integer_timestep(p->time_bin);
+    const float dt = (ti_current - (ti_begin + ti_step / 2)) * timeBase;
+
+    /* Get the total acceleration */
     float a_tot[3] = {p->a_hydro[0], p->a_hydro[1], p->a_hydro[2]};
     if (gp != NULL) {
       a_tot[0] += gp->a_grav[0];
       a_tot[1] += gp->a_grav[1];
       a_tot[2] += gp->a_grav[2];
     }
+
+    /* Extrapolate velocities to current time */
     const float v[3] = {xp->v_full[0] + a_tot[0] * dt,
                         xp->v_full[1] + a_tot[1] * dt,
                         xp->v_full[2] + a_tot[2] * dt};
 
     const float m = hydro_get_mass(p);
-
-    /* if (p->id == ICHECK) */
-    /*   message("bin=%d dti=%lld ti_begin=%lld ti_end=%lld dt=%e v=[%e %e %e]",
-     */
-    /* 	      p->time_bin, dti, ti_begin, ti_end, */
-    /*           dt, v[0], v[1], v[2]); */
+    const double x[3] = {p->x[0], p->x[1], p->x[2]};
 
     /* Collect mass */
     stats.mass += m;
@@ -167,13 +162,12 @@ void stats_collect_part_mapper(void *map_data, int nr_parts, void *extra_data) {
 
     /* Collect energies. */
     stats.E_kin += 0.5f * m * (v[0] * v[0] + v[1] * v[1] + v[2] * v[2]);
+    stats.E_int += m * hydro_get_internal_energy(p, dt);
+    stats.E_rad += cooling_get_radiated_energy(xp);
     stats.E_pot_self += 0.f;
     if (gp != NULL)
       stats.E_pot_ext += m * external_gravity_get_potential_energy(
                                  time, potential, phys_const, gp);
-    stats.E_int += m * hydro_get_internal_energy(p, dt);
-    stats.E_rad += cooling_get_radiated_energy(xp);
-
     /* Collect entropy */
     stats.entropy += m * hydro_get_entropy(p, dt);
   }