Commit 660d8d42 authored by Josh Borrow's avatar Josh Borrow Committed by Matthieu Schaller
Browse files

Move entropy floor until after density drift

parent c154e08c
......@@ -887,23 +887,6 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
/* Predict the internal energy */
p->u += p->u_dt * dt_therm;
internal_energy_ratio *= p->u;
/* Now we can use this to 'update' the value of the smoothed pressure. To
* truly update this variable, we would need another loop over neighbours
* using the new internal energies of everyone, but that's not feasible. */
p->pressure_bar *= internal_energy_ratio;
/* Check against entropy floor */
const float floor_A = entropy_floor(p, cosmo, floor_props);
const float floor_u = gas_internal_energy_from_entropy(p->rho, floor_A);
/* Check against absolute minimum */
const float min_u =
hydro_props->minimal_internal_energy / cosmo->a_factor_internal_energy;
p->u = max(p->u, floor_u);
p->u = max(p->u, min_u);
const float h_inv = 1.f / p->h;
......@@ -927,6 +910,26 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
p->pressure_bar *= expf_exact;
}
/* Check against entropy floor */
const float floor_A = entropy_floor(p, cosmo, floor_props);
const float floor_u = gas_internal_energy_from_entropy(p->rho, floor_A);
/* Check against absolute minimum */
const float min_u =
hydro_props->minimal_internal_energy / cosmo->a_factor_internal_energy;
p->u = max(p->u, floor_u);
p->u = max(p->u, min_u);
/* Now that p->u has been properly bounded, we can use it to apply the
* drift for the pressure */
internal_energy_ratio *= p->u;
/* Now we can use this to 'update' the value of the smoothed pressure. To
* truly update this variable, we would need another loop over neighbours
* using the new internal energies of everyone, but that's not feasible. */
p->pressure_bar *= internal_energy_ratio;
/* Compute the new sound speed */
const float soundspeed =
gas_soundspeed_from_pressure(p->rho, p->pressure_bar);
......
......@@ -825,25 +825,15 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
/* Predict the internal energy */
p->u += p->u_dt * dt_therm;
/* Check against entropy floor */
const float floor_A = entropy_floor(p, cosmo, floor_props);
const float floor_u = gas_internal_energy_from_entropy(p->rho, floor_A);
/* Check against absolute minimum */
const float min_u =
hydro_props->minimal_internal_energy / cosmo->a_factor_internal_energy;
p->u = max(p->u, floor_u);
p->u = max(p->u, min_u);
const float h_inv = 1.f / p->h;
/* Predict smoothing length */
const float w1 = p->force.h_dt * h_inv * dt_drift;
if (fabsf(w1) < 0.2f)
if (fabsf(w1) < 0.2f) {
p->h *= approx_expf(w1); /* 4th order expansion of exp(w) */
else
} else {
p->h *= expf(w1);
}
/* Predict density and weighted pressure */
const float w2 = -hydro_dimension * w1;
......@@ -856,6 +846,17 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
p->rho *= expf_exact;
}
/* Check against entropy floor */
const float floor_A = entropy_floor(p, cosmo, floor_props);
const float floor_u = gas_internal_energy_from_entropy(p->rho, floor_A);
/* Check against absolute minimum */
const float min_u =
hydro_props->minimal_internal_energy / cosmo->a_factor_internal_energy;
p->u = max(p->u, floor_u);
p->u = max(p->u, min_u);
/* Compute the new sound speed */
const float pressure = gas_pressure_from_internal_energy(p->rho, p->u);
const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure);
......
......@@ -718,6 +718,24 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
/* Predict the entropy */
p->entropy += p->entropy_dt * dt_therm;
const float h_inv = 1.f / p->h;
/* Predict smoothing length */
const float w1 = p->force.h_dt * h_inv * dt_drift;
if (fabsf(w1) < 0.2f) {
p->h *= approx_expf(w1); /* 4th order expansion of exp(w) */
} else {
p->h *= expf(w1);
}
/* Predict density */
const float w2 = -hydro_dimension * w1;
if (fabsf(w2) < 0.2f) {
p->rho *= approx_expf(w2); /* 4th order expansion of exp(w) */
} else {
p->rho *= expf(w2);
}
/* Check against entropy floor */
const float floor_A = entropy_floor(p, cosmo, floor_props);
......@@ -732,22 +750,6 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
p->entropy = max(p->entropy, floor_A);
p->entropy = max(p->entropy, min_A);
const float h_inv = 1.f / p->h;
/* Predict smoothing length */
const float w1 = p->force.h_dt * h_inv * dt_drift;
if (fabsf(w1) < 0.2f)
p->h *= approx_expf(w1); /* 4th order expansion of exp(w) */
else
p->h *= expf(w1);
/* Predict density */
const float w2 = -hydro_dimension * w1;
if (fabsf(w2) < 0.2f)
p->rho *= approx_expf(w2); /* 4th order expansion of exp(w) */
else
p->rho *= expf(w2);
/* Re-compute the pressure */
float comoving_pressure = gas_pressure_from_entropy(p->rho, p->entropy);
comoving_pressure =
......
......@@ -705,32 +705,34 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
/* Predict the internal energy */
p->u += p->u_dt * dt_therm;
/* Check against entropy floor */
const float floor_A = entropy_floor(p, cosmo, floor_props);
const float floor_u = gas_internal_energy_from_entropy(p->rho, floor_A);
/* Check against absolute minimum */
const float min_u =
hydro_props->minimal_internal_energy / cosmo->a_factor_internal_energy;
p->u = max(p->u, floor_u);
p->u = max(p->u, min_u);
const float h_inv = 1.f / p->h;
/* Predict smoothing length */
const float w1 = p->force.h_dt * h_inv * dt_drift;
if (fabsf(w1) < 0.2f)
if (fabsf(w1) < 0.2f) {
p->h *= approx_expf(w1); /* 4th order expansion of exp(w) */
else
} else {
p->h *= expf(w1);
}
/* Predict density */
const float w2 = -hydro_dimension * w1;
if (fabsf(w2) < 0.2f)
if (fabsf(w2) < 0.2f) {
p->rho *= approx_expf(w2); /* 4th order expansion of exp(w) */
else
} else {
p->rho *= expf(w2);
}
/* Check against entropy floor */
const float floor_A = entropy_floor(p, cosmo, floor_props);
const float floor_u = gas_internal_energy_from_entropy(p->rho, floor_A);
/* Check against absolute minimum */
const float min_u =
hydro_props->minimal_internal_energy / cosmo->a_factor_internal_energy;
p->u = max(p->u, floor_u);
p->u = max(p->u, min_u);
/* Compute the new pressure */
const float pressure = gas_pressure_from_internal_energy(p->rho, p->u);
......
......@@ -750,23 +750,6 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
/* Predict the internal energy */
p->u += p->u_dt * dt_therm;
internal_energy_ratio *= p->u;
/* Now we can use this to 'update' the value of the smoothed pressure. To
* truly update this variable, we would need another loop over neighbours
* using the new internal energies of everyone, but that's not feasible. */
p->pressure_bar *= internal_energy_ratio;
/* Check against entropy floor */
const float floor_A = entropy_floor(p, cosmo, floor_props);
const float floor_u = gas_internal_energy_from_entropy(p->rho, floor_A);
/* Check against absolute minimum */
const float min_u =
hydro_props->minimal_internal_energy / cosmo->a_factor_internal_energy;
p->u = max(p->u, floor_u);
p->u = max(p->u, min_u);
const float h_inv = 1.f / p->h;
......@@ -790,6 +773,26 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
p->pressure_bar *= expf_exact;
}
/* Check against entropy floor */
const float floor_A = entropy_floor(p, cosmo, floor_props);
const float floor_u = gas_internal_energy_from_entropy(p->rho, floor_A);
/* Check against absolute minimum */
const float min_u =
hydro_props->minimal_internal_energy / cosmo->a_factor_internal_energy;
p->u = max(p->u, floor_u);
p->u = max(p->u, min_u);
/* Now that p->u has been properly bounded, we can use it to apply the
* drift for the pressure */
internal_energy_ratio *= p->u;
/* Now we can use this to 'update' the value of the smoothed pressure. To
* truly update this variable, we would need another loop over neighbours
* using the new internal energies of everyone, but that's not feasible. */
p->pressure_bar *= internal_energy_ratio;
/* Compute the new sound speed */
hydro_update_soundspeed(p, cosmo);
......
......@@ -747,23 +747,6 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
/* Predict the internal energy */
p->u += p->u_dt * dt_therm;
internal_energy_ratio *= p->u;
/* Now we can use this to 'update' the value of the smoothed pressure. To
* truly update this variable, we would need another loop over neighbours
* using the new internal energies of everyone, but that's not feasible. */
p->pressure_bar *= internal_energy_ratio;
/* Check against entropy floor */
const float floor_A = entropy_floor(p, cosmo, floor_props);
const float floor_u = gas_internal_energy_from_entropy(p->rho, floor_A);
/* Check against absolute minimum */
const float min_u =
hydro_props->minimal_internal_energy / cosmo->a_factor_internal_energy;
p->u = max(p->u, floor_u);
p->u = max(p->u, min_u);
const float h_inv = 1.f / p->h;
......@@ -787,6 +770,26 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
p->pressure_bar *= expf_exact;
}
/* Check against entropy floor */
const float floor_A = entropy_floor(p, cosmo, floor_props);
const float floor_u = gas_internal_energy_from_entropy(p->rho, floor_A);
/* Check against absolute minimum */
const float min_u =
hydro_props->minimal_internal_energy / cosmo->a_factor_internal_energy;
p->u = max(p->u, floor_u);
p->u = max(p->u, min_u);
/* Now that p->u has been properly bounded, we can use it to apply the
* drift for the pressure */
internal_energy_ratio *= p->u;
/* Now we can use this to 'update' the value of the smoothed pressure. To
* truly update this variable, we would need another loop over neighbours
* using the new internal energies of everyone, but that's not feasible. */
p->pressure_bar *= internal_energy_ratio;
/* Compute the new sound speed */
const float soundspeed =
gas_soundspeed_from_pressure(p->rho, p->pressure_bar);
......
......@@ -683,17 +683,6 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
/* Predict the entropy */
p->entropy += p->entropy_dt * dt_therm;
/* Check against entropy floor */
const float floor_A = entropy_floor(p, cosmo, floor_props);
/* Check against absolute minimum */
const float min_u = hydro_props->minimal_internal_energy;
const float min_A =
gas_entropy_from_internal_energy(p->rho * cosmo->a3_inv, min_u);
p->entropy = max(p->entropy, floor_A);
p->entropy = max(p->entropy, min_A);
const float h_inv = 1.f / p->h;
/* Predict smoothing length */
......@@ -713,6 +702,17 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
p->rho_bar *= expf(w2);
}
/* Check against entropy floor */
const float floor_A = entropy_floor(p, cosmo, floor_props);
/* Check against absolute minimum */
const float min_u = hydro_props->minimal_internal_energy;
const float min_A =
gas_entropy_from_internal_energy(p->rho * cosmo->a3_inv, min_u);
p->entropy = max(p->entropy, floor_A);
p->entropy = max(p->entropy, min_A);
/* Compute the pressure */
const float pressure = gas_pressure_from_entropy(p->rho_bar, p->entropy);
......
......@@ -901,17 +901,6 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
/* Predict the internal energy */
p->u += p->u_dt * dt_therm;
/* Check against entropy floor */
const float floor_A = entropy_floor(p, cosmo, floor_props);
const float floor_u = gas_internal_energy_from_entropy(p->rho, floor_A);
/* Check against absolute minimum */
const float min_u =
hydro_props->minimal_internal_energy / cosmo->a_factor_internal_energy;
p->u = max(p->u, floor_u);
p->u = max(p->u, min_u);
const float h_inv = 1.f / p->h;
/* Predict smoothing length */
......@@ -932,6 +921,18 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
p->rho *= expf_exact;
}
/* Check against entropy floor - explicitly do this after drifting the
* density as this has a density dependence. */
const float floor_A = entropy_floor(p, cosmo, floor_props);
const float floor_u = gas_internal_energy_from_entropy(p->rho, floor_A);
/* Check against absolute minimum */
const float min_u =
hydro_props->minimal_internal_energy / cosmo->a_factor_internal_energy;
p->u = max(p->u, floor_u);
p->u = max(p->u, min_u);
/* Compute the new sound speed */
const float pressure = gas_pressure_from_internal_energy(p->rho, p->u);
const float pressure_including_floor =
......
Supports Markdown
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