Merge branch 'v_sig_reset' into 'master'

Update the signal velocity in more locations

See merge request !901

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);
 }
 
 /**

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);
 }
 
 /**

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);
 }
 
 /**

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);
 }
 
 /**

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);
 }
 
 /**

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);

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);
 }
 
 /**