Magnetosonic speed in timescale used in updating the viscosity alphas

src/hydro/SPHENIX/hydro.h

@@ -37,6 +37,7 @@
 #include "hydro_properties.h"
 #include "hydro_space.h"
 #include "kernel_hydro.h"
+#include "mhd.h"
 #include "minmax.h"
 #include "pressure_floor.h"
 
@@ -840,7 +841,8 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
       cosmo->a_factor_sound_speed;
 
   const float sound_crossing_time_inverse =
-      soundspeed_physical * kernel_support_physical_inv;
+    mhd_get_magnetosonic_speed(p,1.0f,1.0f) * kernel_support_physical_inv;  
+    // soundspeed_physical * kernel_support_physical_inv;
 
   /* Construct time differential of div.v implicitly following the ANARCHY spec
    */

src/hydro/SPHENIX/hydro_parameters.h

@@ -50,7 +50,7 @@
 /*! Cosmology default beta=3.0.
  * Alpha can be set in the parameter file.
  * Beta is defined as in e.g. Price (2010) Eqn (103) */
-#define const_viscosity_beta 1.5f
+#define const_viscosity_beta 3.0f
 
 /*! The viscosity that the particles are reset to after being hit by a
  * feedback event. This should be set to the same value as the
@@ -62,7 +62,7 @@
 
 /*! The "initial" hydro viscosity, or the fixed value for non-variable
  * schemes. This usually takes the value 0.8. */
-#define hydro_props_default_viscosity_alpha 0.8f
+#define hydro_props_default_viscosity_alpha 0.1f
 
 /*! Minimal value for the viscosity alpha in variable schemes. */
 #define hydro_props_default_viscosity_alpha_min 0.0f

src/mhd/DirectInduction/mhd.h

@@ -411,7 +411,7 @@ __attribute__((always_inline)) INLINE static float mhd_get_psi_over_ch_dt(
   const float par = hydro_props->mhd.par_dedner;
 
   const float div_B = p->mhd_data.divB;
-  const float div_v = a * a * hydro_get_div_v(p) - 3.0f * a * a * H;
+  const float div_v = a * a * p->viscosity.div_v - 3.0f * a * a * H;
   const float psi_over_ch = p->mhd_data.psi_over_ch;
 
   const float hyperbolic_term =
@@ -576,7 +576,7 @@ __attribute__((always_inline)) INLINE static void mhd_convert_quantities(
   p->mhd_data.monopole_beta = hydro_props->mhd.monopole_subtraction;
   /* Set Artificial Difussion */
   p->mhd_data.art_diff_beta = hydro_props->mhd.art_diffusion;
-
+  
   /* Convert B into B/rho */
   p->mhd_data.B_over_rho[0] /= p->rho;
   p->mhd_data.B_over_rho[1] /= p->rho;
@@ -586,7 +586,7 @@ __attribute__((always_inline)) INLINE static void mhd_convert_quantities(
   p->mhd_data.B_over_rho[0] *= pow(cosmo->a, 1.5f * hydro_gamma);
   p->mhd_data.B_over_rho[1] *= pow(cosmo->a, 1.5f * hydro_gamma);
   p->mhd_data.B_over_rho[2] *= pow(cosmo->a, 1.5f * hydro_gamma);
-
+  
   xp->mhd_data.B_over_rho_full[0] = p->mhd_data.B_over_rho[0];
   xp->mhd_data.B_over_rho_full[1] = p->mhd_data.B_over_rho[1];
   xp->mhd_data.B_over_rho_full[2] = p->mhd_data.B_over_rho[2];

[Federico Andrés Stasyszyn]

breaks the cosmo runs. mhd_get.... is in comoving and soundspeed_physical is in physical

[Orestis Karapiperis]

Good catch, thanks