Skip to content
Snippets Groups Projects
Commit be2a599c authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Merge branch 'fix_phantom' into 'master' 

Three fixes to the PHANTOM SPH scheme

Closes #803

See merge request !1607
parents ef0d40bf 20cbe767
No related branches found
No related tags found
12 merge requests!1715Update planetary strength after planetary plus's master rebase,!1693More threapool plotting tweaks,!1668before Mag.Egy in all the flavors,!1663Initial sync to work again,!1662Initial sync from previous months,!1642When searching for more particles in a ghost task we walk up the cell tree,...,!1633When searching for more particles in a ghost task we walk up the cell tree,...,!1631Solving issues with different Hydro Schemes,!1620Mhd canvas,!1619Mhd canvas,!1617Update to from main branch to Canvas to MHD_FS,!1607Three fixes to the PHANTOM SPH scheme
Hide whitespace changes
Inline Side-by-side
src/hydro/Phantom/hydro_iact.h
+ 9
− 11
View file @ be2a599c
...@@ -333,9 +333,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_force( ...@@ -333,9 +333,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
const float omega_ij = min(dvdr_Hubble, 0.f); const float omega_ij = min(dvdr_Hubble, 0.f);
const float mu_ij = fac_mu * r_inv * omega_ij; /* This is 0 or negative */ const float mu_ij = fac_mu * r_inv * omega_ij; /* This is 0 or negative */
/* Compute sound speeds and signal velocity */ /* Compute the signal velocity */
const float v_sig = pi->force.soundspeed + pj->force.soundspeed - const float v_sig = signal_velocity(dx, pi, pj, mu_ij, const_viscosity_beta);
const_viscosity_beta * mu_ij;
/* Variable smoothing length term */ /* Variable smoothing length term */
const float f_ij = 1.f - pi->force.f / mj; const float f_ij = 1.f - pi->force.f / mj;
...@@ -389,7 +388,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_force( ...@@ -389,7 +388,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
const float alpha_diff = 0.5f * (pi->diffusion.alpha + pj->diffusion.alpha); const float alpha_diff = 0.5f * (pi->diffusion.alpha + pj->diffusion.alpha);
/* wi_dx + wj_dx / 2 is F_ij */ /* wi_dx + wj_dx / 2 is F_ij */
const float diff_du_term = alpha_diff * v_diff * (pi->u - pj->u) * 0.5f * const float diff_du_term = alpha_diff * v_diff * (pi->u - pj->u) * 0.5f *
(wi_dr * f_ij / pi->rho + wj_dr * f_ji / pi->rho); (wi_dr * f_ij / rhoi + wj_dr * f_ji / rhoj);
/* Assemble the energy equation term */ /* Assemble the energy equation term */
const float du_dt_i = sph_du_term_i + visc_du_term + diff_du_term; const float du_dt_i = sph_du_term_i + visc_du_term + diff_du_term;
...@@ -400,8 +399,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_force( ...@@ -400,8 +399,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
pj->u_dt += du_dt_j * mi; pj->u_dt += du_dt_j * mi;
/* Get the time derivative for h. */ /* Get the time derivative for h. */
pi->force.h_dt -= mj * dvdr * f_ij * r_inv / rhoj * wi_dr; pi->force.h_dt -= mj * dvdr * r_inv / rhoj * wi_dr;
pj->force.h_dt -= mi * dvdr * f_ji * r_inv / rhoi * wj_dr; pj->force.h_dt -= mi * dvdr * r_inv / rhoi * wj_dr;
} }
/** /**
...@@ -466,9 +465,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force( ...@@ -466,9 +465,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
const float omega_ij = min(dvdr_Hubble, 0.f); const float omega_ij = min(dvdr_Hubble, 0.f);
const float mu_ij = fac_mu * r_inv * omega_ij; /* This is 0 or negative */ const float mu_ij = fac_mu * r_inv * omega_ij; /* This is 0 or negative */
/* Compute sound speeds and signal velocity */ /* Compute the signal velocity */
const float v_sig = pi->force.soundspeed + pj->force.soundspeed - const float v_sig = signal_velocity(dx, pi, pj, mu_ij, const_viscosity_beta);
const_viscosity_beta * mu_ij;
/* Variable smoothing length term */ /* Variable smoothing length term */
const float f_ij = 1.f - pi->force.f / mj; const float f_ij = 1.f - pi->force.f / mj;
...@@ -516,7 +514,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force( ...@@ -516,7 +514,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
const float alpha_diff = 0.5f * (pi->diffusion.alpha + pj->diffusion.alpha); const float alpha_diff = 0.5f * (pi->diffusion.alpha + pj->diffusion.alpha);
/* wi_dx + wj_dx / 2 is F_ij */ /* wi_dx + wj_dx / 2 is F_ij */
const float diff_du_term = alpha_diff * v_diff * (pi->u - pj->u) * 0.5f * const float diff_du_term = alpha_diff * v_diff * (pi->u - pj->u) * 0.5f *
(wi_dr * f_ij / pi->rho + wj_dr * f_ji / pi->rho); (wi_dr * f_ij / rhoi + wj_dr * f_ji / rhoj);
/* Assemble the energy equation term */ /* Assemble the energy equation term */
const float du_dt_i = sph_du_term_i + visc_du_term + diff_du_term; const float du_dt_i = sph_du_term_i + visc_du_term + diff_du_term;
...@@ -525,7 +523,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force( ...@@ -525,7 +523,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
pi->u_dt += du_dt_i * mj; pi->u_dt += du_dt_i * mj;
/* Get the time derivative for h. */ /* Get the time derivative for h. */
pi->force.h_dt -= mj * dvdr * f_ij * r_inv / rhoj * wi_dr; pi->force.h_dt -= mj * dvdr * r_inv / rhoj * wi_dr;
} }
#endif /* SWIFT_PHANTOM_HYDRO_IACT_H */ #endif /* SWIFT_PHANTOM_HYDRO_IACT_H */
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment