From b4243314ef8abf94a878ab32c7901758e9cc2b8c Mon Sep 17 00:00:00 2001
From: boson112358 <13271927606@163.com>
Date: Fri, 1 Nov 2024 11:40:43 +0000
Subject: [PATCH] correct the mu_j term and add scale factor in balsara
---
src/hydro/REMIX/hydro.h | 16 ++++++++++------
src/hydro/REMIX/hydro_visc_difn.h | 2 +-
2 files changed, 11 insertions(+), 7 deletions(-)
diff --git a/src/hydro/REMIX/hydro.h b/src/hydro/REMIX/hydro.h
index 13ebc90561..0a10d8b30a 100644
--- a/src/hydro/REMIX/hydro.h
+++ b/src/hydro/REMIX/hydro.h
@@ -669,12 +669,16 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
const float soundspeed =
gas_soundspeed_from_internal_energy(p->rho_evol, p->u);
- float div_v = p->dv_norm_kernel[0][0] + p->dv_norm_kernel[1][1] +
- p->dv_norm_kernel[2][2];
+ /* Compute balsara switch including scale factor*/
+ const float a_inv2 = cosmo->a2_inv;
+ const float fac_B = cosmo->a_factor_Balsara_eps;
+
+ float div_v = (p->dv_norm_kernel[0][0] + p->dv_norm_kernel[1][1] +
+ p->dv_norm_kernel[2][2]) * a_inv2 + cosmo->H * hydro_dimension;
float curl_v[3];
- curl_v[0] = p->dv_norm_kernel[1][2] - p->dv_norm_kernel[2][1];
- curl_v[1] = p->dv_norm_kernel[2][0] - p->dv_norm_kernel[0][2];
- curl_v[2] = p->dv_norm_kernel[0][1] - p->dv_norm_kernel[1][0];
+ curl_v[0] = (p->dv_norm_kernel[1][2] - p->dv_norm_kernel[2][1]) * a_inv2;
+ curl_v[1] = (p->dv_norm_kernel[2][0] - p->dv_norm_kernel[0][2]) * a_inv2;
+ curl_v[2] = (p->dv_norm_kernel[0][1] - p->dv_norm_kernel[1][0]) * a_inv2;
float mod_curl_v = sqrtf(curl_v[0] * curl_v[0] + curl_v[1] * curl_v[1] +
curl_v[2] * curl_v[2]);
@@ -684,7 +688,7 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
balsara = 0.f;
} else {
balsara = fabsf(div_v) /
- (fabsf(div_v) + mod_curl_v + 0.0001f * soundspeed / p->h);
+ (fabsf(div_v) + mod_curl_v + 0.0001f * soundspeed * fac_B / p->h);
}
/* Compute the pressure */
diff --git a/src/hydro/REMIX/hydro_visc_difn.h b/src/hydro/REMIX/hydro_visc_difn.h
index 2979368c37..27e08cd191 100644
--- a/src/hydro/REMIX/hydro_visc_difn.h
+++ b/src/hydro/REMIX/hydro_visc_difn.h
@@ -305,7 +305,7 @@ __attribute__((always_inline)) INLINE static void hydro_set_Qi_Qj(
float mu_j =
fac_mu * min(0.f, ((vtilde_i[0] - vtilde_j[0]) * dx[0] +
(vtilde_i[1] - vtilde_j[1]) * dx[1] +
- (vtilde_i[2] - vtilde_j[2]) * dx[2]) *
+ (vtilde_i[2] - vtilde_j[2]) * dx[2] + a2_Hubble * r2) *
hj_inv / (r2 * hj_inv * hj_inv + epsilon * epsilon));
const float ci = pi->force.soundspeed;
--
GitLab