From d6ebadaf7cb76804a74ccdb832a3b1eba41eeb19 Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <matthieu.schaller@durham.ac.uk>
Date: Mon, 19 Sep 2016 21:22:31 +0100
Subject: [PATCH] When updating the sound-speed of a particle, also update its
 signal velocity.

---
 src/hydro/Gadget2/hydro.h         | 23 +++++++++++++++++++----
 src/hydro/Minimal/hydro.h         | 23 +++++++++++++++++++----
 src/hydro/PressureEntropy/hydro.h | 25 +++++++++++++++++++++----
 3 files changed, 59 insertions(+), 12 deletions(-)

diff --git a/src/hydro/Gadget2/hydro.h b/src/hydro/Gadget2/hydro.h
index 282f7f4819..0bdf56fe84 100644
--- a/src/hydro/Gadget2/hydro.h
+++ b/src/hydro/Gadget2/hydro.h
@@ -37,6 +37,7 @@
 #include "equation_of_state.h"
 #include "hydro_properties.h"
 #include "kernel_hydro.h"
+#include "minmax.h"
 
 /**
  * @brief Returns the internal energy of a particle
@@ -116,8 +117,9 @@ __attribute__((always_inline)) INLINE static float hydro_get_mass(
  * @brief Modifies the thermal state of a particle to the imposed internal
  * energy
  *
- * This overrides the current state of the particle but does *not* change its
- * time-derivatives
+ * This overwrites the current state of the particle but does *not* change its
+ * time-derivatives. Entropy, pressure, sound-speed and signal velocity will be
+ * updated.
  *
  * @param p The particle
  * @param u The new internal energy
@@ -133,17 +135,24 @@ __attribute__((always_inline)) INLINE static void hydro_set_internal_energy(
   /* Compute the new sound speed */
   const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure);
 
+  /* Update the signal velocity */
+  const float v_sig_old = p->force.v_sig;
+  const float v_sig_new = p->force.v_sig - p->force.soundspeed + soundspeed;
+  const float v_sig = max(v_sig_old, v_sig_new);
+
   const float rho_inv = 1.f / p->rho;
 
   p->force.soundspeed = soundspeed;
   p->force.P_over_rho2 = pressure * rho_inv * rho_inv;
+  p->force.v_sig = v_sig;
 }
 
 /**
  * @brief Modifies the thermal state of a particle to the imposed entropy
  *
- * This overrides the current state of the particle but does *not* change its
- * time-derivatives
+ * This overwrites the current state of the particle but does *not* change its
+ * time-derivatives. Entropy, pressure, sound-speed and signal velocity will be
+ * updated.
  *
  * @param p The particle
  * @param S The new entropy
@@ -159,10 +168,16 @@ __attribute__((always_inline)) INLINE static void hydro_set_entropy(
   /* Compute the new sound speed */
   const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure);
 
+  /* Update the signal velocity */
+  const float v_sig_old = p->force.v_sig;
+  const float v_sig_new = p->force.v_sig - p->force.soundspeed + soundspeed;
+  const float v_sig = max(v_sig_old, v_sig_new);
+
   const float rho_inv = 1.f / p->rho;
 
   p->force.soundspeed = soundspeed;
   p->force.P_over_rho2 = pressure * rho_inv * rho_inv;
+  p->force.v_sig = v_sig;
 }
 
 /**
diff --git a/src/hydro/Minimal/hydro.h b/src/hydro/Minimal/hydro.h
index 42dc2aa3d8..3015f26c6b 100644
--- a/src/hydro/Minimal/hydro.h
+++ b/src/hydro/Minimal/hydro.h
@@ -39,6 +39,7 @@
 #include "equation_of_state.h"
 #include "hydro_properties.h"
 #include "kernel_hydro.h"
+#include "minmax.h"
 
 /**
  * @brief Returns the internal energy of a particle
@@ -124,8 +125,9 @@ __attribute__((always_inline)) INLINE static float hydro_get_mass(
  * @brief Modifies the thermal state of a particle to the imposed internal
  * energy
  *
- * This overrides the current state of the particle but does *not* change its
- * time-derivatives
+ * This overwrites the current state of the particle but does *not* change its
+ * time-derivatives. Internal energy, pressure, sound-speed and signal velocity
+ * will be updated.
  *
  * @param p The particle
  * @param u The new internal energy
@@ -141,15 +143,22 @@ __attribute__((always_inline)) INLINE static void hydro_set_internal_energy(
   /* Compute the new sound speed */
   const float soundspeed = gas_soundspeed_from_internal_energy(p->rho, p->u);
 
+  /* Update the signal velocity */
+  const float v_sig_old = p->force.v_sig;
+  const float v_sig_new = p->force.v_sig - p->force.soundspeed + soundspeed;
+  const float v_sig = max(v_sig_old, v_sig_new);
+
   p->force.soundspeed = soundspeed;
   p->force.pressure = pressure;
+  p->force.v_sig = v_sig;
 }
 
 /**
  * @brief Modifies the thermal state of a particle to the imposed entropy
  *
- * This overrides the current state of the particle but does *not* change its
- * time-derivatives
+ * This overwrites the current state of the particle but does *not* change its
+ * time-derivatives. Internal energy, pressure, sound-speed and signal velocity
+ * will be updated.
  *
  * @param p The particle
  * @param S The new entropy
@@ -165,8 +174,14 @@ __attribute__((always_inline)) INLINE static void hydro_set_entropy(
   /* Compute the new sound speed */
   const float soundspeed = gas_soundspeed_from_internal_energy(p->rho, p->u);
 
+  /* Update the signal velocity */
+  const float v_sig_old = p->force.v_sig;
+  const float v_sig_new = p->force.v_sig - p->force.soundspeed + soundspeed;
+  const float v_sig = max(v_sig_old, v_sig_new);
+
   p->force.soundspeed = soundspeed;
   p->force.pressure = pressure;
+  p->force.v_sig = v_sig;
 }
 
 /**
diff --git a/src/hydro/PressureEntropy/hydro.h b/src/hydro/PressureEntropy/hydro.h
index a5e0d41804..ede16fe518 100644
--- a/src/hydro/PressureEntropy/hydro.h
+++ b/src/hydro/PressureEntropy/hydro.h
@@ -37,6 +37,7 @@
 #include "equation_of_state.h"
 #include "hydro_properties.h"
 #include "kernel_hydro.h"
+#include "minmax.h"
 
 /**
  * @brief Returns the internal energy of a particle
@@ -116,8 +117,9 @@ __attribute__((always_inline)) INLINE static float hydro_get_mass(
  * @brief Modifies the thermal state of a particle to the imposed internal
  * energy
  *
- * This overrides the current state of the particle but does *not* change its
- * time-derivatives
+ * This overwrites the current state of the particle but does *not* change its
+ * time-derivatives. Entropy, pressure, sound-speed and signal velocity will be
+ * updated.
  *
  * @param p The particle
  * @param u The new internal energy
@@ -133,17 +135,25 @@ __attribute__((always_inline)) INLINE static void hydro_set_internal_energy(
 
   /* Compute the sound speed from the pressure*/
   const float soundspeed = gas_soundspeed_from_pressure(p->rho_bar, pressure);
+
+  /* Update the signal velocity */
+  const float v_sig_old = p->force.v_sig;
+  const float v_sig_new = p->force.v_sig - p->force.soundspeed + soundspeed;
+  const float v_sig = max(v_sig_old, v_sig_new);
+
   const float rho_bar_inv = 1.f / p->rho_bar;
 
   p->force.soundspeed = soundspeed;
   p->force.P_over_rho2 = pressure * rho_bar_inv * rho_bar_inv;
+  p->force.v_sig = v_sig;
 }
 
 /**
  * @brief Modifies the thermal state of a particle to the imposed entropy
  *
- * This overrides the current state of the particle but does *not* change its
- * time-derivatives
+ * This overwrites the current state of the particle but does *not* change its
+ * time-derivatives. Entropy, pressure, sound-speed and signal velocity will be
+ * updated.
  *
  * @param p The particle
  * @param S The new entropy
@@ -159,10 +169,17 @@ __attribute__((always_inline)) INLINE static void hydro_set_entropy(
 
   /* Compute the sound speed from the pressure*/
   const float soundspeed = gas_soundspeed_from_pressure(p->rho_bar, pressure);
+
+  /* Update the signal velocity */
+  const float v_sig_old = p->force.v_sig;
+  const float v_sig_new = p->force.v_sig - p->force.soundspeed + soundspeed;
+  const float v_sig = max(v_sig_old, v_sig_new);
+
   const float rho_bar_inv = 1.f / p->rho_bar;
 
   p->force.soundspeed = soundspeed;
   p->force.P_over_rho2 = pressure * rho_bar_inv * rho_bar_inv;
+  p->force.v_sig = v_sig;
 }
 
 /**
-- 
GitLab