From 94a7cafa11ec97c0d4154c6ab27486afa728e4ba Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <schaller@strw.leidenuniv.nl>
Date: Mon, 11 Feb 2019 22:52:59 +0100
Subject: [PATCH] Make sure the EAGLE cooling does not lead to a rate going
below the entropy floor.
---
src/cooling.h | 6 ++++++
src/cooling/EAGLE/cooling.c | 26 +++++++++++++++++++-------
src/cooling/EAGLE/cooling.h | 22 ++++++++++++----------
src/runner.c | 6 ++++--
4 files changed, 41 insertions(+), 19 deletions(-)
diff --git a/src/cooling.h b/src/cooling.h
index c1a78e256f..875ef50544 100644
--- a/src/cooling.h
+++ b/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"
diff --git a/src/cooling/EAGLE/cooling.c b/src/cooling/EAGLE/cooling.c
index d657beeaa4..8f4dc31718 100644
--- a/src/cooling/EAGLE/cooling.c
+++ b/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.
diff --git a/src/cooling/EAGLE/cooling.h b/src/cooling/EAGLE/cooling.h
index a811007136..d95c75e58a 100644
--- a/src/cooling/EAGLE/cooling.h
+++ b/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);
diff --git a/src/runner.c b/src/runner.c
index e15bed84c7..5548bef3dd 100644
--- a/src/runner.c
+++ b/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);
}
}
}
--
GitLab