PU: add update_soundspeed and remove computation in get_comoving_soundspeed

src/hydro/PressureEnergy/hydro.h

@@ -212,21 +212,33 @@ hydro_get_drifted_physical_entropy(const struct part *restrict p,
   return gas_entropy_from_internal_energy(p->rho, p->u);
 }
 
+
 /**
- * @brief Returns the comoving sound speed of a particle
+ * @brief Update the sound speed of a particle
  *
- * @param p The particle of interest
+ * @param p The particle of interest.
+ * @param cosmo The cosmological model.
  */
-__attribute__((always_inline)) INLINE static float
-hydro_get_comoving_soundspeed(const struct part *restrict p) {
+__attribute__((always_inline)) INLINE static void
+hydro_update_soundspeed(struct part *restrict p,
+			const struct cosmology *cosmo) {
 
   /* Compute the sound speed -- see theory section for justification */
   /* IDEAL GAS ONLY -- P-U does not work with generic EoS. */
   const float comoving_pressure =
-      pressure_floor_get_comoving_pressure(p, p->rho, p->pressure_bar);
-  const float square_rooted = sqrtf(hydro_gamma * comoving_pressure / p->rho);
+    pressure_floor_get_comoving_pressure(p, p->pressure_bar, cosmo);
+  p->force.soundspeed = sqrtf(hydro_gamma * comoving_pressure / p->rho);
+}
+
+/**
+ * @brief Returns the comoving sound speed of a particle
+ *
+ * @param p The particle of interest.
+ */
+__attribute__((always_inline)) INLINE static float
+hydro_get_comoving_soundspeed(const struct part *restrict p) {
 
-  return square_rooted;
+  return p->force.soundspeed;
 }
 
 /**
@@ -418,11 +430,8 @@ 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);
-
   /* Update variables. */
-  p->force.soundspeed = soundspeed;
+  hydro_update_soundspeed(p, cosmo);
 }
 
 /**
@@ -631,6 +640,7 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
   const float abs_div_v = fabsf(p->density.div_v);
 
   /* Compute the sound speed -- see theory section for justification */
+  hydro_update_soundspeed(p, cosmo);
   const float soundspeed = hydro_get_comoving_soundspeed(p);
 
   /* Compute the Balsara switch */
@@ -646,11 +656,10 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
 
   /* Get the pressures */
   const float comoving_pressure_with_floor =
-      pressure_floor_get_comoving_pressure(p, p->rho, p->pressure_bar);
+    pressure_floor_get_comoving_pressure(p, p->pressure_bar, cosmo);
 
   /* Update variables. */
   p->force.f = grad_h_term;
-  p->force.soundspeed = soundspeed;
   p->force.balsara = balsara;
   p->force.pressure_bar_with_floor = comoving_pressure_with_floor;
 }
@@ -685,7 +694,8 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration(
  * @param xp The extended data of this particle.
  */
 __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values(
-    struct part *restrict p, const struct xpart *restrict xp) {
+    struct part *restrict p, const struct xpart *restrict xp,
+    const struct cosmology *cosmo) {
 
   /* Re-set the predicted velocities */
   p->v[0] = xp->v_full[0];
@@ -696,9 +706,7 @@ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values(
   p->u = xp->u_full;
 
   /* Compute the sound speed */
-  const float soundspeed = hydro_get_comoving_soundspeed(p);
-
-  p->force.soundspeed = soundspeed;
+  hydro_update_soundspeed(p, cosmo);
 }
 
 /**
@@ -761,13 +769,11 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
   }
 
   /* Compute the new sound speed */
-  const float soundspeed = hydro_get_comoving_soundspeed(p);
-
-  p->force.soundspeed = soundspeed;
+  hydro_update_soundspeed(p, cosmo);
 
   /* update the required variables */
   const float comoving_pressure_with_floor =
-      pressure_floor_get_comoving_pressure(p, p->rho, p->pressure_bar);
+    pressure_floor_get_comoving_pressure(p, p->pressure_bar, cosmo);
   p->force.pressure_bar_with_floor = comoving_pressure_with_floor;
 }