Skip to content

GitLab

Commit 8bf391d4 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

P_over_rho --> P_over_rho2

parent 579822fc
Hide whitespace changes
Inline Side-by-side
src/hydro/Gadget2/hydro.h
View file @ 8bf391d4
......@@ -135,22 +135,28 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
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 dt = (ti_current - (p->ti_begin + p->ti_end) * 0.5f) * timeBase;
const float pressure = (p->entropy + p->force.entropy_dt * dt) * pow_gamma(p->rho);
const float dt = (ti_current - (p->ti_begin + p->ti_end) / 2) * timeBase;
const float pressure =
(p->entropy + p->force.entropy_dt * dt) * pow_gamma(p->rho);
const float irho = 1.f / p->rho;
/* Divide the pressure by the density and density gradient */
const float P_over_rho = pressure / (p->rho * p->rho) * p->rho_dh;
const float P_over_rho2 = pressure * irho * irho * p->rho_dh;
/* Compute the sound speed */
const float soundspeed = sqrtf(hydro_gamma * pressure / p->rho);
const float soundspeed = sqrtf(hydro_gamma * pressure * irho);
/* Compute the Balsara switch */
float balsara = fabsf(p->density.div_v) /
(fabsf(p->density.div_v) + curl_v + 0.0001f * p->force.soundspeed / fac_mu / p->h);
const float balsara =
abs_div_v / (abs_div_v + curl_v + 0.0001f * soundspeed / fac_mu / p->h);
/* Update variables. */
p->force.P_over_rho = P_over_rho;
p->force.P_over_rho2 = P_over_rho2;
p->force.soundspeed = soundspeed;
p->force.balsara = balsara;
}
......@@ -196,16 +202,18 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
/* Drift the pressure */
const float dt_entr = (t1 - (p->ti_begin + p->ti_end) / 2) * timeBase;
const float pressure =
(p->entropy + p->force.entropy_dt * dt_entr) * pow_gamma(p->rho);
(p->entropy + p->force.entropy_dt * dt_entr) * pow_gamma(p->rho);
/* Divide the pressure by the density and density gradient */
const float P_over_rho = pressure / (p->rho * p->rho) * p->rho_dh;
const float irho = 1.f / p->rho;
/* Divide the pressure by the density and density gradient */
const float P_over_rho2 = pressure * irho * irho * p->rho_dh;
/* Compute the new sound speed */
const float soundspeed = sqrtf(hydro_gamma * pressure / p->rho);
const float soundspeed = sqrtf(hydro_gamma * pressure * irho);
/* Update variables */
p->force.P_over_rho = P_over_rho;
p->force.P_over_rho2 = P_over_rho2;
p->force.soundspeed = soundspeed;
}
......@@ -219,7 +227,8 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
__attribute__((always_inline)) INLINE static void hydro_end_force(
struct part *restrict p) {
p->force.entropy_dt *= hydro_gamma_minus_one * pow_minus_gamma_minus_one(p->rho);
p->force.entropy_dt *=
hydro_gamma_minus_one * pow_minus_gamma_minus_one(p->rho);
}
/**
......
src/hydro/Gadget2/hydro_debug.h
View file @ 8bf391d4
......@@ -23,7 +23,8 @@ __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=%d, wcount_dh=%.3e, m=%.3e, dh_drho=%.3e, rho=%.3e, P_over_rho=%.3e, "
"wcount=%d, wcount_dh=%.3e, m=%.3e, dh_drho=%.3e, rho=%.3e, "
"P_over_rho2=%.3e, "
"S=%.3e, "
"dS/dt=%.3e, c=%.3e\n"
"divV=%.3e, rotV=[%.3e,%.3e,%.3e], balsara=%.3e \n "
......@@ -31,7 +32,8 @@ __attribute__((always_inline)) INLINE static void hydro_debug_particle(
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, (int)p->density.wcount, p->density.wcount_dh, p->mass, p->rho_dh,
p->rho, p->force.P_over_rho, p->entropy, p->force.entropy_dt, p->force.soundspeed,
p->density.div_v, p->density.rot_v[0], p->density.rot_v[1],
p->density.rot_v[2], p->force.balsara, p->force.v_sig, p->h_dt, p->ti_begin, p->ti_end);
p->rho, p->force.P_over_rho2, p->entropy, p->force.entropy_dt,
p->force.soundspeed, p->density.div_v, p->density.rot_v[0],
p->density.rot_v[1], p->density.rot_v[2], p->force.balsara,
p->force.v_sig, p->h_dt, p->ti_begin, p->ti_end);
}
src/hydro/Gadget2/hydro_iact.h
View file @ 8bf391d4
......@@ -134,7 +134,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_vec_density(
vj[k].v = vec_set(pj[0]->v[k], pj[1]->v[k], pj[2]->v[k], pj[3]->v[k],
pj[4]->v[k], pj[5]->v[k], pj[6]->v[k], pj[7]->v[k]);
}
/* Get each component of particle separation. (Dx={dx1,dy1,dz1,dx2,dy2,dz2,...,dxn,dyn,dzn})*/
/* Get each component of particle separation.
* (Dx={dx1,dy1,dz1,dx2,dy2,dz2,...,dxn,dyn,dzn})*/
for (k = 0; k < 3; k++)
dx[k].v = vec_set(Dx[0 + k], Dx[3 + k], Dx[6 + k], Dx[9 + k], Dx[12 + k],
Dx[15 + k], Dx[18 + k], Dx[21 + k]);
......@@ -154,7 +155,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_vec_density(
/* Get the radius and inverse radius. */
r2.v = vec_load(R2);
ri.v = vec_rsqrt(r2.v);
/*vec_rsqrt does not have the level of accuracy we need, so an extra term is added below.*/
/*vec_rsqrt does not have the level of accuracy we need, so an extra term is
* added below.*/
ri.v = ri.v - vec_set1(0.5f) * ri.v * (r2.v * ri.v * ri.v - vec_set1(1.0f));
r.v = r2.v * ri.v;
......@@ -221,10 +223,11 @@ __attribute__((always_inline)) INLINE static void runner_iact_vec_density(
#else
error("The Gadget2 serial version of runner_iact_density was called when the vectorised version should have been used.")
error(
"The Gadget2 serial version of runner_iact_density was called when the "
"vectorised version should have been used.")
#endif
}
/**
......@@ -304,7 +307,8 @@ runner_iact_nonsym_vec_density(float *R2, float *Dx, float *Hi, float *Hj,
vj[k].v = vec_set(pj[0]->v[k], pj[1]->v[k], pj[2]->v[k], pj[3]->v[k],
pj[4]->v[k], pj[5]->v[k], pj[6]->v[k], pj[7]->v[k]);
}
/* Get each component of particle separation. (Dx={dx1,dy1,dz1,dx2,dy2,dz2,...,dxn,dyn,dzn})*/
/* Get each component of particle separation.
* (Dx={dx1,dy1,dz1,dx2,dy2,dz2,...,dxn,dyn,dzn})*/
for (k = 0; k < 3; k++)
dx[k].v = vec_set(Dx[0 + k], Dx[3 + k], Dx[6 + k], Dx[9 + k], Dx[12 + k],
Dx[15 + k], Dx[18 + k], Dx[21 + k]);
......@@ -323,7 +327,8 @@ runner_iact_nonsym_vec_density(float *R2, float *Dx, float *Hi, float *Hj,
/* Get the radius and inverse radius. */
r2.v = vec_load(R2);
ri.v = vec_rsqrt(r2.v);
/*vec_rsqrt does not have the level of accuracy we need, so an extra term is added below.*/
/*vec_rsqrt does not have the level of accuracy we need, so an extra term is
* added below.*/
ri.v = ri.v - vec_set1(0.5f) * ri.v * (r2.v * ri.v * ri.v - vec_set1(1.0f));
r.v = r2.v * ri.v;
......@@ -369,10 +374,11 @@ runner_iact_nonsym_vec_density(float *R2, float *Dx, float *Hi, float *Hj,
#else
error("The Gadget2 serial version of runner_iact_nonsym_density was called when the vectorised version should have been used.")
error(
"The Gadget2 serial version of runner_iact_nonsym_density was called "
"when the vectorised version should have been used.")
#endif
}
/**
......@@ -409,8 +415,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
const float wj_dr = hj2_inv * hj2_inv * wj_dx;
/* Compute gradient terms */
const float P_over_rho_i = pi->force.P_over_rho;
const float P_over_rho_j = pj->force.P_over_rho;
const float P_over_rho2_i = pi->force.P_over_rho2;
const float P_over_rho2_j = pj->force.P_over_rho2;
/* Compute sound speeds */
const float ci = pi->force.soundspeed;
......@@ -424,7 +430,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
/* Balsara term */
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);
const float mu_ij = fac_mu * r_inv * omega_ij; /* This is 0 or negative */
......@@ -439,7 +445,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
/* Now, convolve with the kernel */
const float visc_term = 0.5f * visc * (wi_dr + wj_dr) * r_inv;
const float sph_term = (P_over_rho_i * wi_dr + P_over_rho_j * wj_dr) * r_inv;
const float sph_term = (P_over_rho2_i * wi_dr + P_over_rho2_j * wj_dr) * r_inv;
/* Eventually got the acceleration */
const float acc = visc_term + sph_term;
......@@ -511,8 +517,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
const float wj_dr = hj2_inv * hj2_inv * wj_dx;
/* Compute gradient terms */
const float P_over_rho_i = pi->force.P_over_rho;
const float P_over_rho_j = pj->force.P_over_rho;
const float P_over_rho2_i = pi->force.P_over_rho2;
const float P_over_rho2_j = pj->force.P_over_rho2;
/* Compute sound speeds */
const float ci = pi->force.soundspeed;
......@@ -541,7 +547,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
/* Now, convolve with the kernel */
const float visc_term = 0.5f * visc * (wi_dr + wj_dr) * r_inv;
const float sph_term = (P_over_rho_i * wi_dr + P_over_rho_j * wj_dr) * r_inv;
const float sph_term = (P_over_rho2_i * wi_dr + P_over_rho2_j * wj_dr) * r_inv;
/* Eventually got the acceleration */
const float acc = visc_term + sph_term;
......
src/hydro/Gadget2/hydro_part.h
View file @ 8bf391d4
......@@ -25,8 +25,6 @@ struct xpart {
/* Velocity at the last full step. */
float v_full[3];
} __attribute__((aligned(xpart_align)));
......@@ -90,15 +88,15 @@ struct part {
/* Balsara switch */
float balsara;
/* Signal velocity. */
float v_sig;
/* Pressure over density*/
float P_over_rho;
float P_over_rho2;
/* Particle sound speed. */
float soundspeed;
/* Signal velocity. */
float v_sig;
/* Entropy time derivative */
float entropy_dt;
......
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