Make use of the EoS functions everywhere in PE-SPH.

src/const.h

@@ -64,8 +64,8 @@
 
 /* SPH variant to use */
 //#define MINIMAL_SPH
-#define GADGET2_SPH
-//#define HOPKINS_PE_SPH
+//#define GADGET2_SPH
+#define HOPKINS_PE_SPH
 //#define DEFAULT_SPH
 //#define GIZMO_SPH
 

src/hydro/PressureEntropy/hydro.h

@@ -129,12 +129,12 @@ __attribute__((always_inline)) INLINE static void hydro_set_internal_energy(
   p->entropy_one_over_gamma = pow_one_over_gamma(p->entropy);
 
   /* Compute the pressure */
-  const float entropy = p->entropy;
-  const float pressure = gas_pressure_from_entropy(p->rho_bar, entropy);
+  const float pressure = gas_pressure_from_entropy(p->rho_bar, p->entropy);
 
   /* Compute the sound speed from the pressure*/
+  const float soundspeed = gas_soundspeed_from_pressure(p->rho_bar, pressure);
   const float rho_bar_inv = 1.f / p->rho_bar;
-  const float soundspeed = sqrtf(hydro_gamma * pressure * rho_bar_inv);
+
   p->force.soundspeed = soundspeed;
   p->force.P_over_rho2 = pressure * rho_bar_inv * rho_bar_inv;
 }
@@ -155,12 +155,12 @@ __attribute__((always_inline)) INLINE static void hydro_set_entropy(
   p->entropy_one_over_gamma = pow_one_over_gamma(p->entropy);
 
   /* Compute the pressure */
-  const float entropy = p->entropy;
-  const float pressure = gas_pressure_from_entropy(p->rho_bar, entropy);
+  const float pressure = gas_pressure_from_entropy(p->rho_bar, p->entropy);
 
   /* Compute the sound speed from the pressure*/
+  const float soundspeed = gas_soundspeed_from_pressure(p->rho_bar, pressure);
   const float rho_bar_inv = 1.f / p->rho_bar;
-  const float soundspeed = sqrtf(hydro_gamma * pressure * rho_bar_inv);
+
   p->force.soundspeed = soundspeed;
   p->force.P_over_rho2 = pressure * rho_bar_inv * rho_bar_inv;
 }
@@ -216,10 +216,11 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
  */
 __attribute__((always_inline)) INLINE static void hydro_init_part(
     struct part *restrict p) {
-  p->density.wcount = 0.f;
-  p->density.wcount_dh = 0.f;
+
   p->rho = 0.f;
   p->rho_bar = 0.f;
+  p->density.wcount = 0.f;
+  p->density.wcount_dh = 0.f;
   p->density.rho_dh = 0.f;
   p->density.pressure_dh = 0.f;
 
@@ -313,8 +314,11 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
   const float pressure = hydro_get_pressure(p, half_dt);
 
   /* Compute the sound speed from the pressure*/
+  const float soundspeed = gas_soundspeed_from_pressure(p->rho_bar, pressure);
+
+  /* Divide the pressure by the density squared to get the SPH term */
   const float rho_bar_inv = 1.f / p->rho_bar;
-  const float soundspeed = sqrtf(hydro_gamma * pressure * rho_bar_inv);
+  const float P_over_rho2 = pressure * rho_bar_inv * rho_bar_inv;
 
   /* Compute the Balsara switch */
   const float balsara =
@@ -325,7 +329,7 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
 
   /* Update variables. */
   p->force.soundspeed = soundspeed;
-  p->force.P_over_rho2 = pressure * rho_bar_inv * rho_bar_inv;
+  p->force.P_over_rho2 = P_over_rho2;
   p->force.balsara = balsara;
   p->force.f = grad_h_term;
 }
@@ -394,14 +398,17 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
   /* Compute the pressure */
   const float pressure = gas_pressure_from_entropy(p->rho_bar, entropy);
 
-  /* Compute the sound speed from the pressure*/
+  /* Compute the new sound speed */
+  const float soundspeed = gas_soundspeed_from_pressure(p->rho_bar, pressure);
+
+  /* Divide the pressure by the density squared to get the SPH term */
   const float rho_bar_inv = 1.f / p->rho_bar;
-  const float soundspeed = sqrtf(hydro_gamma * pressure * rho_bar_inv);
+  const float P_over_rho2 = pressure * rho_bar_inv * rho_bar_inv;
 
   /* Update the variables */
   p->entropy_one_over_gamma = pow_one_over_gamma(entropy);
   p->force.soundspeed = soundspeed;
-  p->force.P_over_rho2 = pressure * rho_bar_inv * rho_bar_inv;
+  p->force.P_over_rho2 = P_over_rho2;
 }
 
 /**
@@ -446,13 +453,16 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
   /* Compute the pressure */
   const float pressure = gas_pressure_from_entropy(p->rho_bar, p->entropy);
 
-  /* Compute the sound speed from the pressure*/
+  /* Compute the new sound speed */
+  const float soundspeed = gas_soundspeed_from_pressure(p->rho_bar, pressure);
+
+  /* Divide the pressure by the density squared to get the SPH term */
   const float rho_bar_inv = 1.f / p->rho_bar;
-  const float soundspeed = sqrtf(hydro_gamma * pressure * rho_bar_inv);
+  const float P_over_rho2 = pressure * rho_bar_inv * rho_bar_inv;
 
   p->entropy_one_over_gamma = pow_one_over_gamma(p->entropy);
   p->force.soundspeed = soundspeed;
-  p->force.P_over_rho2 = pressure * rho_bar_inv * rho_bar_inv;
+  p->force.P_over_rho2 = P_over_rho2;
 }
 
 /**
@@ -470,14 +480,17 @@ __attribute__((always_inline)) INLINE static void hydro_convert_quantities(
   p->entropy_one_over_gamma = pow_one_over_gamma(p->entropy);
 
   /* Compute the pressure */
-  const float entropy = p->entropy;
-  const float pressure = gas_pressure_from_entropy(p->rho_bar, entropy);
+  const float pressure = gas_pressure_from_entropy(p->rho_bar, p->entropy);
 
-  /* Compute the sound speed from the pressure*/
+  /* Compute the sound speed */
+  const float soundspeed = gas_soundspeed_from_pressure(p->rho_bar, pressure);
+
+  /* Divide the pressure by the density squared to get the SPH term */
   const float rho_bar_inv = 1.f / p->rho_bar;
-  const float soundspeed = sqrtf(hydro_gamma * pressure * rho_bar_inv);
+  const float P_over_rho2 = pressure * rho_bar_inv * rho_bar_inv;
+
   p->force.soundspeed = soundspeed;
-  p->force.P_over_rho2 = pressure * rho_bar_inv * rho_bar_inv;
+  p->force.P_over_rho2 = P_over_rho2;
 }
 
 #endif /* SWIFT_PRESSURE_ENTROPY_HYDRO_H */