Commit 8508aff1 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Impose the minimal energy in the kick operation.

parent 7394924f
......@@ -504,7 +504,8 @@ __attribute__((always_inline)) INLINE static void hydro_end_force(
* @param dt_therm The time-step for this kick (for thermodynamic quantities)
*/
__attribute__((always_inline)) INLINE static void hydro_kick_extra(
struct part *restrict p, struct xpart *restrict xp, float dt_therm) {}
struct part *restrict p, struct xpart *restrict xp, float dt_therm,
const struct cosmology *cosmo, const struct hydro_props *hydro_props) {}
/**
* @brief Converts hydro quantity of a particle at the start of a run
......
......@@ -525,19 +525,29 @@ __attribute__((always_inline)) INLINE static void hydro_end_force(
* @param p The particle to act upon
* @param xp The particle extended data to act upon
* @param dt_therm The time-step for this kick (for thermodynamic quantities)
* @param cosmo The cosmological model.
* @param hydro_props The constants used in the scheme
*/
__attribute__((always_inline)) INLINE static void hydro_kick_extra(
struct part *restrict p, struct xpart *restrict xp, float dt_therm) {
struct part *restrict p, struct xpart *restrict xp, float dt_therm,
const struct cosmology *cosmo, const struct hydro_props *hydro_props) {
/* Do not decrease the entropy by more than a factor of 2 */
if (dt_therm > 0. && p->entropy_dt * dt_therm < -0.5f * xp->entropy_full) {
/* message("Warning! Limiting entropy_dt. Possible cooling error.\n
* entropy_full = %g \n entropy_dt * dt =%g \n", */
/* xp->entropy_full,p->entropy_dt * dt); */
p->entropy_dt = -0.5f * xp->entropy_full / dt_therm;
}
xp->entropy_full += p->entropy_dt * dt_therm;
/* Apply the minimal energy limit */
const float density = p->rho * cosmo->a3_inv;
const float min_energy = hydro_props->minimal_internal_energy;
const float min_entropy =
gas_entropy_from_internal_energy(density, min_energy);
if (xp->entropy_full < min_entropy) {
xp->entropy_full = min_entropy;
p->entropy_dt = 0.f;
}
/* Compute the pressure */
const float pressure = gas_pressure_from_entropy(p->rho, xp->entropy_full);
......
......@@ -611,7 +611,8 @@ __attribute__((always_inline)) INLINE static void hydro_end_force(
* @param half_dt Half the physical time step.
*/
__attribute__((always_inline)) INLINE static void hydro_kick_extra(
struct part* p, struct xpart* xp, float dt) {
struct part* p, struct xpart* xp, float dt, const struct cosmology* cosmo,
const struct hydro_props* hydro_props) {
float a_grav[3];
......@@ -628,6 +629,14 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
p->conserved.energy += p->conserved.flux.energy * dt;
#endif
/* Apply the minimal energy limit */
const float min_energy =
hydro_props->minimal_internal_energy * cosmo->a_factor_internal_energy;
if (p->conserved.energy < min_energy) {
p->conserved.energy = min_energy;
p->conserved.flux.energy = 0.f;
}
gizmo_check_physical_quantity("mass", p->conserved.mass);
gizmo_check_physical_quantity("energy", p->conserved.energy);
......
......@@ -522,9 +522,12 @@ __attribute__((always_inline)) INLINE static void hydro_end_force(
* @param p The particle to act upon.
* @param xp The particle extended data to act upon.
* @param dt_therm The time-step for this kick (for thermodynamic quantities).
* @param cosmo The cosmological model.
* @param hydro_props The constants used in the scheme
*/
__attribute__((always_inline)) INLINE static void hydro_kick_extra(
struct part *restrict p, struct xpart *restrict xp, float dt_therm) {
struct part *restrict p, struct xpart *restrict xp, float dt_therm,
const struct cosmology *cosmo, const struct hydro_props *hydro_props) {
/* Do not decrease the energy by more than a factor of 2*/
if (dt_therm > 0. && p->u_dt * dt_therm < -0.5f * xp->u_full) {
......@@ -532,6 +535,14 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
}
xp->u_full += p->u_dt * dt_therm;
/* Apply the minimal energy limit */
const float min_energy =
hydro_props->minimal_internal_energy * cosmo->a_factor_internal_energy;
if (xp->u_full < min_energy) {
xp->u_full = min_energy;
p->u_dt = 0.f;
}
/* Compute the pressure */
const float pressure = gas_pressure_from_internal_energy(p->rho, xp->u_full);
......
......@@ -536,9 +536,12 @@ __attribute__((always_inline)) INLINE static void hydro_end_force(
* @param p The particle to act upon
* @param xp The particle extended data to act upon
* @param dt_therm The time-step for this kick (for thermodynamic quantities)
* @param cosmo The cosmological model.
* @param hydro_props The constants used in the scheme
*/
__attribute__((always_inline)) INLINE static void hydro_kick_extra(
struct part *restrict p, struct xpart *restrict xp, float dt_therm) {
struct part *restrict p, struct xpart *restrict xp, float dt_therm,
const struct cosmology *cosmo, const struct hydro_props *hydro_props) {
/* Do not decrease the entropy (temperature) by more than a factor of 2*/
if (dt_therm > 0. && p->entropy_dt * dt_therm < -0.5f * xp->entropy_full) {
......@@ -546,6 +549,16 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
}
xp->entropy_full += p->entropy_dt * dt_therm;
/* Apply the minimal energy limit */
const float density = p->rho_bar * cosmo->a3_inv;
const float min_energy = hydro_props->minimal_internal_energy;
const float min_entropy =
gas_entropy_from_internal_energy(density, min_energy);
if (xp->entropy_full < min_entropy) {
xp->entropy_full = min_entropy;
p->entropy_dt = 0.f;
}
/* Compute the pressure */
const float pressure =
gas_pressure_from_entropy(p->rho_bar, xp->entropy_full);
......
......@@ -68,13 +68,16 @@ __attribute__((always_inline)) INLINE static void kick_gpart(
* @param dt_kick_hydro The kick time-step for hydro accelerations.
* @param dt_kick_grav The kick time-step for gravity accelerations.
* @param dt_kick_therm The kick time-step for changes in thermal state.
* @param cosmo The cosmological model.
* @param hydro_props The constants used in the scheme
* @param ti_start The starting (integer) time of the kick (for debugging
* checks).
* @param ti_end The ending (integer) time of the kick (for debugging checks).
*/
__attribute__((always_inline)) INLINE static void kick_part(
struct part *restrict p, struct xpart *restrict xp, double dt_kick_hydro,
double dt_kick_grav, double dt_kick_therm, integertime_t ti_start,
double dt_kick_grav, double dt_kick_therm, const struct cosmology *cosmo,
const struct hydro_props *hydro_props, integertime_t ti_start,
integertime_t ti_end) {
#ifdef SWIFT_DEBUG_CHECKS
......@@ -107,7 +110,7 @@ __attribute__((always_inline)) INLINE static void kick_part(
}
/* Extra kick work */
hydro_kick_extra(p, xp, dt_kick_therm);
hydro_kick_extra(p, xp, dt_kick_therm, cosmo, hydro_props);
if (p->gpart != NULL) gravity_kick_extra(p->gpart, dt_kick_grav);
}
......
......@@ -981,6 +981,7 @@ void runner_do_kick1(struct runner *r, struct cell *c, int timer) {
const struct engine *e = r->e;
const struct cosmology *cosmo = e->cosmology;
const struct hydro_props *hydro_props = e->hydro_properties;
const int with_cosmology = (e->policy & engine_policy_cosmology);
struct part *restrict parts = c->parts;
struct xpart *restrict xparts = c->xparts;
......@@ -1044,8 +1045,8 @@ void runner_do_kick1(struct runner *r, struct cell *c, int timer) {
}
/* do the kick */
kick_part(p, xp, dt_kick_hydro, dt_kick_grav, dt_kick_therm, ti_begin,
ti_begin + ti_step / 2);
kick_part(p, xp, dt_kick_hydro, dt_kick_grav, dt_kick_therm, cosmo,
hydro_props, ti_begin, ti_begin + ti_step / 2);
/* Update the accelerations to be used in the drift for hydro */
if (p->gpart != NULL) {
......@@ -1150,6 +1151,7 @@ void runner_do_kick2(struct runner *r, struct cell *c, int timer) {
const struct engine *e = r->e;
const struct cosmology *cosmo = e->cosmology;
const struct hydro_props *hydro_props = e->hydro_properties;
const int with_cosmology = (e->policy & engine_policy_cosmology);
const int count = c->count;
const int gcount = c->gcount;
......@@ -1209,8 +1211,8 @@ void runner_do_kick2(struct runner *r, struct cell *c, int timer) {
}
/* Finish the time-step with a second half-kick */
kick_part(p, xp, dt_kick_hydro, dt_kick_grav, dt_kick_therm,
ti_begin + ti_step / 2, ti_begin + ti_step);
kick_part(p, xp, dt_kick_hydro, dt_kick_grav, dt_kick_therm, cosmo,
hydro_props, ti_begin + ti_step / 2, ti_begin + ti_step);
#ifdef SWIFT_DEBUG_CHECKS
/* Check that kick and the drift are synchronized */
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment