diff --git a/doc/Doxyfile.in b/doc/Doxyfile.in
index 0658dd27e043b11c7480eaf8b01c0c40809e215c..f444d05000bd079a1554c945df4b01bd8ca04d93 100644
--- a/doc/Doxyfile.in
+++ b/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
diff --git a/src/cooling/const_du/cooling.h b/src/cooling/const_du/cooling.h
index 634723f5a5a9af56aefc5fee1d409f583b243eba..70f704ae457c5ec71335ed6cba9758b53457cbb6 100644
--- a/src/cooling/const_du/cooling.h
+++ b/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;
}
diff --git a/src/cooling/const_lambda/cooling.h b/src/cooling/const_lambda/cooling.h
index 3fd583b55a31534996e23e4dce71bcae9394d7de..5b1cbf3954306d30e04001f84f1f900e89189d2b 100644
--- a/src/cooling/const_lambda/cooling.h
+++ b/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);
diff --git a/src/timestep.h b/src/timestep.h
index d92f88d06451892ce47db4b9468db9714bb52baa..6f6ef74f2bce60f32d643f26ebca377730385438 100644
--- a/src/timestep.h
+++ b/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 =