diff --git a/src/const.h b/src/const.h
index 85d64f5f3bf97b3f6ca6df8f230e0e0990687d34..928835c5ae0ac70d641243b0d183c2c0652cdc8a 100644
--- a/src/const.h
+++ b/src/const.h
@@ -68,7 +68,9 @@
#define GIZMO_UNPHYSICAL_RESCUE
/* Show a warning message if an unphysical value was reset (only works if
GIZMO_UNPHYSICAL_RESCUE is also selected). */
+#ifdef SWIFT_DEBUG_CHECKS
#define GIZMO_UNPHYSICAL_WARNING
+#endif
/* Parameters that control how GIZMO handles pathological particle
configurations. */
diff --git a/src/hydro/Gizmo/hydro_gradients.h b/src/hydro/Gizmo/hydro_gradients.h
index 4a2f603e03321e40b744e5e787f26c849d6e2a8a..cd08e19efb76848b21d483fe5360e32f209eb65d 100644
--- a/src/hydro/Gizmo/hydro_gradients.h
+++ b/src/hydro/Gizmo/hydro_gradients.h
@@ -93,17 +93,16 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_finalize(
*/
__attribute__((always_inline)) INLINE static void hydro_gradients_predict(
struct part* restrict pi, struct part* restrict pj, float hi, float hj,
- const float* dx, float r, float* xij_i, float* Wi, float* Wj) {
+ const float* dx, float r, const float* xij_i, float* Wi, float* Wj) {
- float dWi[5], dWj[5];
- float xij_j[3];
/* perform gradient reconstruction in space and time */
/* Compute interface position (relative to pj, since we don't need the actual
* position) eqn. (8) */
const float xfac = hj / (hi + hj);
- for (int k = 0; k < 3; k++) xij_j[k] = xfac * dx[k];
+ const float xij_j[3] = {xfac * dx[0], xfac * dx[1], xfac * dx[2]};
+ float dWi[5];
dWi[0] = pi->primitives.gradients.rho[0] * xij_i[0] +
pi->primitives.gradients.rho[1] * xij_i[1] +
pi->primitives.gradients.rho[2] * xij_i[2];
@@ -120,6 +119,7 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_predict(
pi->primitives.gradients.P[1] * xij_i[1] +
pi->primitives.gradients.P[2] * xij_i[2];
+ float dWj[5];
dWj[0] = pj->primitives.gradients.rho[0] * xij_j[0] +
pj->primitives.gradients.rho[1] * xij_j[1] +
pj->primitives.gradients.rho[2] * xij_j[2];
diff --git a/src/hydro/Gizmo/hydro_iact.h b/src/hydro/Gizmo/hydro_iact.h
index 7017056f2e5a784b2eed6273f13bc7d1752617fc..2274d0fa75d92ffa814c321b019803f83421bbc4 100644
--- a/src/hydro/Gizmo/hydro_iact.h
+++ b/src/hydro/Gizmo/hydro_iact.h
@@ -51,15 +51,14 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
float r2, const float *dx, float hi, float hj, struct part *restrict pi,
struct part *restrict pj, float a, float H) {
- float r = sqrtf(r2);
- float xi, xj;
- float h_inv;
float wi, wj, wi_dx, wj_dx;
- int k, l;
+
+ /* Get r and h inverse. */
+ const float r = sqrtf(r2);
/* Compute density of pi. */
- h_inv = 1.f / hi;
- xi = r * h_inv;
+ const float hi_inv = 1.f / hi;
+ const float xi = r * hi_inv;
kernel_deval(xi, &wi, &wi_dx);
pi->density.wcount += wi;
@@ -67,16 +66,17 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
/* these are eqns. (1) and (2) in the summary */
pi->geometry.volume += wi;
- for (k = 0; k < 3; k++)
- for (l = 0; l < 3; l++) pi->geometry.matrix_E[k][l] += dx[k] * dx[l] * wi;
+ for (int k = 0; k < 3; k++)
+ for (int l = 0; l < 3; l++)
+ pi->geometry.matrix_E[k][l] += dx[k] * dx[l] * wi;
pi->geometry.centroid[0] -= dx[0] * wi;
pi->geometry.centroid[1] -= dx[1] * wi;
pi->geometry.centroid[2] -= dx[2] * wi;
/* Compute density of pj. */
- h_inv = 1.f / hj;
- xj = r * h_inv;
+ const float hj_inv = 1.f / hj;
+ const float xj = r * hj_inv;
kernel_deval(xj, &wj, &wj_dx);
pj->density.wcount += wj;
@@ -84,8 +84,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
/* these are eqns. (1) and (2) in the summary */
pj->geometry.volume += wj;
- for (k = 0; k < 3; k++)
- for (l = 0; l < 3; l++) pj->geometry.matrix_E[k][l] += dx[k] * dx[l] * wj;
+ for (int k = 0; k < 3; k++)
+ for (int l = 0; l < 3; l++)
+ pj->geometry.matrix_E[k][l] += dx[k] * dx[l] * wj;
pj->geometry.centroid[0] += dx[0] * wj;
pj->geometry.centroid[1] += dx[1] * wj;
@@ -117,17 +118,13 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
float r2, const float *dx, float hi, float hj, struct part *restrict pi,
const struct part *restrict pj, float a, float H) {
- float r;
- float xi;
- float h_inv;
float wi, wi_dx;
- int k, l;
- /* Get r and r inverse. */
- r = sqrtf(r2);
+ /* Get r and h inverse. */
+ const float r = sqrtf(r2);
- h_inv = 1.f / hi;
- xi = r * h_inv;
+ const float hi_inv = 1.f / hi;
+ const float xi = r * hi_inv;
kernel_deval(xi, &wi, &wi_dx);
pi->density.wcount += wi;
@@ -135,8 +132,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
/* these are eqns. (1) and (2) in the summary */
pi->geometry.volume += wi;
- for (k = 0; k < 3; k++)
- for (l = 0; l < 3; l++) pi->geometry.matrix_E[k][l] += dx[k] * dx[l] * wi;
+ for (int k = 0; k < 3; k++)
+ for (int l = 0; l < 3; l++)
+ pi->geometry.matrix_E[k][l] += dx[k] * dx[l] * wi;
pi->geometry.centroid[0] -= dx[0] * wi;
pi->geometry.centroid[1] -= dx[1] * wi;
@@ -222,34 +220,24 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
float r2, const float *dx, float hi, float hj, struct part *restrict pi,
struct part *restrict pj, int mode, float a, float H) {
- float r = sqrtf(r2);
- float xi, xj;
- float hi_inv, hi_inv_dim;
- float hj_inv, hj_inv_dim;
- float wi, wj, wi_dx, wj_dx;
- int k, l;
- float A[3];
- float Anorm;
- float Bi[3][3];
- float Bj[3][3];
- float Vi, Vj;
- float xij_i[3], xfac, xijdotdx;
- float vmax, dvdotdx;
- float vi[3], vj[3], vij[3];
- float Wi[5], Wj[5];
- float n_unit[3];
+ const float r_inv = 1.f / sqrtf(r2);
+ const float r = r2 * r_inv;
/* Initialize local variables */
- for (k = 0; k < 3; k++) {
- for (l = 0; l < 3; l++) {
+ float Bi[3][3];
+ float Bj[3][3];
+ float vi[3], vj[3];
+ 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];
}
vi[k] = pi->v[k]; /* particle velocities */
vj[k] = pj->v[k];
}
- Vi = pi->geometry.volume;
- Vj = pj->geometry.volume;
+ const float Vi = pi->geometry.volume;
+ const float Vj = pj->geometry.volume;
+ float Wi[5], Wj[5];
Wi[0] = pi->primitives.rho;
Wi[1] = pi->primitives.v[0];
Wi[2] = pi->primitives.v[1];
@@ -262,56 +250,63 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
Wj[4] = pj->primitives.P;
/* calculate the maximal signal velocity */
+ float vmax;
if (Wi[0] > 0.0f && Wj[0] > 0.0f) {
- vmax =
- sqrtf(hydro_gamma * Wi[4] / Wi[0]) + sqrtf(hydro_gamma * Wj[4] / Wj[0]);
- } else {
+ const float ci = gas_soundspeed_from_pressure(Wi[0], Wi[4]);
+ const float cj = gas_soundspeed_from_pressure(Wj[0], Wj[4]);
+ vmax = ci + cj;
+ } else
vmax = 0.0f;
- }
- dvdotdx = (Wi[1] - Wj[1]) * dx[0] + (Wi[2] - Wj[2]) * dx[1] +
- (Wi[3] - Wj[3]) * dx[2];
- dvdotdx = min(dvdotdx, (vi[0] - vj[0]) * dx[0] + (vi[1] - vj[1]) * dx[1] +
- (vi[2] - vj[2]) * dx[2]);
- if (dvdotdx < 0.) {
- /* the magical factor 3 also appears in Gadget2 */
- vmax -= 3.f * dvdotdx / r;
- }
+
+ float dvdr = (pi->v[0] - pj->v[0]) * dx[0] + (pi->v[1] - pj->v[1]) * dx[1] +
+ (pi->v[2] - pj->v[2]) * dx[2];
+
+ /* Velocity on the axis linking the particles */
+ float dvdotdx = (Wi[1] - Wj[1]) * dx[0] + (Wi[2] - Wj[2]) * dx[1] +
+ (Wi[3] - Wj[3]) * dx[2];
+ dvdotdx = min(dvdotdx, dvdr);
+
+ /* We only care about this velocity for particles moving towards each others
+ */
+ dvdotdx = min(dvdotdx, 0.f);
+
+ /* Get the signal velocity */
+ /* the magical factor 3 also appears in Gadget2 */
+ vmax -= 3.f * dvdotdx * r_inv;
+
+ /* Store the signal velocity */
pi->timestepvars.vmax = max(pi->timestepvars.vmax, vmax);
- if (mode == 1) {
- pj->timestepvars.vmax = max(pj->timestepvars.vmax, vmax);
- }
+ if (mode == 1) pj->timestepvars.vmax = max(pj->timestepvars.vmax, vmax);
/* Compute kernel of pi. */
- hi_inv = 1.f / hi;
- hi_inv_dim = pow_dimension(hi_inv);
- xi = r * hi_inv;
+ float wi, wi_dx;
+ const float hi_inv = 1.f / hi;
+ const float hi_inv_dim = pow_dimension(hi_inv);
+ const float xi = r * hi_inv;
kernel_deval(xi, &wi, &wi_dx);
/* Compute kernel of pj. */
- hj_inv = 1.f / hj;
- hj_inv_dim = pow_dimension(hj_inv);
- xj = r * hj_inv;
+ float wj, wj_dx;
+ const float hj_inv = 1.f / hj;
+ const float hj_inv_dim = pow_dimension(hj_inv);
+ const float xj = r * hj_inv;
kernel_deval(xj, &wj, &wj_dx);
/* Compute h_dt. We are going to use an SPH-like estimate of div_v for that */
- float dvdr = (pi->v[0] - pj->v[0]) * dx[0] + (pi->v[1] - pj->v[1]) * dx[1] +
- (pi->v[2] - pj->v[2]) * dx[2];
- float ri = 1.0f / r;
- float hidp1 = pow_dimension_plus_one(hi_inv);
- float hjdp1 = pow_dimension_plus_one(hj_inv);
- float wi_dr = hidp1 * wi_dx;
- float wj_dr = hjdp1 * wj_dx;
- dvdr *= ri;
- if (pj->primitives.rho > 0.) {
+ const float hidp1 = pow_dimension_plus_one(hi_inv);
+ const float hjdp1 = pow_dimension_plus_one(hj_inv);
+ const float wi_dr = hidp1 * wi_dx;
+ const float wj_dr = hjdp1 * wj_dx;
+ dvdr *= r_inv;
+ if (pj->primitives.rho > 0.)
pi->force.h_dt -= pj->conserved.mass * dvdr / pj->primitives.rho * wi_dr;
- }
- if (mode == 1 && pi->primitives.rho > 0.) {
+ if (mode == 1 && pi->primitives.rho > 0.)
pj->force.h_dt -= pi->conserved.mass * dvdr / pi->primitives.rho * wj_dr;
- }
- /* Compute area */
+ /* Compute (square of) area */
/* eqn. (7) */
- Anorm = 0.0f;
+ float Anorm2 = 0.0f;
+ float A[3];
if (pi->density.wcorr > const_gizmo_min_wcorr &&
pj->density.wcorr > const_gizmo_min_wcorr) {
/* in principle, we use Vi and Vj as weights for the left and right
@@ -327,30 +322,31 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
Xj = Xi;
}
#endif
- for (k = 0; k < 3; k++) {
+ for (int k = 0; k < 3; k++) {
/* we add a minus sign since dx is pi->x - pj->x */
A[k] = -Xi * (Bi[k][0] * dx[0] + Bi[k][1] * dx[1] + Bi[k][2] * dx[2]) *
wi * hi_inv_dim -
Xj * (Bj[k][0] * dx[0] + Bj[k][1] * dx[1] + Bj[k][2] * dx[2]) *
wj * hj_inv_dim;
- Anorm += A[k] * A[k];
+ Anorm2 += A[k] * A[k];
}
} else {
/* ill condition gradient matrix: revert to SPH face area */
- Anorm = -(hidp1 * Vi * Vi * wi_dx + hjdp1 * Vj * Vj * wj_dx) * ri;
+ const float Anorm =
+ -(hidp1 * Vi * Vi * wi_dx + hjdp1 * Vj * Vj * wj_dx) * r_inv;
A[0] = -Anorm * dx[0];
A[1] = -Anorm * dx[1];
A[2] = -Anorm * dx[2];
- Anorm *= Anorm * r2;
+ Anorm2 = Anorm * Anorm * r2;
}
- if (Anorm == 0.f) {
- /* if the interface has no area, nothing happens and we return */
- /* continuing results in dividing by zero and NaN's... */
- return;
- }
+ /* if the interface has no area, nothing happens and we return */
+ /* continuing results in dividing by zero and NaN's... */
+ if (Anorm2 == 0.f) return;
- Anorm = sqrtf(Anorm);
+ /* Compute the area */
+ const float Anorm_inv = 1. / sqrtf(Anorm2);
+ const float Anorm = Anorm2 * Anorm_inv;
#ifdef SWIFT_DEBUG_CHECKS
/* For stability reasons, we do require A and dx to have opposite
@@ -358,29 +354,32 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
always points from particle i to particle j, as it would in a real
moving-mesh code). If not, our scheme is no longer upwind and hence can
become unstable. */
- float dA_dot_dx = A[0] * dx[0] + A[1] * dx[1] + A[2] * dx[2];
+ const float dA_dot_dx = A[0] * dx[0] + A[1] * dx[1] + A[2] * dx[2];
/* In GIZMO, Phil Hopkins reverts to an SPH integration scheme if this
happens. We curently just ignore this case and display a message. */
const float rdim = pow_dimension(r);
- if (dA_dot_dx > 1.e-6 * rdim) {
+ if (dA_dot_dx > 1.e-6f * rdim) {
message("Ill conditioned gradient matrix (%g %g %g %g %g)!", dA_dot_dx,
Anorm, Vi, Vj, r);
}
#endif
/* compute the normal vector of the interface */
- for (k = 0; k < 3; k++) n_unit[k] = A[k] / Anorm;
+ const float n_unit[3] = {A[0] * Anorm_inv, A[1] * Anorm_inv,
+ A[2] * Anorm_inv};
/* Compute interface position (relative to pi, since we don't need the actual
- * position) */
- /* eqn. (8) */
- xfac = hi / (hi + hj);
- for (k = 0; k < 3; k++) xij_i[k] = -xfac * dx[k];
+ * position) eqn. (8) */
+ const float xfac = -hi / (hi + hj);
+ const float xij_i[3] = {xfac * dx[0], xfac * dx[1], xfac * dx[2]};
/* Compute interface velocity */
/* eqn. (9) */
- xijdotdx = xij_i[0] * dx[0] + xij_i[1] * dx[1] + xij_i[2] * dx[2];
- for (k = 0; k < 3; k++) vij[k] = vi[k] + (vi[k] - vj[k]) * xijdotdx / r2;
+ float xijdotdx = xij_i[0] * dx[0] + xij_i[1] * dx[1] + xij_i[2] * dx[2];
+ xijdotdx *= r_inv * r_inv;
+ const float vij[3] = {vi[0] + (vi[0] - vj[0]) * xijdotdx,
+ vi[1] + (vi[1] - vj[1]) * xijdotdx,
+ vi[2] + (vi[2] - vj[2]) * xijdotdx};
/* complete calculation of position of interface */
/* NOTE: dx is not necessarily just pi->x - pj->x but can also contain
@@ -388,8 +387,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
we have to use xij w.r.t. the actual particle.
=> we need a separate xij for pi and pj... */
/* tldr: we do not need the code below, but we do need the same code as above
- but then
- with i and j swapped */
+ but then with i and j swapped */
// for ( k = 0 ; k < 3 ; k++ )
// xij[k] += pi->x[k];
@@ -418,10 +416,10 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
totflux[4] *= Anorm;
/* Store mass flux */
- float mflux = totflux[0];
- pi->gravity.mflux[0] += mflux * dx[0];
- pi->gravity.mflux[1] += mflux * dx[1];
- pi->gravity.mflux[2] += mflux * dx[2];
+ const float mflux_i = totflux[0];
+ pi->gravity.mflux[0] += mflux_i * dx[0];
+ pi->gravity.mflux[1] += mflux_i * dx[1];
+ pi->gravity.mflux[2] += mflux_i * dx[2];
/* Update conserved variables */
/* eqn. (16) */
@@ -432,13 +430,13 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
pi->conserved.flux.energy -= totflux[4];
#ifndef GIZMO_TOTAL_ENERGY
- float ekin = 0.5f * (pi->primitives.v[0] * pi->primitives.v[0] +
- pi->primitives.v[1] * pi->primitives.v[1] +
- pi->primitives.v[2] * pi->primitives.v[2]);
+ const float ekin_i = 0.5f * (pi->primitives.v[0] * pi->primitives.v[0] +
+ pi->primitives.v[1] * pi->primitives.v[1] +
+ pi->primitives.v[2] * pi->primitives.v[2]);
pi->conserved.flux.energy += totflux[1] * pi->primitives.v[0];
pi->conserved.flux.energy += totflux[2] * pi->primitives.v[1];
pi->conserved.flux.energy += totflux[3] * pi->primitives.v[2];
- pi->conserved.flux.energy -= totflux[0] * ekin;
+ pi->conserved.flux.energy -= totflux[0] * ekin_i;
#endif
/* Note that this used to be much more complicated in early implementations of
@@ -448,10 +446,10 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
* for the flux over the entire time step. */
if (mode == 1) {
/* Store mass flux */
- mflux = totflux[0];
- pj->gravity.mflux[0] -= mflux * dx[0];
- pj->gravity.mflux[1] -= mflux * dx[1];
- pj->gravity.mflux[2] -= mflux * dx[2];
+ const float mflux_j = totflux[0];
+ pj->gravity.mflux[0] -= mflux_j * dx[0];
+ pj->gravity.mflux[1] -= mflux_j * dx[1];
+ pj->gravity.mflux[2] -= mflux_j * dx[2];
pj->conserved.flux.mass += totflux[0];
pj->conserved.flux.momentum[0] += totflux[1];
@@ -460,13 +458,13 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
pj->conserved.flux.energy += totflux[4];
#ifndef GIZMO_TOTAL_ENERGY
- ekin = 0.5f * (pj->primitives.v[0] * pj->primitives.v[0] +
- pj->primitives.v[1] * pj->primitives.v[1] +
- pj->primitives.v[2] * pj->primitives.v[2]);
+ const float ekin_j = 0.5f * (pj->primitives.v[0] * pj->primitives.v[0] +
+ pj->primitives.v[1] * pj->primitives.v[1] +
+ pj->primitives.v[2] * pj->primitives.v[2]);
pj->conserved.flux.energy -= totflux[1] * pj->primitives.v[0];
pj->conserved.flux.energy -= totflux[2] * pj->primitives.v[1];
pj->conserved.flux.energy -= totflux[3] * pj->primitives.v[2];
- pj->conserved.flux.energy += totflux[0] * ekin;
+ pj->conserved.flux.energy += totflux[0] * ekin_j;
#endif
}
}
diff --git a/src/hydro/Gizmo/hydro_slope_limiters_face.h b/src/hydro/Gizmo/hydro_slope_limiters_face.h
index 5a955217ee9e5b9efe876b1be2dc78d7f0262fb1..b5ebcc5eab8ad515c791686e1c72b3eafe1db2be 100644
--- a/src/hydro/Gizmo/hydro_slope_limiters_face.h
+++ b/src/hydro/Gizmo/hydro_slope_limiters_face.h
@@ -88,7 +88,7 @@ hydro_slope_limit_face_quantity(float phi_i, float phi_j, float phi_mid0,
* @param r Distance between particle i and particle j.
*/
__attribute__((always_inline)) INLINE static void hydro_slope_limit_face(
- float *Wi, float *Wj, float *dWi, float *dWj, float *xij_i, float *xij_j,
+ float *Wi, float *Wj, float *dWi, float *dWj, const float *xij_i, const float *xij_j,
float r) {
const float xij_i_norm =
diff --git a/src/riemann/riemann_exact.h b/src/riemann/riemann_exact.h
index 4a7b40ae78c4c3f1a880539367172894fe298995..9081eb31060ee20ef0dada0621b9bc71fad4db08 100644
--- a/src/riemann/riemann_exact.h
+++ b/src/riemann/riemann_exact.h
@@ -47,14 +47,14 @@
* @param W The left or right state vector
* @param a The left or right sound speed
*/
-__attribute__((always_inline)) INLINE static float riemann_fb(float p, float* W,
+__attribute__((always_inline)) INLINE static float riemann_fb(float p,
+ const float* W,
float a) {
- float fval = 0.;
- float A, B;
+ float fval;
if (p > W[4]) {
- A = hydro_two_over_gamma_plus_one / W[0];
- B = hydro_gamma_minus_one_over_gamma_plus_one * W[4];
+ const float A = hydro_two_over_gamma_plus_one / W[0];
+ const float B = hydro_gamma_minus_one_over_gamma_plus_one * W[4];
fval = (p - W[4]) * sqrtf(A / (p + B));
} else {
fval = hydro_two_over_gamma_minus_one * a *
@@ -75,7 +75,8 @@ __attribute__((always_inline)) INLINE static float riemann_fb(float p, float* W,
* @param aR The right sound speed
*/
__attribute__((always_inline)) INLINE static float riemann_f(
- float p, float* WL, float* WR, float vL, float vR, float aL, float aR) {
+ float p, const float* WL, const float* WR, float vL, float vR, float aL,
+ float aR) {
return riemann_fb(p, WL, aL) + riemann_fb(p, WR, aR) + (vR - vL);
}
@@ -87,15 +88,13 @@ __attribute__((always_inline)) INLINE static float riemann_f(
* @param W The left or right state vector
* @param a The left or right sound speed
*/
-__attribute__((always_inline)) INLINE static float riemann_fprimeb(float p,
- float* W,
- float a) {
+__attribute__((always_inline)) INLINE static float riemann_fprimeb(
+ float p, const float* W, float a) {
- float fval = 0.;
- float A, B;
+ float fval;
if (p > W[4]) {
- A = hydro_two_over_gamma_plus_one / W[0];
- B = hydro_gamma_minus_one_over_gamma_plus_one * W[4];
+ const float A = hydro_two_over_gamma_plus_one / W[0];
+ const float B = hydro_gamma_minus_one_over_gamma_plus_one * W[4];
fval = (1.0f - 0.5f * (p - W[4]) / (B + p)) * sqrtf(A / (p + B));
} else {
fval = 1.0f / W[0] / a * pow_minus_gamma_plus_one_over_two_gamma(p / W[4]);
@@ -113,7 +112,7 @@ __attribute__((always_inline)) INLINE static float riemann_fprimeb(float p,
* @param aR The right sound speed
*/
__attribute__((always_inline)) INLINE static float riemann_fprime(
- float p, float* WL, float* WR, float aL, float aR) {
+ float p, const float* WL, const float* WR, float aL, float aR) {
return riemann_fprimeb(p, WL, aL) + riemann_fprimeb(p, WR, aR);
}
@@ -125,11 +124,10 @@ __attribute__((always_inline)) INLINE static float riemann_fprime(
* @param W The left or right state vector
*/
__attribute__((always_inline)) INLINE static float riemann_gb(float p,
- float* W) {
+ const float* W) {
- float A, B;
- A = hydro_two_over_gamma_plus_one / W[0];
- B = hydro_gamma_minus_one_over_gamma_plus_one * W[4];
+ const float A = hydro_two_over_gamma_plus_one / W[0];
+ const float B = hydro_gamma_minus_one_over_gamma_plus_one * W[4];
return sqrtf(A / (p + B));
}
@@ -147,7 +145,7 @@ __attribute__((always_inline)) INLINE static float riemann_gb(float p,
* @param aR The right sound speed
*/
__attribute__((always_inline)) INLINE static float riemann_guess_p(
- float* WL, float* WR, float vL, float vR, float aL, float aR) {
+ const float* WL, const float* WR, float vL, float vR, float aL, float aR) {
float pguess, pmin, pmax, qmax;
float ppv;
@@ -199,8 +197,8 @@ __attribute__((always_inline)) INLINE static float riemann_guess_p(
*/
__attribute__((always_inline)) INLINE static float riemann_solve_brent(
float lower_limit, float upper_limit, float lowf, float upf,
- float error_tol, float* WL, float* WR, float vL, float vR, float aL,
- float aR) {
+ float error_tol, const float* WL, const float* WR, float vL, float vR,
+ float aL, float aR) {
float a, b, c, d, s;
float fa, fb, fc, fs;
@@ -306,7 +304,7 @@ __attribute__((always_inline)) INLINE static float riemann_solve_brent(
* @param n_unit Normal vector of the interface
*/
__attribute__((always_inline)) INLINE static void riemann_solver_solve(
- float* WL, float* WR, float* Whalf, float* n_unit) {
+ const float* WL, const float* WR, float* Whalf, const float* n_unit) {
/* velocity of the left and right state in a frame aligned with n_unit */
float vL, vR, vhalf;
@@ -510,7 +508,8 @@ __attribute__((always_inline)) INLINE static void riemann_solver_solve(
}
__attribute__((always_inline)) INLINE static void riemann_solve_for_flux(
- float* Wi, float* Wj, float* n_unit, float* vij, float* totflux) {
+ const float* Wi, const float* Wj, const float* n_unit, const float* vij,
+ float* totflux) {
float Whalf[5];
float flux[5][3];
diff --git a/src/riemann/riemann_hllc.h b/src/riemann/riemann_hllc.h
index 1cdc4e809cf4cfc6865a1e0cbb75f709056f7c64..2578db44aa0fabb9f2d238abea497a2f0a6be71f 100644
--- a/src/riemann/riemann_hllc.h
+++ b/src/riemann/riemann_hllc.h
@@ -38,7 +38,8 @@
#endif
__attribute__((always_inline)) INLINE static void riemann_solve_for_flux(
- float *WL, float *WR, float *n, float *vij, float *totflux) {
+ const float *WL, const float *WR, const float *n, const float *vij,
+ float *totflux) {
/* Handle pure vacuum */
if (!WL[0] && !WR[0]) {
diff --git a/src/riemann/riemann_trrs.h b/src/riemann/riemann_trrs.h
index 44fbe6ffb8fb1f557045b6e466c87e7d8b591aab..e5605ab7d09700b8a79792cb552d6ef52569c7b5 100644
--- a/src/riemann/riemann_trrs.h
+++ b/src/riemann/riemann_trrs.h
@@ -52,7 +52,7 @@
* @param n_unit Normal vector of the interface
*/
__attribute__((always_inline)) INLINE static void riemann_solver_solve(
- float* WL, float* WR, float* Whalf, float* n_unit) {
+ const float* WL, const float* WR, float* Whalf, const float* n_unit) {
float aL, aR;
float PLR;
float vL, vR;
@@ -160,7 +160,8 @@ __attribute__((always_inline)) INLINE static void riemann_solver_solve(
}
__attribute__((always_inline)) INLINE static void riemann_solve_for_flux(
- float* Wi, float* Wj, float* n_unit, float* vij, float* totflux) {
+ const float* Wi, const float* Wj, const float* n_unit, const float* vij,
+ float* totflux) {
float Whalf[5];
float flux[5][3];
diff --git a/src/riemann/riemann_vacuum.h b/src/riemann/riemann_vacuum.h
index ace1b672e42e797cd2125bc27d8d35232ecd45dc..c1144a2557889df20ece7dea6918af804ee03232 100644
--- a/src/riemann/riemann_vacuum.h
+++ b/src/riemann/riemann_vacuum.h
@@ -24,7 +24,7 @@
* @brief Check if the given input states are vacuum or will generate vacuum
*/
__attribute__((always_inline)) INLINE static int riemann_is_vacuum(
- float* WL, float* WR, float vL, float vR, float aL, float aR) {
+ const float* WL, const float* WR, float vL, float vR, float aL, float aR) {
/* vacuum */
if (!WL[0] || !WR[0]) return 1;
@@ -53,8 +53,8 @@ __attribute__((always_inline)) INLINE static int riemann_is_vacuum(
* @param n_unit Normal vector of the interface
*/
__attribute__((always_inline)) INLINE static void riemann_solve_vacuum(
- float* WL, float* WR, float vL, float vR, float aL, float aR, float* Whalf,
- float* n_unit) {
+ const float* WL, const float* WR, float vL, float vR, float aL, float aR,
+ float* Whalf, const float* n_unit) {
float SL, SR;
float vhalf;
@@ -202,8 +202,8 @@ __attribute__((always_inline)) INLINE static void riemann_solve_vacuum(
* @brief Solve the vacuum Riemann problem and return the fluxes
*/
__attribute__((always_inline)) INLINE static void riemann_solve_vacuum_flux(
- float* WL, float* WR, float vL, float vR, float aL, float aR, float* n_unit,
- float* vij, float* totflux) {
+ const float* WL, const float* WR, float vL, float vR, float aL, float aR,
+ const float* n_unit, const float* vij, float* totflux) {
float Whalf[5];
float flux[5][3];