diff --git a/src/hydro/SPHENIX/hydro.h b/src/hydro/SPHENIX/hydro.h index c8de61ed33c102ebc2b4c5576528e78ed3585b47..ed158d7d10fb720d8c638b1bd76c759741ce11e4 100644 --- a/src/hydro/SPHENIX/hydro.h +++ b/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 */ diff --git a/src/hydro/SPHENIX/hydro_parameters.h b/src/hydro/SPHENIX/hydro_parameters.h index d77693a50e2b33c4e5917badbd937af50da4892e..06b77e9eecaebccea5fc444e21772025d961af87 100644 --- a/src/hydro/SPHENIX/hydro_parameters.h +++ b/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 */