diff --git a/src/hydro/AnarchyDU/hydro.h b/src/hydro/AnarchyDU/hydro.h
index a8b84c46aca9a117efe8173c492bbfd678106a5f..559570650b7086539181e24213dcb30b00bbdb16 100644
--- a/src/hydro/AnarchyDU/hydro.h
+++ b/src/hydro/AnarchyDU/hydro.h
@@ -398,17 +398,23 @@ hydro_set_drifted_physical_internal_energy(struct part *p,
const struct cosmology *cosmo,
const float u) {
+ /* There is no need to use the floor here as this function is called in the
+ * feedback, so the new value of the internal energy should be strictly
+ * higher than the old value. */
+
p->u = u / cosmo->a_factor_internal_energy;
/* Now recompute the extra quantities */
/* Compute the sound speed */
- const float soundspeed = hydro_get_comoving_soundspeed(p);
- const float pressure = hydro_get_comoving_pressure(p);
+ const float pressure = gas_pressure_from_internal_energy(p->rho, p->u);
+ const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure);
/* Update variables. */
p->force.soundspeed = soundspeed;
p->force.pressure = pressure;
+
+ p->viscosity.v_sig = max(p->viscosity.v_sig, 2.f * soundspeed);
}
/**
@@ -842,6 +848,9 @@ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values(
p->force.pressure = pressure;
p->force.soundspeed = soundspeed;
+
+ /* Update the signal velocity, if we need to. */
+ p->viscosity.v_sig = max(p->viscosity.v_sig, 2.f * soundspeed);
}
/**
@@ -903,10 +912,13 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
/* Compute the new sound speed */
const float pressure = gas_pressure_from_internal_energy(p->rho, p->u);
- const float soundspeed = hydro_get_comoving_soundspeed(p);
+ const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure);
p->force.pressure = pressure;
p->force.soundspeed = soundspeed;
+
+ /* Update signal velocity if we need to */
+ p->viscosity.v_sig = max(p->viscosity.v_sig, 2.f * soundspeed);
}
/**
diff --git a/src/hydro/AnarchyPU/hydro.h b/src/hydro/AnarchyPU/hydro.h
index 806d4aed204c71d3b9265d1bbea1098fdb83ebb1..04882ee12ef6e8d764ef6dfcd92bd3b10b4a7f2f 100644
--- a/src/hydro/AnarchyPU/hydro.h
+++ b/src/hydro/AnarchyPU/hydro.h
@@ -219,9 +219,7 @@ hydro_get_comoving_soundspeed(const struct part *restrict p) {
/* Compute the sound speed -- see theory section for justification */
/* IDEAL GAS ONLY -- P-U does not work with generic EoS. */
- const float square_rooted = sqrtf(hydro_gamma * p->pressure_bar / p->rho);
-
- return square_rooted;
+ return gas_soundspeed_from_pressure(p->rho, p->pressure_bar);
}
/**
@@ -408,15 +406,28 @@ hydro_set_drifted_physical_internal_energy(struct part *p,
const struct cosmology *cosmo,
const float u) {
+ /* Store ratio of new internal energy to old internal energy, as we use this
+ * in the drifting of the pressure. */
+ float internal_energy_ratio = 1.f / p->u;
+
+ /* Update the internal energy */
p->u = u / cosmo->a_factor_internal_energy;
+ 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;
/* Now recompute the extra quantities */
/* Compute the sound speed */
- const float soundspeed = hydro_get_comoving_soundspeed(p);
+ const float soundspeed = gas_soundspeed_from_pressure(p->rho; p->pressure_bar);
/* Update variables. */
p->force.soundspeed = soundspeed;
+
+ p->viscosity.v_sig = max(p->viscosity.v_sig, 2.f * soundspeed);
}
/**
@@ -467,10 +478,7 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
const float dt_cfl = 2.f * kernel_gamma * CFL_condition * cosmo->a * p->h /
(cosmo->a_factor_sound_speed * p->viscosity.v_sig);
- const float dt_u_change =
- (p->u_dt != 0.0f) ? fabsf(const_max_u_change * p->u / p->u_dt) : FLT_MAX;
-
- return fminf(dt_cfl, dt_u_change);
+ return dt_cfl;
}
/**
@@ -871,8 +879,18 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
const struct hydro_props *hydro_props,
const struct entropy_floor_properties *floor_props) {
+ /* Store ratio of new internal energy to old internal energy, as we use this
+ * in the drifting of the pressure. */
+ float internal_energy_ratio = 1.f / p->u;
+
/* 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);
@@ -908,9 +926,11 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
}
/* Compute the new sound speed */
- const float soundspeed = hydro_get_comoving_soundspeed(p);
-
+ const float soundspeed = gas_soundspeed_from_pressure(p->rho, p->pressure_bar);
p->force.soundspeed = soundspeed;
+
+ /* Update the signal velocity */
+ p->viscosity.v_sig = max(p->viscosity.v_sig, 2.f * soundspeed);
}
/**
diff --git a/src/hydro/Default/hydro.h b/src/hydro/Default/hydro.h
index 469fee137ee624ab649718f25c490f61e35cfb59..5a62770d9c4ddc1381ee9b390c06d0e8e607d32b 100644
--- a/src/hydro/Default/hydro.h
+++ b/src/hydro/Default/hydro.h
@@ -403,12 +403,14 @@ hydro_set_drifted_physical_internal_energy(struct part *p,
/* Now recompute the extra quantities */
/* Compute the sound speed */
- const float soundspeed = hydro_get_comoving_soundspeed(p);
- const float pressure = hydro_get_comoving_pressure(p);
+ const float pressure = gas_pressure_from_internal_energy(p->rho, p->u);
+ const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure);
/* Update variables. */
p->force.soundspeed = soundspeed;
p->force.pressure = pressure;
+
+ p->viscosity.v_sig = max(p->viscosity.v_sig, 2.f * soundspeed);
}
/**
@@ -782,10 +784,12 @@ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values(
/* Compute the sound speed */
const float pressure = gas_pressure_from_internal_energy(p->rho, p->u);
- const float soundspeed = hydro_get_comoving_soundspeed(p);
+ const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure);
p->force.pressure = pressure;
p->force.soundspeed = soundspeed;
+
+ p->viscosity.v_sig = max(p->viscosity.v_sig, 2.f * soundspeed);
}
/**
@@ -847,10 +851,12 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
/* Compute the new sound speed */
const float pressure = gas_pressure_from_internal_energy(p->rho, p->u);
- const float soundspeed = hydro_get_comoving_soundspeed(p);
+ const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure);
p->force.pressure = pressure;
p->force.soundspeed = soundspeed;
+
+ p->viscosity.v_sig = max(p->viscosity.v_sig, 2.f * soundspeed);
}
/**
diff --git a/src/hydro/Gadget2/hydro.h b/src/hydro/Gadget2/hydro.h
index 3fa6f19cc145890feacbd7284368d5378654bf38..00fc59e70c2f421f39a647d104e1369e6c6f626b 100644
--- a/src/hydro/Gadget2/hydro.h
+++ b/src/hydro/Gadget2/hydro.h
@@ -403,6 +403,8 @@ hydro_set_drifted_physical_internal_energy(struct part *p,
/* Update variables. */
p->force.P_over_rho2 = P_over_rho2;
p->force.soundspeed = soundspeed;
+
+ p->force.v_sig = max(p->force.v_sig, 2.f * soundspeed);
}
/**
@@ -650,7 +652,7 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration(
p->force.h_dt = 0.0f;
/* Reset maximal signal velocity */
- p->force.v_sig = p->force.soundspeed;
+ p->force.v_sig = 2.f * p->force.soundspeed;
}
/**
@@ -757,6 +759,8 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
/* Update variables */
p->force.soundspeed = soundspeed;
p->force.P_over_rho2 = P_over_rho2;
+
+ p->force.v_sig = max(p->force.v_sig, 2.f * soundspeed);
}
/**
diff --git a/src/hydro/Minimal/hydro.h b/src/hydro/Minimal/hydro.h
index 3d7f43579033afd8f5c29e765b31fee145d9c590..cdfbdbf3bcf3790b8cba84ce74084d3885cf75bf 100644
--- a/src/hydro/Minimal/hydro.h
+++ b/src/hydro/Minimal/hydro.h
@@ -390,10 +390,14 @@ hydro_set_drifted_physical_internal_energy(struct part *p,
/* Now recompute the extra quantities */
/* Compute the sound speed */
- const float soundspeed = hydro_get_comoving_soundspeed(p);
+ const float pressure = gas_pressure_from_internal_energy(p->rho, p->u);
+ const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure);
/* Update variables. */
+ p->force.pressure = pressure;
p->force.soundspeed = soundspeed;
+
+ p->force.v_sig = max(p->force.v_sig, 2.f * soundspeed);
}
/**
@@ -633,7 +637,7 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration(
/* Reset the time derivatives. */
p->u_dt = 0.0f;
p->force.h_dt = 0.0f;
- p->force.v_sig = p->force.soundspeed;
+ p->force.v_sig = 2.f * p->force.soundspeed;
}
/**
@@ -665,6 +669,8 @@ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values(
/* Update variables */
p->force.pressure = pressure;
p->force.soundspeed = soundspeed;
+
+ p->force.v_sig = max(p->force.v_sig, 2.f * soundspeed);
}
/**
@@ -728,6 +734,8 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
p->force.pressure = pressure;
p->force.soundspeed = soundspeed;
+
+ p->force.v_sig = max(p->force.v_sig, 2.f * soundspeed);
}
/**
diff --git a/src/hydro/PressureEnergy/hydro.h b/src/hydro/PressureEnergy/hydro.h
index e49445f59edd7a3a09860bfffe81a56d6a05abd4..97a768bbab5959c5995dde9070c819ffa9f724ca 100644
--- a/src/hydro/PressureEnergy/hydro.h
+++ b/src/hydro/PressureEnergy/hydro.h
@@ -226,6 +226,9 @@ __attribute__((always_inline)) INLINE static void hydro_update_soundspeed(
const float comoving_pressure =
pressure_floor_get_comoving_pressure(p, p->pressure_bar, cosmo);
p->force.soundspeed = gas_soundspeed_from_pressure(p->rho, comoving_pressure);
+
+ /* Also update the signal velocity; this could be a huge change! */
+ p->force.v_sig = max(p->force.v_sig, 2.f * p->force.soundspeed);
}
/**
@@ -424,9 +427,18 @@ hydro_set_drifted_physical_internal_energy(struct part *p,
const struct cosmology *cosmo,
const float u) {
+ /* Store ratio of new internal energy to old internal energy, as we use this
+ * in the drifting of the pressure. */
+ float internal_energy_ratio = 1.f / p->u;
+
+ /* Update the internal energy */
p->u = u / cosmo->a_factor_internal_energy;
+ internal_energy_ratio *= p->u;
- /* Now recompute the extra quantities */
+ /* 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;
/* Update variables. */
hydro_update_soundspeed(p, cosmo);
@@ -476,10 +488,7 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
const float dt_cfl = 2.f * kernel_gamma * CFL_condition * cosmo->a * p->h /
(cosmo->a_factor_sound_speed * p->force.v_sig);
- const float dt_u_change =
- (p->u_dt != 0.0f) ? fabsf(const_max_u_change * p->u / p->u_dt) : FLT_MAX;
-
- return fminf(dt_cfl, dt_u_change);
+ return dt_cfl;
}
/**
@@ -677,7 +686,7 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration(
/* Reset the time derivatives. */
p->u_dt = 0.0f;
p->force.h_dt = 0.0f;
- p->force.v_sig = p->force.soundspeed;
+ p->force.v_sig = 2.f * p->force.soundspeed;
}
/**
@@ -726,9 +735,19 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
float dt_therm, const struct cosmology *cosmo,
const struct hydro_props *hydro_props,
const struct entropy_floor_properties *floor_props) {
+
+ /* Store ratio of new internal energy to old internal energy, as we use this
+ * in the drifting of the pressure. */
+ float internal_energy_ratio = 1.f / p->u;
/* 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);
diff --git a/src/hydro/PressureEnergyMorrisMonaghanAV/hydro.h b/src/hydro/PressureEnergyMorrisMonaghanAV/hydro.h
index 3f4578b7bd05ec292d59535474a56c635e4bfb9c..4f1b560d8e5b63553a9698cc5b75792ce6f3ee9a 100644
--- a/src/hydro/PressureEnergyMorrisMonaghanAV/hydro.h
+++ b/src/hydro/PressureEnergyMorrisMonaghanAV/hydro.h
@@ -164,9 +164,8 @@ hydro_get_comoving_soundspeed(const struct part *restrict p) {
/* Compute the sound speed -- see theory section for justification */
/* IDEAL GAS ONLY -- P-U does not work with generic EoS. */
- const float square_rooted = sqrtf(hydro_gamma * p->pressure_bar / p->rho);
- return square_rooted;
+ return gas_soundspeed_from_pressure(p->rho, p->pressure_bar);
}
/**
@@ -407,15 +406,28 @@ hydro_set_drifted_physical_internal_energy(struct part *p,
const struct cosmology *cosmo,
const float u) {
+ /* Store ratio of new internal energy to old internal energy, as we use this
+ * in the drifting of the pressure. */
+ float internal_energy_ratio = 1.f / p->u;
+
+ /* Update the internal energy */
p->u = u / cosmo->a_factor_internal_energy;
+ 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;
/* Now recompute the extra quantities */
/* Compute the sound speed */
- const float soundspeed = hydro_get_comoving_soundspeed(p);
+ const float soundspeed =
+ gas_soundspeed_from_pressure(p->rho, p->pressure_bar);
/* Update variables. */
p->force.soundspeed = soundspeed;
+ p->viscosity.v_sig = max(p->viscosity.v_sig, 2.f * soundspeed);
}
/**
@@ -463,10 +475,7 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
const float dt_cfl = 2.f * kernel_gamma * CFL_condition * cosmo->a * p->h /
(cosmo->a_factor_sound_speed * p->force.v_sig);
- const float dt_u_change =
- (p->u_dt != 0.0f) ? fabsf(const_max_u_change * p->u / p->u_dt) : FLT_MAX;
-
- return fminf(dt_cfl, dt_u_change);
+ return dt_cfl;
}
/**
@@ -680,7 +689,7 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration(
/* Reset the time derivatives. */
p->u_dt = 0.0f;
p->force.h_dt = 0.0f;
- p->force.v_sig = p->force.soundspeed;
+ p->force.v_sig = 2.f * p->force.soundspeed;
}
/**
@@ -732,8 +741,18 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
const struct hydro_props *hydro_props,
const struct entropy_floor_properties *floor_props) {
+ /* Store ratio of new internal energy to old internal energy, as we use this
+ * in the drifting of the pressure. */
+ float internal_energy_ratio = 1.f / p->u;
+
/* 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);
@@ -769,9 +788,10 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
}
/* Compute the new sound speed */
- const float soundspeed = hydro_get_comoving_soundspeed(p);
+ const float soundspeed = gas_soundspeed_from_pressure(p->rho, p->pressure_bar);
p->force.soundspeed = soundspeed;
+ p->viscosity.v_sig = max(p->viscosity.v_sig, 2.f * soundspeed);
}
/**