From 8508aff1d3941c9395d3ae320b47420a3b2c65bb Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <matthieu.schaller@durham.ac.uk>
Date: Sat, 17 Mar 2018 19:21:37 +0100
Subject: [PATCH] Impose the minimal energy in the kick operation.
---
src/hydro/Default/hydro.h | 3 ++-
src/hydro/Gadget2/hydro.h | 18 ++++++++++++++----
src/hydro/Gizmo/hydro.h | 11 ++++++++++-
src/hydro/Minimal/hydro.h | 13 ++++++++++++-
src/hydro/PressureEntropy/hydro.h | 15 ++++++++++++++-
src/kick.h | 7 +++++--
src/runner.c | 10 ++++++----
7 files changed, 63 insertions(+), 14 deletions(-)
diff --git a/src/hydro/Default/hydro.h b/src/hydro/Default/hydro.h
index d541ee7774..7da0054562 100644
--- a/src/hydro/Default/hydro.h
+++ b/src/hydro/Default/hydro.h
@@ -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
diff --git a/src/hydro/Gadget2/hydro.h b/src/hydro/Gadget2/hydro.h
index ad60465951..bc06a24e2a 100644
--- a/src/hydro/Gadget2/hydro.h
+++ b/src/hydro/Gadget2/hydro.h
@@ -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);
diff --git a/src/hydro/Gizmo/hydro.h b/src/hydro/Gizmo/hydro.h
index 66cfac1be5..bc90a6790f 100644
--- a/src/hydro/Gizmo/hydro.h
+++ b/src/hydro/Gizmo/hydro.h
@@ -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);
diff --git a/src/hydro/Minimal/hydro.h b/src/hydro/Minimal/hydro.h
index 1070df105a..3f9d99683b 100644
--- a/src/hydro/Minimal/hydro.h
+++ b/src/hydro/Minimal/hydro.h
@@ -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);
diff --git a/src/hydro/PressureEntropy/hydro.h b/src/hydro/PressureEntropy/hydro.h
index 78cc9c1217..87d46c6d43 100644
--- a/src/hydro/PressureEntropy/hydro.h
+++ b/src/hydro/PressureEntropy/hydro.h
@@ -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);
diff --git a/src/kick.h b/src/kick.h
index 9b1f4f1811..9d10f1e78d 100644
--- a/src/kick.h
+++ b/src/kick.h
@@ -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);
}
diff --git a/src/runner.c b/src/runner.c
index 1b8f9a2fc8..35996a3ba9 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -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 */
--
GitLab