gitlab is upgraded to version 13, please report any issues and enjoy

Commit a6760b65 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Merge branch 'change_f_ij_sphenix' into 'master'

Changed f_ij factors in SPHENIX to respect unequal m

See merge request !1105
parents f96a961b 5aa1e810
......@@ -603,17 +603,21 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_gradient(
const float balsara =
abs_div_v / (abs_div_v + curl_v + 0.0001f * soundspeed * fac_B / p->h);
/* Compute the "grad h" term */
const float rho_inv = 1.f / p->rho;
float rho_dh = p->density.rho_dh;
/* Compute the "grad h" term - Note here that we have \tilde{x}
* as 1 as we use the local number density to find neighbours. This
* introduces a j-component that is considered in the force loop,
* meaning that this cached grad_h_term gives:
*
* f_ij = 1.f - grad_h_term_i / m_j */
const float common_factor = p->h * hydro_dimension_inv / p->density.wcount;
float grad_h_term = common_factor * p->density.rho_dh /
(1.f + common_factor * p->density.wcount_dh);
/* Ignore changing-kernel effects when h ~= h_max */
if (p->h > 0.9999f * hydro_props->h_max) {
rho_dh = 0.f;
grad_h_term = 0.f;
}
const float grad_h_term =
1.f / (1.f + hydro_dimension_inv * p->h * rho_dh * rho_inv);
/* Update variables. */
p->force.f = grad_h_term;
p->force.pressure = pressure_including_floor;
......
......@@ -366,6 +366,10 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
const float v_sig = pi->force.soundspeed + pj->force.soundspeed -
const_viscosity_beta * mu_ij;
/* Variable smoothing length term */
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;
......@@ -377,11 +381,12 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
-0.25f * 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;
const float visc_acc_term =
0.5f * visc * (wi_dr * f_ij + wj_dr * f_ji) * r_inv;
/* Compute gradient terms */
const float P_over_rho2_i = pressurei / (rhoi * rhoi) * pi->force.f;
const float P_over_rho2_j = pressurej / (rhoj * rhoj) * pj->force.f;
const float P_over_rho2_i = pressurei / (rhoi * rhoi) * f_ij;
const float P_over_rho2_j = pressurej / (rhoj * rhoj) * f_ji;
/* SPH acceleration term */
const float sph_acc_term =
......@@ -419,7 +424,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
fabsf(fac_mu * r_inv * dvdr_Hubble));
/* wi_dx + wj_dx / 2 is F_ij */
const float diff_du_term =
v_diff * (pi->u - pj->u) * (wi_dr / rhoi + wj_dr / rhoj);
v_diff * (pi->u - pj->u) * (f_ij * wi_dr / rhoi + f_ji * wj_dr / rhoj);
/* Assemble the energy equation term */
const float du_dt_i = sph_du_term_i + visc_du_term + diff_du_term;
......@@ -464,7 +469,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
const float r_inv = 1.0f / r;
/* Recover some data */
// const float mi = pi->mass;
const float mi = pi->mass;
const float mj = pj->mass;
const float rhoi = pi->rho;
......@@ -505,6 +510,10 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
const float v_sig = pi->force.soundspeed + pj->force.soundspeed -
const_viscosity_beta * mu_ij;
/* Variable smoothing length term */
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;
......@@ -516,11 +525,12 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
-0.25f * 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;
const float visc_acc_term =
0.5f * visc * (wi_dr * f_ij + wj_dr * f_ji) * r_inv;
/* Compute gradient terms */
const float P_over_rho2_i = pressurei / (rhoi * rhoi) * pi->force.f;
const float P_over_rho2_j = pressurej / (rhoj * rhoj) * pj->force.f;
const float P_over_rho2_i = pressurei / (rhoi * rhoi) * f_ij;
const float P_over_rho2_j = pressurej / (rhoj * rhoj) * f_ji;
/* SPH acceleration term */
const float sph_acc_term =
......@@ -553,7 +563,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
fabsf(fac_mu * r_inv * dvdr_Hubble));
/* wi_dx + wj_dx / 2 is F_ij */
const float diff_du_term =
v_diff * (pi->u - pj->u) * (wi_dr / rhoi + wj_dr / rhoj);
v_diff * (pi->u - pj->u) * (f_ij * wi_dr / rhoi + f_ji * wj_dr / rhoj);
/* Assemble the energy equation term */
const float du_dt_i = sph_du_term_i + visc_du_term + diff_du_term;
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment