diff --git a/src/const.h b/src/const.h
index 7a8bcbf6a82c3ce8474d5e449bbc23cf8f9b2ac5..a02e0b4437213e3032ee9c5b57bc8f487b220223 100644
--- a/src/const.h
+++ b/src/const.h
@@ -63,8 +63,8 @@
 //#define WENDLAND_C6_KERNEL
 
 /* SPH variant to use */
-//#define MINIMAL_SPH
-#define GADGET2_SPH
+#define MINIMAL_SPH
+//#define GADGET2_SPH
 //#define HOPKINS_PE_SPH
 //#define DEFAULT_SPH
 //#define GIZMO_SPH
diff --git a/src/hydro/Minimal/hydro.h b/src/hydro/Minimal/hydro.h
index 16939d49ac033ab587fd983581bd70b8e941920a..d4dfc68fd35f235806c32cd0a39603371ddd2233 100644
--- a/src/hydro/Minimal/hydro.h
+++ b/src/hydro/Minimal/hydro.h
@@ -141,6 +141,12 @@ __attribute__((always_inline)) INLINE static void hydro_set_internal_energy(
 
   /* Compute the pressure */
   const float pressure = gas_pressure_from_internal_energy(p->rho, p->u);
+
+  /* Compute the sound speed from the pressure*/
+  const float rho_inv = 1.f / p->rho;
+  const float soundspeed = sqrtf(hydro_gamma * pressure * rho_inv);
+
+  p->force.soundspeed = soundspeed;
   p->force.pressure = pressure;
 }
 
@@ -160,6 +166,12 @@ __attribute__((always_inline)) INLINE static void hydro_set_entropy(
 
   /* Compute the pressure */
   const float pressure = gas_pressure_from_internal_energy(p->rho, p->u);
+
+  /* Compute the sound speed from the pressure*/
+  const float rho_inv = 1.f / p->rho;
+  const float soundspeed = sqrtf(hydro_gamma * pressure * rho_inv);
+
+  p->force.soundspeed = soundspeed;
   p->force.pressure = pressure;
 }
 
@@ -277,10 +289,15 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
     struct part *restrict p, struct xpart *restrict xp, int ti_current,
     double timeBase) {
 
+  /* Compute the pressure */
   const float half_dt = (ti_current - (p->ti_begin + p->ti_end) / 2) * timeBase;
   const float pressure = hydro_get_pressure(p, half_dt);
+
   const float rho_inv = 1.f / p->rho;
 
+  /* Compute the sound speed */
+  const float soundspeed = sqrtf(hydro_gamma * pressure * rho_inv);
+
   /* Compute the "grad h" term */
   const float grad_h_term =
       1.f / (1.f + hydro_dimension_inv * p->h * p->density.rho_dh * rho_inv);
@@ -288,6 +305,7 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
   /* Update variables. */
   p->force.f = grad_h_term;
   p->force.pressure = pressure;
+  p->force.soundspeed = soundspeed;
 }
 
 /**
@@ -347,7 +365,15 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
 
   /* Drift the pressure */
   const float dt_entr = (t1 - (p->ti_begin + p->ti_end) / 2) * timeBase;
-  p->force.pressure = hydro_get_pressure(p, dt_entr);
+  const float pressure = hydro_get_pressure(p, dt_entr);
+
+  const float rho_inv = 1.f / p->rho;
+
+  /* Compute the new sound speed */
+  const float soundspeed = sqrtf(hydro_gamma * pressure * rho_inv);
+
+  p->force.pressure = pressure;
+  p->force.soundspeed = soundspeed;
 }
 
 /**
@@ -392,6 +418,12 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
 
   /* Compute the pressure */
   const float pressure = gas_pressure_from_internal_energy(p->rho, p->u);
+
+  /* Compute the sound speed from the pressure*/
+  const float rho_inv = 1.f / p->rho;
+  const float soundspeed = sqrtf(hydro_gamma * pressure * rho_inv);
+
+  p->force.soundspeed = soundspeed;
   p->force.pressure = pressure;
 }
 
diff --git a/src/hydro/Minimal/hydro_iact.h b/src/hydro/Minimal/hydro_iact.h
index 3e0770df6adf96eb7ac3c62393d70609d33dee18..169947b99e92d9bd1b0870d502a49e311820ff81 100644
--- a/src/hydro/Minimal/hydro_iact.h
+++ b/src/hydro/Minimal/hydro_iact.h
@@ -165,8 +165,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
   const float mu_ij = fac_mu * r_inv * omega_ij; /* This is 0 or negative */
 
   /* Compute sound speeds and signal velocity */
-  const float ci = sqrtf(hydro_gamma * pressurei / rhoi);
-  const float cj = sqrtf(hydro_gamma * pressurej / rhoj);
+  const float ci = pi->force.soundspeed;
+  const float cj = pj->force.soundspeed;
   const float v_sig = ci + cj - 3.f * mu_ij;
 
   /* Construct the full viscosity term */
@@ -276,8 +276,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
   const float mu_ij = fac_mu * r_inv * omega_ij; /* This is 0 or negative */
 
   /* Compute sound speeds and signal velocity */
-  const float ci = sqrtf(hydro_gamma * pressurei / rhoi);
-  const float cj = sqrtf(hydro_gamma * pressurej / rhoj);
+  const float ci = pi->force.soundspeed;
+  const float cj = pj->force.soundspeed;
   const float v_sig = ci + cj - 3.f * mu_ij;
 
   /* Construct the full viscosity term */
diff --git a/src/hydro/Minimal/hydro_part.h b/src/hydro/Minimal/hydro_part.h
index cc8086ae641938a6d1a7614576104ffd76b529f1..8542177278998d5e0b830dc164988611549ef24d 100644
--- a/src/hydro/Minimal/hydro_part.h
+++ b/src/hydro/Minimal/hydro_part.h
@@ -131,6 +131,9 @@ struct part {
       /*! Particle pressure. */
       float pressure;
 
+      /*! Particle soundspeed. */
+      float soundspeed;
+
       /*! Particle signal velocity */
       float v_sig;