Moved Balsara switch to shock indicator - Helps considerably with some MHD shock tubes

src/hydro/SPHENIX/hydro.h

@@ -757,30 +757,6 @@ __attribute__((always_inline)) INLINE static void hydro_end_gradient(
   /* Include the extra factors in the del^2 u */
 
   p->diffusion.laplace_u *= 2.f * h_inv_dim_plus_one;
-
-  /* Compute Balsara switch */
-  const float fac_B = cosmo->a_factor_Balsara_eps;
-
-  /* Compute the norm of the curl */
-  const float curl_v = sqrtf(p->viscosity.rot_v[0] * p->viscosity.rot_v[0] +
-                             p->viscosity.rot_v[1] * p->viscosity.rot_v[1] +
-                             p->viscosity.rot_v[2] * p->viscosity.rot_v[2]);
-
-  /* Compute the norm of div v */
-  const float abs_div_v = fabsf(p->viscosity.div_v);
-
-  /* Compute the sound speed  */
-  const float pressure = hydro_get_comoving_pressure(p);
-  const float pressure_including_floor =
-      pressure_floor_get_comoving_pressure(p, pressure_floor, pressure, cosmo);
-  const float soundspeed =
-      gas_soundspeed_from_pressure(p->rho, pressure_including_floor);
-
-  /* Compute the Balsara switch */
-  const float balsara =
-      abs_div_v / (abs_div_v + curl_v + 0.0001f * soundspeed * fac_B / p->h);
-
-  p->force.balsara = balsara;
   
 #ifdef SWIFT_HYDRO_DENSITY_CHECKS
   p->n_gradient += kernel_root;
@@ -865,13 +841,39 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
   const float soundspeed_physical =
       gas_soundspeed_from_pressure(p->rho, pressure_including_floor) *
       cosmo->a_factor_sound_speed;
-
+  const float fac_B = cosmo->a_factor_Balsara_eps;
+  
   const float sound_crossing_time_inverse =
       soundspeed_physical * kernel_support_physical_inv;
 
+  /* Compute the norm of the curl */
+  const float curl_v = sqrtf(p->viscosity.rot_v[0] * p->viscosity.rot_v[0] +
+                             p->viscosity.rot_v[1] * p->viscosity.rot_v[1] +
+                             p->viscosity.rot_v[2] * p->viscosity.rot_v[2]);
+
+  /* Compute the norm of div v */
+  const float abs_div_v = fabsf(p->viscosity.div_v);
+
+  /* Compute the sound speed  */
+  const float pressure = hydro_get_comoving_pressure(p);
+  const float pressure_including_floor =
+      pressure_floor_get_comoving_pressure(p, pressure_floor, pressure, cosmo);
+  const float soundspeed =
+      gas_soundspeed_from_pressure(p->rho, pressure_including_floor);
+
+  /* Get the squares of the quantities necessary for the Balsara-like switch */
+  const float fac_B_2 = fac_B * fac_B;
+  const float curl_v_2 = curl_v * curl_v;
+  const float abs_div_v_2 = abs_div_v * abs_div_v;
+  const float soundspeed_2 = soundspeed * soundspeed;
+  const float h_2 = p->h * p->h;
+  
+  /* Compute the Balsara-like switch (Price et a. (2018), eq. 47; a simplified version of eq. 18 in Cullen & Dehnen (2012)) */
+  const float balsara =
+      abs_div_v_2 / (abs_div_v_2 + curl_v_2 + 0.0001f * soundspeed_2 * fac_B_2 / h_2);
+  
   /* Construct time differential of div.v implicitly following the ANARCHY spec
    */
-
   const float div_v_dt =
       dt_alpha == 0.f
           ? 0.f
@@ -882,7 +884,7 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
   /* Source term is only activated if flow is converging (i.e. in the pre-
    * shock region) */
   const float S = p->viscosity.div_v < 0.f
-                      ? kernel_support_physical * kernel_support_physical *
+                      ? balsara * kernel_support_physical * kernel_support_physical *
                             max(0.f, -1.f * div_v_dt)
                       : 0.f;
 

src/hydro/SPHENIX/hydro_csds.h

@@ -77,7 +77,7 @@ INLINE static void *csds_hydro_convert_secondary(const struct part *p,
   const float secondary[7] = {
       hydro_get_comoving_entropy(p, xp),
       hydro_get_comoving_pressure(p),
-      p->viscosity.alpha * p->force.balsara,
+      p->viscosity.alpha,
       p->diffusion.alpha,
       p->diffusion.laplace_u,
       p->viscosity.div_v,

src/hydro/SPHENIX/hydro_iact.h

@@ -428,15 +428,11 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
   const float f_ij = 1.f - pi->force.f / mj;
   const float f_ji = 1.f - pj->force.f / mi;
 
-  /* Balsara term */
-  const float balsara_i = pi->force.balsara;
-  const float balsara_j = pj->force.balsara;
-
   /* Construct the full viscosity term */
   const float rho_ij = rhoi + rhoj;
   const float alpha = pi->viscosity.alpha + pj->viscosity.alpha;
   const float visc =
-      -0.25f * alpha * v_sig * mu_ij * (balsara_i + balsara_j) / rho_ij;
+      -0.5f * alpha * v_sig * mu_ij / rho_ij;
 
   /* Convolve with the kernel */
   const float visc_acc_term =
@@ -580,15 +576,11 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
   const float f_ij = 1.f - pi->force.f / mj;
   const float f_ji = 1.f - pj->force.f / mi;
 
-  /* Balsara term */
-  const float balsara_i = pi->force.balsara;
-  const float balsara_j = pj->force.balsara;
-
   /* Construct the full viscosity term */
   const float rho_ij = rhoi + rhoj;
   const float alpha = pi->viscosity.alpha + pj->viscosity.alpha;
   const float visc =
-      -0.25f * alpha * v_sig * mu_ij * (balsara_i + balsara_j) / rho_ij;
+      -0.5f * alpha * v_sig * mu_ij / rho_ij;
 
   /* Convolve with the kernel */
   const float visc_acc_term =

src/hydro/SPHENIX/hydro_io.h

@@ -171,7 +171,7 @@ INLINE static void convert_part_softening(const struct engine* e,
 INLINE static void convert_viscosity(const struct engine* e,
                                      const struct part* p,
                                      const struct xpart* xp, float* ret) {
-  ret[0] = p->viscosity.alpha * p->force.balsara;
+  ret[0] = p->viscosity.alpha;
 }
 
 INLINE static void convert_diffusion(const struct engine* e,
@@ -236,8 +236,7 @@ INLINE static void hydro_write_particles(const struct part* parts,
   list[9] = io_make_output_field_convert_part(
       "ViscosityParameters", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f, parts, xparts,
       convert_viscosity,
-      "Visosity coefficient (alpha_visc) of the particles, multiplied by the "
-      "balsara switch");
+      "Visosity coefficient (alpha_visc) of the particles");
 
   list[10] = io_make_output_field_convert_part(
       "DiffusionParameters", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f, parts, xparts,

src/hydro/SPHENIX/hydro_part.h

@@ -208,15 +208,11 @@ struct part {
     /*! Time derivative of smoothing length  */
     float h_dt;
 
-    /*! Balsara switch */
-    float balsara;
-
     /*! Maximal alpha (viscosity) over neighbours */
     float alpha_visc_max_ngb;
 
   } force;
-  //};
-
+  
   /*! Additional data used for adaptive softening */
   struct adaptive_softening_part_data adaptive_softening_data;