Cleanup rho dh and use of the EoS functions

A few more rather unimportant cosmetic changes:

• moved rho_dh to density.rho_dh everywhere.
• added a force.f term to contain the transformed value of rho_dh.
• When computing the sound speed, I am now using the functions defined in equation_of_state.h everywhere.

The particles are the same size in memory as before despite the change. I had to update the tolerances for the tests as rho_dh changed meaning.

src/hydro/Gadget2/hydro.h

@@ -126,6 +126,17 @@ __attribute__((always_inline)) INLINE static void hydro_set_internal_energy(
     struct part *restrict p, float u) {
 
   p->entropy = gas_entropy_from_internal_energy(p->rho, u);
+
+  /* Compute the new pressure */
+  const float pressure = gas_pressure_from_internal_energy(p->rho, u);
+
+  /* Compute the new sound speed */
+  const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure);
+
+  const float rho_inv = 1.f / p->rho;
+
+  p->force.soundspeed = soundspeed;
+  p->force.P_over_rho2 = pressure * rho_inv * rho_inv;
 }
 
 /**
@@ -141,6 +152,17 @@ __attribute__((always_inline)) INLINE static void hydro_set_entropy(
     struct part *restrict p, float S) {
 
   p->entropy = S;
+
+  /* Compute the pressure */
+  const float pressure = gas_pressure_from_entropy(p->rho, p->entropy);
+
+  /* Compute the new sound speed */
+  const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure);
+
+  const float rho_inv = 1.f / p->rho;
+
+  p->force.soundspeed = soundspeed;
+  p->force.P_over_rho2 = pressure * rho_inv * rho_inv;
 }
 
 /**
@@ -195,7 +217,7 @@ __attribute__((always_inline)) INLINE static void hydro_init_part(
   p->density.wcount = 0.f;
   p->density.wcount_dh = 0.f;
   p->rho = 0.f;
-  p->rho_dh = 0.f;
+  p->density.rho_dh = 0.f;
   p->density.div_v = 0.f;
   p->density.rot_v[0] = 0.f;
   p->density.rot_v[1] = 0.f;
@@ -222,20 +244,17 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
 
   /* Final operation on the density (add self-contribution). */
   p->rho += p->mass * kernel_root;
-  p->rho_dh -= hydro_dimension * p->mass * kernel_root;
+  p->density.rho_dh -= hydro_dimension * p->mass * kernel_root;
   p->density.wcount += kernel_root;
 
   /* Finish the calculation by inserting the missing h-factors */
   p->rho *= h_inv_dim;
-  p->rho_dh *= h_inv_dim_plus_one;
+  p->density.rho_dh *= h_inv_dim_plus_one;
   p->density.wcount *= kernel_norm;
   p->density.wcount_dh *= h_inv * kernel_gamma * kernel_norm;
 
   const float rho_inv = 1.f / p->rho;
 
-  /* Compute the derivative term */
-  p->rho_dh = 1.f / (1.f + hydro_dimension_inv * p->h * p->rho_dh * rho_inv);
-
   /* Finish calculation of the velocity curl components */
   p->density.rot_v[0] *= h_inv_dim_plus_one * rho_inv;
   p->density.rot_v[1] *= h_inv_dim_plus_one * rho_inv;
@@ -273,19 +292,23 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
   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;
-
-  /* Divide the pressure by the density and density gradient */
-  const float P_over_rho2 = pressure * rho_inv * rho_inv * p->rho_dh;
-
   /* Compute the sound speed */
-  const float soundspeed = sqrtf(hydro_gamma * pressure * rho_inv);
+  const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure);
+
+  /* Divide the pressure by the density squared to get the SPH term */
+  const float rho_inv = 1.f / p->rho;
+  const float P_over_rho2 = pressure * rho_inv * rho_inv;
 
   /* Compute the Balsara switch */
   const float balsara =
       abs_div_v / (abs_div_v + curl_v + 0.0001f * soundspeed / fac_mu / p->h);
 
+  /* 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);
+
   /* Update variables. */
+  p->force.f = grad_h_term;
   p->force.P_over_rho2 = P_over_rho2;
   p->force.soundspeed = soundspeed;
   p->force.balsara = balsara;
@@ -349,17 +372,16 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
   const float dt_entr = (t1 - (p->ti_begin + p->ti_end) / 2) * timeBase;
   const float pressure = hydro_get_pressure(p, dt_entr);
 
-  const float rho_inv = 1.f / p->rho;
-
-  /* Divide the pressure by the density and density gradient */
-  const float P_over_rho2 = pressure * rho_inv * rho_inv * p->rho_dh;
-
   /* Compute the new sound speed */
-  const float soundspeed = sqrtf(hydro_gamma * pressure * rho_inv);
+  const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure);
+
+  /* Divide the pressure by the density squared to get the SPH term */
+  const float rho_inv = 1.f / p->rho;
+  const float P_over_rho2 = pressure * rho_inv * rho_inv;
 
   /* Update variables */
-  p->force.P_over_rho2 = P_over_rho2;
   p->force.soundspeed = soundspeed;
+  p->force.P_over_rho2 = P_over_rho2;
 }
 
 /**
@@ -400,6 +422,19 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
   /* Do not 'overcool' when timestep increases */
   if (p->entropy + p->entropy_dt * half_dt < 0.5f * p->entropy)
     p->entropy_dt = -0.5f * p->entropy / half_dt;
+
+  /* Compute the pressure */
+  const float pressure = gas_pressure_from_entropy(p->rho, p->entropy);
+
+  /* Compute the new sound speed */
+  const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure);
+
+  /* Divide the pressure by the density squared to get the SPH term */
+  const float rho_inv = 1.f / p->rho;
+  const float P_over_rho2 = pressure * rho_inv * rho_inv;
+
+  p->force.soundspeed = soundspeed;
+  p->force.P_over_rho2 = P_over_rho2;
 }
 
 /**
@@ -414,6 +449,19 @@ __attribute__((always_inline)) INLINE static void hydro_convert_quantities(
 
   /* We read u in the entropy field. We now get S from u */
   p->entropy = gas_entropy_from_internal_energy(p->rho, p->entropy);
+
+  /* Compute the pressure */
+  const float pressure = gas_pressure_from_entropy(p->rho, p->entropy);
+
+  /* Compute the sound speed */
+  const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure);
+
+  /* Divide the pressure by the density squared to get the SPH term */
+  const float rho_inv = 1.f / p->rho;
+  const float P_over_rho2 = pressure * rho_inv * rho_inv;
+
+  p->force.soundspeed = soundspeed;
+  p->force.P_over_rho2 = P_over_rho2;
 }
 
 #endif /* SWIFT_GADGET2_HYDRO_H */