diff --git a/src/hydro/SPHENIX/hydro.h b/src/hydro/SPHENIX/hydro.h index 8a6a55f16bffbd4ad89cfa3d365914368b441efa..dc6826c49e738bf7456aa7c0ab624b3dd27b8b5b 100644 --- a/src/hydro/SPHENIX/hydro.h +++ b/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; diff --git a/src/hydro/SPHENIX/hydro_csds.h b/src/hydro/SPHENIX/hydro_csds.h index e593ab97b026b0fd7ca9cd9ceb5737e3ff562b61..8d81c99c3009a638cf88293fff7aebdb45f91322 100644 --- a/src/hydro/SPHENIX/hydro_csds.h +++ b/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, diff --git a/src/hydro/SPHENIX/hydro_iact.h b/src/hydro/SPHENIX/hydro_iact.h index 285a92c9eac9e76c333239cdf41c13bd73daab16..4a96f22c5123d2f23b7e13bf482f89d242929ba7 100644 --- a/src/hydro/SPHENIX/hydro_iact.h +++ b/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 = diff --git a/src/hydro/SPHENIX/hydro_io.h b/src/hydro/SPHENIX/hydro_io.h index d123f6d1f703bdb45020f04e34b3a6c85b55d417..211ef38e9968a2737d3aa3d743b7701c27922fa5 100644 --- a/src/hydro/SPHENIX/hydro_io.h +++ b/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, diff --git a/src/hydro/SPHENIX/hydro_part.h b/src/hydro/SPHENIX/hydro_part.h index af582a6cef9b4a371bc501866008db5a1f1bcc9f..859176e34823793d2e54c2646e494e5c71351bb5 100644 --- a/src/hydro/SPHENIX/hydro_part.h +++ b/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;