diff --git a/src/hydro/PressureEnergy/hydro.h b/src/hydro/PressureEnergy/hydro.h
index cba76d2851440f0d34d9c565884d0f509a83c4e1..263b04913a740485755f4da776b64b4ed80562f2 100644
--- a/src/hydro/PressureEnergy/hydro.h
+++ b/src/hydro/PressureEnergy/hydro.h
@@ -312,6 +312,11 @@ __attribute__((always_inline)) INLINE static void hydro_init_part(
p->density.rho_dh = 0.f;
p->pressure_bar = 0.f;
p->density.pressure_bar_dh = 0.f;
+
+ p->density.div_v = 0.f;
+ p->density.rot_v[0] = 0.f;
+ p->density.rot_v[1] = 0.f;
+ p->density.rot_v[2] = 0.f;
}
/**
@@ -401,12 +406,26 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
struct part *restrict p, struct xpart *restrict xp,
const struct cosmology *cosmo) {
+ const float fac_mu = cosmo->a_factor_mu;
+
+ /* Compute the norm of the curl */
+ const float curl_v = sqrtf(p->density.rot_v[0] * p->density.rot_v[0] +
+ p->density.rot_v[1] * p->density.rot_v[1] +
+ p->density.rot_v[2] * p->density.rot_v[2]);
+
+ /* Compute the norm of div v */
+ const float abs_div_v = fabsf(p->density.div_v);
+
/* Compute the pressure */
const float pressure = gas_pressure_from_internal_energy(p->rho, p->u);
/* Compute the sound speed */
const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure);
+ /* Compute the Balsara switch */
+ const float balsara =
+ abs_div_v / (abs_div_v + curl_v + 0.0001f * soundspeed * fac_mu / p->h);
+
/* Compute the "grad h" term */
const float common_factor = p->h / (hydro_dimension * p->density.wcount);
const float grad_h_term =
@@ -417,6 +436,7 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
p->force.f = grad_h_term;
p->force.pressure = pressure;
p->force.soundspeed = soundspeed;
+ p->force.balsara = balsara;
}
/**
diff --git a/src/hydro/PressureEnergy/hydro_iact.h b/src/hydro/PressureEnergy/hydro_iact.h
index f2fbe67d05ad7037d80caa243721703adf324655..8ab6da92f22412759f62c815c660a18cea61c8c5 100644
--- a/src/hydro/PressureEnergy/hydro_iact.h
+++ b/src/hydro/PressureEnergy/hydro_iact.h
@@ -52,6 +52,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
struct part* pj, float a, float H) {
float wi, wj, wi_dx, wj_dx;
+ float dv[3], curlvr[3];
const float r = sqrtf(r2);
@@ -84,6 +85,33 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
mi * pi->u * (hydro_dimension * wj + uj * wj_dx);
pj->density.wcount += wj;
pj->density.wcount_dh -= (hydro_dimension * wj + uj * wj_dx);
+
+ /* Now we need to compute the div terms */
+ const float r_inv = 1.f / r;
+ const float faci = mj * wi_dx * r_inv;
+ const float facj = mi * wj_dx * r_inv;
+
+ /* Compute dv dot r */
+ dv[0] = pi->v[0] - pj->v[0];
+ dv[1] = pi->v[1] - pj->v[1];
+ dv[2] = pi->v[2] - pj->v[2];
+ const float dvdr = dv[0] * dx[0] + dv[1] * dx[1] + dv[2] * dx[2];
+
+ pi->density.div_v -= faci * dvdr;
+ pj->density.div_v -= facj * dvdr;
+
+ /* Compute dv cross r */
+ curlvr[0] = dv[1] * dx[2] - dv[2] * dx[1];
+ curlvr[1] = dv[2] * dx[0] - dv[0] * dx[2];
+ curlvr[2] = dv[0] * dx[1] - dv[1] * dx[0];
+
+ pi->density.rot_v[0] += faci * curlvr[0];
+ pi->density.rot_v[1] += faci * curlvr[1];
+ pi->density.rot_v[2] += faci * curlvr[2];
+
+ pj->density.rot_v[0] += facj * curlvr[0];
+ pj->density.rot_v[1] += facj * curlvr[1];
+ pj->density.rot_v[2] += facj * curlvr[2];
}
/**
@@ -103,6 +131,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
const struct part* pj, float a, float H) {
float wi, wi_dx;
+ float dv[3], curlvr[3];
/* Get the masses. */
const float mj = pj->mass;
@@ -121,6 +150,26 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
mj * pj->u * (hydro_dimension * wi + ui * wi_dx);
pi->density.wcount += wi;
pi->density.wcount_dh -= (hydro_dimension * wi + ui * wi_dx);
+
+ const float r_inv = 1.f / r;
+ const float faci = mj * wi_dx * r_inv;
+
+ /* Compute dv dot r */
+ dv[0] = pi->v[0] - pj->v[0];
+ dv[1] = pi->v[1] - pj->v[1];
+ dv[2] = pi->v[2] - pj->v[2];
+ const float dvdr = dv[0] * dx[0] + dv[1] * dx[1] + dv[2] * dx[2];
+
+ pi->density.div_v -= faci * dvdr;
+
+ /* Compute dv cross r */
+ curlvr[0] = dv[1] * dx[2] - dv[2] * dx[1];
+ curlvr[1] = dv[2] * dx[0] - dv[0] * dx[2];
+ curlvr[2] = dv[0] * dx[1] - dv[1] * dx[0];
+
+ pi->density.rot_v[0] += faci * curlvr[0];
+ pi->density.rot_v[1] += faci * curlvr[1];
+ pi->density.rot_v[2] += faci * curlvr[2];
}
/**
@@ -189,9 +238,14 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
const float cj = pj->force.soundspeed;
const float v_sig = ci + cj - 3.f * mu_ij;
+ /* 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 = 0.5f * (rhoi + rhoj);
- const float visc = -0.5f * const_viscosity_alpha * v_sig * mu_ij / rho_ij;
+ const float visc = -0.25f * const_viscosity_alpha * v_sig * mu_ij *
+ (balsara_i + balsara_j) / rho_ij;
/* Convolve with the kernel */
const float visc_acc_term = 0.5f * visc * (wi_dr + wj_dr) * r_inv;
@@ -229,7 +283,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
const float du_dt_i = sph_du_term_i + visc_du_term;
const float du_dt_j = sph_du_term_j + visc_du_term;
- /* Internal energy time derivatibe */
+ /* Internal energy time derivative */
pi->u_dt += du_dt_i * mj;
pj->u_dt += du_dt_j * mi;
@@ -309,9 +363,14 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
const float cj = pj->force.soundspeed;
const float v_sig = ci + cj - 3.f * mu_ij;
+ /* 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 = 0.5f * (rhoi + rhoj);
- const float visc = -0.5f * const_viscosity_alpha * v_sig * mu_ij / rho_ij;
+ const float visc = -0.25f * const_viscosity_alpha * v_sig * mu_ij *
+ (balsara_i + balsara_j) / rho_ij;
/* Convolve with the kernel */
const float visc_acc_term = 0.5f * visc * (wi_dr + wj_dr) * r_inv;
diff --git a/src/hydro/PressureEnergy/hydro_part.h b/src/hydro/PressureEnergy/hydro_part.h
index 30f66abe264960f39eca7bc00e7b1097db230f40..378ac1ffbaba516ef79a4560444229d1af7b143c 100644
--- a/src/hydro/PressureEnergy/hydro_part.h
+++ b/src/hydro/PressureEnergy/hydro_part.h
@@ -160,6 +160,8 @@ struct part {
/*! Time derivative of smoothing length */
float h_dt;
+ /*! Balsara switch */
+ float balsara;
} force;
};