Modify the particle getters to distinguish between kicked and drifted thermal quantities.

doc/Doxyfile.in

@@ -761,11 +761,13 @@ WARN_LOGFILE           =
 
 INPUT                  =  @top_srcdir@ @top_srcdir@/src @top_srcdir@/tests @top_srcdir@/examples
 INPUT		       += @top_srcdir@/src/hydro/Minimal
+INPUT		       += @top_srcdir@/src/hydro/Gadget2
 INPUT		       += @top_srcdir@/src/gravity/Default
 INPUT		       += @top_srcdir@/src/stars/Default
 INPUT		       += @top_srcdir@/src/riemann
 INPUT		       += @top_srcdir@/src/potential/point_mass
 INPUT		       += @top_srcdir@/src/equation_of_state/ideal_gas
+INPUT		       += @top_srcdir@/src/cooling/const_lambda
 INPUT		       += @top_srcdir@/src/cooling/EAGLE
 INPUT		       += @top_srcdir@/src/chemistry/EAGLE
 

examples/SmallCosmoVolume_cooling/plotTempEvolution.py

@@ -24,7 +24,7 @@ k_in_J_K = 1.38064852e-23
 mH_in_kg = 1.6737236e-27
 
 # Number of snapshots generated
-n_snapshots = 160
+n_snapshots = 200
 
 import matplotlib
 matplotlib.use("Agg")

src/cooling/const_lambda/cooling.h

@@ -98,7 +98,7 @@ __attribute__((always_inline)) INLINE static void cooling_cool_part(
   const float u_floor = hydro_props->minimal_internal_energy;
 
   /* Current energy */
-  const float u_old = hydro_get_physical_internal_energy2(p, xp, cosmo);
+  const float u_old = hydro_get_physical_internal_energy(p, xp, cosmo);
 
   /* Current du_dt in physical coordinates */
   const float hydro_du_dt = hydro_get_physical_internal_energy_dt(p, cosmo);
@@ -151,13 +151,15 @@ __attribute__((always_inline)) INLINE static void cooling_cool_part(
  * @param hydro_props The properties of the hydro scheme.
  * @param us The internal system of units.
  * @param p Pointer to the particle data.
+ * @param xp Pointer to the extended data of the particle.
  */
 __attribute__((always_inline)) INLINE static float cooling_timestep(
     const struct cooling_function_data* restrict cooling,
     const struct phys_const* restrict phys_const,
     const struct cosmology* restrict cosmo,
     const struct unit_system* restrict us,
-    const struct hydro_props* hydro_props, const struct part* restrict p) {
+    const struct hydro_props* hydro_props, const struct part* restrict p,
+    const struct xpart* restrict xp) {
 
   /* Start with the case where there is no limit */
   if (cooling->cooling_tstep_mult == FLT_MAX) {
@@ -165,7 +167,7 @@ __attribute__((always_inline)) INLINE static float cooling_timestep(
   }
 
   /* Get current internal energy and cooling rate */
-  const float u = hydro_get_physical_internal_energy(p, cosmo);
+  const float u = hydro_get_physical_internal_energy(p, xp, cosmo);
   const double cooling_du_dt_cgs =
       cooling_rate_cgs(cosmo, hydro_props, cooling, p);
 
@@ -174,7 +176,7 @@ __attribute__((always_inline)) INLINE static float cooling_timestep(
       cooling_du_dt_cgs * cooling->conv_factor_energy_rate_from_cgs;
 
   /* If we are close to (or below) the limit, we ignore the condition */
-  if (u < 1.01f * hydro_props->minimal_internal_energy)
+  if (u < 1.01f * hydro_props->minimal_internal_energy || cooling_du_dt == 0.f)
     return FLT_MAX;
   else
     return cooling->cooling_tstep_mult * u / fabsf(cooling_du_dt);

src/cooling/const_lambda/cooling_io.h

@@ -32,7 +32,8 @@
 
 /**
  * @brief Writes the current model of SPH to the file
- * @param h_grpsph The HDF5 group in which to write
+ * @param h_grp The HDF5 group in which to write
+ * @param cooling the parameters of the cooling function.
  */
 __attribute__((always_inline)) INLINE static void cooling_write_flavour(
     hid_t h_grp, const struct cooling_function_data* cooling) {

src/hydro/Gadget2/hydro.h

@@ -42,42 +42,60 @@
 #include "minmax.h"
 
 /**
- * @brief Returns the comoving internal energy of a particle
+ * @brief Returns the comoving internal energy of a particle at the last
+ * time the particle was kicked.
  *
  * @param p The particle of interest
+ * @param xp The extended data of the particle of interest.
  */
 __attribute__((always_inline)) INLINE static float
-hydro_get_comoving_internal_energy(const struct part *restrict p) {
+hydro_get_comoving_internal_energy(const struct part *restrict p,
+                                   const struct xpart *restrict xp) {
 
-  return gas_internal_energy_from_entropy(p->rho, p->entropy);
+  return gas_internal_energy_from_entropy(p->rho, xp->entropy_full);
 }
 
 /**
- * @brief Returns the physical internal energy of a particle
+ * @brief Returns the physical internal energy of a particle at the last
+ * time the particle was kicked.
  *
  * @param p The particle of interest.
+ * @param xp The extended data of the particle of interest.
  * @param cosmo The cosmological model.
  */
 __attribute__((always_inline)) INLINE static float
 hydro_get_physical_internal_energy(const struct part *restrict p,
+                                   const struct xpart *restrict xp,
                                    const struct cosmology *cosmo) {
 
-  return gas_internal_energy_from_entropy(p->rho * cosmo->a3_inv, p->entropy);
+  return gas_internal_energy_from_entropy(p->rho * cosmo->a3_inv,
+                                          xp->entropy_full);
 }
 
 /**
- * @brief Returns the physical internal energy of a particle
+ * @brief Returns the comoving internal energy of a particle drifted to the
+ * current time.
+ *
+ * @param p The particle of interest
+ */
+__attribute__((always_inline)) INLINE static float
+hydro_get_drifted_comoving_internal_energy(const struct part *restrict p) {
+
+  return gas_internal_energy_from_entropy(p->rho, p->entropy);
+}
+
+/**
+ * @brief Returns the physical internal energy of a particle drifted to the
+ * current time.
  *
  * @param p The particle of interest.
  * @param cosmo The cosmological model.
  */
 __attribute__((always_inline)) INLINE static float
-hydro_get_physical_internal_energy2(const struct part *restrict p,
-                                    const struct xpart *restrict xp,
-                                    const struct cosmology *cosmo) {
+hydro_get_drifted_physical_internal_energy(const struct part *restrict p,
+                                           const struct cosmology *cosmo) {
 
-  return gas_internal_energy_from_entropy(p->rho * cosmo->a3_inv,
-                                          xp->entropy_full);
+  return gas_internal_energy_from_entropy(p->rho * cosmo->a3_inv, p->entropy);
 }
 
 /**
@@ -94,7 +112,8 @@ __attribute__((always_inline)) INLINE static float hydro_get_comoving_pressure(
 /**
  * @brief Returns the physical pressure of a particle
  *
- * @param p The particle of interest
+ * @param p The particle of interest.
+ * @param cosmo The cosmological model.
  */
 __attribute__((always_inline)) INLINE static float hydro_get_physical_pressure(
     const struct part *restrict p, const struct cosmology *cosmo) {
@@ -103,24 +122,57 @@ __attribute__((always_inline)) INLINE static float hydro_get_physical_pressure(
 }
 
 /**
- * @brief Returns the comoving entropy of a particle
+ * @brief Returns the comoving entropy of a particle at the last
+ * time the particle was kicked.
  *
  * @param p The particle of interest.
+ * @param xp The extended data of the particle of interest.
  */
 __attribute__((always_inline)) INLINE static float hydro_get_comoving_entropy(
-    const struct part *restrict p) {
+    const struct part *restrict p, const struct xpart *restrict xp) {
 
-  return p->entropy;
+  return xp->entropy_full;
 }
 
 /**
- * @brief Returns the physical entropy of a particle
+ * @brief Returns the physical entropy of a particl at the last
+ * time the particle was kicked.
  *
  * @param p The particle of interest.
+ * @param xp The extended data of the particle of interest.
  * @param cosmo The cosmological model.
  */
 __attribute__((always_inline)) INLINE static float hydro_get_physical_entropy(
-    const struct part *restrict p, const struct cosmology *cosmo) {
+    const struct part *restrict p, const struct xpart *restrict xp,
+    const struct cosmology *cosmo) {
+
+  /* Note: no cosmological conversion required here with our choice of
+   * coordinates. */
+  return xp->entropy_full;
+}
+
+/**
+ * @brief Returns the comoving entropy of a particle drifted to the
+ * current time.
+ *
+ * @param p The particle of interest.
+ */
+__attribute__((always_inline)) INLINE static float
+hydro_get_drifted_comoving_entropy(const struct part *restrict p) {
+
+  return p->entropy;
+}
+
+/**
+ * @brief Returns the physical entropy of a particle drifted to the
+ * current time.
+ *
+ * @param p The particle of interest.
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static float
+hydro_get_drifted_physical_entropy(const struct part *restrict p,
+                                   const struct cosmology *cosmo) {
 
   /* Note: no cosmological conversion required here with our choice of
    * coordinates. */
@@ -581,6 +633,9 @@ __attribute__((always_inline)) INLINE static void hydro_end_force(
  * @param p The particle to act upon
  * @param xp The particle extended data to act upon
  * @param dt_therm The time-step for this kick (for thermodynamic quantities)
+ * @param dt_grav The time-step for this kick (for gravity forces)
+ * @param dt_hydro The time-step for this kick (for hydro forces)
+ * @param dt_kick_corr The time-step for this kick (for correction of the kick)
  * @param cosmo The cosmological model.
  * @param hydro_props The constants used in the scheme
  */

src/hydro/Gadget2/hydro_io.h

@@ -59,7 +59,7 @@ INLINE static void hydro_read_particles(struct part* parts,
 INLINE static void convert_part_u(const struct engine* e, const struct part* p,
                                   const struct xpart* xp, float* ret) {
 
-  ret[0] = hydro_get_comoving_internal_energy(p);
+  ret[0] = hydro_get_comoving_internal_energy(p, xp);
 }
 
 INLINE static void convert_part_P(const struct engine* e, const struct part* p,
@@ -132,6 +132,7 @@ INLINE static void convert_part_potential(const struct engine* e,
  * @brief Specifies which particle fields to write to a dataset
  *
  * @param parts The particle array.
+ * @param xparts The extended particle data array.
  * @param list The list of i/o properties to write.
  * @param num_fields The number of i/o fields to write.
  */

src/statistics.c

@@ -166,8 +166,8 @@ void stats_collect_part_mapper(void *map_data, int nr_parts, void *extra_data) {
     hydro_get_drifted_velocities(p, xp, dt_kick_hydro, dt_kick_grav, v);
     const double x[3] = {p->x[0], p->x[1], p->x[2]};
     const float m = hydro_get_mass(p);
-    const float entropy = hydro_get_physical_entropy(p, cosmo);
-    const float u_inter = hydro_get_physical_internal_energy(p, cosmo);
+    const float entropy = hydro_get_drifted_physical_entropy(p, cosmo);
+    const float u_inter = hydro_get_drifted_physical_internal_energy(p, cosmo);
 
     /* Collect mass */
     stats.mass += m;

src/timestep.h

@@ -128,7 +128,7 @@ __attribute__((always_inline)) INLINE static integertime_t get_part_timestep(
   if (e->policy & engine_policy_cooling)
     new_dt_cooling =
         cooling_timestep(e->cooling_func, e->physical_constants, e->cosmology,
-                         e->internal_units, e->hydro_properties, p);
+                         e->internal_units, e->hydro_properties, p, xp);
 
   /* Compute the next timestep (gravity condition) */
   float new_dt_grav = FLT_MAX, new_dt_self_grav = FLT_MAX,