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

Added the Balsara viscosity switch to the P-E implementation

parent 3674b337
No related branches found
No related tags found
1 merge request!252Implementation Pressure entropy SPH
Hide whitespace changes
Inline Side-by-side
Showing
with 150 additions and 108 deletions
examples/SedovBlast_1D/sedov.yml
+ 1
− 1
View file @ aefe3c71
...@@ -21,7 +21,7 @@ Snapshots: ...@@ -21,7 +21,7 @@ Snapshots:
# Parameters governing the conserved quantities statistics # Parameters governing the conserved quantities statistics
Statistics: Statistics:
delta_time: 1e-3 # Time between statistics output delta_time: 1e-5 # Time between statistics output
# Parameters for the hydrodynamics scheme # Parameters for the hydrodynamics scheme
SPH: SPH:
......
src/const.h
+ 2
− 2
View file @ aefe3c71
...@@ -41,8 +41,8 @@ ...@@ -41,8 +41,8 @@
/* Dimensionality of the problem */ /* Dimensionality of the problem */
//#define HYDRO_DIMENSION_3D //#define HYDRO_DIMENSION_3D
#define HYDRO_DIMENSION_2D //#define HYDRO_DIMENSION_2D
//#define HYDRO_DIMENSION_1D #define HYDRO_DIMENSION_1D
/* Hydrodynamical adiabatic index. */ /* Hydrodynamical adiabatic index. */
#define HYDRO_GAMMA_5_3 #define HYDRO_GAMMA_5_3
......
src/engine.c
+ 1
− 0
View file @ aefe3c71
...@@ -2659,6 +2659,7 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs) { ...@@ -2659,6 +2659,7 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs) {
/* Apply some conversions (e.g. internal energy -> entropy) */ /* Apply some conversions (e.g. internal energy -> entropy) */
if (!flag_entropy_ICs) space_map_cells_pre(s, 0, cell_convert_hydro, NULL); 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); clocks_gettime(&time2);
......
src/hydro/PressureEntropy/hydro.h
+ 65
− 21
View file @ aefe3c71
...@@ -27,7 +27,8 @@ ...@@ -27,7 +27,8 @@
* The thermal variable is the entropy (S) and the entropy is smoothed over * The thermal variable is the entropy (S) and the entropy is smoothed over
* contact discontinuities to prevent spurious surface tension. * 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" #include "adiabatic_index.h"
...@@ -90,7 +91,7 @@ __attribute__((always_inline)) INLINE static float hydro_get_soundspeed( ...@@ -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 * @param p The particle of interest
*/ */
...@@ -125,6 +126,17 @@ __attribute__((always_inline)) INLINE static void hydro_set_internal_energy( ...@@ -125,6 +126,17 @@ __attribute__((always_inline)) INLINE static void hydro_set_internal_energy(
struct part *restrict p, float u) { struct part *restrict p, float u) {
p->entropy = gas_entropy_from_internal_energy(p->rho_bar, 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( ...@@ -140,6 +152,17 @@ __attribute__((always_inline)) INLINE static void hydro_set_entropy(
struct part *restrict p, float S) { struct part *restrict p, float S) {
p->entropy = 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( ...@@ -176,6 +199,8 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
p->ti_begin = 0; p->ti_begin = 0;
p->ti_end = 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[0] = p->v[0];
xp->v_full[1] = p->v[1]; xp->v_full[1] = p->v[1];
xp->v_full[2] = p->v[2]; xp->v_full[2] = p->v[2];
...@@ -194,15 +219,14 @@ __attribute__((always_inline)) INLINE static void hydro_init_part( ...@@ -194,15 +219,14 @@ __attribute__((always_inline)) INLINE static void hydro_init_part(
p->density.wcount = 0.f; p->density.wcount = 0.f;
p->density.wcount_dh = 0.f; p->density.wcount_dh = 0.f;
p->rho = 0.f; p->rho = 0.f;
p->rho_dh = 0.f;
p->rho_bar = 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.div_v = 0.f; */ p->density.rot_v[0] = 0.f;
/* p->density.rot_v[0] = 0.f; */ p->density.rot_v[1] = 0.f;
/* p->density.rot_v[1] = 0.f; */ p->density.rot_v[2] = 0.f;
/* p->density.rot_v[2] = 0.f; */
} }
/** /**
...@@ -225,17 +249,17 @@ __attribute__((always_inline)) INLINE static void hydro_end_density( ...@@ -225,17 +249,17 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
/* Final operation on the density (add self-contribution). */ /* Final operation on the density (add self-contribution). */
p->rho += p->mass * kernel_root; 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->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; hydro_dimension * p->mass * p->entropy_one_over_gamma * kernel_root;
p->density.wcount += kernel_root; p->density.wcount += kernel_root;
/* Finish the calculation by inserting the missing h-factors */ /* Finish the calculation by inserting the missing h-factors */
p->rho *= h_inv_dim; p->rho *= h_inv_dim;
p->rho_dh *= h_inv_dim;
p->rho_bar *= 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 *= kernel_norm;
p->density.wcount_dh *= h_inv * kernel_gamma * 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( ...@@ -246,17 +270,18 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
p->rho_bar *= entropy_minus_one_over_gamma; p->rho_bar *= entropy_minus_one_over_gamma;
/* Compute the derivative term */ /* Compute the derivative term */
p->rho_dh = 1.f / (1.f + hydro_dimension_inv * p->h * p->rho_dh * rho_inv); p->density.rho_dh =
p->pressure_dh *= 1.f / (1.f + hydro_dimension_inv * h * p->density.rho_dh * rho_inv);
p->h * rho_inv * hydro_dimension_inv * entropy_minus_one_over_gamma; p->density.pressure_dh *=
h * rho_inv * hydro_dimension_inv * entropy_minus_one_over_gamma;
/* Finish calculation of the velocity curl components */ /* Finish calculation of the velocity curl components */
// p->density.rot_v[0] *= 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 * irho; p->density.rot_v[1] *= h_inv_dim_plus_one * rho_inv;
// p->density.rot_v[2] *= h_inv_dim_plus_one * irho; p->density.rot_v[2] *= h_inv_dim_plus_one * rho_inv;
/* Finish calculation of the velocity divergence */ /* 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( ...@@ -273,6 +298,16 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
struct part *restrict p, struct xpart *restrict xp, int ti_current, struct part *restrict p, struct xpart *restrict xp, int ti_current,
double timeBase) { 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 */ /* Compute the pressure */
const float half_dt = (ti_current - (p->ti_begin + p->ti_end) / 2) * timeBase; const float half_dt = (ti_current - (p->ti_begin + p->ti_end) / 2) * timeBase;
const float pressure = hydro_get_pressure(p, half_dt); const float pressure = hydro_get_pressure(p, half_dt);
...@@ -281,9 +316,18 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force( ...@@ -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 rho_bar_inv = 1.f / p->rho_bar;
const float soundspeed = sqrtf(hydro_gamma * pressure * rho_bar_inv); 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. */ /* Update variables. */
p->force.soundspeed = soundspeed; p->force.soundspeed = soundspeed;
p->force.P_over_rho2 = pressure * rho_bar_inv * rho_bar_inv; 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( ...@@ -422,7 +466,7 @@ __attribute__((always_inline)) INLINE static void hydro_convert_quantities(
struct part *restrict p) { struct part *restrict p) {
/* We read u in the entropy field. We now get S from u */ /* 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); p->entropy_one_over_gamma = pow_one_over_gamma(p->entropy);
/* Compute the pressure */ /* Compute the pressure */
......
src/hydro/PressureEntropy/hydro_debug.h
+ 9
− 7
View file @ aefe3c71
...@@ -26,7 +26,8 @@ ...@@ -26,7 +26,8 @@
* The thermal variable is the entropy (S) and the entropy is smoothed over * The thermal variable is the entropy (S) and the entropy is smoothed over
* contact discontinuities to prevent spurious surface tension. * 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( __attribute__((always_inline)) INLINE static void hydro_debug_particle(
...@@ -35,14 +36,15 @@ __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], " "x=[%.3e,%.3e,%.3e], "
"v=[%.3e,%.3e,%.3e],v_full=[%.3e,%.3e,%.3e] \n a=[%.3e,%.3e,%.3e],\n " "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, " "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" "rho_bar=%.3e, P=%.3e, dP_dh=%.3e, P_over_rho2=%.3e, S=%.3e, S^1/g=%.3e, "
"v_sig=%e dh/dt=%.3e t_begin=%d, t_end=%d\n", "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], 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], 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->h, p->density.wcount, p->density.wcount_dh, p->mass, p->density.rho_dh,
p->rho, hydro_get_pressure(p, 0.), 0.f, p->rho, p->rho_bar, hydro_get_pressure(p, 0.), p->density.pressure_dh,
/*p->force.P_over_rho2,*/ p->entropy, p->entropy_dt, p->force.soundspeed, p->force.P_over_rho2, p->entropy, p->entropy_one_over_gamma,
p->force.v_sig, p->force.h_dt, p->ti_begin, p->ti_end); 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 */ #endif /* SWIFT_PRESSURE_ENTROPY_HYDRO_DEBUG_H */
src/hydro/PressureEntropy/hydro_iact.h
+ 58
− 62
View file @ aefe3c71
...@@ -26,7 +26,8 @@ ...@@ -26,7 +26,8 @@
* The thermal variable is the entropy (S) and the entropy is smoothed over * The thermal variable is the entropy (S) and the entropy is smoothed over
* contact discontinuities to prevent spurious surface tension. * 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( ...@@ -37,19 +38,15 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
float wi, wi_dx; float wi, wi_dx;
float wj, wj_dx; float wj, wj_dx;
/* float dv[3], curlvr[3]; */ float dv[3], curlvr[3];
/* Get the masses. */ /* Get the masses. */
const float mi = pi->mass; const float mi = pi->mass;
const float mj = pj->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. */ /* Get r and r inverse. */
const float r = sqrtf(r2); 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 */ /* Compute the kernel function for pi */
const float hi_inv = 1.f / hi; const float hi_inv = 1.f / hi;
...@@ -58,7 +55,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_density( ...@@ -58,7 +55,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
/* Compute contribution to the density */ /* Compute contribution to the density */
pi->rho += mj * wi; 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 */ /* Compute contribution to the number of neighbours */
pi->density.wcount += wi; pi->density.wcount += wi;
...@@ -66,7 +63,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_density( ...@@ -66,7 +63,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
/* Compute contribution to the weighted density */ /* Compute contribution to the weighted density */
pi->rho_bar += mj * pj->entropy_one_over_gamma * wi; 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); mj * pj->entropy_one_over_gamma * (hydro_dimension * wi + ui * wi_dx);
/* Compute the kernel function for pj */ /* Compute the kernel function for pj */
...@@ -76,7 +73,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_density( ...@@ -76,7 +73,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
/* Compute contribution to the density */ /* Compute contribution to the density */
pj->rho += mi * wj; 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 */ /* Compute contribution to the number of neighbours */
pj->density.wcount += wj; pj->density.wcount += wj;
...@@ -84,33 +81,33 @@ __attribute__((always_inline)) INLINE static void runner_iact_density( ...@@ -84,33 +81,33 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
/* Compute contribution to the weighted density */ /* Compute contribution to the weighted density */
pj->rho_bar += mi * pi->entropy_one_over_gamma * wj; 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); mi * pi->entropy_one_over_gamma * (hydro_dimension * wj + uj * wj_dx);
/* const float faci = mj * wi_dx * r_inv; */ const float faci = mj * wi_dx * r_inv;
/* const float facj = mi * wj_dx * r_inv; */ const float facj = mi * wj_dx * r_inv;
/* /\* Compute dv dot r *\/ */ /* Compute dv dot r */
/* dv[0] = pi->v[0] - pj->v[0]; */ dv[0] = pi->v[0] - pj->v[0];
/* dv[1] = pi->v[1] - pj->v[1]; */ dv[1] = pi->v[1] - pj->v[1];
/* dv[2] = pi->v[2] - pj->v[2]; */ dv[2] = pi->v[2] - pj->v[2];
/* const float dvdr = dv[0] * dx[0] + dv[1] * dx[1] + dv[2] * dx[2]; */ const float dvdr = dv[0] * dx[0] + dv[1] * dx[1] + dv[2] * dx[2];
/* pi->density.div_v -= faci * dvdr; */ pi->density.div_v -= faci * dvdr;
/* pj->density.div_v -= facj * dvdr; */ pj->density.div_v -= facj * dvdr;
/* /\* Compute dv cross r *\/ */ /* Compute dv cross r */
/* curlvr[0] = dv[1] * dx[2] - dv[2] * dx[1]; */ curlvr[0] = dv[1] * dx[2] - dv[2] * dx[1];
/* curlvr[1] = dv[2] * dx[0] - dv[0] * dx[2]; */ curlvr[1] = dv[2] * dx[0] - dv[0] * dx[2];
/* curlvr[2] = dv[0] * dx[1] - dv[1] * dx[0]; */ curlvr[2] = dv[0] * dx[1] - dv[1] * dx[0];
/* pi->density.rot_v[0] += faci * curlvr[0]; */ pi->density.rot_v[0] += faci * curlvr[0];
/* pi->density.rot_v[1] += faci * curlvr[1]; */ pi->density.rot_v[1] += faci * curlvr[1];
/* pi->density.rot_v[2] += faci * curlvr[2]; */ pi->density.rot_v[2] += faci * curlvr[2];
/* pj->density.rot_v[0] += facj * curlvr[0]; */ pj->density.rot_v[0] += facj * curlvr[0];
/* pj->density.rot_v[1] += facj * curlvr[1]; */ pj->density.rot_v[1] += facj * curlvr[1];
/* pj->density.rot_v[2] += facj * curlvr[2]; */ pj->density.rot_v[2] += facj * curlvr[2];
} }
/** /**
...@@ -130,17 +127,14 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density( ...@@ -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 r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {
float wi, wi_dx; float wi, wi_dx;
/* float dv[3], curlvr[3]; */ float dv[3], curlvr[3];
/* Get the masses. */ /* Get the masses. */
const float mj = pj->mass; const float mj = pj->mass;
/* Get the entropies. */
// const float entropy_j = pj->entropy;
/* Get r and r inverse. */ /* Get r and r inverse. */
const float r = sqrtf(r2); const float r = sqrtf(r2);
/* const float ri = 1.0f / r; */ const float ri = 1.0f / r;
/* Compute the kernel function */ /* Compute the kernel function */
const float h_inv = 1.0f / hi; const float h_inv = 1.0f / hi;
...@@ -149,34 +143,34 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density( ...@@ -149,34 +143,34 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
/* Compute contribution to the density */ /* Compute contribution to the density */
pi->rho += mj * wi; 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 */ /* Compute contribution to the number of neighbours */
pi->density.wcount += wi; 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 */ /* Compute contribution to the weighted density */
pi->rho_bar += mj * pj->entropy_one_over_gamma * wi; 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); 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 *\/ */ /* Compute dv dot r */
/* dv[0] = pi->v[0] - pj->v[0]; */ dv[0] = pi->v[0] - pj->v[0];
/* dv[1] = pi->v[1] - pj->v[1]; */ dv[1] = pi->v[1] - pj->v[1];
/* dv[2] = pi->v[2] - pj->v[2]; */ dv[2] = pi->v[2] - pj->v[2];
/* const float dvdr = dv[0] * dx[0] + dv[1] * dx[1] + dv[2] * dx[2]; */ const float dvdr = dv[0] * dx[0] + dv[1] * dx[1] + dv[2] * dx[2];
/* pi->density.div_v -= fac * dvdr; */ pi->density.div_v -= fac * dvdr;
/* /\* Compute dv cross r *\/ */ /* Compute dv cross r */
/* curlvr[0] = dv[1] * dx[2] - dv[2] * dx[1]; */ curlvr[0] = dv[1] * dx[2] - dv[2] * dx[1];
/* curlvr[1] = dv[2] * dx[0] - dv[0] * dx[2]; */ curlvr[1] = dv[2] * dx[0] - dv[0] * dx[2];
/* curlvr[2] = dv[0] * dx[1] - dv[1] * dx[0]; */ curlvr[2] = dv[0] * dx[1] - dv[1] * dx[0];
/* pi->density.rot_v[0] += fac * curlvr[0]; */ pi->density.rot_v[0] += fac * curlvr[0];
/* pi->density.rot_v[1] += fac * curlvr[1]; */ pi->density.rot_v[1] += fac * curlvr[1];
/* pi->density.rot_v[2] += fac * curlvr[2]; */ pi->density.rot_v[2] += fac * curlvr[2];
} }
/** /**
...@@ -223,8 +217,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_force( ...@@ -223,8 +217,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
const float wj_dr = hjd_inv * wj_dx; const float wj_dr = hjd_inv * wj_dx;
/* Compute gradient terms */ /* Compute gradient terms */
const float f_i = pi->rho_dh * pi->pressure_dh; const float f_i = pi->force.f;
const float f_j = pj->rho_dh * pj->pressure_dh; const float f_j = pj->force.f;
/* Compute Pressure terms */ /* Compute Pressure terms */
const float P_over_rho2_i = pi->force.P_over_rho2; const float P_over_rho2_i = pi->force.P_over_rho2;
...@@ -244,8 +238,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_force( ...@@ -244,8 +238,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
(pi->v[2] - pj->v[2]) * dx[2]; (pi->v[2] - pj->v[2]) * dx[2];
/* Balsara term */ /* Balsara term */
/* const float balsara_i = pi->force.balsara; */ const float balsara_i = pi->force.balsara;
/* const float balsara_j = pj->force.balsara; */ const float balsara_j = pj->force.balsara;
/* Are the particles moving towards each others ? */ /* Are the particles moving towards each others ? */
const float omega_ij = fminf(dvdr, 0.f); const float omega_ij = fminf(dvdr, 0.f);
...@@ -256,7 +250,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_force( ...@@ -256,7 +250,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
/* Now construct the full viscosity term */ /* Now construct the full viscosity term */
const float rho_ij = 0.5f * (rhoi + rhoj); 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 */ /* Now, convolve with the kernel */
const float visc_term = 0.5f * visc * (wi_dr + wj_dr); 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( ...@@ -332,8 +327,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
const float wj_dr = hjd_inv * wj_dx; const float wj_dr = hjd_inv * wj_dx;
/* Compute gradient terms */ /* Compute gradient terms */
const float f_i = pi->rho_dh * pi->pressure_dh; const float f_i = pi->force.f;
const float f_j = pj->rho_dh * pj->pressure_dh; const float f_j = pj->force.f;
/* Compute Pressure terms */ /* Compute Pressure terms */
const float P_over_rho2_i = pi->force.P_over_rho2; 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( ...@@ -353,8 +348,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
(pi->v[2] - pj->v[2]) * dx[2]; (pi->v[2] - pj->v[2]) * dx[2];
/* Balsara term */ /* Balsara term */
/* const float balsara_i = pi->force.balsara; */ const float balsara_i = pi->force.balsara;
/* const float balsara_j = pj->force.balsara; */ const float balsara_j = pj->force.balsara;
/* Are the particles moving towards each others ? */ /* Are the particles moving towards each others ? */
const float omega_ij = fminf(dvdr, 0.f); const float omega_ij = fminf(dvdr, 0.f);
...@@ -365,7 +360,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force( ...@@ -365,7 +360,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
/* Now construct the full viscosity term */ /* Now construct the full viscosity term */
const float rho_ij = 0.5f * (rhoi + rhoj); 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 */ /* Now, convolve with the kernel */
const float visc_term = 0.5f * visc * (wi_dr + wj_dr); const float visc_term = 0.5f * visc * (wi_dr + wj_dr);
......
src/hydro/PressureEntropy/hydro_io.h
+ 2
− 1
View file @ aefe3c71
...@@ -26,7 +26,8 @@ ...@@ -26,7 +26,8 @@
* The thermal variable is the entropy (S) and the entropy is smoothed over * The thermal variable is the entropy (S) and the entropy is smoothed over
* contact discontinuities to prevent spurious surface tension. * 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" #include "adiabatic_index.h"
......
src/hydro/PressureEntropy/hydro_part.h
+ 12
− 14
View file @ aefe3c71
...@@ -26,7 +26,8 @@ ...@@ -26,7 +26,8 @@
* The thermal variable is the entropy (S) and the entropy is smoothed over * The thermal variable is the entropy (S) and the entropy is smoothed over
* contact discontinuities to prevent spurious surface tension. * 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" #include "cooling_struct.h"
...@@ -72,15 +73,6 @@ struct part { ...@@ -72,15 +73,6 @@ struct part {
/*! Particle density. */ /*! Particle density. */
float rho; 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 */ /*! Particle weighted density */
float rho_bar; float rho_bar;
...@@ -103,11 +95,14 @@ struct part { ...@@ -103,11 +95,14 @@ struct part {
/*! Number of neighbours spatial derivative. */ /*! Number of neighbours spatial derivative. */
float wcount_dh; float wcount_dh;
/*! Derivative of particle weighted pressure with h. */ /*! Derivative of density with respect to h */
// float weightedPressure_dh; float rho_dh;
/*! Derivative of pressure with respect to h */
float pressure_dh;
/*! Particle velocity curl. */ /*! Particle velocity curl. */
// float rot_v[3]; float rot_v[3];
/*! Particle velocity divergence. */ /*! Particle velocity divergence. */
float div_v; float div_v;
...@@ -116,8 +111,11 @@ struct part { ...@@ -116,8 +111,11 @@ struct part {
struct { struct {
/*! Balsara switch */
float balsara;
/*! "Grad h" term */ /*! "Grad h" term */
// float f_ij; float f;
/*! Pressure term */ /*! Pressure term */
float P_over_rho2; float P_over_rho2;
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment