diff --git a/src/hydro/PressureEntropy/hydro.h b/src/hydro/PressureEntropy/hydro.h
index 8c063596efd3be97ebb4da6b6879ac06122bd357..8ef993e5798d25a5ef25e7b72cd9d9ddb3716689 100644
--- a/src/hydro/PressureEntropy/hydro.h
+++ b/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 */
diff --git a/src/hydro/PressureEntropy/hydro_debug.h b/src/hydro/PressureEntropy/hydro_debug.h
index 486543793515795092e7cc97fe7b567b8230be3b..1192d9653be2be29c6e4ec14d5fcb1e6946a3759 100644
--- a/src/hydro/PressureEntropy/hydro_debug.h
+++ b/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 */
diff --git a/src/hydro/PressureEntropy/hydro_part.h b/src/hydro/PressureEntropy/hydro_part.h
index cac585ff79bae737f0e5c09860a38536cbf3a38c..b6e496918fa0e7989a8bddcfc5e8ea6b332c338e 100644
--- a/src/hydro/PressureEntropy/hydro_part.h
+++ b/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;