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

Reduce the number of divisions in the SPH gradient calculation of gismo.

parent 0444cb29
......@@ -97,16 +97,12 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_predict(
float dWi[5], dWj[5];
float xij_j[3];
int k;
float xfac;
/* perform gradient reconstruction in space and time */
/* space */
/* Compute interface position (relative to pj, since we don't need the actual
* position) */
/* eqn. (8) */
xfac = hj / (hi + hj);
for (k = 0; k < 3; k++) xij_j[k] = xfac * dx[k];
* position) eqn. (8) */
const float xfac = hj / (hi + hj);
for (int k = 0; k < 3; k++) xij_j[k] = xfac * dx[k];
dWi[0] = pi->primitives.gradients.rho[0] * xij_i[0] +
pi->primitives.gradients.rho[1] * xij_i[1] +
......@@ -140,6 +136,7 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_predict(
pj->primitives.gradients.P[1] * xij_j[1] +
pj->primitives.gradients.P[2] * xij_j[2];
/* Apply the slope limiter at this interface */
hydro_slope_limit_face(Wi, Wj, dWi, dWj, xij_i, xij_j, r);
Wi[0] += dWi[0];
......@@ -154,6 +151,7 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_predict(
Wj[3] += dWj[3];
Wj[4] += dWj[4];
/* Check that we don't have problematic densities or pressures */
gizmo_check_physical_quantity("density", Wi[0]);
gizmo_check_physical_quantity("pressure", Wi[4]);
gizmo_check_physical_quantity("density", Wj[0]);
......
......@@ -65,18 +65,17 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_collect(
float r2, const float *dx, float hi, float hj, struct part *restrict pi,
struct part *restrict pj) {
float r = sqrtf(r2);
float xi, xj;
float hi_inv, hj_inv;
const float r_inv = 1.f / sqrtf(r2);
const float r = r2 * r_inv;
float wi, wj, wi_dx, wj_dx;
int k, l;
float Bi[3][3];
float Bj[3][3];
float Wi[5], Wj[5];
/* Initialize local variables */
for (k = 0; k < 3; k++) {
for (l = 0; l < 3; l++) {
for (int k = 0; k < 3; k++) {
for (int l = 0; l < 3; l++) {
Bi[k][l] = pi->geometry.matrix_E[k][l];
Bj[k][l] = pj->geometry.matrix_E[k][l];
}
......@@ -93,8 +92,8 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_collect(
Wj[4] = pj->primitives.P;
/* Compute kernel of pi. */
hi_inv = 1.f / hi;
xi = r * hi_inv;
const float hi_inv = 1.f / hi;
const float xi = r * hi_inv;
kernel_deval(xi, &wi, &wi_dx);
if (pi->density.wcorr > const_gizmo_min_wcorr) {
......@@ -153,46 +152,46 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_collect(
/* The gradient matrix was not well-behaved, switch to SPH gradients */
pi->primitives.gradients.rho[0] -=
wi_dx * dx[0] * (pi->primitives.rho - pj->primitives.rho) / r;
wi_dx * dx[0] * (pi->primitives.rho - pj->primitives.rho) * r_inv;
pi->primitives.gradients.rho[1] -=
wi_dx * dx[1] * (pi->primitives.rho - pj->primitives.rho) / r;
wi_dx * dx[1] * (pi->primitives.rho - pj->primitives.rho) * r_inv;
pi->primitives.gradients.rho[2] -=
wi_dx * dx[2] * (pi->primitives.rho - pj->primitives.rho) / r;
wi_dx * dx[2] * (pi->primitives.rho - pj->primitives.rho) * r_inv;
pi->primitives.gradients.v[0][0] -=
wi_dx * dx[0] * (pi->primitives.v[0] - pj->primitives.v[0]) / r;
wi_dx * dx[0] * (pi->primitives.v[0] - pj->primitives.v[0]) * r_inv;
pi->primitives.gradients.v[0][1] -=
wi_dx * dx[1] * (pi->primitives.v[0] - pj->primitives.v[0]) / r;
wi_dx * dx[1] * (pi->primitives.v[0] - pj->primitives.v[0]) * r_inv;
pi->primitives.gradients.v[0][2] -=
wi_dx * dx[2] * (pi->primitives.v[0] - pj->primitives.v[0]) / r;
wi_dx * dx[2] * (pi->primitives.v[0] - pj->primitives.v[0]) * r_inv;
pi->primitives.gradients.v[1][0] -=
wi_dx * dx[0] * (pi->primitives.v[1] - pj->primitives.v[1]) / r;
wi_dx * dx[0] * (pi->primitives.v[1] - pj->primitives.v[1]) * r_inv;
pi->primitives.gradients.v[1][1] -=
wi_dx * dx[1] * (pi->primitives.v[1] - pj->primitives.v[1]) / r;
wi_dx * dx[1] * (pi->primitives.v[1] - pj->primitives.v[1]) * r_inv;
pi->primitives.gradients.v[1][2] -=
wi_dx * dx[2] * (pi->primitives.v[1] - pj->primitives.v[1]) / r;
wi_dx * dx[2] * (pi->primitives.v[1] - pj->primitives.v[1]) * r_inv;
pi->primitives.gradients.v[2][0] -=
wi_dx * dx[0] * (pi->primitives.v[2] - pj->primitives.v[2]) / r;
wi_dx * dx[0] * (pi->primitives.v[2] - pj->primitives.v[2]) * r_inv;
pi->primitives.gradients.v[2][1] -=
wi_dx * dx[1] * (pi->primitives.v[2] - pj->primitives.v[2]) / r;
wi_dx * dx[1] * (pi->primitives.v[2] - pj->primitives.v[2]) * r_inv;
pi->primitives.gradients.v[2][2] -=
wi_dx * dx[2] * (pi->primitives.v[2] - pj->primitives.v[2]) / r;
wi_dx * dx[2] * (pi->primitives.v[2] - pj->primitives.v[2]) * r_inv;
pi->primitives.gradients.P[0] -=
wi_dx * dx[0] * (pi->primitives.P - pj->primitives.P) / r;
wi_dx * dx[0] * (pi->primitives.P - pj->primitives.P) * r_inv;
pi->primitives.gradients.P[1] -=
wi_dx * dx[1] * (pi->primitives.P - pj->primitives.P) / r;
wi_dx * dx[1] * (pi->primitives.P - pj->primitives.P) * r_inv;
pi->primitives.gradients.P[2] -=
wi_dx * dx[2] * (pi->primitives.P - pj->primitives.P) / r;
wi_dx * dx[2] * (pi->primitives.P - pj->primitives.P) * r_inv;
}
hydro_slope_limit_cell_collect(pi, pj, r);
/* Compute kernel of pj. */
hj_inv = 1.f / hj;
xj = r * hj_inv;
const float hj_inv = 1.f / hj;
const float xj = r * hj_inv;
kernel_deval(xj, &wj, &wj_dx);
if (pj->density.wcorr > const_gizmo_min_wcorr) {
......@@ -252,38 +251,38 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_collect(
/* SPH gradients */
pj->primitives.gradients.rho[0] -=
wj_dx * dx[0] * (pi->primitives.rho - pj->primitives.rho) / r;
wj_dx * dx[0] * (pi->primitives.rho - pj->primitives.rho) * r_inv;
pj->primitives.gradients.rho[1] -=
wj_dx * dx[1] * (pi->primitives.rho - pj->primitives.rho) / r;
wj_dx * dx[1] * (pi->primitives.rho - pj->primitives.rho) * r_inv;
pj->primitives.gradients.rho[2] -=
wj_dx * dx[2] * (pi->primitives.rho - pj->primitives.rho) / r;
wj_dx * dx[2] * (pi->primitives.rho - pj->primitives.rho) * r_inv;
pj->primitives.gradients.v[0][0] -=
wj_dx * dx[0] * (pi->primitives.v[0] - pj->primitives.v[0]) / r;
wj_dx * dx[0] * (pi->primitives.v[0] - pj->primitives.v[0]) * r_inv;
pj->primitives.gradients.v[0][1] -=
wj_dx * dx[1] * (pi->primitives.v[0] - pj->primitives.v[0]) / r;
wj_dx * dx[1] * (pi->primitives.v[0] - pj->primitives.v[0]) * r_inv;
pj->primitives.gradients.v[0][2] -=
wj_dx * dx[2] * (pi->primitives.v[0] - pj->primitives.v[0]) / r;
wj_dx * dx[2] * (pi->primitives.v[0] - pj->primitives.v[0]) * r_inv;
pj->primitives.gradients.v[1][0] -=
wj_dx * dx[0] * (pi->primitives.v[1] - pj->primitives.v[1]) / r;
wj_dx * dx[0] * (pi->primitives.v[1] - pj->primitives.v[1]) * r_inv;
pj->primitives.gradients.v[1][1] -=
wj_dx * dx[1] * (pi->primitives.v[1] - pj->primitives.v[1]) / r;
wj_dx * dx[1] * (pi->primitives.v[1] - pj->primitives.v[1]) * r_inv;
pj->primitives.gradients.v[1][2] -=
wj_dx * dx[2] * (pi->primitives.v[1] - pj->primitives.v[1]) / r;
wj_dx * dx[2] * (pi->primitives.v[1] - pj->primitives.v[1]) * r_inv;
pj->primitives.gradients.v[2][0] -=
wj_dx * dx[0] * (pi->primitives.v[2] - pj->primitives.v[2]) / r;
wj_dx * dx[0] * (pi->primitives.v[2] - pj->primitives.v[2]) * r_inv;
pj->primitives.gradients.v[2][1] -=
wj_dx * dx[1] * (pi->primitives.v[2] - pj->primitives.v[2]) / r;
wj_dx * dx[1] * (pi->primitives.v[2] - pj->primitives.v[2]) * r_inv;
pj->primitives.gradients.v[2][2] -=
wj_dx * dx[2] * (pi->primitives.v[2] - pj->primitives.v[2]) / r;
wj_dx * dx[2] * (pi->primitives.v[2] - pj->primitives.v[2]) * r_inv;
pj->primitives.gradients.P[0] -=
wj_dx * dx[0] * (pi->primitives.P - pj->primitives.P) / r;
wj_dx * dx[0] * (pi->primitives.P - pj->primitives.P) * r_inv;
pj->primitives.gradients.P[1] -=
wj_dx * dx[1] * (pi->primitives.P - pj->primitives.P) / r;
wj_dx * dx[1] * (pi->primitives.P - pj->primitives.P) * r_inv;
pj->primitives.gradients.P[2] -=
wj_dx * dx[2] * (pi->primitives.P - pj->primitives.P) / r;
wj_dx * dx[2] * (pi->primitives.P - pj->primitives.P) * r_inv;
}
hydro_slope_limit_cell_collect(pj, pi, r);
......@@ -304,17 +303,15 @@ hydro_gradients_nonsym_collect(float r2, const float *dx, float hi, float hj,
struct part *restrict pi,
struct part *restrict pj) {
float r = sqrtf(r2);
float xi;
float hi_inv;
float wi, wi_dx;
int k, l;
const float r_inv = 1.f / sqrtf(r2);
const float r = r2 * r_inv;
float Bi[3][3];
float Wi[5], Wj[5];
/* Initialize local variables */
for (k = 0; k < 3; k++) {
for (l = 0; l < 3; l++) {
for (int k = 0; k < 3; k++) {
for (int l = 0; l < 3; l++) {
Bi[k][l] = pi->geometry.matrix_E[k][l];
}
}
......@@ -330,8 +327,9 @@ hydro_gradients_nonsym_collect(float r2, const float *dx, float hi, float hj,
Wj[4] = pj->primitives.P;
/* Compute kernel of pi. */
hi_inv = 1.f / hi;
xi = r * hi_inv;
float wi, wi_dx;
const float hi_inv = 1.f / hi;
const float xi = r * hi_inv;
kernel_deval(xi, &wi, &wi_dx);
if (pi->density.wcorr > const_gizmo_min_wcorr) {
......@@ -390,38 +388,38 @@ hydro_gradients_nonsym_collect(float r2, const float *dx, float hi, float hj,
/* Gradient matrix is not well-behaved, switch to SPH gradients */
pi->primitives.gradients.rho[0] -=
wi_dx * dx[0] * (pi->primitives.rho - pj->primitives.rho) / r;
wi_dx * dx[0] * (pi->primitives.rho - pj->primitives.rho) * r_inv;
pi->primitives.gradients.rho[1] -=
wi_dx * dx[1] * (pi->primitives.rho - pj->primitives.rho) / r;
wi_dx * dx[1] * (pi->primitives.rho - pj->primitives.rho) * r_inv;
pi->primitives.gradients.rho[2] -=
wi_dx * dx[2] * (pi->primitives.rho - pj->primitives.rho) / r;
wi_dx * dx[2] * (pi->primitives.rho - pj->primitives.rho) * r_inv;
pi->primitives.gradients.v[0][0] -=
wi_dx * dx[0] * (pi->primitives.v[0] - pj->primitives.v[0]) / r;
wi_dx * dx[0] * (pi->primitives.v[0] - pj->primitives.v[0]) * r_inv;
pi->primitives.gradients.v[0][1] -=
wi_dx * dx[1] * (pi->primitives.v[0] - pj->primitives.v[0]) / r;
wi_dx * dx[1] * (pi->primitives.v[0] - pj->primitives.v[0]) * r_inv;
pi->primitives.gradients.v[0][2] -=
wi_dx * dx[2] * (pi->primitives.v[0] - pj->primitives.v[0]) / r;
wi_dx * dx[2] * (pi->primitives.v[0] - pj->primitives.v[0]) * r_inv;
pi->primitives.gradients.v[1][0] -=
wi_dx * dx[0] * (pi->primitives.v[1] - pj->primitives.v[1]) / r;
wi_dx * dx[0] * (pi->primitives.v[1] - pj->primitives.v[1]) * r_inv;
pi->primitives.gradients.v[1][1] -=
wi_dx * dx[1] * (pi->primitives.v[1] - pj->primitives.v[1]) / r;
wi_dx * dx[1] * (pi->primitives.v[1] - pj->primitives.v[1]) * r_inv;
pi->primitives.gradients.v[1][2] -=
wi_dx * dx[2] * (pi->primitives.v[1] - pj->primitives.v[1]) / r;
wi_dx * dx[2] * (pi->primitives.v[1] - pj->primitives.v[1]) * r_inv;
pi->primitives.gradients.v[2][0] -=
wi_dx * dx[0] * (pi->primitives.v[2] - pj->primitives.v[2]) / r;
wi_dx * dx[0] * (pi->primitives.v[2] - pj->primitives.v[2]) * r_inv;
pi->primitives.gradients.v[2][1] -=
wi_dx * dx[1] * (pi->primitives.v[2] - pj->primitives.v[2]) / r;
wi_dx * dx[1] * (pi->primitives.v[2] - pj->primitives.v[2]) * r_inv;
pi->primitives.gradients.v[2][2] -=
wi_dx * dx[2] * (pi->primitives.v[2] - pj->primitives.v[2]) / r;
wi_dx * dx[2] * (pi->primitives.v[2] - pj->primitives.v[2]) * r_inv;
pi->primitives.gradients.P[0] -=
wi_dx * dx[0] * (pi->primitives.P - pj->primitives.P) / r;
wi_dx * dx[0] * (pi->primitives.P - pj->primitives.P) * r_inv;
pi->primitives.gradients.P[1] -=
wi_dx * dx[1] * (pi->primitives.P - pj->primitives.P) / r;
wi_dx * dx[1] * (pi->primitives.P - pj->primitives.P) * r_inv;
pi->primitives.gradients.P[2] -=
wi_dx * dx[2] * (pi->primitives.P - pj->primitives.P) / r;
wi_dx * dx[2] * (pi->primitives.P - pj->primitives.P) * r_inv;
}
hydro_slope_limit_cell_collect(pi, pj, r);
......@@ -435,12 +433,12 @@ hydro_gradients_nonsym_collect(float r2, const float *dx, float hi, float hj,
__attribute__((always_inline)) INLINE static void hydro_gradients_finalize(
struct part *p) {
float h, ih;
/* add kernel normalization to gradients */
h = p->h;
ih = 1.0f / h;
const float ihdim = pow_dimension(ih);
const float volume = p->geometry.volume;
const float h = p->h;
const float h_inv = 1.0f / h;
const float ihdim = pow_dimension(h_inv);
const float ihdimp1 = pow_dimension_plus_one(h_inv);
if (p->density.wcorr > const_gizmo_min_wcorr) {
p->primitives.gradients.rho[0] *= ihdim;
......@@ -462,9 +460,6 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_finalize(
p->primitives.gradients.P[2] *= ihdim;
} else {
const float ihdimp1 = pow_dimension_plus_one(ih);
float volume = p->geometry.volume;
/* finalize gradients by multiplying with volume */
p->primitives.gradients.rho[0] *= ihdimp1 * volume;
......
......@@ -64,90 +64,91 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_collect(
float r2, const float *dx, float hi, float hj, struct part *restrict pi,
struct part *restrict pj) {
float wi, wi_dx, xi, hi_inv;
float wj, wj_dx, xj, hj_inv;
float r = sqrtf(r2);
const float r_inv = 1.f / sqrtf(r2);
const float r = r2 * r_inv;
hi_inv = 1.0f / hi;
xi = r * hi_inv;
float wi, wi_dx;
const float hi_inv = 1.0f / hi;
const float xi = r * hi_inv;
kernel_deval(xi, &wi, &wi_dx);
/* very basic gradient estimate */
pi->primitives.gradients.rho[0] -=
wi_dx * dx[0] * (pi->primitives.rho - pj->primitives.rho) / r;
wi_dx * dx[0] * (pi->primitives.rho - pj->primitives.rho) * r_inv;
pi->primitives.gradients.rho[1] -=
wi_dx * dx[1] * (pi->primitives.rho - pj->primitives.rho) / r;
wi_dx * dx[1] * (pi->primitives.rho - pj->primitives.rho) * r_inv;
pi->primitives.gradients.rho[2] -=
wi_dx * dx[2] * (pi->primitives.rho - pj->primitives.rho) / r;
wi_dx * dx[2] * (pi->primitives.rho - pj->primitives.rho) * r_inv;
pi->primitives.gradients.v[0][0] -=
wi_dx * dx[0] * (pi->primitives.v[0] - pj->primitives.v[0]) / r;
wi_dx * dx[0] * (pi->primitives.v[0] - pj->primitives.v[0]) * r_inv;
pi->primitives.gradients.v[0][1] -=
wi_dx * dx[1] * (pi->primitives.v[0] - pj->primitives.v[0]) / r;
wi_dx * dx[1] * (pi->primitives.v[0] - pj->primitives.v[0]) * r_inv;
pi->primitives.gradients.v[0][2] -=
wi_dx * dx[2] * (pi->primitives.v[0] - pj->primitives.v[0]) / r;
wi_dx * dx[2] * (pi->primitives.v[0] - pj->primitives.v[0]) * r_inv;
pi->primitives.gradients.v[1][0] -=
wi_dx * dx[0] * (pi->primitives.v[1] - pj->primitives.v[1]) / r;
wi_dx * dx[0] * (pi->primitives.v[1] - pj->primitives.v[1]) * r_inv;
pi->primitives.gradients.v[1][1] -=
wi_dx * dx[1] * (pi->primitives.v[1] - pj->primitives.v[1]) / r;
wi_dx * dx[1] * (pi->primitives.v[1] - pj->primitives.v[1]) * r_inv;
pi->primitives.gradients.v[1][2] -=
wi_dx * dx[2] * (pi->primitives.v[1] - pj->primitives.v[1]) / r;
wi_dx * dx[2] * (pi->primitives.v[1] - pj->primitives.v[1]) * r_inv;
pi->primitives.gradients.v[2][0] -=
wi_dx * dx[0] * (pi->primitives.v[2] - pj->primitives.v[2]) / r;
wi_dx * dx[0] * (pi->primitives.v[2] - pj->primitives.v[2]) * r_inv;
pi->primitives.gradients.v[2][1] -=
wi_dx * dx[1] * (pi->primitives.v[2] - pj->primitives.v[2]) / r;
wi_dx * dx[1] * (pi->primitives.v[2] - pj->primitives.v[2]) * r_inv;
pi->primitives.gradients.v[2][2] -=
wi_dx * dx[2] * (pi->primitives.v[2] - pj->primitives.v[2]) / r;
wi_dx * dx[2] * (pi->primitives.v[2] - pj->primitives.v[2]) * r_inv;
pi->primitives.gradients.P[0] -=
wi_dx * dx[0] * (pi->primitives.P - pj->primitives.P) / r;
wi_dx * dx[0] * (pi->primitives.P - pj->primitives.P) * r_inv;
pi->primitives.gradients.P[1] -=
wi_dx * dx[1] * (pi->primitives.P - pj->primitives.P) / r;
wi_dx * dx[1] * (pi->primitives.P - pj->primitives.P) * r_inv;
pi->primitives.gradients.P[2] -=
wi_dx * dx[2] * (pi->primitives.P - pj->primitives.P) / r;
wi_dx * dx[2] * (pi->primitives.P - pj->primitives.P) * r_inv;
hydro_slope_limit_cell_collect(pi, pj, r);
hj_inv = 1.0f / hj;
xj = r * hj_inv;
float wj, wj_dx;
const float hj_inv = 1.0f / hj;
const float xj = r * hj_inv;
kernel_deval(xj, &wj, &wj_dx);
/* signs are the same as before, since we swap i and j twice */
pj->primitives.gradients.rho[0] -=
wj_dx * dx[0] * (pi->primitives.rho - pj->primitives.rho) / r;
wj_dx * dx[0] * (pi->primitives.rho - pj->primitives.rho) * r_inv;
pj->primitives.gradients.rho[1] -=
wj_dx * dx[1] * (pi->primitives.rho - pj->primitives.rho) / r;
wj_dx * dx[1] * (pi->primitives.rho - pj->primitives.rho) * r_inv;
pj->primitives.gradients.rho[2] -=
wj_dx * dx[2] * (pi->primitives.rho - pj->primitives.rho) / r;
wj_dx * dx[2] * (pi->primitives.rho - pj->primitives.rho) * r_inv;
pj->primitives.gradients.v[0][0] -=
wj_dx * dx[0] * (pi->primitives.v[0] - pj->primitives.v[0]) / r;
wj_dx * dx[0] * (pi->primitives.v[0] - pj->primitives.v[0]) * r_inv;
pj->primitives.gradients.v[0][1] -=
wj_dx * dx[1] * (pi->primitives.v[0] - pj->primitives.v[0]) / r;
wj_dx * dx[1] * (pi->primitives.v[0] - pj->primitives.v[0]) * r_inv;
pj->primitives.gradients.v[0][2] -=
wj_dx * dx[2] * (pi->primitives.v[0] - pj->primitives.v[0]) / r;
wj_dx * dx[2] * (pi->primitives.v[0] - pj->primitives.v[0]) * r_inv;
pj->primitives.gradients.v[1][0] -=
wj_dx * dx[0] * (pi->primitives.v[1] - pj->primitives.v[1]) / r;
wj_dx * dx[0] * (pi->primitives.v[1] - pj->primitives.v[1]) * r_inv;
pj->primitives.gradients.v[1][1] -=
wj_dx * dx[1] * (pi->primitives.v[1] - pj->primitives.v[1]) / r;
wj_dx * dx[1] * (pi->primitives.v[1] - pj->primitives.v[1]) * r_inv;
pj->primitives.gradients.v[1][2] -=
wj_dx * dx[2] * (pi->primitives.v[1] - pj->primitives.v[1]) / r;
wj_dx * dx[2] * (pi->primitives.v[1] - pj->primitives.v[1]) * r_inv;
pj->primitives.gradients.v[2][0] -=
wj_dx * dx[0] * (pi->primitives.v[2] - pj->primitives.v[2]) / r;
wj_dx * dx[0] * (pi->primitives.v[2] - pj->primitives.v[2]) * r_inv;
pj->primitives.gradients.v[2][1] -=
wj_dx * dx[1] * (pi->primitives.v[2] - pj->primitives.v[2]) / r;
wj_dx * dx[1] * (pi->primitives.v[2] - pj->primitives.v[2]) * r_inv;
pj->primitives.gradients.v[2][2] -=
wj_dx * dx[2] * (pi->primitives.v[2] - pj->primitives.v[2]) / r;
wj_dx * dx[2] * (pi->primitives.v[2] - pj->primitives.v[2]) * r_inv;
pj->primitives.gradients.P[0] -=
wj_dx * dx[0] * (pi->primitives.P - pj->primitives.P) / r;
wj_dx * dx[0] * (pi->primitives.P - pj->primitives.P) * r_inv;
pj->primitives.gradients.P[1] -=
wj_dx * dx[1] * (pi->primitives.P - pj->primitives.P) / r;
wj_dx * dx[1] * (pi->primitives.P - pj->primitives.P) * r_inv;
pj->primitives.gradients.P[2] -=
wj_dx * dx[2] * (pi->primitives.P - pj->primitives.P) / r;
wj_dx * dx[2] * (pi->primitives.P - pj->primitives.P) * r_inv;
hydro_slope_limit_cell_collect(pj, pi, r);
}
......@@ -168,48 +169,49 @@ hydro_gradients_nonsym_collect(float r2, const float *dx, float hi, float hj,
struct part *restrict pi,
struct part *restrict pj) {
float wi, wi_dx, xi, hi_inv;
float r = sqrtf(r2);
const float r_inv = 1.f / sqrtf(r2);
const float r = r2 * r_inv;
hi_inv = 1.0f / hi;
xi = r * hi_inv;
float wi, wi_dx;
const float hi_inv = 1.0f / hi;
const float xi = r * hi_inv;
kernel_deval(xi, &wi, &wi_dx);
/* very basic gradient estimate */
pi->primitives.gradients.rho[0] -=
wi_dx * dx[0] * (pi->primitives.rho - pj->primitives.rho) / r;
wi_dx * dx[0] * (pi->primitives.rho - pj->primitives.rho) * r_inv;
pi->primitives.gradients.rho[1] -=
wi_dx * dx[1] * (pi->primitives.rho - pj->primitives.rho) / r;
wi_dx * dx[1] * (pi->primitives.rho - pj->primitives.rho) * r_inv;
pi->primitives.gradients.rho[2] -=
wi_dx * dx[2] * (pi->primitives.rho - pj->primitives.rho) / r;
wi_dx * dx[2] * (pi->primitives.rho - pj->primitives.rho) * r_inv;
pi->primitives.gradients.v[0][0] -=
wi_dx * dx[0] * (pi->primitives.v[0] - pj->primitives.v[0]) / r;
wi_dx * dx[0] * (pi->primitives.v[0] - pj->primitives.v[0]) * r_inv;
pi->primitives.gradients.v[0][1] -=
wi_dx * dx[1] * (pi->primitives.v[0] - pj->primitives.v[0]) / r;
wi_dx * dx[1] * (pi->primitives.v[0] - pj->primitives.v[0]) * r_inv;
pi->primitives.gradients.v[0][2] -=
wi_dx * dx[2] * (pi->primitives.v[0] - pj->primitives.v[0]) / r;
wi_dx * dx[2] * (pi->primitives.v[0] - pj->primitives.v[0]) * r_inv;
pi->primitives.gradients.v[1][0] -=
wi_dx * dx[0] * (pi->primitives.v[1] - pj->primitives.v[1]) / r;
wi_dx * dx[0] * (pi->primitives.v[1] - pj->primitives.v[1]) * r_inv;
pi->primitives.gradients.v[1][1] -=
wi_dx * dx[1] * (pi->primitives.v[1] - pj->primitives.v[1]) / r;
wi_dx * dx[1] * (pi->primitives.v[1] - pj->primitives.v[1]) * r_inv;
pi->primitives.gradients.v[1][2] -=
wi_dx * dx[2] * (pi->primitives.v[1] - pj->primitives.v[1]) / r;
wi_dx * dx[2] * (pi->primitives.v[1] - pj->primitives.v[1]) * r_inv;
pi->primitives.gradients.v[2][0] -=
wi_dx * dx[0] * (pi->primitives.v[2] - pj->primitives.v[2]) / r;
wi_dx * dx[0] * (pi->primitives.v[2] - pj->primitives.v[2]) * r_inv;
pi->primitives.gradients.v[2][1] -=
wi_dx * dx[1] * (pi->primitives.v[2] - pj->primitives.v[2]) / r;
wi_dx * dx[1] * (pi->primitives.v[2] - pj->primitives.v[2]) * r_inv;
pi->primitives.gradients.v[2][2] -=
wi_dx * dx[2] * (pi->primitives.v[2] - pj->primitives.v[2]) / r;
wi_dx * dx[2] * (pi->primitives.v[2] - pj->primitives.v[2]) * r_inv;
pi->primitives.gradients.P[0] -=
wi_dx * dx[0] * (pi->primitives.P - pj->primitives.P) / r;
wi_dx * dx[0] * (pi->primitives.P - pj->primitives.P) * r_inv;
pi->primitives.gradients.P[1] -=
wi_dx * dx[1] * (pi->primitives.P - pj->primitives.P) / r;
wi_dx * dx[1] * (pi->primitives.P - pj->primitives.P) * r_inv;
pi->primitives.gradients.P[2] -=
wi_dx * dx[2] * (pi->primitives.P - pj->primitives.P) / r;
wi_dx * dx[2] * (pi->primitives.P - pj->primitives.P) * r_inv;
hydro_slope_limit_cell_collect(pi, pj, r);
}
......@@ -225,8 +227,7 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_finalize(
const float h = p->h;
const float ih = 1.0f / h;
const float ihdimp1 = pow_dimension_plus_one(ih);
float volume = p->geometry.volume;
const float volume = p->geometry.volume;
/* finalize gradients by multiplying with volume */
p->primitives.gradients.rho[0] *= ihdimp1 * volume;
......
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