diff --git a/src/cooling.h b/src/cooling.h index c1a78e256fdd77fcb1f5cde074f843bd16d412ec..875ef5054491f783d526e7c8e2caf3e005c8a5a0 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 d657beeaa4f36d1efd3c5cc407f8c1d689037693..8f4dc31718e54ba56120ce823b1d25b12d3ffb67 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 a811007136168646a98df0c0da7b43227de082f8..d95c75e58aecfd8fb4816f1e50c7a3f379a08e51 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 e15bed84c75bca8a33ab30cea929fd973af2fb23..5548bef3dd280f4f407f08eaff909918c4d36bcd 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); } } }