Use the cooling time-step condition when computing dt. Better documentation....

Use the cooling time-step condition when computing dt. Better documentation. More uniform interface accross cooling functions.

doc/Doxyfile.in

@@ -763,7 +763,7 @@ INPUT                  =  @top_srcdir@ @top_srcdir@/src @top_srcdir@/tests @top_
 INPUT		       += @top_srcdir@/src/hydro/Minimal
 INPUT		       += @top_srcdir@/src/gravity/Default
 INPUT		       += @top_srcdir@/src/riemann
-INPUT		       += @top_srcdir@/src/cooling/const_lambda
+INPUT		       += @top_srcdir@/src/cooling/const_du
 
 # This tag can be used to specify the character encoding of the source files
 # that doxygen parses. Internally doxygen uses the UTF-8 encoding. Doxygen uses

src/cooling/const_du/cooling.h

@@ -22,11 +22,12 @@
 #define SWIFT_COOLING_CONST_DU_H
 
 /**
- * @file src/cooling/const/cooling.h
+ * @file src/cooling/const_du/cooling.h
  * @brief Routines related to the "constant cooling" cooling function.
  *
  * This is the simplest possible cooling function. A constant cooling rate with
- * a minimal energy floor is applied.
+ * a minimal energy floor is applied. Should be used as a template for more
+ * realistic functions.
  */
 
 /* Some standard headers. */
@@ -59,6 +60,9 @@ struct cooling_data {
 /**
  * @brief Apply the cooling function to a particle.
  *
+ * In this simple example we just apply the const cooling rate
+ * and check that we don't go below the given floor.
+ *
  * @param phys_const The physical constants in internal units.
  * @param us The internal system of units.
  * @param cooling The #cooling_data used in the run.
@@ -66,8 +70,10 @@ struct cooling_data {
  * @param dt The time-step of this particle.
  */
 __attribute__((always_inline)) INLINE static void cooling_cool_part(
-    const struct phys_const* const phys_const, const struct UnitSystem* us,
-    const struct cooling_data* cooling, struct part* p, float dt) {
+    const struct phys_const* restrict phys_const,
+    const struct UnitSystem* restrict us,
+    const struct cooling_data* restrict cooling, struct part* restrict p,
+    float dt) {
 
   /* Get current internal energy (dt=0) */
   const float u_old = hydro_get_internal_energy(p, 0.f);
@@ -91,17 +97,22 @@ __attribute__((always_inline)) INLINE static void cooling_cool_part(
 /**
  * @brief Computes the cooling time-step.
  *
+ * In this simple example, we return \f$ \alpha \frac{u}{du/dt} \f$.
+ * This is used to compute the time-step of the particle. Cooling functions
+ * that are solved implicitly can simply return FLT_MAX here.
+ *
  * @param cooling The #cooling_data used in the run.
  * @param phys_const The physical constants in internal units.
+ * @param us The internal system of units.
  * @param p Pointer to the particle data.
  */
 __attribute__((always_inline)) INLINE static double cooling_timestep(
-    const struct cooling_data* cooling,
-    const struct phys_const* const phys_const, const struct part* const p) {
+    const struct cooling_data* restrict cooling,
+    const struct phys_const* restrict phys_const,
+    const struct UnitSystem* restrict us, const struct part* restrict p) {
 
   const float cooling_rate = cooling->cooling_rate;
-  const float internal_energy =
-      hydro_get_internal_energy(p, 0);  // dt = 0 because using current entropy
+  const float internal_energy = hydro_get_internal_energy(p, 0);
   return cooling->cooling_tstep_mult * internal_energy / cooling_rate;
 }
 

src/cooling/const_lambda/cooling.h

@@ -106,8 +106,10 @@ __attribute__((always_inline)) INLINE static float cooling_rate(
  * @param dt The time-step of this particle.
  */
 __attribute__((always_inline)) INLINE static void cooling_cool_part(
-    const struct phys_const* const phys_const, const struct UnitSystem* us,
-    const struct cooling_data* cooling, struct part* p, float dt) {
+    const struct phys_const* restrict phys_const,
+    const struct UnitSystem* restrict us,
+    const struct cooling_data* restrict cooling, struct part* restrict p,
+    float dt) {
 
   /* Get current internal energy (dt=0) */
   const float u_old = hydro_get_internal_energy(p, 0.f);
@@ -146,9 +148,9 @@ __attribute__((always_inline)) INLINE static void cooling_cool_part(
  * @param p Pointer to the particle data.
  */
 __attribute__((always_inline)) INLINE static float cooling_timestep(
-    const struct cooling_data* cooling,
-    const struct phys_const* const phys_const, const struct UnitSystem* us,
-    const struct part* const p) {
+    const struct cooling_data* restrict cooling,
+    const struct phys_const* restrict phys_const,
+    const struct UnitSystem* restrict us, const struct part* restrict p) {
 
   /* Get du_dt */
   const float du_dt = cooling_rate(phys_const, us, cooling, p);

src/timestep.h

@@ -74,12 +74,11 @@ __attribute__((always_inline)) INLINE static int get_gpart_timestep(
   /*     gravity_compute_timestep_self(e->physical_constants, gp); */
   const float new_dt_self = FLT_MAX;  // MATTHIEU
 
-  float new_dt =
-      (new_dt_external < new_dt_self) ? new_dt_external : new_dt_self;
+  float new_dt = min(new_dt_external, new_dt_self);
 
   /* Limit timestep within the allowed range */
-  new_dt = (new_dt < e->dt_max) ? new_dt : e->dt_max;
-  new_dt = (new_dt > e->dt_min) ? new_dt : e->dt_min;
+  new_dt = min(new_dt, e->dt_max);
+  new_dt = max(new_dt, e->dt_min);
 
   /* Convert to integer time */
   const int new_dti =
@@ -102,6 +101,12 @@ __attribute__((always_inline)) INLINE static int get_part_timestep(
   /* Compute the next timestep (hydro condition) */
   const float new_dt_hydro = hydro_compute_timestep(p, xp, e->hydro_properties);
 
+  /* Compute the next timestep (cooling condition) */
+  float new_dt_cooling = FLT_MAX;
+  if (e->policy & engine_policy_cooling)
+    new_dt_cooling = cooling_timestep(e->cooling_data, e->physical_constants,
+                                      e->internalUnits, p);
+
   /* Compute the next timestep (gravity condition) */
   float new_dt_grav = FLT_MAX;
   if (p->gpart != NULL) {
@@ -112,12 +117,11 @@ __attribute__((always_inline)) INLINE static int get_part_timestep(
     /*     gravity_compute_timestep_self(e->physical_constants, p->gpart); */
     const float new_dt_self = FLT_MAX;  // MATTHIEU
 
-    new_dt_grav =
-        (new_dt_external < new_dt_self) ? new_dt_external : new_dt_self;
+    new_dt_grav = min(new_dt_external, new_dt_self);
   }
 
   /* Final time-step is minimum of hydro and gravity */
-  float new_dt = (new_dt_hydro < new_dt_grav) ? new_dt_hydro : new_dt_grav;
+  float new_dt = min(min(new_dt_hydro, new_dt_cooling), new_dt_grav);
 
   /* Limit change in h */
   const float dt_h_change =
@@ -125,11 +129,11 @@ __attribute__((always_inline)) INLINE static int get_part_timestep(
           ? fabsf(e->hydro_properties->log_max_h_change * p->h / p->force.h_dt)
           : FLT_MAX;
 
-  new_dt = (new_dt < dt_h_change) ? new_dt : dt_h_change;
+  new_dt = min(new_dt, dt_h_change);
 
   /* Limit timestep within the allowed range */
-  new_dt = (new_dt < e->dt_max) ? new_dt : e->dt_max;
-  new_dt = (new_dt > e->dt_min) ? new_dt : e->dt_min;
+  new_dt = min(new_dt, e->dt_max);
+  new_dt = max(new_dt, e->dt_min);
 
   /* Convert to integer time */
   const int new_dti =