Merge branch 'sphenix_minor_updates' into 'master'

Added updates for sphenix scheme (minor parameter changes)

See merge request !1079

src/hydro/SPHENIX/hydro.h

@@ -746,26 +746,33 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
 
   /* Construct the source term for the AV; if shock detected this is _positive_
    * as div_v_dt should be _negative_ before the shock hits */
-  const float S = kernel_support_physical * kernel_support_physical *
-                  max(0.f, -1.f * div_v_dt);
-  /* 0.25 factor comes from our definition of v_sig (sum of soundspeeds rather
-   * than mean). */
-  /* Note this is v_sig_physical squared, not comoving */
-  const float v_sig_square = 0.25 * v_sig_physical * v_sig_physical;
+  /* 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 *
+                            max(0.f, -1.f * div_v_dt)
+                      : 0.f;
+
+  /* We want the decay to occur based on the thermodynamic properties
+   * of the gas - hence use the soundspeed instead of the signal velocity */
+  const float soundspeed_square = soundspeed_physical * soundspeed_physical;
 
   /* Calculate the current appropriate value of the AV based on the above */
   const float alpha_loc =
-      hydro_props->viscosity.alpha_max * S / (v_sig_square + S);
+      hydro_props->viscosity.alpha_max * S / (soundspeed_square + S);
 
   if (alpha_loc > p->viscosity.alpha) {
     /* Reset the value of alpha to the appropriate value */
     p->viscosity.alpha = alpha_loc;
   } else {
     /* Integrate the alpha forward in time to decay back to alpha = alpha_loc */
-    p->viscosity.alpha =
-        alpha_loc + (p->viscosity.alpha - alpha_loc) *
-                        expf(-dt_alpha * sound_crossing_time_inverse *
-                             hydro_props->viscosity.length);
+    /* This type of integration is stable w.r.t. different time-step lengths
+     * (Price 2018) */
+    const float timescale_ratio =
+        dt_alpha * sound_crossing_time_inverse * hydro_props->viscosity.length;
+
+    p->viscosity.alpha += alpha_loc * timescale_ratio;
+    p->viscosity.alpha /= (1.f + timescale_ratio);
   }
 
   /* Check that we did not hit the minimum */

src/hydro/SPHENIX/hydro_parameters.h

@@ -71,7 +71,7 @@
 #define hydro_props_default_viscosity_alpha_max 2.0f
 
 /*! Decay length for the viscosity scheme. This is scheme dependent. */
-#define hydro_props_default_viscosity_length 0.25f
+#define hydro_props_default_viscosity_length 0.05f
 
 /* Diffusion parameters -- FIXED -- MUST BE DEFINED AT COMPILE-TIME */