Preliminary version of Hopkins' PE-SPH. Works with fixdt. Need to think about...

Preliminary version of Hopkins' PE-SPH. Works with fixdt. Need to think about the drift and grad h terms.

src/adiabatic_index.h

@@ -149,4 +149,91 @@ __attribute__((always_inline)) INLINE static float pow_minus_gamma_minus_one(
 #endif
 }
 
+/**
+ * @brief Returns one over the argument to the power given by one over the
+ * adiabatic index.
+ *
+ * Computes \f$x^{\frac{1}{\gamma}}\f$.
+ */
+__attribute__((always_inline)) INLINE static float pow_one_over_gamma(float x) {
+
+#if defined(HYDRO_GAMMA_5_3)
+
+  return powf(x, 0.6f); /* x^(3/5) */
+
+#elif defined(HYDRO_GAMMA_4_3)
+
+  return powf(x, 0.75f); /* x^(3/4) */
+
+#elif defined(HYDRO_GAMMA_2_1)
+
+  return sqrtf(x); /* x^(1/2) */
+
+#else
+
+  error("The adiabatic index is not defined !");
+  return 0.f;
+
+#endif
+}
+
+/**
+ * @brief Returns the argument to the power given by two over the adiabatic
+ * index.
+ *
+ * Computes \f$x^{\frac{2}{\gamma}}\f$.
+ */
+__attribute__((always_inline)) INLINE static float pow_two_over_gamma(float x) {
+
+#if defined(HYDRO_GAMMA_5_3)
+
+  return powf(x, 1.2f); /* x^(6/5) */
+
+#elif defined(HYDRO_GAMMA_4_3)
+
+  const float sqrt = sqrtf(x);
+  return sqrt * sqrt * sqrt;
+
+#elif defined(HYDRO_GAMMA_2_1)
+
+  return x; /* x^(2/2) */
+
+#else
+
+  error("The adiabatic index is not defined !");
+  return 0.f;
+
+#endif
+}
+
+/**
+ * @brief Returns the argument to the power one minus two over the
+ * adiabatic index.
+ *
+ * Computes \f$x^{1 - \frac{2}{\gamma}}\f$.
+ */
+__attribute__((always_inline)) INLINE static float pow_one_minus_two_over_gamma(
+    float x) {
+
+#if defined(HYDRO_GAMMA_5_3)
+
+  return powf(x, -0.2f); /* x^(-1/5) */
+
+#elif defined(HYDRO_GAMMA_4_3)
+
+  const float sqrt = sqrtf(x);
+  return 1.f / sqrt;
+
+#elif defined(HYDRO_GAMMA_2_1)
+
+  return 1.f; /* x^0 */
+
+#else
+
+  error("The adiabatic index is not defined !");
+  return 0.f;
+
+#endif
+}
+
 #endif /* SWIFT_ADIABATIC_INDEX_H */

src/const.h

@@ -37,8 +37,8 @@
 #define const_max_u_change 0.1f
 
 /* Dimensionality of the problem */
-#define HYDRO_DIMENSION_3D
-//#define HYDRO_DIMENSION_2D
+//#define HYDRO_DIMENSION_3D
+#define HYDRO_DIMENSION_2D
 //#define HYDRO_DIMENSION_1D
 
 /* Hydrodynamical adiabatic index. */

src/hydro.h

@@ -37,7 +37,7 @@
 #elif defined(HOPKINS_PE_SPH)
 #include "./hydro/PressureEntropy/hydro.h"
 #include "./hydro/PressureEntropy/hydro_iact.h"
-#define SPH_IMPLEMENTATION "Pressure-Entropy formulation of SPH (Hopkins 2013)"
+#define SPH_IMPLEMENTATION "Pressure-Entropy SPH (Hopkins 2013)"
 #elif defined(DEFAULT_SPH)
 #include "./hydro/Default/hydro.h"
 #include "./hydro/Default/hydro_iact.h"

src/hydro/PressureEntropy/hydro.h

@@ -30,6 +30,7 @@
  */
 
 #include "adiabatic_index.h"
+#include "approx_math.h"
 #include "dimension.h"
 #include "equation_of_state.h"
 #include "hydro_properties.h"
@@ -170,11 +171,13 @@ __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.div_v = 0.f;
-  p->density.rot_v[0] = 0.f;
-  p->density.rot_v[1] = 0.f;
-  p->density.rot_v[2] = 0.f;
+  p->weightedPressure = 0.f;
+  p->density.weightedPressure_dh = 0.f;
+  // p->rho_dh = 0.f;
+  /* p->density.div_v = 0.f; */
+  /* p->density.rot_v[0] = 0.f; */
+  /* p->density.rot_v[1] = 0.f; */
+  /* p->density.rot_v[2] = 0.f; */
 }
 
 /**
@@ -191,33 +194,40 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
 
   /* Some smoothing length multiples. */
   const float h = p->h;
-  const float h_inv = 1.0f / h;                       /* 1/h */
-  const float h_inv_dim = pow_dimension(h_inv);       /* 1/h^d */
-  const float h_inv_dim_plus_one = h_inv_dim * h_inv; /* 1/h^(d+1) */
+  const float h_inv = 1.0f / h;                 /* 1/h */
+  const float h_inv_dim = pow_dimension(h_inv); /* 1/h^d */
+  // const float h_inv_dim_plus_one = h_inv_dim * h_inv; /* 1/h^(d+1) */
 
   /* 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.wcount += kernel_root;
+  // p->density.wcount_dh -= hydro_dimension * kernel_root;
+  p->weightedPressure += p->mass * pow_one_over_gamma(p->entropy) * kernel_root;
+  p->density.weightedPressure_dh -=
+      hydro_dimension * p->mass * pow_one_over_gamma(p->entropy) * 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.wcount *= kernel_norm;
   p->density.wcount_dh *= h_inv * kernel_gamma * kernel_norm;
+  p->weightedPressure *= h_inv_dim;
+  p->density.weightedPressure_dh *= h_inv * kernel_gamma * kernel_norm;
 
-  const float irho = 1.f / p->rho;
+  /* Final operation on the weighted pressure */
+  p->weightedPressure = pow_gamma(p->weightedPressure);
+
+  /* const float irho = 1.f / p->rho; */
 
   /* Compute the derivative term */
-  p->rho_dh = 1.f / (1.f + hydro_dimension_inv * p->h * p->rho_dh * irho);
+  // p->rho_dh = 1.f / (1.f + hydro_dimension_inv * p->h * p->rho_dh * irho);
 
   /* Finish calculation of the velocity curl components */
-  p->density.rot_v[0] *= h_inv_dim_plus_one * irho;
-  p->density.rot_v[1] *= h_inv_dim_plus_one * irho;
-  p->density.rot_v[2] *= h_inv_dim_plus_one * irho;
+  // p->density.rot_v[0] *= h_inv_dim_plus_one * irho;
+  // p->density.rot_v[1] *= h_inv_dim_plus_one * irho;
+  // p->density.rot_v[2] *= h_inv_dim_plus_one * irho;
 
   /* Finish calculation of the velocity divergence */
-  p->density.div_v *= h_inv_dim_plus_one * irho;
+  // p->density.div_v *= h_inv_dim_plus_one * irho;
 }
 
 /**
@@ -234,36 +244,27 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
     struct part *restrict p, struct xpart *restrict xp, int ti_current,
     double timeBase) {
 
-  const float fac_mu = 1.f; /* Will change with cosmological integration */
-
-  /* Compute the norm of the curl */
-  const float curl_v = sqrtf(p->density.rot_v[0] * p->density.rot_v[0] +
-                             p->density.rot_v[1] * p->density.rot_v[1] +
-                             p->density.rot_v[2] * p->density.rot_v[2]);
-
-  /* Compute the norm of div v */
-  const float abs_div_v = fabsf(p->density.div_v);
-
   /* 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);
 
+  /* Compute the sound speed from the actual pressure*/
+  const float pressure = hydro_get_pressure(p, half_dt);
   const float irho = 1.f / p->rho;
-
-  /* Divide the pressure by the density and density gradient */
-  const float P_over_rho2 = pressure * irho * irho * p->rho_dh;
-
-  /* Compute the sound speed */
   const float soundspeed = sqrtf(hydro_gamma * pressure * irho);
 
-  /* Compute the Balsara switch */
-  const float balsara =
-      abs_div_v / (abs_div_v + curl_v + 0.0001f * soundspeed / fac_mu / p->h);
+  /* Compute the "i" part of the f_ij term */
+  const float number_density = p->density.wcount / kernel_norm;
+  const float factor = p->h / (hydro_dimension * number_density);
+  const float f_ij =
+      factor * p->density.weightedPressure_dh / (1.f + factor * number_density);
+
+  /* Comput the pressure term */
+  const float pressure_term = pow_one_minus_two_over_gamma(p->weightedPressure);
 
   /* Update variables. */
-  p->force.P_over_rho2 = P_over_rho2;
   p->force.soundspeed = soundspeed;
-  p->force.balsara = balsara;
+  p->force.f_ij = f_ij;
+  p->force.pressure_term = pressure_term;
 }
 
 /**
@@ -309,15 +310,20 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
 
   const float irho = 1.f / p->rho;
 
-  /* Divide the pressure by the density and density gradient */
-  const float P_over_rho2 = pressure * irho * irho * p->rho_dh;
-
   /* Compute the new sound speed */
   const float soundspeed = sqrtf(hydro_gamma * pressure * irho);
 
+  const float w1 = -hydro_dimension * p->force.h_dt * dt_entr / p->h;
+  if (fabsf(w1) < 0.2f)
+    p->weightedPressure *= approx_expf(w1); /* 4th order expansion of exp(w) */
+  else
+    p->weightedPressure *= expf(w1);
+
+  const float pressure_term = pow_one_minus_two_over_gamma(p->weightedPressure);
+
   /* Update variables */
-  p->force.P_over_rho2 = P_over_rho2;
   p->force.soundspeed = soundspeed;
+  p->force.pressure_term = pressure_term;
 }
 
 /**

src/hydro/PressureEntropy/hydro_debug.h

@@ -36,14 +36,12 @@ __attribute__((always_inline)) INLINE static void hydro_debug_particle(
       "v=[%.3e,%.3e,%.3e],v_full=[%.3e,%.3e,%.3e] \n a=[%.3e,%.3e,%.3e],\n "
       "h=%.3e, wcount=%.3f, wcount_dh=%.3e, m=%.3e, dh_drho=%.3e, rho=%.3e, "
       "P=%.3e, P_over_rho2=%.3e, S=%.3e, dS/dt=%.3e, c=%.3e\n"
-      "divV=%.3e, rotV=[%.3e,%.3e,%.3e], balsara=%.3e \n "
       "v_sig=%e dh/dt=%.3e t_begin=%d, t_end=%d\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, p->density.wcount, p->density.wcount_dh, p->mass, p->rho_dh, p->rho,
-      hydro_get_pressure(p, 0.), p->force.P_over_rho2, p->entropy,
-      p->entropy_dt, p->force.soundspeed, p->density.div_v, p->density.rot_v[0],
-      p->density.rot_v[1], p->density.rot_v[2], p->force.balsara,
+      p->h, p->density.wcount, p->density.wcount_dh, p->mass, 0.f /*p->rho_dh*/,
+      p->rho, hydro_get_pressure(p, 0.), 0.f,
+      /*p->force.P_over_rho2,*/ p->entropy, p->entropy_dt, p->force.soundspeed,
       p->force.v_sig, p->force.h_dt, p->ti_begin, p->ti_end);
 }
 

src/hydro/PressureEntropy/hydro_iact.h

@@ -37,15 +37,19 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
 
   float wi, wi_dx;
   float wj, wj_dx;
-  float dv[3], curlvr[3];
+  /* float dv[3], curlvr[3]; */
 
   /* Get the masses. */
   const float mi = pi->mass;
   const float mj = pj->mass;
 
+  /* Get the entropies. */
+  const float entropy_i = pi->entropy;
+  const float entropy_j = pj->entropy;
+
   /* Get r and r inverse. */
   const float r = sqrtf(r2);
-  const float r_inv = 1.0f / r;
+  /* const float r_inv = 1.0f / r; */
 
   /* Compute the kernel function for pi */
   const float hi_inv = 1.f / hi;
@@ -54,11 +58,15 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
 
   /* Compute contribution to the density */
   pi->rho += mj * wi;
-  pi->rho_dh -= mj * (hydro_dimension * wi + ui * wi_dx);
+
+  /* Weighted pressure */
+  pi->weightedPressure += mj * pow_one_over_gamma(pj->entropy) * wi;
+  pi->density.weightedPressure_dh -=
+      mj * pow_one_over_gamma(entropy_j) * (hydro_dimension * wi + ui * wi_dx);
 
   /* Compute contribution to the number of neighbours */
   pi->density.wcount += wi;
-  pi->density.wcount_dh -= ui * wi_dx;
+  pi->density.wcount_dh -= ui * wi_dx;  //(hydro_dimension * wi + ui * wi_dx);
 
   /* Compute the kernel function for pj */
   const float hj_inv = 1.f / hj;
@@ -67,36 +75,41 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
 
   /* Compute contribution to the density */
   pj->rho += mi * wj;
-  pj->rho_dh -= mi * (hydro_dimension * wj + uj * wj_dx);
+  // pj->rho_dh -= mi * (hydro_dimension * wj + uj * wj_dx);
+
+  /* Weighted pressure */
+  pj->weightedPressure += mi * pow_one_over_gamma(pi->entropy) * wj;
+  pj->density.weightedPressure_dh -=
+      mi * pow_one_over_gamma(entropy_i) * (hydro_dimension * wj + uj * wj_dx);
 
   /* Compute contribution to the number of neighbours */
   pj->density.wcount += wj;
-  pj->density.wcount_dh -= uj * wj_dx;
+  pj->density.wcount_dh -= uj * wj_dx;  //(hydro_dimension * wj + uj * wj_dx);
 
-  const float faci = mj * wi_dx * r_inv;
-  const float facj = mi * wj_dx * r_inv;
+  /* const float faci = mj * wi_dx * r_inv; */
+  /* const float facj = mi * wj_dx * r_inv; */
 
-  /* Compute dv dot r */
-  dv[0] = pi->v[0] - pj->v[0];
-  dv[1] = pi->v[1] - pj->v[1];
-  dv[2] = pi->v[2] - pj->v[2];
-  const float dvdr = dv[0] * dx[0] + dv[1] * dx[1] + dv[2] * dx[2];
+  /* /\* Compute dv dot r *\/ */
+  /* dv[0] = pi->v[0] - pj->v[0]; */
+  /* dv[1] = pi->v[1] - pj->v[1]; */
+  /* dv[2] = pi->v[2] - pj->v[2]; */
+  /* const float dvdr = dv[0] * dx[0] + dv[1] * dx[1] + dv[2] * dx[2]; */
 
-  pi->density.div_v -= faci * dvdr;
-  pj->density.div_v -= facj * dvdr;
+  /* pi->density.div_v -= faci * dvdr; */
+  /* pj->density.div_v -= facj * dvdr; */
 
-  /* Compute dv cross r */
-  curlvr[0] = dv[1] * dx[2] - dv[2] * dx[1];
-  curlvr[1] = dv[2] * dx[0] - dv[0] * dx[2];
-  curlvr[2] = dv[0] * dx[1] - dv[1] * dx[0];
+  /* /\* Compute dv cross r *\/ */
+  /* curlvr[0] = dv[1] * dx[2] - dv[2] * dx[1]; */
+  /* curlvr[1] = dv[2] * dx[0] - dv[0] * dx[2]; */
+  /* curlvr[2] = dv[0] * dx[1] - dv[1] * dx[0]; */
 
-  pi->density.rot_v[0] += faci * curlvr[0];
-  pi->density.rot_v[1] += faci * curlvr[1];
-  pi->density.rot_v[2] += faci * curlvr[2];
+  /* pi->density.rot_v[0] += faci * curlvr[0]; */
+  /* pi->density.rot_v[1] += faci * curlvr[1]; */
+  /* pi->density.rot_v[2] += faci * curlvr[2]; */
 
-  pj->density.rot_v[0] += facj * curlvr[0];
-  pj->density.rot_v[1] += facj * curlvr[1];
-  pj->density.rot_v[2] += facj * curlvr[2];
+  /* pj->density.rot_v[0] += facj * curlvr[0]; */
+  /* pj->density.rot_v[1] += facj * curlvr[1]; */
+  /* pj->density.rot_v[2] += facj * curlvr[2]; */
 }
 
 /**
@@ -116,14 +129,17 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
     float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {
 
   float wi, wi_dx;
-  float dv[3], curlvr[3];
+  /* float dv[3], curlvr[3]; */
 
   /* Get the masses. */
   const float mj = pj->mass;
 
+  /* Get the entropies. */
+  const float entropy_j = pj->entropy;
+
   /* Get r and r inverse. */
   const float r = sqrtf(r2);
-  const float ri = 1.0f / r;
+  /* const float ri = 1.0f / r; */
 
   /* Compute the kernel function */
   const float h_inv = 1.0f / hi;
@@ -132,29 +148,33 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
 
   /* Compute contribution to the density */
   pi->rho += mj * wi;
-  pi->rho_dh -= mj * (hydro_dimension * wi + u * wi_dx);
+
+  /* Weighted pressure */
+  pi->weightedPressure += mj * pow_one_over_gamma(pj->entropy) * wi;
+  pi->density.weightedPressure_dh -=
+      mj * pow_one_over_gamma(entropy_j) * (hydro_dimension * wi + u * wi_dx);
 
   /* Compute contribution to the number of neighbours */
   pi->density.wcount += wi;
-  pi->density.wcount_dh -= u * wi_dx;
+  pi->density.wcount_dh -= u * wi_dx;  //(hydro_dimension * wi + u * wi_dx);
 
-  const float fac = mj * wi_dx * ri;
+  /* const float fac = mj * wi_dx * ri; */
 
-  /* Compute dv dot r */
-  dv[0] = pi->v[0] - pj->v[0];
-  dv[1] = pi->v[1] - pj->v[1];
-  dv[2] = pi->v[2] - pj->v[2];
-  const float dvdr = dv[0] * dx[0] + dv[1] * dx[1] + dv[2] * dx[2];
-  pi->density.div_v -= fac * dvdr;
+  /* /\* Compute dv dot r *\/ */
+  /* dv[0] = pi->v[0] - pj->v[0]; */
+  /* dv[1] = pi->v[1] - pj->v[1]; */
+  /* dv[2] = pi->v[2] - pj->v[2]; */
+  /* const float dvdr = dv[0] * dx[0] + dv[1] * dx[1] + dv[2] * dx[2]; */
+  /* pi->density.div_v -= fac * dvdr; */
 
-  /* Compute dv cross r */
-  curlvr[0] = dv[1] * dx[2] - dv[2] * dx[1];
-  curlvr[1] = dv[2] * dx[0] - dv[0] * dx[2];
-  curlvr[2] = dv[0] * dx[1] - dv[1] * dx[0];
+  /* /\* Compute dv cross r *\/ */
+  /* curlvr[0] = dv[1] * dx[2] - dv[2] * dx[1]; */
+  /* curlvr[1] = dv[2] * dx[0] - dv[0] * dx[2]; */
+  /* curlvr[2] = dv[0] * dx[1] - dv[1] * dx[0]; */
 
-  pi->density.rot_v[0] += fac * curlvr[0];
-  pi->density.rot_v[1] += fac * curlvr[1];
-  pi->density.rot_v[2] += fac * curlvr[2];
+  /* pi->density.rot_v[0] += fac * curlvr[0]; */
+  /* pi->density.rot_v[1] += fac * curlvr[1]; */
+  /* pi->density.rot_v[2] += fac * curlvr[2]; */
 }
 
 /**
@@ -185,6 +205,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
   const float mj = pj->mass;
   const float rhoi = pi->rho;
   const float rhoj = pj->rho;
+  const float entropy_i = pi->entropy;
+  const float entropy_j = pj->entropy;
+  const float entropy_product = pow_one_over_gamma(entropy_i * entropy_j);
 
   /* Get the kernel for hi. */
   const float hi_inv = 1.0f / hi;
@@ -201,8 +224,10 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
   const float wj_dr = hjd_inv * wj_dx;
 
   /* Compute gradient terms */
-  const float P_over_rho2_i = pi->force.P_over_rho2;
-  const float P_over_rho2_j = pj->force.P_over_rho2;
+  const float f_ij =
+      1.f;  // - pi->force.f_ij / (mj * pow_one_over_gamma(entropy_j));
+  const float f_ji =
+      1.f;  // - pj->force.f_ij / (mi * pow_one_over_gamma(entropy_i));
 
   /* Compute sound speeds */
   const float ci = pi->force.soundspeed;
@@ -214,8 +239,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
                      (pi->v[2] - pj->v[2]) * dx[2];
 
   /* Balsara term */
-  const float balsara_i = pi->force.balsara;
-  const float balsara_j = pj->force.balsara;
+  /* const float balsara_i = pi->force.balsara; */
+  /* const float balsara_j = pj->force.balsara; */
 
   /* Are the particles moving towards each others ? */
   const float omega_ij = fminf(dvdr, 0.f);
@@ -226,13 +251,22 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
 
   /* Now construct the full viscosity term */
   const float rho_ij = 0.5f * (rhoi + rhoj);
-  const float visc = -0.25f * const_viscosity_alpha * v_sig * mu_ij *
-                     (balsara_i + balsara_j) / rho_ij;
+  const float visc = -0.5f * const_viscosity_alpha * v_sig * mu_ij / rho_ij;
 
   /* Now, convolve with the kernel */
   const float visc_term = 0.5f * visc * (wi_dr + wj_dr) * r_inv;
-  const float sph_term =
-      (P_over_rho2_i * wi_dr + P_over_rho2_j * wj_dr) * r_inv;
+  const float sph_term = entropy_product *
+                         (f_ij * pi->force.pressure_term * wi_dr +
+                          f_ji * pj->force.pressure_term * wj_dr) *
+                         r_inv;
+
+  /* if(pi->entropy != pj->entropy) { */
+  /* message("Si = %f Pi = %f, Pi_term = %f", pi->entropy, pi->weightedPressure,
+   * pi->force.pressure_term); */
+  /* message("Sj = %f Pj = %f, Pj_term = %f", pj->entropy, pj->weightedPressure,
+   * pj->force.pressure_term); */
+  /* error("done"); */
+  /* } */
 
   /* Eventually got the acceleration */
   const float acc = visc_term + sph_term;
@@ -287,6 +321,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
   const float mj = pj->mass;
   const float rhoi = pi->rho;
   const float rhoj = pj->rho;
+  const float entropy_i = pi->entropy;
+  const float entropy_j = pj->entropy;
+  const float entropy_product = pow_one_over_gamma(entropy_i * entropy_j);
 
   /* Get the kernel for hi. */
   const float hi_inv = 1.0f / hi;
@@ -303,8 +340,10 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
   const float wj_dr = hjd_inv * wj_dx;
 
   /* Compute gradient terms */
-  const float P_over_rho2_i = pi->force.P_over_rho2;
-  const float P_over_rho2_j = pj->force.P_over_rho2;
+  const float f_ij =
+      1.f;  // - pi->force.f_ij / (mj * pow_one_over_gamma(entropy_j));
+  const float f_ji =
+      1.f;  // - pj->force.f_ij / (mi * pow_one_over_gamma(entropy_i));
 
   /* Compute sound speeds */
   const float ci = pi->force.soundspeed;
@@ -316,8 +355,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
                      (pi->v[2] - pj->v[2]) * dx[2];
 
   /* Balsara term */
-  const float balsara_i = pi->force.balsara;
-  const float balsara_j = pj->force.balsara;
+  /* const float balsara_i = pi->force.balsara; */
+  /* const float balsara_j = pj->force.balsara; */
 
   /* Are the particles moving towards each others ? */
   const float omega_ij = fminf(dvdr, 0.f);
@@ -328,13 +367,14 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
 
   /* Now construct the full viscosity term */
   const float rho_ij = 0.5f * (rhoi + rhoj);
-  const float visc = -0.25f * const_viscosity_alpha * v_sig * mu_ij *
-                     (balsara_i + balsara_j) / rho_ij;
+  const float visc = -0.5f * const_viscosity_alpha * v_sig * mu_ij / rho_ij;
 
   /* Now, convolve with the kernel */
   const float visc_term = 0.5f * visc * (wi_dr + wj_dr) * r_inv;
-  const float sph_term =
-      (P_over_rho2_i * wi_dr + P_over_rho2_j * wj_dr) * r_inv;
+  const float sph_term = entropy_product *
+                         (f_ij * pi->force.pressure_term * wi_dr +
+                          f_ji * pj->force.pressure_term * wj_dr) *
+                         r_inv;
 
   /* Eventually got the acceleration */
   const float acc = visc_term + sph_term;

src/hydro/PressureEntropy/hydro_io.h

@@ -55,8 +55,12 @@ void hydro_read_particles(struct part* parts, struct io_props* list,
                                 parts, mass);
   list[3] = io_make_input_field("SmoothingLength", FLOAT, 1, COMPULSORY,
                                 UNIT_CONV_LENGTH, parts, h);
-  list[4] = io_make_input_field("InternalEnergy", FLOAT, 1, COMPULSORY,
-                                UNIT_CONV_ENERGY_PER_UNIT_MASS, parts, entropy);
+  /* list[4] = io_make_input_field("InternalEnergy", FLOAT, 1, COMPULSORY, */
+  /*                               UNIT_CONV_ENERGY_PER_UNIT_MASS, parts,
+   * entropy); */
+  list[4] =
+      io_make_input_field("InternalEnergy", FLOAT, 1, COMPULSORY,
+                          UNIT_CONV_ENTROPY_PER_UNIT_MASS, parts, entropy);
   list[5] = io_make_input_field("ParticleIDs", ULONGLONG, 1, COMPULSORY,
                                 UNIT_CONV_NO_UNITS, parts, id);
   list[6] = io_make_input_field("Accelerations", FLOAT, 3, OPTIONAL,
@@ -85,7 +89,7 @@ float convert_P(struct engine* e, struct part* p) {
 void hydro_write_particles(struct part* parts, struct io_props* list,
                            int* num_fields) {
 
-  *num_fields = 10;
+  *num_fields = 11;
 
   /* List what we want to write */
   list[0] = io_make_output_field("Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH,
@@ -105,10 +109,12 @@ void hydro_write_particles(struct part* parts, struct io_props* list,
                                  UNIT_CONV_ACCELERATION, parts, a_hydro);
   list[7] =
       io_make_output_field("Density", FLOAT, 1, UNIT_CONV_DENSITY, parts, rho);
-  list[8] = io_make_output_field("Entropy", FLOAT, 1,
-                                 UNIT_CONV_ENTROPY_PER_UNIT_MASS, parts, rho);
+  list[8] = io_make_output_field(
+      "Entropy", FLOAT, 1, UNIT_CONV_ENTROPY_PER_UNIT_MASS, parts, entropy);
   list[9] = io_make_output_field_convert_part(
       "Pressure", FLOAT, 1, UNIT_CONV_PRESSURE, parts, rho, convert_P);
+  list[10] = io_make_output_field("WeightedPressure", FLOAT, 1,
+                                  UNIT_CONV_PRESSURE, parts, weightedPressure);
 }