Soundspeed gets pre-computed before the force loop

src/hydro/Gadget2/hydro.h

@@ -126,6 +126,9 @@ INLINE static void hydro_prepare_force(struct part* p, struct xpart* xp, float t
   const float dt = time - 0.5f * (p->t_begin + p->t_end);
   p->force.pressure =
       (p->entropy + p->entropy_dt * dt) * powf(p->rho, const_hydro_gamma);
+
+  /* Compute the sound speed */
+  p->force.soundspeed = sqrtf(const_hydro_gamma * p->force.pressure / p->rho);
 }
 
 /**
@@ -164,12 +167,13 @@ __attribute__((always_inline))
 __attribute__((always_inline)) INLINE static void hydro_predict_extra(
     struct part* p, struct xpart* xp, float t0, float t1) {
 
-  // p->rho *= expf(-p->div_v * dt);
-  // p->h *= expf(0.33333333f * p->div_v * dt)
-
+  /* Drift the pressure */
   const float dt_entr = t1 - 0.5f * (p->t_begin + p->t_end);
   p->force.pressure =
       (p->entropy + p->entropy_dt * dt_entr) * powf(p->rho, const_hydro_gamma);
+
+  /* Compute the new sound speed */
+  p->force.soundspeed = sqrtf(const_hydro_gamma * p->force.pressure / p->rho);
 }
 
 /**

src/hydro/Gadget2/hydro_debug.h

@@ -24,13 +24,14 @@ __attribute__((always_inline))
       "v=[%.3e,%.3e,%.3e],v_full=[%.3e,%.3e,%.3e] \n a=[%.3e,%.3e,%.3e],\n "
       "h=%.3e, "
       "wcount=%d, m=%.3e, dh_drho=%.3e, rho=%.3e, P=%.3e, S=%.3e, "
-      "dS/dt=%.3e,\n"
+      "dS/dt=%.3e, c=%.3e\n"
       "divV=%.3e, curlV=%.3e, rotV=[%.3e,%.3e,%.3e]  \n "
       "v_sig=%e dh/dt=%.3e t_begin=%.3e, t_end=%.3e\n",
       p->x[0], p->x[1], p->x[2], p->v[0], p->v[1], p->v[2], xp->v_full[0],
       xp->v_full[1], xp->v_full[2], p->a_hydro[0], p->a_hydro[1], p->a_hydro[2],
       p->h, (int)p->density.wcount, p->mass, p->rho_dh, p->rho, p->force.pressure,
-      p->entropy, p->entropy_dt, p->div_v, p->force.curl_v, p->density.rot_v[0],
+      p->entropy, p->entropy_dt, p->force.soundspeed,
+      p->div_v, p->force.curl_v, p->density.rot_v[0],
       p->density.rot_v[1], p->density.rot_v[2], p->force.v_sig, p->force.h_dt,
       p->t_begin, p->t_end);
 }

src/hydro/Gadget2/hydro_iact.h

@@ -223,8 +223,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
   const float P_over_rho_j = pressurej / (rhoj * rhoj) * pj->rho_dh;
 
   /* Compute sound speeds */
-  const float ci = sqrtf(const_hydro_gamma * pressurei / rhoi);
-  const float cj = sqrtf(const_hydro_gamma * pressurej / rhoj);
+  const float ci = pi->force.soundspeed;
+  const float cj = pj->force.soundspeed;
   float v_sig = ci + cj;
 
   /* Compute dv dot r. */
@@ -334,8 +334,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
   const float P_over_rho_j = pressurej / (rhoj * rhoj) * pj->rho_dh;
 
   /* Compute sound speeds */
-  const float ci = sqrtf(const_hydro_gamma * pressurei / rhoi);
-  const float cj = sqrtf(const_hydro_gamma * pressurej / rhoj);
+  const float ci = pi->force.soundspeed;
+  const float cj = pj->force.soundspeed;
   float v_sig = ci + cj;
 
   /* Compute dv dot r. */

src/hydro/Gadget2/hydro_part.h

@@ -95,6 +95,9 @@ struct part {
       /* Particle pressure */
       float pressure;
 
+      /* Particle sound speed */
+      float soundspeed;
+      
     } force;
     
   };