Updated PressureEntropy-SPH to the new scheme

src/hydro/PressureEntropy/hydro.h

@@ -202,27 +202,6 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
   return dt_cfl;
 }
 
-/**
- * @brief Initialises the particles for the first time
- *
- * This function is called only once just after the ICs have been
- * read in to do some conversions.
- *
- * @param p The particle to act upon
- * @param xp The extended particle data to act upon
- */
-__attribute__((always_inline)) INLINE static void hydro_first_init_part(
-    struct part *restrict p, struct xpart *restrict xp) {
-
-  p->ti_begin = 0;
-  p->ti_end = 0;
-  p->rho_bar = 0.f;
-  p->entropy_one_over_gamma = pow_one_over_gamma(p->entropy);
-  xp->v_full[0] = p->v[0];
-  xp->v_full[1] = p->v[1];
-  xp->v_full[2] = p->v[2];
-}
-
 /**
  * @brief Prepares a particle for the density calculation.
  *
@@ -302,12 +281,9 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
  *
  * @param p The particle to act upon
  * @param xp The extended particle data to act upon
- * @param ti_current The current time (on the timeline)
- * @param timeBase The minimal time-step size
  */
 __attribute__((always_inline)) INLINE static void hydro_prepare_force(
-    struct part *restrict p, struct xpart *restrict xp, int ti_current,
-    double timeBase) {
+    struct part *restrict p, struct xpart *restrict xp) {
 
   const float fac_mu = 1.f; /* Will change with cosmological integration */
 
@@ -320,9 +296,7 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
   const float abs_div_v = fabsf(p->density.div_v);
 
   /* Compute the pressure */
-  const float half_dt = (ti_current - (p->ti_begin + p->ti_end) / 2) * timeBase;
-  const float entropy = hydro_get_entropy(p, half_dt);
-  const float pressure = gas_pressure_from_entropy(p->rho_bar, entropy);
+  const float pressure = gas_pressure_from_entropy(p->rho_bar, p->entropy);
 
   /* Compute the sound speed from the pressure*/
   const float soundspeed = gas_soundspeed_from_pressure(p->rho_bar, pressure);
@@ -375,19 +349,34 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration(
   p->force.v_sig = 0.0f;
 }
 
+/**
+ * @brief Sets the values to be predicted in the drifts to their values at a
+ * kick time
+ *
+ * @param p The particle.
+ * @param xp The extended data of this particle.
+ */
+__attribute__((always_inline)) INLINE static void hydro_reset_predicted_values(
+    struct part *restrict p, const struct xpart *restrict xp) {
+
+  /* Re-set the predicted velocities */
+  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;
+}
+
 /**
  * @brief Predict additional particle fields forward in time when drifting
  *
  * @param p The particle
  * @param xp The extended data of the particle
  * @param dt The drift time-step.
- * @param t0 The time at the start of the drift (on the timeline).
- * @param t1 The time at the end of the drift (on the timeline).
- * @param timeBase The minimal time-step size
  */
 __attribute__((always_inline)) INLINE static void hydro_predict_extra(
-    struct part *restrict p, const struct xpart *restrict xp, float dt, int t0,
-    int t1, double timeBase) {
+    struct part *restrict p, const struct xpart *restrict xp, float dt) {
 
   const float h_inv = 1.f / p->h;
 
@@ -408,12 +397,11 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
     p->rho_bar *= expf(w2);
   }
 
-  /* Drift the entropy */
-  const float dt_entr = (t1 - (p->ti_begin + p->ti_end) / 2) * timeBase;
-  const float entropy = hydro_get_entropy(p, dt_entr);
+  /* Predict the entropy */
+  p->entropy += p->entropy_dt * dt;
 
   /* Compute the pressure */
-  const float pressure = gas_pressure_from_entropy(p->rho_bar, entropy);
+  const float pressure = gas_pressure_from_entropy(p->rho_bar, p->entropy);
 
   /* Compute the new sound speed */
   const float soundspeed = gas_soundspeed_from_pressure(p->rho_bar, pressure);
@@ -423,7 +411,7 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
   const float P_over_rho2 = pressure * rho_bar_inv * rho_bar_inv;
 
   /* Update the variables */
-  p->entropy_one_over_gamma = pow_one_over_gamma(entropy);
+  p->entropy_one_over_gamma = pow_one_over_gamma(p->entropy);
   p->force.soundspeed = soundspeed;
   p->force.P_over_rho2 = P_over_rho2;
 }
@@ -453,22 +441,18 @@ __attribute__((always_inline)) INLINE static void hydro_end_force(
  * @param half_dt The half time-step for this kick
  */
 __attribute__((always_inline)) INLINE static void hydro_kick_extra(
-    struct part *restrict p, struct xpart *restrict xp, float dt,
-    float half_dt) {
+    struct part *restrict p, struct xpart *restrict xp, 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;
+  if (entropy_change > -0.5f * xp->entropy_full)
+    xp->entropy_full += entropy_change;
   else
-    p->entropy *= 0.5f;
-
-  /* Do not 'overcool' when timestep increases */
-  if (p->entropy + p->entropy_dt * half_dt < 0.5f * p->entropy)
-    p->entropy_dt = -0.5f * p->entropy / half_dt;
+    xp->entropy_full *= 0.5f;
 
   /* Compute the pressure */
-  const float pressure = gas_pressure_from_entropy(p->rho_bar, p->entropy);
+  const float pressure =
+      gas_pressure_from_entropy(p->rho_bar, xp->entropy_full);
 
   /* Compute the new sound speed */
   const float soundspeed = gas_soundspeed_from_pressure(p->rho_bar, pressure);
@@ -490,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_bar, p->entropy);
+  xp->entropy_full = gas_entropy_from_internal_energy(p->rho_bar, p->entropy);
+  p->entropy = xp->entropy_full;
   p->entropy_one_over_gamma = pow_one_over_gamma(p->entropy);
 
   /* Compute the pressure */
@@ -510,4 +495,27 @@ __attribute__((always_inline)) INLINE static void hydro_convert_quantities(
   p->force.P_over_rho2 = P_over_rho2;
 }
 
+/**
+ * @brief Initialises the particles for the first time
+ *
+ * This function is called only once just after the ICs have been
+ * read in to do some conversions.
+ *
+ * @param p The particle to act upon
+ * @param xp The extended particle data to act upon
+ */
+__attribute__((always_inline)) INLINE static void hydro_first_init_part(
+    struct part *restrict p, struct xpart *restrict xp) {
+
+  p->time_bin = 0;
+  p->rho_bar = 0.f;
+  p->entropy_one_over_gamma = pow_one_over_gamma(p->entropy);
+  xp->v_full[0] = p->v[0];
+  xp->v_full[1] = p->v[1];
+  xp->v_full[2] = p->v[2];
+
+  hydro_reset_acceleration(p);
+  hydro_init_part(p);
+}
+
 #endif /* SWIFT_PRESSURE_ENTROPY_HYDRO_H */

src/hydro/PressureEntropy/hydro_debug.h

@@ -29,7 +29,6 @@
  * Follows eqautions (19), (21) and (22) of Hopkins, P., MNRAS, 2013,
  * Volume 428, Issue 4, pp. 2840-2856 with a simple Balsara viscosity term.
  */
-
 __attribute__((always_inline)) INLINE static void hydro_debug_particle(
     const struct part* p, const struct xpart* xp) {
   printf(
@@ -37,14 +36,14 @@ __attribute__((always_inline)) INLINE static void hydro_debug_particle(
       "v=[%.3e,%.3e,%.3e],v_full=[%.3e,%.3e,%.3e] \n a=[%.3e,%.3e,%.3e],\n "
       "h=%.3e, wcount=%.3f, wcount_dh=%.3e, m=%.3e, dh_drho=%.3e, rho=%.3e, "
       "rho_bar=%.3e, P=%.3e, dP_dh=%.3e, P_over_rho2=%.3e, S=%.3e, S^1/g=%.3e, "
-      "dS/dt=%.3e,\nc=%.3e v_sig=%e dh/dt=%.3e t_begin=%d, t_end=%d\n",
+      "dS/dt=%.3e,\nc=%.3e v_sig=%e dh/dt=%.3e time_bin=%d\n",
       p->x[0], p->x[1], p->x[2], p->v[0], p->v[1], p->v[2], xp->v_full[0],
       xp->v_full[1], xp->v_full[2], p->a_hydro[0], p->a_hydro[1], p->a_hydro[2],
       p->h, p->density.wcount, p->density.wcount_dh, p->mass, p->density.rho_dh,
       p->rho, p->rho_bar, hydro_get_pressure(p, 0.), p->density.pressure_dh,
       p->force.P_over_rho2, p->entropy, p->entropy_one_over_gamma,
       p->entropy_dt, p->force.soundspeed, p->force.v_sig, p->force.h_dt,
-      p->ti_begin, p->ti_end);
+      p->time_bin);
 }
 
 #endif /* SWIFT_PRESSURE_ENTROPY_HYDRO_DEBUG_H */

src/hydro/PressureEntropy/hydro_part.h

@@ -41,6 +41,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;
 
@@ -49,6 +52,12 @@ struct xpart {
 /* Data of a single particle. */
 struct part {
 
+  /*! Particle ID. */
+  long long id;
+
+  /*! Pointer to corresponding gravity part. */
+  struct gpart* gpart;
+
   /*! Particle position. */
   double x[3];
 
@@ -64,12 +73,6 @@ struct part {
   /*! Particle mass. */
   float mass;
 
-  /*! Particle time of beginning of time-step. */
-  int ti_begin;
-
-  /*! Particle time of end of time-step. */
-  int ti_end;
-
   /*! Particle density. */
   float rho;
 
@@ -132,11 +135,18 @@ struct part {
     } force;
   };
 
-  /*! Particle ID. */
-  long long id;
+  /* Time-step length */
+  timebin_t time_bin;
 
-  /*! Pointer to corresponding gravity part. */
-  struct gpart* gpart;
+#ifdef SWIFT_DEBUG_CHECKS
+
+  /* Time of the last drift */
+  integertime_t ti_drift;
+
+  /* Time of the last kick */
+  integertime_t ti_kick;
+
+#endif
 
 } SWIFT_STRUCT_ALIGN;