Make sure the EAGLE cooling does not lead to a rate going below the entropy floor.

src/cooling.h

@@ -27,6 +27,12 @@
 /* Config parameters. */
 #include "../config.h"
 
+/* Local includes */
+#include "parser.h"
+#include "physical_constants.h"
+#include "restart.h"
+#include "units.h"
+
 /* Import the right cooling definition */
 #if defined(COOLING_NONE)
 #include "./cooling/none/cooling.h"

src/cooling/EAGLE/cooling.c

@@ -36,6 +36,7 @@
 #include "cooling_rates.h"
 #include "cooling_struct.h"
 #include "cooling_tables.h"
+#include "entropy_floor.h"
 #include "error.h"
 #include "exp10.h"
 #include "hydro.h"
@@ -450,11 +451,12 @@ INLINE static double bisection_iter(
  * @param dt The cooling time-step of this particle.
  * @param dt_therm The hydro time-step of this particle.
  */
-void cooling_cool_part(const struct phys_const *restrict phys_const,
-                       const struct unit_system *restrict us,
-                       const struct cosmology *restrict cosmo,
-                       const struct hydro_props *restrict hydro_properties,
-                       const struct cooling_function_data *restrict cooling,
+void cooling_cool_part(const struct phys_const *phys_const,
+                       const struct unit_system *us,
+                       const struct cosmology *cosmo,
+                       const struct hydro_props *hydro_properties,
+                       const struct entropy_floor_properties *floor_props,
+                       const struct cooling_function_data *cooling,
                        struct part *restrict p, struct xpart *restrict xp,
                        const float dt, const float dt_therm) {
 
@@ -596,11 +598,21 @@ void cooling_cool_part(const struct phys_const *restrict phys_const,
 
   /* We now need to check that we are not going to go below any of the limits */
 
+  /* Limit imposed by the entropy floor */
+  const double A_floor = entropy_floor(p, cosmo, floor_props);
+  const double u_floor = gas_internal_energy_from_entropy(p->rho, A_floor);
+
+  /* Absolute minimum */
+  const double u_minimal = hydro_properties->minimal_internal_energy;
+
+  /* Largest of both limits */
+  const double u_limit = max(u_minimal, u_floor);
+
   /* First, check whether we may end up below the minimal energy after
    * this step 1/2 kick + another 1/2 kick that could potentially be for
    * a time-step twice as big. We hence check for 1.5 delta_u. */
-  if (u_start + 1.5 * delta_u < hydro_properties->minimal_internal_energy) {
-    delta_u = (hydro_properties->minimal_internal_energy - u_start) / 1.5;
+  if (u_start + 1.5 * delta_u < u_limit) {
+    delta_u = (u_limit - u_start) / 1.5;
   }
 
   /* Second, check whether the energy used in the prediction could get negative.

src/cooling/EAGLE/cooling.h

@@ -26,20 +26,22 @@
 
 /* Local includes. */
 #include "cooling_struct.h"
-#include "cosmology.h"
-#include "hydro_properties.h"
-#include "part.h"
-#include "physical_constants.h"
-#include "units.h"
+
+struct part;
+struct xpart;
+struct cosmology;
+struct hydro_props;
+struct entropy_floor_properties;
 
 void cooling_update(const struct cosmology *cosmo,
                     struct cooling_function_data *cooling);
 
-void cooling_cool_part(const struct phys_const *restrict phys_const,
-                       const struct unit_system *restrict us,
-                       const struct cosmology *restrict cosmo,
-                       const struct hydro_props *restrict hydro_properties,
-                       const struct cooling_function_data *restrict cooling,
+void cooling_cool_part(const struct phys_const *phys_const,
+                       const struct unit_system *us,
+                       const struct cosmology *cosmo,
+                       const struct hydro_props *hydro_properties,
+                       const struct entropy_floor_properties *floor_props,
+                       const struct cooling_function_data *cooling,
                        struct part *restrict p, struct xpart *restrict xp,
                        const float dt, const float dt_therm);
 

src/runner.c

@@ -416,6 +416,7 @@ void runner_do_cooling(struct runner *r, struct cell *c, int timer) {
   const struct phys_const *constants = e->physical_constants;
   const struct unit_system *us = e->internal_units;
   const struct hydro_props *hydro_props = e->hydro_properties;
+  const struct entropy_floor_properties *entropy_floor_props = e->entropy_floor;
   const double time_base = e->time_base;
   const integertime_t ti_current = e->ti_current;
   struct part *restrict parts = c->hydro.parts;
@@ -459,8 +460,9 @@ void runner_do_cooling(struct runner *r, struct cell *c, int timer) {
         }
 
         /* Let's cool ! */
-        cooling_cool_part(constants, us, cosmo, hydro_props, cooling_func, p,
-                          xp, dt_cool, dt_therm);
+        cooling_cool_part(constants, us, cosmo, hydro_props,
+                          entropy_floor_props, cooling_func, p, xp, dt_cool,
+                          dt_therm);
       }
     }
   }