From aefe3c7154feb2f4a832ea01c951a444982b92c6 Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <matthieu.schaller@durham.ac.uk>
Date: Tue, 13 Sep 2016 12:23:50 +0100
Subject: [PATCH] Added the Balsara viscosity switch to the P-E implementation
---
examples/SedovBlast_1D/sedov.yml | 2 +-
src/const.h | 4 +-
src/engine.c | 1 +
src/hydro/PressureEntropy/hydro.h | 86 ++++++++++++-----
src/hydro/PressureEntropy/hydro_debug.h | 16 ++--
src/hydro/PressureEntropy/hydro_iact.h | 120 ++++++++++++------------
src/hydro/PressureEntropy/hydro_io.h | 3 +-
src/hydro/PressureEntropy/hydro_part.h | 26 +++--
8 files changed, 150 insertions(+), 108 deletions(-)
diff --git a/examples/SedovBlast_1D/sedov.yml b/examples/SedovBlast_1D/sedov.yml
index 6f519835d2..1ecfeb3245 100644
--- a/examples/SedovBlast_1D/sedov.yml
+++ b/examples/SedovBlast_1D/sedov.yml
@@ -21,7 +21,7 @@ Snapshots:
# Parameters governing the conserved quantities statistics
Statistics:
- delta_time: 1e-3 # Time between statistics output
+ delta_time: 1e-5 # Time between statistics output
# Parameters for the hydrodynamics scheme
SPH:
diff --git a/src/const.h b/src/const.h
index 2fc1068d41..4bdcc3aff1 100644
--- a/src/const.h
+++ b/src/const.h
@@ -41,8 +41,8 @@
/* Dimensionality of the problem */
//#define HYDRO_DIMENSION_3D
-#define HYDRO_DIMENSION_2D
-//#define HYDRO_DIMENSION_1D
+//#define HYDRO_DIMENSION_2D
+#define HYDRO_DIMENSION_1D
/* Hydrodynamical adiabatic index. */
#define HYDRO_GAMMA_5_3
diff --git a/src/engine.c b/src/engine.c
index 7adc94b427..4a6460d8cd 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -2659,6 +2659,7 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs) {
/* Apply some conversions (e.g. internal energy -> entropy) */
if (!flag_entropy_ICs) space_map_cells_pre(s, 0, cell_convert_hydro, NULL);
+ if (1) engine_launch(e, e->nr_threads, mask, submask);
clocks_gettime(&time2);
diff --git a/src/hydro/PressureEntropy/hydro.h b/src/hydro/PressureEntropy/hydro.h
index b2a56cf260..33cf5267eb 100644
--- a/src/hydro/PressureEntropy/hydro.h
+++ b/src/hydro/PressureEntropy/hydro.h
@@ -27,7 +27,8 @@
* The thermal variable is the entropy (S) and the entropy is smoothed over
* contact discontinuities to prevent spurious surface tension.
*
- * Follows Hopkins, P., MNRAS, 2013, Volume 428, Issue 4, pp. 2840-2856
+ * Follows eqautions (19), (21) and (22) of Hopkins, P., MNRAS, 2013,
+ * Volume 428, Issue 4, pp. 2840-2856 with a simple Balsara viscosity term.
*/
#include "adiabatic_index.h"
@@ -90,7 +91,7 @@ __attribute__((always_inline)) INLINE static float hydro_get_soundspeed(
}
/**
- * @brief Returns the density of a particle
+ * @brief Returns the physical density of a particle
*
* @param p The particle of interest
*/
@@ -125,6 +126,17 @@ __attribute__((always_inline)) INLINE static void hydro_set_internal_energy(
struct part *restrict p, float u) {
p->entropy = gas_entropy_from_internal_energy(p->rho_bar, u);
+ p->entropy_one_over_gamma = pow_one_over_gamma(p->entropy);
+
+ /* Compute the pressure */
+ const float entropy = p->entropy;
+ const float pressure = gas_pressure_from_entropy(p->rho_bar, entropy);
+
+ /* Compute the sound speed from the pressure*/
+ const float rho_bar_inv = 1.f / p->rho_bar;
+ const float soundspeed = sqrtf(hydro_gamma * pressure * rho_bar_inv);
+ p->force.soundspeed = soundspeed;
+ p->force.P_over_rho2 = pressure * rho_bar_inv * rho_bar_inv;
}
/**
@@ -140,6 +152,17 @@ __attribute__((always_inline)) INLINE static void hydro_set_entropy(
struct part *restrict p, float S) {
p->entropy = S;
+ p->entropy_one_over_gamma = pow_one_over_gamma(p->entropy);
+
+ /* Compute the pressure */
+ const float entropy = p->entropy;
+ const float pressure = gas_pressure_from_entropy(p->rho_bar, entropy);
+
+ /* Compute the sound speed from the pressure*/
+ const float rho_bar_inv = 1.f / p->rho_bar;
+ const float soundspeed = sqrtf(hydro_gamma * pressure * rho_bar_inv);
+ p->force.soundspeed = soundspeed;
+ p->force.P_over_rho2 = pressure * rho_bar_inv * rho_bar_inv;
}
/**
@@ -176,6 +199,8 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
p->ti_begin = 0;
p->ti_end = 0;
+ p->rho_bar = 0.f;
+ p->entropy_one_over_gamma = pow_one_over_gamma(p->entropy);
xp->v_full[0] = p->v[0];
xp->v_full[1] = p->v[1];
xp->v_full[2] = p->v[2];
@@ -194,15 +219,14 @@ __attribute__((always_inline)) INLINE static void hydro_init_part(
p->density.wcount = 0.f;
p->density.wcount_dh = 0.f;
p->rho = 0.f;
- p->rho_dh = 0.f;
p->rho_bar = 0.f;
- p->pressure_dh = 0.f;
+ p->density.rho_dh = 0.f;
+ p->density.pressure_dh = 0.f;
- // p->density.weightedPressure_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; */
+ 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;
}
/**
@@ -225,17 +249,17 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
/* Final operation on the density (add self-contribution). */
p->rho += p->mass * kernel_root;
- p->rho_dh -= hydro_dimension * p->mass * kernel_root;
p->rho_bar += p->mass * p->entropy_one_over_gamma * kernel_root;
- p->pressure_dh -=
+ p->density.rho_dh -= hydro_dimension * p->mass * kernel_root;
+ p->density.pressure_dh -=
hydro_dimension * p->mass * p->entropy_one_over_gamma * kernel_root;
p->density.wcount += kernel_root;
/* Finish the calculation by inserting the missing h-factors */
p->rho *= h_inv_dim;
- p->rho_dh *= h_inv_dim;
p->rho_bar *= h_inv_dim;
- p->pressure_dh *= h_inv_dim_plus_one;
+ p->density.rho_dh *= h_inv_dim_plus_one;
+ p->density.pressure_dh *= h_inv_dim_plus_one;
p->density.wcount *= kernel_norm;
p->density.wcount_dh *= h_inv * kernel_gamma * kernel_norm;
@@ -246,17 +270,18 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
p->rho_bar *= entropy_minus_one_over_gamma;
/* Compute the derivative term */
- p->rho_dh = 1.f / (1.f + hydro_dimension_inv * p->h * p->rho_dh * rho_inv);
- p->pressure_dh *=
- p->h * rho_inv * hydro_dimension_inv * entropy_minus_one_over_gamma;
+ p->density.rho_dh =
+ 1.f / (1.f + hydro_dimension_inv * h * p->density.rho_dh * rho_inv);
+ p->density.pressure_dh *=
+ h * rho_inv * hydro_dimension_inv * entropy_minus_one_over_gamma;
/* Finish calculation of the velocity curl components */
- // p->density.rot_v[0] *= h_inv_dim_plus_one * irho;
- // p->density.rot_v[1] *= h_inv_dim_plus_one * irho;
- // p->density.rot_v[2] *= h_inv_dim_plus_one * irho;
+ p->density.rot_v[0] *= h_inv_dim_plus_one * rho_inv;
+ p->density.rot_v[1] *= h_inv_dim_plus_one * rho_inv;
+ p->density.rot_v[2] *= h_inv_dim_plus_one * rho_inv;
/* Finish calculation of the velocity divergence */
- // p->density.div_v *= h_inv_dim_plus_one * irho;
+ p->density.div_v *= h_inv_dim_plus_one * rho_inv;
}
/**
@@ -273,6 +298,16 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
struct part *restrict p, struct xpart *restrict xp, int ti_current,
double timeBase) {
+ const float fac_mu = 1.f; /* Will change with cosmological integration */
+
+ /* 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 half_dt = (ti_current - (p->ti_begin + p->ti_end) / 2) * timeBase;
const float pressure = hydro_get_pressure(p, half_dt);
@@ -281,9 +316,18 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
const float rho_bar_inv = 1.f / p->rho_bar;
const float soundspeed = sqrtf(hydro_gamma * pressure * rho_bar_inv);
+ /* Compute the Balsara switch */
+ const float balsara =
+ abs_div_v / (abs_div_v + curl_v + 0.0001f * soundspeed / fac_mu / p->h);
+
+ /* Compute "grad h" term */
+ const float grad_h_term = p->density.rho_dh * p->density.pressure_dh;
+
/* Update variables. */
p->force.soundspeed = soundspeed;
p->force.P_over_rho2 = pressure * rho_bar_inv * rho_bar_inv;
+ p->force.balsara = balsara;
+ p->force.f = grad_h_term;
}
/**
@@ -422,7 +466,7 @@ __attribute__((always_inline)) INLINE static void hydro_convert_quantities(
struct part *restrict p) {
/* We read u in the entropy field. We now get S from u */
- p->entropy = gas_entropy_from_internal_energy(p->rho, p->entropy);
+ p->entropy = gas_entropy_from_internal_energy(p->rho_bar, p->entropy);
p->entropy_one_over_gamma = pow_one_over_gamma(p->entropy);
/* Compute the pressure */
diff --git a/src/hydro/PressureEntropy/hydro_debug.h b/src/hydro/PressureEntropy/hydro_debug.h
index daa0d2a7fd..4865437935 100644
--- a/src/hydro/PressureEntropy/hydro_debug.h
+++ b/src/hydro/PressureEntropy/hydro_debug.h
@@ -26,7 +26,8 @@
* The thermal variable is the entropy (S) and the entropy is smoothed over
* contact discontinuities to prevent spurious surface tension.
*
- * Follows Hopkins, P., MNRAS, 2013, Volume 428, Issue 4, pp. 2840-2856
+ * Follows eqautions (19), (21) and (22) of Hopkins, P., MNRAS, 2013,
+ * Volume 428, Issue 4, pp. 2840-2856 with a simple Balsara viscosity term.
*/
__attribute__((always_inline)) INLINE static void hydro_debug_particle(
@@ -35,14 +36,15 @@ __attribute__((always_inline)) INLINE static void hydro_debug_particle(
"x=[%.3e,%.3e,%.3e], "
"v=[%.3e,%.3e,%.3e],v_full=[%.3e,%.3e,%.3e] \n a=[%.3e,%.3e,%.3e],\n "
"h=%.3e, wcount=%.3f, wcount_dh=%.3e, m=%.3e, dh_drho=%.3e, rho=%.3e, "
- "P=%.3e, P_over_rho2=%.3e, S=%.3e, dS/dt=%.3e, c=%.3e\n"
- "v_sig=%e dh/dt=%.3e t_begin=%d, t_end=%d\n",
+ "rho_bar=%.3e, P=%.3e, dP_dh=%.3e, P_over_rho2=%.3e, S=%.3e, S^1/g=%.3e, "
+ "dS/dt=%.3e,\nc=%.3e v_sig=%e dh/dt=%.3e t_begin=%d, t_end=%d\n",
p->x[0], p->x[1], p->x[2], p->v[0], p->v[1], p->v[2], xp->v_full[0],
xp->v_full[1], xp->v_full[2], p->a_hydro[0], p->a_hydro[1], p->a_hydro[2],
- p->h, p->density.wcount, p->density.wcount_dh, p->mass, 0.f /*p->rho_dh*/,
- p->rho, hydro_get_pressure(p, 0.), 0.f,
- /*p->force.P_over_rho2,*/ p->entropy, p->entropy_dt, p->force.soundspeed,
- p->force.v_sig, p->force.h_dt, p->ti_begin, p->ti_end);
+ p->h, p->density.wcount, p->density.wcount_dh, p->mass, p->density.rho_dh,
+ p->rho, p->rho_bar, hydro_get_pressure(p, 0.), p->density.pressure_dh,
+ p->force.P_over_rho2, p->entropy, p->entropy_one_over_gamma,
+ p->entropy_dt, p->force.soundspeed, p->force.v_sig, p->force.h_dt,
+ p->ti_begin, p->ti_end);
}
#endif /* SWIFT_PRESSURE_ENTROPY_HYDRO_DEBUG_H */
diff --git a/src/hydro/PressureEntropy/hydro_iact.h b/src/hydro/PressureEntropy/hydro_iact.h
index 6068652093..18e22233f6 100644
--- a/src/hydro/PressureEntropy/hydro_iact.h
+++ b/src/hydro/PressureEntropy/hydro_iact.h
@@ -26,7 +26,8 @@
* The thermal variable is the entropy (S) and the entropy is smoothed over
* contact discontinuities to prevent spurious surface tension.
*
- * Follows Hopkins, P., MNRAS, 2013, Volume 428, Issue 4, pp. 2840-2856
+ * Follows eqautions (19), (21) and (22) of Hopkins, P., MNRAS, 2013,
+ * Volume 428, Issue 4, pp. 2840-2856 with a simple Balsara viscosity term.
*/
/**
@@ -37,19 +38,15 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
float wi, wi_dx;
float wj, wj_dx;
- /* float dv[3], curlvr[3]; */
+ float dv[3], curlvr[3];
/* Get the masses. */
const float mi = pi->mass;
const float mj = pj->mass;
- /* /\* Get the entropies. *\/ */
- /* const float entropy_i = pi->entropy; */
- /* const float entropy_j = pj->entropy; */
-
/* Get r and r inverse. */
const float r = sqrtf(r2);
- /* const float r_inv = 1.0f / r; */
+ const float r_inv = 1.0f / r;
/* Compute the kernel function for pi */
const float hi_inv = 1.f / hi;
@@ -58,7 +55,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
/* Compute contribution to the density */
pi->rho += mj * wi;
- pi->rho_dh -= mj * (hydro_dimension * wi + ui * wi_dx);
+ pi->density.rho_dh -= mj * (hydro_dimension * wi + ui * wi_dx);
/* Compute contribution to the number of neighbours */
pi->density.wcount += wi;
@@ -66,7 +63,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
/* Compute contribution to the weighted density */
pi->rho_bar += mj * pj->entropy_one_over_gamma * wi;
- pi->pressure_dh -=
+ pi->density.pressure_dh -=
mj * pj->entropy_one_over_gamma * (hydro_dimension * wi + ui * wi_dx);
/* Compute the kernel function for pj */
@@ -76,7 +73,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
/* Compute contribution to the density */
pj->rho += mi * wj;
- pj->rho_dh -= mi * (hydro_dimension * wj + uj * wj_dx);
+ pj->density.rho_dh -= mi * (hydro_dimension * wj + uj * wj_dx);
/* Compute contribution to the number of neighbours */
pj->density.wcount += wj;
@@ -84,33 +81,33 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
/* Compute contribution to the weighted density */
pj->rho_bar += mi * pi->entropy_one_over_gamma * wj;
- pj->pressure_dh -=
+ pj->density.pressure_dh -=
mi * pi->entropy_one_over_gamma * (hydro_dimension * wj + uj * wj_dx);
- /* const float faci = mj * wi_dx * r_inv; */
- /* const float facj = mi * wj_dx * r_inv; */
+ 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]; */
+ /* 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; */
+ 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]; */
+ /* 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]; */
+ 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]; */
+ pj->density.rot_v[0] += facj * curlvr[0];
+ pj->density.rot_v[1] += facj * curlvr[1];
+ pj->density.rot_v[2] += facj * curlvr[2];
}
/**
@@ -130,17 +127,14 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {
float wi, wi_dx;
- /* float dv[3], curlvr[3]; */
+ float dv[3], curlvr[3];
/* Get the masses. */
const float mj = pj->mass;
- /* Get the entropies. */
- // const float entropy_j = pj->entropy;
-
/* Get r and r inverse. */
const float r = sqrtf(r2);
- /* const float ri = 1.0f / r; */
+ const float ri = 1.0f / r;
/* Compute the kernel function */
const float h_inv = 1.0f / hi;
@@ -149,34 +143,34 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
/* Compute contribution to the density */
pi->rho += mj * wi;
- pi->rho_dh -= mj * (hydro_dimension * wi + ui * wi_dx);
+ pi->density.rho_dh -= mj * (hydro_dimension * wi + ui * wi_dx);
/* Compute contribution to the number of neighbours */
pi->density.wcount += wi;
- pi->density.wcount_dh -= ui * wi_dx; //(hydro_dimension * wi + u * wi_dx);
+ pi->density.wcount_dh -= ui * wi_dx;
/* Compute contribution to the weighted density */
pi->rho_bar += mj * pj->entropy_one_over_gamma * wi;
- pi->pressure_dh -=
+ pi->density.pressure_dh -=
mj * pj->entropy_one_over_gamma * (hydro_dimension * wi + ui * wi_dx);
- /* const float fac = mj * wi_dx * ri; */
+ const float fac = mj * wi_dx * ri;
- /* /\* 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 -= fac * dvdr; */
+ /* 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 -= fac * 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]; */
+ /* 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] += fac * curlvr[0]; */
- /* pi->density.rot_v[1] += fac * curlvr[1]; */
- /* pi->density.rot_v[2] += fac * curlvr[2]; */
+ pi->density.rot_v[0] += fac * curlvr[0];
+ pi->density.rot_v[1] += fac * curlvr[1];
+ pi->density.rot_v[2] += fac * curlvr[2];
}
/**
@@ -223,8 +217,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
const float wj_dr = hjd_inv * wj_dx;
/* Compute gradient terms */
- const float f_i = pi->rho_dh * pi->pressure_dh;
- const float f_j = pj->rho_dh * pj->pressure_dh;
+ const float f_i = pi->force.f;
+ const float f_j = pj->force.f;
/* Compute Pressure terms */
const float P_over_rho2_i = pi->force.P_over_rho2;
@@ -244,8 +238,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
(pi->v[2] - pj->v[2]) * dx[2];
/* Balsara term */
- /* const float balsara_i = pi->force.balsara; */
- /* const float balsara_j = pj->force.balsara; */
+ const float balsara_i = pi->force.balsara;
+ const float balsara_j = pj->force.balsara;
/* Are the particles moving towards each others ? */
const float omega_ij = fminf(dvdr, 0.f);
@@ -256,7 +250,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
/* Now 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;
/* Now, convolve with the kernel */
const float visc_term = 0.5f * visc * (wi_dr + wj_dr);
@@ -332,8 +327,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
const float wj_dr = hjd_inv * wj_dx;
/* Compute gradient terms */
- const float f_i = pi->rho_dh * pi->pressure_dh;
- const float f_j = pj->rho_dh * pj->pressure_dh;
+ const float f_i = pi->force.f;
+ const float f_j = pj->force.f;
/* Compute Pressure terms */
const float P_over_rho2_i = pi->force.P_over_rho2;
@@ -353,8 +348,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
(pi->v[2] - pj->v[2]) * dx[2];
/* Balsara term */
- /* const float balsara_i = pi->force.balsara; */
- /* const float balsara_j = pj->force.balsara; */
+ const float balsara_i = pi->force.balsara;
+ const float balsara_j = pj->force.balsara;
/* Are the particles moving towards each others ? */
const float omega_ij = fminf(dvdr, 0.f);
@@ -365,7 +360,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
/* Now 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;
/* Now, convolve with the kernel */
const float visc_term = 0.5f * visc * (wi_dr + wj_dr);
diff --git a/src/hydro/PressureEntropy/hydro_io.h b/src/hydro/PressureEntropy/hydro_io.h
index c57e8eb684..3263b6b7a9 100644
--- a/src/hydro/PressureEntropy/hydro_io.h
+++ b/src/hydro/PressureEntropy/hydro_io.h
@@ -26,7 +26,8 @@
* The thermal variable is the entropy (S) and the entropy is smoothed over
* contact discontinuities to prevent spurious surface tension.
*
- * Follows Hopkins, P., MNRAS, 2013, Volume 428, Issue 4, pp. 2840-2856
+ * Follows eqautions (19), (21) and (22) of Hopkins, P., MNRAS, 2013,
+ * Volume 428, Issue 4, pp. 2840-2856 with a simple Balsara viscosity term.
*/
#include "adiabatic_index.h"
diff --git a/src/hydro/PressureEntropy/hydro_part.h b/src/hydro/PressureEntropy/hydro_part.h
index 1ef61c4264..cac585ff79 100644
--- a/src/hydro/PressureEntropy/hydro_part.h
+++ b/src/hydro/PressureEntropy/hydro_part.h
@@ -26,7 +26,8 @@
* The thermal variable is the entropy (S) and the entropy is smoothed over
* contact discontinuities to prevent spurious surface tension.
*
- * Follows Hopkins, P., MNRAS, 2013, Volume 428, Issue 4, pp. 2840-2856
+ * Follows eqautions (19), (21) and (22) of Hopkins, P., MNRAS, 2013,
+ * Volume 428, Issue 4, pp. 2840-2856 with a simple Balsara viscosity term.
*/
#include "cooling_struct.h"
@@ -72,15 +73,6 @@ struct part {
/*! Particle density. */
float rho;
- /*! Derivative of density with respect to h */
- float rho_dh;
-
- /*! Particle pressure. */
- float pressure;
-
- /*! Derivative of pressure with respect to h */
- float pressure_dh;
-
/*! Particle weighted density */
float rho_bar;
@@ -103,11 +95,14 @@ struct part {
/*! Number of neighbours spatial derivative. */
float wcount_dh;
- /*! Derivative of particle weighted pressure with h. */
- // float weightedPressure_dh;
+ /*! Derivative of density with respect to h */
+ float rho_dh;
+
+ /*! Derivative of pressure with respect to h */
+ float pressure_dh;
/*! Particle velocity curl. */
- // float rot_v[3];
+ float rot_v[3];
/*! Particle velocity divergence. */
float div_v;
@@ -116,8 +111,11 @@ struct part {
struct {
+ /*! Balsara switch */
+ float balsara;
+
/*! "Grad h" term */
- // float f_ij;
+ float f;
/*! Pressure term */
float P_over_rho2;
--
GitLab