diff --git a/src/Makefile.am b/src/Makefile.am
index a7bc4774ea3fb162a2b3af1ea825607659e34c6a..1bd80b669396d4b4597f6149acafe6b1a4165f43 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -65,7 +65,7 @@ AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h inline.h kernel_hydro.h kernel_gravity.h \
kernel_long_gravity.h vector.h cache.h runner_doiact.h runner_doiact_vec.h runner_doiact_grav.h runner_doiact_fft.h \
runner_doiact_nosort.h units.h intrinsics.h minmax.h kick.h timestep.h drift.h adiabatic_index.h io_properties.h \
- dimension.h part_type.h periodic.h memswap.h dump.h logger.h \
+ dimension.h part_type.h periodic.h memswap.h dump.h logger.h sign.h \
gravity.h gravity_io.h gravity_cache.h \
gravity/Default/gravity.h gravity/Default/gravity_iact.h gravity/Default/gravity_io.h \
gravity/Default/gravity_debug.h gravity/Default/gravity_part.h \
@@ -110,6 +110,7 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h
hydro/Shadowswift/voronoi_cell.h \
riemann/riemann_hllc.h riemann/riemann_trrs.h \
riemann/riemann_exact.h riemann/riemann_vacuum.h \
+ riemann/riemann_checks.h \
stars.h stars_io.h \
stars/Default/star.h stars/Default/star_iact.h stars/Default/star_io.h \
stars/Default/star_debug.h stars/Default/star_part.h \
diff --git a/src/adiabatic_index.h b/src/adiabatic_index.h
index e7798f7be6e28536191475f3c940c6f84a92bdad..6a157a41c54ee8fa912d263b84578a7dcd0d4cc0 100644
--- a/src/adiabatic_index.h
+++ b/src/adiabatic_index.h
@@ -42,6 +42,7 @@
#define hydro_gamma 1.66666666666666667f
#define hydro_gamma_minus_one 0.66666666666666667f
+#define hydro_gamma_plus_one 2.66666666666666667f
#define hydro_one_over_gamma_minus_one 1.5f
#define hydro_gamma_plus_one_over_two_gamma 0.8f
#define hydro_gamma_minus_one_over_two_gamma 0.2f
@@ -56,6 +57,7 @@
#define hydro_gamma 1.4f
#define hydro_gamma_minus_one 0.4f
+#define hydro_gamma_plus_one 2.4f
#define hydro_one_over_gamma_minus_one 2.5f
#define hydro_gamma_plus_one_over_two_gamma 0.857142857f
#define hydro_gamma_minus_one_over_two_gamma 0.142857143f
@@ -70,6 +72,7 @@
#define hydro_gamma 1.33333333333333333f
#define hydro_gamma_minus_one 0.33333333333333333f
+#define hydro_gamma_plus_one 2.33333333333333333f
#define hydro_one_over_gamma_minus_one 3.f
#define hydro_gamma_plus_one_over_two_gamma 0.875f
#define hydro_gamma_minus_one_over_two_gamma 0.125f
@@ -84,6 +87,7 @@
#define hydro_gamma 2.f
#define hydro_gamma_minus_one 1.f
+#define hydro_gamma_plus_one 3.f
#define hydro_one_over_gamma_minus_one 1.f
#define hydro_gamma_plus_one_over_two_gamma 0.75f
#define hydro_gamma_minus_one_over_two_gamma 0.25f
diff --git a/src/cell.c b/src/cell.c
index ea036c97081e625c4e4ca2bb1d06cd9ea8251bd3..a1f57edcca3ff1fb42a93cdf6e222d48a00797ac 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -2405,6 +2405,15 @@ void cell_drift_part(struct cell *c, const struct engine *e, int force) {
drift_part(p, xp, dt_drift, dt_kick_hydro, dt_kick_grav, dt_therm,
ti_old_part, ti_current);
+#ifdef SWIFT_DEBUG_CHECKS
+ /* Make sure the particle does not drift by more than a box length. */
+ if (fabsf(xp->v_full[0] * dt_drift) > e->s->dim[0] ||
+ fabsf(xp->v_full[1] * dt_drift) > e->s->dim[1] ||
+ fabsf(xp->v_full[2] * dt_drift) > e->s->dim[2]) {
+ error("Particle drifts by more than a box length!");
+ }
+#endif
+
/* Limit h to within the allowed range */
p->h = min(p->h, hydro_h_max);
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.h b/src/hydro/Gizmo/hydro.h
index bc90a6790fa8123b19badf63761a6dd1378197f6..079e66669eb4e3b574742b09ca6b7e80d5386fb7 100644
--- a/src/hydro/Gizmo/hydro.h
+++ b/src/hydro/Gizmo/hydro.h
@@ -362,8 +362,9 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
#endif
/* sanity checks */
- gizmo_check_physical_quantity("density", p->primitives.rho);
- gizmo_check_physical_quantity("pressure", p->primitives.P);
+ gizmo_check_physical_quantities("density", "pressure", p->primitives.rho,
+ p->primitives.v[0], p->primitives.v[1],
+ p->primitives.v[2], p->primitives.P);
#ifdef GIZMO_LLOYD_ITERATION
/* overwrite primitive variables to make sure they still have safe values */
@@ -563,7 +564,6 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
p->primitives.v[2] +=
p->conserved.flux.momentum[2] * dt_drift / p->conserved.mass;
-// MATTHIEU: Bert is this correct?
#if !defined(EOS_ISOTHERMAL_GAS)
const float u = p->conserved.energy + p->conserved.flux.energy * dt_therm;
p->primitives.P =
@@ -576,6 +576,10 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
error("Zero or negative smoothing length (%g)!", p->h);
}
#endif
+
+ gizmo_check_physical_quantities("density", "pressure", p->primitives.rho,
+ p->primitives.v[0], p->primitives.v[1],
+ p->primitives.v[2], p->primitives.P);
}
/**
@@ -632,13 +636,14 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
/* Apply the minimal energy limit */
const float min_energy =
hydro_props->minimal_internal_energy * cosmo->a_factor_internal_energy;
- if (p->conserved.energy < min_energy) {
- p->conserved.energy = min_energy;
+ if (p->conserved.energy < min_energy * p->conserved.mass) {
+ p->conserved.energy = min_energy * p->conserved.mass;
p->conserved.flux.energy = 0.f;
}
- gizmo_check_physical_quantity("mass", p->conserved.mass);
- gizmo_check_physical_quantity("energy", p->conserved.energy);
+ gizmo_check_physical_quantities(
+ "mass", "energy", p->conserved.mass, p->conserved.momentum[0],
+ p->conserved.momentum[1], p->conserved.momentum[2], p->conserved.energy);
#ifdef SWIFT_DEBUG_CHECKS
/* Note that this check will only have effect if no GIZMO_UNPHYSICAL option
@@ -675,6 +680,10 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
p->conserved.momentum[0] += dt * p->conserved.mass * a_grav[0];
p->conserved.momentum[1] += dt * p->conserved.mass * a_grav[1];
p->conserved.momentum[2] += dt * p->conserved.mass * a_grav[2];
+
+ p->conserved.energy += dt * (p->gravity.mflux[0] * a_grav[0] +
+ p->gravity.mflux[1] * a_grav[1] +
+ p->gravity.mflux[2] * a_grav[2]);
}
hydro_velocities_set(p, xp);
diff --git a/src/hydro/Gizmo/hydro_gradients.h b/src/hydro/Gizmo/hydro_gradients.h
index cd53c1dfce670abe58e2aa108f28b1146b3f1013..c5d95e4dab4f574dd20e8355aacdd176530ec63d 100644
--- a/src/hydro/Gizmo/hydro_gradients.h
+++ b/src/hydro/Gizmo/hydro_gradients.h
@@ -93,21 +93,15 @@ __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) {
-
- float dWi[5], dWj[5];
- float xij_j[3];
- int k;
- float xfac;
+ const float* dx, float r, const float* xij_i, float* Wi, float* Wj) {
/* 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);
+ 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];
@@ -124,6 +118,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];
@@ -140,6 +135,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,10 +150,10 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_predict(
Wj[3] += dWj[3];
Wj[4] += dWj[4];
- gizmo_check_physical_quantity("density", Wi[0]);
- gizmo_check_physical_quantity("pressure", Wi[4]);
- gizmo_check_physical_quantity("density", Wj[0]);
- gizmo_check_physical_quantity("pressure", Wj[4]);
+ gizmo_check_physical_quantities("density", "pressure", Wi[0], Wi[1], Wi[2],
+ Wi[3], Wi[4]);
+ gizmo_check_physical_quantities("density", "pressure", Wj[0], Wj[1], Wj[2],
+ Wj[3], Wj[4]);
}
#endif /* SWIFT_HYDRO_GIZMO_GRADIENTS_H */
diff --git a/src/hydro/Gizmo/hydro_gradients_gizmo.h b/src/hydro/Gizmo/hydro_gradients_gizmo.h
index 0140279745c8a2520479175ecf6e5e24d3fba1d6..9f23f82e870cffb82a73964707bd28fa68ed00d7 100644
--- a/src/hydro/Gizmo/hydro_gradients_gizmo.h
+++ b/src/hydro/Gizmo/hydro_gradients_gizmo.h
@@ -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.0 / 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.0 / 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.0 / 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;
diff --git a/src/hydro/Gizmo/hydro_gradients_sph.h b/src/hydro/Gizmo/hydro_gradients_sph.h
index b3e3eae5b7ca11b305085058960904a4a3877845..fbaa056443df56691f1414f2d661ba9edc459734 100644
--- a/src/hydro/Gizmo/hydro_gradients_sph.h
+++ b/src/hydro/Gizmo/hydro_gradients_sph.h
@@ -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;
diff --git a/src/hydro/Gizmo/hydro_iact.h b/src/hydro/Gizmo/hydro_iact.h
index 08c2480f4722d8704c250e5089786e93a5116da4..0c888e4159c6da6eb07eb0afbf758efce1c697ee 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.0 / 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.0 / 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.0 / 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. * 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.0 / 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.0 / 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
@@ -322,35 +317,36 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
float Xi = Vi;
float Xj = Vj;
#ifdef GIZMO_VOLUME_CORRECTION
- if (fabsf(Vi - Vj) / fminf(Vi, Vj) > 1.5 * hydro_dimension) {
+ if (fabsf(Vi - Vj) / min(Vi, Vj) > 1.5f * hydro_dimension) {
Xi = (Vi * hj + Vj * hi) / (hi + hj);
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.) {
- /* 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_cell.h b/src/hydro/Gizmo/hydro_slope_limiters_cell.h
index ccab0a5d00e8c13a2529efb55e49c87360a436ef..599be632b64e3b06ab1e0339a6ad4b560b2abece 100644
--- a/src/hydro/Gizmo/hydro_slope_limiters_cell.h
+++ b/src/hydro/Gizmo/hydro_slope_limiters_cell.h
@@ -57,29 +57,29 @@ hydro_slope_limit_cell_collect(struct part* pi, struct part* pj, float r) {
/* basic slope limiter: collect the maximal and the minimal value for the
* primitive variables among the ngbs */
pi->primitives.limiter.rho[0] =
- fmin(pj->primitives.rho, pi->primitives.limiter.rho[0]);
+ min(pj->primitives.rho, pi->primitives.limiter.rho[0]);
pi->primitives.limiter.rho[1] =
- fmax(pj->primitives.rho, pi->primitives.limiter.rho[1]);
+ max(pj->primitives.rho, pi->primitives.limiter.rho[1]);
pi->primitives.limiter.v[0][0] =
- fmin(pj->primitives.v[0], pi->primitives.limiter.v[0][0]);
+ min(pj->primitives.v[0], pi->primitives.limiter.v[0][0]);
pi->primitives.limiter.v[0][1] =
- fmax(pj->primitives.v[0], pi->primitives.limiter.v[0][1]);
+ max(pj->primitives.v[0], pi->primitives.limiter.v[0][1]);
pi->primitives.limiter.v[1][0] =
- fmin(pj->primitives.v[1], pi->primitives.limiter.v[1][0]);
+ min(pj->primitives.v[1], pi->primitives.limiter.v[1][0]);
pi->primitives.limiter.v[1][1] =
- fmax(pj->primitives.v[1], pi->primitives.limiter.v[1][1]);
+ max(pj->primitives.v[1], pi->primitives.limiter.v[1][1]);
pi->primitives.limiter.v[2][0] =
- fmin(pj->primitives.v[2], pi->primitives.limiter.v[2][0]);
+ min(pj->primitives.v[2], pi->primitives.limiter.v[2][0]);
pi->primitives.limiter.v[2][1] =
- fmax(pj->primitives.v[2], pi->primitives.limiter.v[2][1]);
+ max(pj->primitives.v[2], pi->primitives.limiter.v[2][1]);
pi->primitives.limiter.P[0] =
- fmin(pj->primitives.P, pi->primitives.limiter.P[0]);
+ min(pj->primitives.P, pi->primitives.limiter.P[0]);
pi->primitives.limiter.P[1] =
- fmax(pj->primitives.P, pi->primitives.limiter.P[1]);
+ max(pj->primitives.P, pi->primitives.limiter.P[1]);
- pi->primitives.limiter.maxr = fmax(r, pi->primitives.limiter.maxr);
+ pi->primitives.limiter.maxr = max(r, pi->primitives.limiter.maxr);
}
/**
@@ -90,8 +90,7 @@ hydro_slope_limit_cell_collect(struct part* pi, struct part* pj, float r) {
__attribute__((always_inline)) INLINE static void hydro_slope_limit_cell(
struct part* p) {
- float gradrho[3], gradv[3][3], gradP[3];
- float gradtrue, gradmax, gradmin, alpha;
+ float gradtrue, gradrho[3], gradv[3][3], gradP[3];
gradrho[0] = p->primitives.gradients.rho[0];
gradrho[1] = p->primitives.gradients.rho[1];
@@ -117,9 +116,11 @@ __attribute__((always_inline)) INLINE static void hydro_slope_limit_cell(
gradrho[2] * gradrho[2]);
if (gradtrue) {
gradtrue *= p->primitives.limiter.maxr;
- gradmax = p->primitives.limiter.rho[1] - p->primitives.rho;
- gradmin = p->primitives.rho - p->primitives.limiter.rho[0];
- alpha = fmin(1.0f, fmin(gradmax / gradtrue, gradmin / gradtrue));
+ const float gradmax = p->primitives.limiter.rho[1] - p->primitives.rho;
+ const float gradmin = p->primitives.rho - p->primitives.limiter.rho[0];
+ const float gradtrue_inv = 1.f / gradtrue;
+ const float alpha =
+ min3(1.0f, gradmax * gradtrue_inv, gradmin * gradtrue_inv);
p->primitives.gradients.rho[0] *= alpha;
p->primitives.gradients.rho[1] *= alpha;
p->primitives.gradients.rho[2] *= alpha;
@@ -129,9 +130,11 @@ __attribute__((always_inline)) INLINE static void hydro_slope_limit_cell(
gradv[0][2] * gradv[0][2]);
if (gradtrue) {
gradtrue *= p->primitives.limiter.maxr;
- gradmax = p->primitives.limiter.v[0][1] - p->primitives.v[0];
- gradmin = p->primitives.v[0] - p->primitives.limiter.v[0][0];
- alpha = fmin(1.0f, fmin(gradmax / gradtrue, gradmin / gradtrue));
+ const float gradmax = p->primitives.limiter.v[0][1] - p->primitives.v[0];
+ const float gradmin = p->primitives.v[0] - p->primitives.limiter.v[0][0];
+ const float gradtrue_inv = 1.f / gradtrue;
+ const float alpha =
+ min3(1.0f, gradmax * gradtrue_inv, gradmin * gradtrue_inv);
p->primitives.gradients.v[0][0] *= alpha;
p->primitives.gradients.v[0][1] *= alpha;
p->primitives.gradients.v[0][2] *= alpha;
@@ -141,9 +144,11 @@ __attribute__((always_inline)) INLINE static void hydro_slope_limit_cell(
gradv[1][2] * gradv[1][2]);
if (gradtrue) {
gradtrue *= p->primitives.limiter.maxr;
- gradmax = p->primitives.limiter.v[1][1] - p->primitives.v[1];
- gradmin = p->primitives.v[1] - p->primitives.limiter.v[1][0];
- alpha = fmin(1.0f, fmin(gradmax / gradtrue, gradmin / gradtrue));
+ const float gradmax = p->primitives.limiter.v[1][1] - p->primitives.v[1];
+ const float gradmin = p->primitives.v[1] - p->primitives.limiter.v[1][0];
+ const float gradtrue_inv = 1.f / gradtrue;
+ const float alpha =
+ min3(1.0f, gradmax * gradtrue_inv, gradmin * gradtrue_inv);
p->primitives.gradients.v[1][0] *= alpha;
p->primitives.gradients.v[1][1] *= alpha;
p->primitives.gradients.v[1][2] *= alpha;
@@ -153,9 +158,11 @@ __attribute__((always_inline)) INLINE static void hydro_slope_limit_cell(
gradv[2][2] * gradv[2][2]);
if (gradtrue) {
gradtrue *= p->primitives.limiter.maxr;
- gradmax = p->primitives.limiter.v[2][1] - p->primitives.v[2];
- gradmin = p->primitives.v[2] - p->primitives.limiter.v[2][0];
- alpha = fmin(1.0f, fmin(gradmax / gradtrue, gradmin / gradtrue));
+ const float gradmax = p->primitives.limiter.v[2][1] - p->primitives.v[2];
+ const float gradmin = p->primitives.v[2] - p->primitives.limiter.v[2][0];
+ const float gradtrue_inv = 1.f / gradtrue;
+ const float alpha =
+ min3(1.0f, gradmax * gradtrue_inv, gradmin * gradtrue_inv);
p->primitives.gradients.v[2][0] *= alpha;
p->primitives.gradients.v[2][1] *= alpha;
p->primitives.gradients.v[2][2] *= alpha;
@@ -165,9 +172,11 @@ __attribute__((always_inline)) INLINE static void hydro_slope_limit_cell(
sqrtf(gradP[0] * gradP[0] + gradP[1] * gradP[1] + gradP[2] * gradP[2]);
if (gradtrue) {
gradtrue *= p->primitives.limiter.maxr;
- gradmax = p->primitives.limiter.P[1] - p->primitives.P;
- gradmin = p->primitives.P - p->primitives.limiter.P[0];
- alpha = fmin(1.0f, fmin(gradmax / gradtrue, gradmin / gradtrue));
+ const float gradmax = p->primitives.limiter.P[1] - p->primitives.P;
+ const float gradmin = p->primitives.P - p->primitives.limiter.P[0];
+ const float gradtrue_inv = 1.f / gradtrue;
+ const float alpha =
+ min3(1.0f, gradmax * gradtrue_inv, gradmin * gradtrue_inv);
p->primitives.gradients.P[0] *= alpha;
p->primitives.gradients.P[1] *= alpha;
p->primitives.gradients.P[2] *= alpha;
diff --git a/src/hydro/Gizmo/hydro_slope_limiters_face.h b/src/hydro/Gizmo/hydro_slope_limiters_face.h
index 0d96dc4d434213008bd56104674c9bc6c60437b5..11e600d8feb877d5fe771684e8173782b3ff5923 100644
--- a/src/hydro/Gizmo/hydro_slope_limiters_face.h
+++ b/src/hydro/Gizmo/hydro_slope_limiters_face.h
@@ -32,52 +32,43 @@
#ifndef SWIFT_GIZMO_SLOPE_LIMITER_FACE_H
#define SWIFT_GIZMO_SLOPE_LIMITER_FACE_H
+#include <float.h>
+
+#include "sign.h"
+
__attribute__((always_inline)) INLINE static float
hydro_slope_limit_face_quantity(float phi_i, float phi_j, float phi_mid0,
- float xij_norm, float r) {
+ float xij_norm, float r_inv) {
- float delta1, delta2, phimin, phimax, phibar, phiplus, phiminus, phi_mid;
const float psi1 = 0.5f;
const float psi2 = 0.25f;
- if (phi_i == phi_j) {
- return 0.0f;
- }
+ const float delta1 = psi1 * fabsf(phi_i - phi_j);
+ const float delta2 = psi2 * fabsf(phi_i - phi_j);
- delta1 = psi1 * fabs(phi_i - phi_j);
- delta2 = psi2 * fabs(phi_i - phi_j);
+ const float phimin = min(phi_i, phi_j);
+ const float phimax = max(phi_i, phi_j);
- phimin = fmin(phi_i, phi_j);
- phimax = fmax(phi_i, phi_j);
+ const float phibar = phi_i + xij_norm * r_inv * (phi_j - phi_i);
- phibar = phi_i + xij_norm / r * (phi_j - phi_i);
+ float phiplus, phiminus, phi_mid;
- /* if sign(phimax+delta1) == sign(phimax) */
- if ((phimax + delta1) * phimax > 0.0f) {
+ if (same_signf(phimax + delta1, phimax))
phiplus = phimax + delta1;
- } else {
- if (phimax != 0.) {
- phiplus = phimax / (1.0f + delta1 / fabs(phimax));
- } else {
- phiplus = 0.;
- }
- }
+ else
+ phiplus = phimax / (1.0f + delta1 / (fabsf(phimax) + FLT_MIN));
- /* if sign(phimin-delta1) == sign(phimin) */
- if ((phimin - delta1) * phimin > 0.0f) {
+ if (same_signf(phimin - delta1, phimin))
phiminus = phimin - delta1;
- } else {
- if (phimin != 0.) {
- phiminus = phimin / (1.0f + delta1 / fabs(phimin));
- } else {
- phiminus = 0.;
- }
- }
+ else
+ phiminus = phimin / (1.0f + delta1 / (fabsf(phimin) + FLT_MIN));
if (phi_i < phi_j) {
- phi_mid = fmax(phiminus, fmin(phibar + delta2, phi_mid0));
+ const float temp = min(phibar + delta2, phi_mid0);
+ phi_mid = max(phiminus, temp);
} else {
- phi_mid = fmin(phiplus, fmax(phibar - delta2, phi_mid0));
+ const float temp = max(phibar - delta2, phi_mid0);
+ phi_mid = min(phiplus, temp);
}
return phi_mid - phi_i;
@@ -97,38 +88,38 @@ 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 r) {
+ float *Wi, float *Wj, float *dWi, float *dWj, const float *xij_i,
+ const float *xij_j, float r) {
- float xij_i_norm, xij_j_norm;
-
- xij_i_norm =
+ const float xij_i_norm =
sqrtf(xij_i[0] * xij_i[0] + xij_i[1] * xij_i[1] + xij_i[2] * xij_i[2]);
- xij_j_norm =
+ const float xij_j_norm =
sqrtf(xij_j[0] * xij_j[0] + xij_j[1] * xij_j[1] + xij_j[2] * xij_j[2]);
+ const float r_inv = 1.f / r;
+
dWi[0] = hydro_slope_limit_face_quantity(Wi[0], Wj[0], Wi[0] + dWi[0],
- xij_i_norm, r);
+ xij_i_norm, r_inv);
dWi[1] = hydro_slope_limit_face_quantity(Wi[1], Wj[1], Wi[1] + dWi[1],
- xij_i_norm, r);
+ xij_i_norm, r_inv);
dWi[2] = hydro_slope_limit_face_quantity(Wi[2], Wj[2], Wi[2] + dWi[2],
- xij_i_norm, r);
+ xij_i_norm, r_inv);
dWi[3] = hydro_slope_limit_face_quantity(Wi[3], Wj[3], Wi[3] + dWi[3],
- xij_i_norm, r);
+ xij_i_norm, r_inv);
dWi[4] = hydro_slope_limit_face_quantity(Wi[4], Wj[4], Wi[4] + dWi[4],
- xij_i_norm, r);
+ xij_i_norm, r_inv);
dWj[0] = hydro_slope_limit_face_quantity(Wj[0], Wi[0], Wj[0] + dWj[0],
- xij_j_norm, r);
+ xij_j_norm, r_inv);
dWj[1] = hydro_slope_limit_face_quantity(Wj[1], Wi[1], Wj[1] + dWj[1],
- xij_j_norm, r);
+ xij_j_norm, r_inv);
dWj[2] = hydro_slope_limit_face_quantity(Wj[2], Wi[2], Wj[2] + dWj[2],
- xij_j_norm, r);
+ xij_j_norm, r_inv);
dWj[3] = hydro_slope_limit_face_quantity(Wj[3], Wi[3], Wj[3] + dWj[3],
- xij_j_norm, r);
+ xij_j_norm, r_inv);
dWj[4] = hydro_slope_limit_face_quantity(Wj[4], Wi[4], Wj[4] + dWj[4],
- xij_j_norm, r);
+ xij_j_norm, r_inv);
}
#endif /* SWIFT_GIZMO_SLOPE_LIMITER_FACE_H */
diff --git a/src/hydro/Gizmo/hydro_unphysical.h b/src/hydro/Gizmo/hydro_unphysical.h
index f1e64a281c13a4e31eb8871139d8f28c95c26928..81c644e6cce00e0d871aafc39140e420e042a926 100644
--- a/src/hydro/Gizmo/hydro_unphysical.h
+++ b/src/hydro/Gizmo/hydro_unphysical.h
@@ -41,14 +41,28 @@
#endif
#define gizmo_check_physical_quantity(name, quantity) \
- if (quantity < 0.) { \
+ if (quantity < 0.f) { \
gizmo_unphysical_message(name, quantity); \
- quantity = 0.; \
+ quantity = 0.f; \
+ }
+
+#define gizmo_check_physical_quantities( \
+ mass_name, energy_name, mass, momentum_x, momentum_y, momentum_z, energy) \
+ gizmo_check_physical_quantity(mass_name, mass); \
+ gizmo_check_physical_quantity(energy_name, energy); \
+ /* now check for vacuum and make sure we have a real vacuum */ \
+ if (mass == 0.f || energy == 0.f) { \
+ mass = 0.f; \
+ momentum_x = 0.f; \
+ momentum_y = 0.f; \
+ momentum_z = 0.f; \
+ energy = 0.f; \
}
#else // defined(GIZMO_UNPHYSICAL_ERROR) || defined(GIZMO_UNPHYSICAL_RESCUE)
-#define gizmo_check_physical_quantity(name, quantity)
+#define gizmo_check_physical_quantities( \
+ mass_name, energy_name, mass, momentum_x, momentum_y, momentum_z, energy)
#endif // defined(GIZMO_UNPHYSICAL_ERROR) || defined(GIZMO_UNPHYSICAL_RESCUE)
diff --git a/src/hydro/Gizmo/hydro_velocities.h b/src/hydro/Gizmo/hydro_velocities.h
index 2f7bc0d54d6130c253f72be34463078a8931725d..ea30d5a6c9c74d34ffd73fa6ab941640f37b02e4 100644
--- a/src/hydro/Gizmo/hydro_velocities.h
+++ b/src/hydro/Gizmo/hydro_velocities.h
@@ -30,9 +30,9 @@ __attribute__((always_inline)) INLINE static void hydro_velocities_init(
struct part* restrict p, struct xpart* restrict xp) {
#ifdef GIZMO_FIX_PARTICLES
- p->v[0] = 0.;
- p->v[1] = 0.;
- p->v[2] = 0.;
+ p->v[0] = 0.f;
+ p->v[1] = 0.f;
+ p->v[2] = 0.f;
#else
p->v[0] = p->primitives.v[0];
p->v[1] = p->primitives.v[1];
@@ -87,17 +87,20 @@ __attribute__((always_inline)) INLINE static void hydro_velocities_set(
/* We first set the particle velocity. */
#ifdef GIZMO_FIX_PARTICLES
- p->v[0] = 0.;
- p->v[1] = 0.;
- p->v[2] = 0.;
+ p->v[0] = 0.f;
+ p->v[1] = 0.f;
+ p->v[2] = 0.f;
#else // GIZMO_FIX_PARTICLES
- if (p->conserved.mass > 0. && p->primitives.rho > 0.) {
+ if (p->conserved.mass > 0.f && p->primitives.rho > 0.f) {
+
+ const float inverse_mass = 1.f / p->conserved.mass;
+
/* Normal case: set particle velocity to fluid velocity. */
- p->v[0] = p->conserved.momentum[0] / p->conserved.mass;
- p->v[1] = p->conserved.momentum[1] / p->conserved.mass;
- p->v[2] = p->conserved.momentum[2] / p->conserved.mass;
+ p->v[0] = p->conserved.momentum[0] * inverse_mass;
+ p->v[1] = p->conserved.momentum[1] * inverse_mass;
+ p->v[2] = p->conserved.momentum[2] * inverse_mass;
#ifdef GIZMO_STEER_MOTION
@@ -112,18 +115,17 @@ __attribute__((always_inline)) INLINE static void hydro_velocities_set(
ds[2] = p->geometry.centroid[2];
const float d = sqrtf(ds[0] * ds[0] + ds[1] * ds[1] + ds[2] * ds[2]);
const float R = get_radius_dimension_sphere(p->geometry.volume);
- const float eta = 0.25;
+ const float eta = 0.25f;
const float etaR = eta * R;
- const float xi = 1.;
+ const float xi = 1.f;
const float soundspeed =
sqrtf(hydro_gamma * p->primitives.P / p->primitives.rho);
/* We only apply the correction if the offset between centroid and position
- is
- too large. */
- if (d > 0.9 * etaR) {
+ is too large. */
+ if (d > 0.9f * etaR) {
float fac = xi * soundspeed / d;
- if (d < 1.1 * etaR) {
- fac *= 5. * (d - 0.9 * etaR) / etaR;
+ if (d < 1.1f * etaR) {
+ fac *= 5.f * (d - 0.9f * etaR) / etaR;
}
p->v[0] -= ds[0] * fac;
p->v[1] -= ds[1] * fac;
@@ -133,9 +135,9 @@ __attribute__((always_inline)) INLINE static void hydro_velocities_set(
#endif // GIZMO_STEER_MOTION
} else {
/* Vacuum particles have no fluid velocity. */
- p->v[0] = 0.;
- p->v[1] = 0.;
- p->v[2] = 0.;
+ p->v[0] = 0.f;
+ p->v[1] = 0.f;
+ p->v[2] = 0.f;
}
#endif // GIZMO_FIX_PARTICLES
diff --git a/src/riemann/riemann_checks.h b/src/riemann/riemann_checks.h
new file mode 100644
index 0000000000000000000000000000000000000000..2a158f2f84b3c383d75750ca9600f73cbf280400
--- /dev/null
+++ b/src/riemann/riemann_checks.h
@@ -0,0 +1,207 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2018 Bert Vandenbroucke (bert.vandenbroucke@ugent.be)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+
+#ifndef SWIFT_RIEMANN_CHECKS_H
+#define SWIFT_RIEMANN_CHECKS_H
+
+#include "error.h"
+
+#ifdef SWIFT_DEBUG_CHECKS
+
+/**
+ * @brief Check if the given state vector is physically valid.
+ *
+ * A valid state vector has 5 finite elements (not NaN), and a positive density
+ * and pressure (first and fifth element).
+ *
+ * @param W State vector.
+ * @return 0 if the state vector is valid, 1 otherwise.
+ */
+__attribute__((always_inline)) INLINE static int riemann_check_state(
+ const float *W) {
+
+ int errorFlag = 0;
+
+ /* check the density: should be finite and positive */
+ if (W[0] != W[0]) {
+ message("NaN density!");
+ errorFlag = 1;
+ }
+ if (W[0] < 0.f) {
+ message("Negative density!");
+ errorFlag = 1;
+ }
+
+ /* check the velocities: should be finite */
+ if (W[1] != W[1]) {
+ message("NaN x velocity!");
+ errorFlag = 1;
+ }
+ if (W[2] != W[2]) {
+ message("NaN y velocity!");
+ errorFlag = 1;
+ }
+ if (W[3] != W[3]) {
+ message("NaN z velocity!");
+ errorFlag = 1;
+ }
+
+ /* check the pressure: should be positive and finite */
+ if (W[4] != W[4]) {
+ message("NaN pressure!");
+ errorFlag = 1;
+ }
+ if (W[4] < 0.f) {
+ message("Negative pressure!");
+ errorFlag = 1;
+ }
+
+ return errorFlag;
+}
+
+/**
+ * @brief Check that the given vector is physically valid.
+ *
+ * A valid vector has 3 finite elements (not NaN).
+ *
+ * @param x Vector to check.
+ * @return 0 if the vector is valid, 1 otherwise.
+ */
+__attribute__((always_inline)) INLINE static int riemann_check_vector(
+ const float *x) {
+
+ int errorFlag = 0;
+
+ /* check that all components are finite */
+ if (x[0] != x[0]) {
+ message("NaN x component!");
+ errorFlag = 1;
+ }
+ if (x[1] != x[1]) {
+ message("NaN y component!");
+ errorFlag = 1;
+ }
+ if (x[2] != x[2]) {
+ message("NaN z component!");
+ errorFlag = 1;
+ }
+
+ return errorFlag;
+}
+
+/**
+ * @brief Check that the given input for a Riemann solver is physically valid.
+ *
+ * Physically valid input consists of a valid left and right state, and valid
+ * surface normal and surface velocity vectors.
+ * If no valid input is provided, an error is thrown.
+ *
+ * @param WL Left state vector.
+ * @param WR Right state vector.
+ * @param n Surface normal vector.
+ * @param vij Surface velocity vector.
+ */
+__attribute__((always_inline)) INLINE static void riemann_check_input(
+ const float *WL, const float *WR, const float *n, const float *vij) {
+
+ int errorFlag = 0;
+
+ if (riemann_check_state(WL)) {
+ message("Invalid left state!");
+ errorFlag = 1;
+ }
+
+ if (riemann_check_state(WR)) {
+ message("Invalid right state!");
+ errorFlag = 1;
+ }
+
+ if (riemann_check_vector(n)) {
+ message("Invalid surface normal vector!");
+ errorFlag = 1;
+ }
+
+ if (riemann_check_vector(vij)) {
+ message("Invalid face velocity vector!");
+ errorFlag = 1;
+ }
+
+ if (errorFlag) {
+ message("WL: %g %g %g %g %g", WL[0], WL[1], WL[2], WL[3], WL[4]);
+ message("WR: %g %g %g %g %g", WR[0], WR[1], WR[2], WR[3], WR[4]);
+ message("n: %g %g %g", n[0], n[1], n[2]);
+ message("vij: %g %g %g", vij[0], vij[1], vij[2]);
+ error("Invalid Riemann solver input!");
+ }
+}
+
+/**
+ * @brief Check that the given output from a Riemann solver is physically valid.
+ *
+ * Physically valid output consists of 5 finite (not NaN) flux components.
+ * If no valid output is provided, an error is thrown.
+ *
+ * @param WL Left state vector.
+ * @param WR Right state vector.
+ * @param n Surface normal vector.
+ * @param vij Surface velocity vector.
+ * @param totflux Riemann solver flux result.
+ */
+__attribute__((always_inline)) INLINE static void riemann_check_output(
+ const float *WL, const float *WR, const float *n, const float *vij,
+ const float *totflux) {
+
+ int errorFlag = 0;
+
+ /* check that all the fluxes are finite */
+ if (totflux[0] != totflux[0]) {
+ message("NaN mass flux!");
+ errorFlag = 1;
+ }
+ if (totflux[1] != totflux[1]) {
+ message("NaN x momentum flux!");
+ errorFlag = 1;
+ }
+ if (totflux[2] != totflux[2]) {
+ message("NaN y momentum flux!");
+ errorFlag = 1;
+ }
+ if (totflux[3] != totflux[3]) {
+ message("NaN z momentum flux!");
+ errorFlag = 1;
+ }
+ if (totflux[4] != totflux[4]) {
+ message("NaN energy flux!");
+ errorFlag = 1;
+ }
+
+ if (errorFlag) {
+ message("WL: %g %g %g %g %g", WL[0], WL[1], WL[2], WL[3], WL[4]);
+ message("WR: %g %g %g %g %g", WR[0], WR[1], WR[2], WR[3], WR[4]);
+ message("n: %g %g %g", n[0], n[1], n[2]);
+ message("vij: %g %g %g", vij[0], vij[1], vij[2]);
+ message("totflux: %g %g %g %g %g", totflux[0], totflux[1], totflux[2],
+ totflux[3], totflux[4]);
+ error("Invalid Riemann solver output!");
+ }
+}
+
+#endif
+
+#endif /* SWIFT_RIEMANN_CHECKS_H */
diff --git a/src/riemann/riemann_exact.h b/src/riemann/riemann_exact.h
index 4a7b40ae78c4c3f1a880539367172894fe298995..0c6c1c86273b42f521102e6e91d642d3dff30b27 100644
--- a/src/riemann/riemann_exact.h
+++ b/src/riemann/riemann_exact.h
@@ -1,6 +1,7 @@
/*******************************************************************************
* This file is part of SWIFT.
- * Coypright (c) 2015 Bert Vandenbroucke (bert.vandenbroucke@ugent.be)
+ * Copyright (c) 2015 Bert Vandenbroucke (bert.vandenbroucke@ugent.be)
+ * 2018 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as published
@@ -38,6 +39,7 @@
#include "adiabatic_index.h"
#include "error.h"
#include "minmax.h"
+#include "riemann_checks.h"
#include "riemann_vacuum.h"
/**
@@ -47,14 +49,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 +77,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 +90,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 +114,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 +126,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 +147,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 +199,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 +306,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;
@@ -320,30 +320,6 @@ __attribute__((always_inline)) INLINE static void riemann_solver_solve(
float SHR, STR;
float pdpL, SL;
float SHL, STL;
- int errorFlag = 0;
-
- /* sanity checks */
- if (WL[0] != WL[0] || WL[4] != WL[4]) {
- printf("NaN WL!\n");
- errorFlag = 1;
- }
- if (WR[0] != WR[0] || WR[4] != WR[4]) {
- printf("NaN WR!\n");
- errorFlag = 1;
- }
- if (WL[0] < 0.0f || WL[4] < 0.0f) {
- printf("Negative WL!\n");
- errorFlag = 1;
- }
- if (WR[0] < 0.0f || WR[4] < 0.0f) {
- printf("Negative WR!\n");
- errorFlag = 1;
- }
- if (errorFlag) {
- printf("WL: %g %g %g %g %g\n", WL[0], WL[1], WL[2], WL[3], WL[4]);
- printf("WR: %g %g %g %g %g\n", WR[0], WR[1], WR[2], WR[3], WR[4]);
- error("Riemman solver input error!\n");
- }
/* calculate velocities in interface frame */
vL = WL[1] * n_unit[0] + WL[2] * n_unit[1] + WL[3] * n_unit[2];
@@ -510,7 +486,12 @@ __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) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+ riemann_check_input(Wi, Wj, n_unit, vij);
+#endif
float Whalf[5];
float flux[5][3];
@@ -556,6 +537,10 @@ __attribute__((always_inline)) INLINE static void riemann_solve_for_flux(
flux[3][0] * n_unit[0] + flux[3][1] * n_unit[1] + flux[3][2] * n_unit[2];
totflux[4] =
flux[4][0] * n_unit[0] + flux[4][1] * n_unit[1] + flux[4][2] * n_unit[2];
+
+#ifdef SWIFT_DEBUG_CHECKS
+ riemann_check_output(Wi, Wj, n_unit, vij, totflux);
+#endif
}
#endif /* SWIFT_RIEMANN_EXACT_H */
diff --git a/src/riemann/riemann_exact_isothermal.h b/src/riemann/riemann_exact_isothermal.h
deleted file mode 100644
index 26fb649f02ddd62c102e31c9191c20f9003c0133..0000000000000000000000000000000000000000
--- a/src/riemann/riemann_exact_isothermal.h
+++ /dev/null
@@ -1,466 +0,0 @@
-/*******************************************************************************
- * This file is part of SWIFT.
- * Coypright (c) 2016 Bert Vandenbroucke (bert.vandenbroucke@ugent.be)
- *
- * This program is free software: you can redistribute it and/or modify
- * it under the terms of the GNU Lesser General Public License as published
- * by the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public License
- * along with this program. If not, see <http://www.gnu.org/licenses/>.
- *
- ******************************************************************************/
-
-#ifndef SWIFT_RIEMANN_EXACT_ISOTHERMAL_H
-#define SWIFT_RIEMANN_EXACT_ISOTHERMAL_H
-
-#include <float.h>
-#include "adiabatic_index.h"
-#include "minmax.h"
-#include "riemann_vacuum.h"
-
-#define const_isothermal_soundspeed \
- sqrtf(hydro_gamma_minus_one* const_isothermal_internal_energy)
-
-/**
- * @brief Relative difference between the middle state velocity and the left or
- * right state velocity used in the middle state density iteration.
- *
- * @param rho Current estimate of the middle state density.
- * @param W Left or right state vector.
- * @return Density dependent part of the middle state velocity.
- */
-__attribute__((always_inline)) INLINE static float riemann_fb(float rho,
- float* W) {
- if (rho < W[0]) {
- return const_isothermal_soundspeed * logf(rho / W[0]);
- } else {
- return const_isothermal_soundspeed *
- (sqrtf(rho / W[0]) - sqrtf(W[0] / rho));
- }
-}
-
-/**
- * @brief Derivative w.r.t. rho of the function riemann_fb.
- *
- * @param rho Current estimate of the middle state density.
- * @param W Left or right state vector.
- * @return Derivative of riemann_fb.
- */
-__attribute__((always_inline)) INLINE static float riemann_fprimeb(float rho,
- float* W) {
- if (rho < W[0]) {
- return const_isothermal_soundspeed * W[0] / rho;
- } else {
- return 0.5 * const_isothermal_soundspeed *
- (sqrtf(rho / W[0]) + sqrtf(W[0] / rho)) / rho;
- }
-}
-
-/**
- * @brief Difference between the left and right middle state velocity estimates.
- *
- * Since the middle state velocity takes on a constant value, we want to get
- * this difference as close to zero as possible.
- *
- * @param rho Current estimate of the middle state density.
- * @param WL Left state vector.
- * @param WR Right state vector.
- * @param vL Left state velocity along the interface normal.
- * @param vR Right state velocity along the interface normal.
- * @return Difference between the left and right middle state velocity
- * estimates.
- */
-__attribute__((always_inline)) INLINE static float riemann_f(
- float rho, float* WL, float* WR, float vL, float vR) {
- return riemann_fb(rho, WR) + riemann_fb(rho, WL) + vR - vL;
-}
-
-/**
- * @brief Derivative of riemann_f w.r.t. rho.
- *
- * @param rho Current estimate of the middle state density.
- * @param WL Left state vector.
- * @param WR Right state vector.
- * @return Derivative of riemann_f.
- */
-__attribute__((always_inline)) INLINE static float riemann_fprime(float rho,
- float* WL,
- float* WR) {
- return riemann_fprimeb(rho, WL) + riemann_fprimeb(rho, WR);
-}
-
-/**
- * @brief Get a good first guess for the middle state density.
- *
- * @param WL The left state vector
- * @param WR The right state vector
- * @param vL The left velocity along the interface normal
- * @param vR The right velocity along the interface normal
- */
-__attribute__((always_inline)) INLINE static float riemann_guess_rho(float* WL,
- float* WR,
- float vL,
- float vR) {
-
- /* Currently three possibilities and not really an algorithm to decide which
- one to choose: */
- /* just the average */
- // return 0.5f * (WL[0] + WR[0]);
-
- /* two rarefaction approximation */
- return sqrtf(WL[0] * WR[0] * expf((vL - vR) / const_isothermal_soundspeed));
-
- /* linearized primitive variable approximation */
- return 0.25f * (WL[0] + WR[0]) * (vL - vR) / const_isothermal_soundspeed +
- 0.5f * (WL[0] + WR[0]);
-}
-
-/**
- * @brief Find the zeropoint of riemann_f(rho) using Brent's method.
- *
- * @param lower_limit Lower limit for the method (riemann_f(lower_limit) < 0)
- * @param upper_limit Upper limit for the method (riemann_f(upper_limit) > 0)
- * @param lowf Value of riemann_f(lower_limit).
- * @param upf Value of riemann_f(upper_limit).
- * @param error_tol Tolerance used to decide if the solution is converged
- * @param WL Left state vector
- * @param WR Right state vector
- * @param vL The left velocity along the interface normal
- * @param vR The right velocity along the interface normal
- */
-__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 a, b, c, d, s;
- float fa, fb, fc, fs;
- float tmp, tmp2;
- int mflag;
- int i;
-
- a = lower_limit;
- b = upper_limit;
- c = 0.0f;
- d = FLT_MAX;
-
- fa = lowf;
- fb = upf;
-
- fc = 0.0f;
- s = 0.0f;
- fs = 0.0f;
-
- /* if f(a) f(b) >= 0 then error-exit */
- if (fa * fb >= 0.0f) {
- error(
- "Brent's method called with equal sign function values!\n"
- "f(%g) = %g, f(%g) = %g\n",
- a, fa, b, fb);
- /* return NaN */
- return 0.0f / 0.0f;
- }
-
- /* if |f(a)| < |f(b)| then swap (a,b) */
- if (fabs(fa) < fabs(fb)) {
- tmp = a;
- a = b;
- b = tmp;
- tmp = fa;
- fa = fb;
- fb = tmp;
- }
-
- c = a;
- fc = fa;
- mflag = 1;
- i = 0;
-
- while (!(fb == 0.0f) && (fabs(a - b) > error_tol * 0.5f * (a + b))) {
- if ((fa != fc) && (fb != fc)) /* Inverse quadratic interpolation */
- s = a * fb * fc / (fa - fb) / (fa - fc) +
- b * fa * fc / (fb - fa) / (fb - fc) +
- c * fa * fb / (fc - fa) / (fc - fb);
- else
- /* Secant Rule */
- s = b - fb * (b - a) / (fb - fa);
-
- tmp2 = 0.25f * (3.0f * a + b);
- if (!(((s > tmp2) && (s < b)) || ((s < tmp2) && (s > b))) ||
- (mflag && (fabs(s - b) >= (0.5f * fabs(b - c)))) ||
- (!mflag && (fabs(s - b) >= (0.5f * fabs(c - d)))) ||
- (mflag && (fabs(b - c) < error_tol)) ||
- (!mflag && (fabs(c - d) < error_tol))) {
- s = 0.5f * (a + b);
- mflag = 1;
- } else {
- mflag = 0;
- }
- fs = riemann_f(s, WL, WR, vL, vR);
- d = c;
- c = b;
- fc = fb;
- if (fa * fs < 0.) {
- b = s;
- fb = fs;
- } else {
- a = s;
- fa = fs;
- }
-
- /* if |f(a)| < |f(b)| then swap (a,b) */
- if (fabs(fa) < fabs(fb)) {
- tmp = a;
- a = b;
- b = tmp;
- tmp = fa;
- fa = fb;
- fb = tmp;
- }
- i++;
- }
- return b;
-}
-
-/**
- * @brief Solve the Riemann problem between the given left and right state and
- * along the given interface normal
- *
- * @param WL The left state vector
- * @param WR The right state vector
- * @param Whalf Empty state vector in which the result will be stored
- * @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) {
-
- /* velocity of the left and right state in a frame aligned with n_unit */
- float vL, vR, vhalf;
- /* variables used for finding rhostar */
- float rho, rhoguess, frho, frhoguess;
- /* variables used for sampling the solution */
- float u, S, SH, ST;
-
- int errorFlag = 0;
-
- /* sanity checks */
- if (WL[0] != WL[0]) {
- printf("NaN WL!\n");
- errorFlag = 1;
- }
- if (WR[0] != WR[0]) {
- printf("NaN WR!\n");
- errorFlag = 1;
- }
- if (WL[0] < 0.0f) {
- printf("Negative WL!\n");
- errorFlag = 1;
- }
- if (WR[0] < 0.0f) {
- printf("Negative WR!\n");
- errorFlag = 1;
- }
- if (errorFlag) {
- printf("WL: %g %g %g %g %g\n", WL[0], WL[1], WL[2], WL[3], WL[4]);
- printf("WR: %g %g %g %g %g\n", WR[0], WR[1], WR[2], WR[3], WR[4]);
- error("Riemman solver input error!\n");
- }
-
- /* calculate velocities in interface frame */
- vL = WL[1] * n_unit[0] + WL[2] * n_unit[1] + WL[3] * n_unit[2];
- vR = WR[1] * n_unit[0] + WR[2] * n_unit[1] + WR[3] * n_unit[2];
-
- if (WL[0] == 0. || WR[0] == 0.) {
- error(
- "One of the states is vacuum, the isothermal solver cannot solve "
- "this!");
- }
-
- rho = 0.;
- /* obtain a first guess for p */
- rhoguess = riemann_guess_rho(WL, WR, vL, vR);
-
-#ifdef SWIFT_DEBUG_CHECKS
- if (rhoguess <= 0.) {
- error("Zero or negative initial density guess.");
- }
-#endif
-
- /* we know that the value of riemann_f for rho=0 is negative (it is -inf). */
- frho = -1.;
- frhoguess = riemann_f(rhoguess, WL, WR, vL, vR);
- /* ok, rhostar is close to 0, better use Brent's method... */
- /* we use Newton-Raphson until we find a suitable interval */
- if (frho * frhoguess >= 0.0f) {
- /* Newton-Raphson until convergence or until suitable interval is found
- to use Brent's method */
- unsigned int counter = 0;
- while (fabs(rho - rhoguess) > 5.e-7f * (rho + rhoguess) &&
- frhoguess < 0.0f) {
- rho = rhoguess;
- rhoguess = rhoguess - frhoguess / riemann_fprime(rhoguess, WL, WR);
- frhoguess = riemann_f(rhoguess, WL, WR, vL, vR);
- counter++;
- if (counter > 1000) {
- error(
- "Stuck in Newton-Raphson (rho: %g, rhoguess: %g, frhoguess: %g, "
- "fprime: %g, rho-rhoguess: %g, WL: %g %g %g, WR: %g %g %g)!\n",
- rho, rhoguess, frhoguess, riemann_fprime(rhoguess, WL, WR),
- (rho - rhoguess), WL[0], vL, WL[4], WR[0], vR, WR[4]);
- }
- }
- }
- /* As soon as there is a suitable interval: use Brent's method */
- if (1.e6 * fabs(rho - rhoguess) > 0.5f * (rho + rhoguess) &&
- frhoguess > 0.0f) {
- rho = 0.0f;
- frho = -1.;
- /* use Brent's method to find the zeropoint */
- rho = riemann_solve_brent(rho, rhoguess, frho, frhoguess, 1.e-6, WL, WR, vL,
- vR);
- } else {
- rho = rhoguess;
- }
-
- /* calculate the middle state velocity */
- u = 0.5f * (vL - riemann_fb(rho, WL) + vR + riemann_fb(rho, WR));
-
- /* sample the solution */
- if (u > 0.0f) {
- /* left state */
- Whalf[1] = WL[1];
- Whalf[2] = WL[2];
- Whalf[3] = WL[3];
- if (WL[0] < rho) {
- /* left shock wave */
- S = vL - const_isothermal_soundspeed * sqrtf(rho / WL[0]);
- if (S >= 0.) {
- /* to the left of the shock */
- Whalf[0] = WL[0];
- vhalf = 0.0f;
- } else {
- /* to the right of the shock */
- Whalf[0] = rho;
- vhalf = u - vL;
- }
- } else {
- /* left rarefaction wave */
- SH = vL - const_isothermal_soundspeed;
- ST = u - const_isothermal_soundspeed;
- if (SH > 0.) {
- /* to the left of the rarefaction */
- Whalf[0] = WL[0];
- vhalf = 0.0f;
- } else if (ST > 0.0f) {
- /* inside the rarefaction */
- Whalf[0] = WL[0] * expf(vL / const_isothermal_soundspeed - 1.0f);
- vhalf = const_isothermal_soundspeed - vL;
- } else {
- /* to the right of the rarefaction */
- Whalf[0] = rho;
- vhalf = u - vL;
- }
- }
- } else {
- /* right state */
- Whalf[1] = WR[1];
- Whalf[2] = WR[2];
- Whalf[3] = WR[3];
- if (WR[0] < rho) {
- /* right shock wave */
- S = vR + const_isothermal_soundspeed * sqrtf(rho / WR[0]);
- if (S > 0.0f) {
- /* to the left of the shock wave: middle state */
- Whalf[0] = rho;
- vhalf = u - vR;
- } else {
- /* to the right of the shock wave: right state */
- Whalf[0] = WR[0];
- vhalf = 0.0f;
- }
- } else {
- /* right rarefaction wave */
- SH = vR + const_isothermal_soundspeed;
- ST = u + const_isothermal_soundspeed;
- if (ST > 0.0f) {
- /* to the left of rarefaction: middle state */
- Whalf[0] = rho;
- vhalf = u - vR;
- } else if (SH > 0.0f) {
- /* inside rarefaction */
- Whalf[0] = WR[0] * expf(-vR / const_isothermal_soundspeed - 1.0f);
- vhalf = -const_isothermal_soundspeed - vR;
- } else {
- /* to the right of rarefaction: right state */
- Whalf[0] = WR[0];
- vhalf = 0.0f;
- }
- }
- }
-
- /* add the velocity solution along the interface normal to the velocities */
- Whalf[1] += vhalf * n_unit[0];
- Whalf[2] += vhalf * n_unit[1];
- Whalf[3] += vhalf * n_unit[2];
-
- /* the pressure is completely irrelevant in this case */
- Whalf[4] =
- Whalf[0] * const_isothermal_soundspeed * const_isothermal_soundspeed;
-}
-
-__attribute__((always_inline)) INLINE static void riemann_solve_for_flux(
- float* Wi, float* Wj, float* n_unit, float* vij, float* totflux) {
-
- float Whalf[5];
- float flux[5][3];
- float vtot[3];
- float rhoe;
-
- riemann_solver_solve(Wi, Wj, Whalf, n_unit);
-
- flux[0][0] = Whalf[0] * Whalf[1];
- flux[0][1] = Whalf[0] * Whalf[2];
- flux[0][2] = Whalf[0] * Whalf[3];
-
- vtot[0] = Whalf[1] + vij[0];
- vtot[1] = Whalf[2] + vij[1];
- vtot[2] = Whalf[3] + vij[2];
- flux[1][0] = Whalf[0] * vtot[0] * Whalf[1] + Whalf[4];
- flux[1][1] = Whalf[0] * vtot[0] * Whalf[2];
- flux[1][2] = Whalf[0] * vtot[0] * Whalf[3];
- flux[2][0] = Whalf[0] * vtot[1] * Whalf[1];
- flux[2][1] = Whalf[0] * vtot[1] * Whalf[2] + Whalf[4];
- flux[2][2] = Whalf[0] * vtot[1] * Whalf[3];
- flux[3][0] = Whalf[0] * vtot[2] * Whalf[1];
- flux[3][1] = Whalf[0] * vtot[2] * Whalf[2];
- flux[3][2] = Whalf[0] * vtot[2] * Whalf[3] + Whalf[4];
-
- /* eqn. (15) */
- /* F_P = \rho e ( \vec{v} - \vec{v_{ij}} ) + P \vec{v} */
- /* \rho e = P / (\gamma-1) + 1/2 \rho \vec{v}^2 */
- rhoe = Whalf[4] / hydro_gamma_minus_one +
- 0.5f * Whalf[0] *
- (vtot[0] * vtot[0] + vtot[1] * vtot[1] + vtot[2] * vtot[2]);
- flux[4][0] = rhoe * Whalf[1] + Whalf[4] * vtot[0];
- flux[4][1] = rhoe * Whalf[2] + Whalf[4] * vtot[1];
- flux[4][2] = rhoe * Whalf[3] + Whalf[4] * vtot[2];
-
- totflux[0] =
- flux[0][0] * n_unit[0] + flux[0][1] * n_unit[1] + flux[0][2] * n_unit[2];
- totflux[1] =
- flux[1][0] * n_unit[0] + flux[1][1] * n_unit[1] + flux[1][2] * n_unit[2];
- totflux[2] =
- flux[2][0] * n_unit[0] + flux[2][1] * n_unit[1] + flux[2][2] * n_unit[2];
- totflux[3] =
- flux[3][0] * n_unit[0] + flux[3][1] * n_unit[1] + flux[3][2] * n_unit[2];
- totflux[4] =
- flux[4][0] * n_unit[0] + flux[4][1] * n_unit[1] + flux[4][2] * n_unit[2];
-}
-
-#endif /* SWIFT_RIEMANN_EXACT_ISOTHERMAL_H */
diff --git a/src/riemann/riemann_hllc.h b/src/riemann/riemann_hllc.h
index e9a32cb93689e4beccc2a53831c56fcb612eac54..147c8eda53a7b9cdc69c64889060d6e7696b41b9 100644
--- a/src/riemann/riemann_hllc.h
+++ b/src/riemann/riemann_hllc.h
@@ -1,6 +1,7 @@
/*******************************************************************************
* This file is part of SWIFT.
- * Coypright (c) 2015 Bert Vandenbroucke (bert.vandenbroucke@ugent.be)
+ * Copyright (c) 2015 Bert Vandenbroucke (bert.vandenbroucke@ugent.be)
+ * 2018 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as published
@@ -30,6 +31,7 @@
#include "adiabatic_index.h"
#include "error.h"
#include "minmax.h"
+#include "riemann_checks.h"
#include "riemann_vacuum.h"
#ifndef EOS_IDEAL_GAS
@@ -38,28 +40,28 @@
#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) {
- float uL, uR, aL, aR;
- float rhobar, abar, pPVRS, pstar, qL, qR, SL, SR, Sstar;
- float v2, eL, eR;
- float UstarL[5], UstarR[5];
+#ifdef SWIFT_DEBUG_CHECKS
+ riemann_check_input(WL, WR, n, vij);
+#endif
/* Handle pure vacuum */
if (!WL[0] && !WR[0]) {
- totflux[0] = 0.;
- totflux[1] = 0.;
- totflux[2] = 0.;
- totflux[3] = 0.;
- totflux[4] = 0.;
+ totflux[0] = 0.f;
+ totflux[1] = 0.f;
+ totflux[2] = 0.f;
+ totflux[3] = 0.f;
+ totflux[4] = 0.f;
return;
}
/* STEP 0: obtain velocity in interface frame */
- uL = WL[1] * n[0] + WL[2] * n[1] + WL[3] * n[2];
- uR = WR[1] * n[0] + WR[2] * n[1] + WR[3] * n[2];
- aL = sqrtf(hydro_gamma * WL[4] / WL[0]);
- aR = sqrtf(hydro_gamma * WR[4] / WR[0]);
+ const float uL = WL[1] * n[0] + WL[2] * n[1] + WL[3] * n[2];
+ const float uR = WR[1] * n[0] + WR[2] * n[1] + WR[3] * n[2];
+ const float aL = sqrtf(hydro_gamma * WL[4] / WL[0]);
+ const float aR = sqrtf(hydro_gamma * WR[4] / WR[0]);
/* Handle vacuum: vacuum does not require iteration and is always exact */
if (riemann_is_vacuum(WL, WR, uL, uR, aL, aR)) {
@@ -68,30 +70,33 @@ __attribute__((always_inline)) INLINE static void riemann_solve_for_flux(
}
/* STEP 1: pressure estimate */
- rhobar = 0.5 * (WL[0] + WR[0]);
- abar = 0.5 * (aL + aR);
- pPVRS = 0.5 * (WL[4] + WR[4]) - 0.5 * (uR - uL) * rhobar * abar;
- pstar = max(0., pPVRS);
+ const float rhobar = 0.5f * (WL[0] + WR[0]);
+ const float abar = 0.5f * (aL + aR);
+ const float pPVRS = 0.5f * (WL[4] + WR[4]) - 0.5f * (uR - uL) * rhobar * abar;
+ const float pstar = max(0.f, pPVRS);
/* STEP 2: wave speed estimates
all these speeds are along the interface normal, since uL and uR are */
- qL = 1.;
- if (pstar > WL[4]) {
- qL = sqrtf(1. +
- 0.5 * (hydro_gamma + 1.) / hydro_gamma * (pstar / WL[4] - 1.));
+ float qL = 1.f;
+ if (pstar > WL[4] && WL[4] > 0.f) {
+ qL = sqrtf(1.f +
+ 0.5f * hydro_gamma_plus_one * hydro_one_over_gamma *
+ (pstar / WL[4] - 1.f));
}
- qR = 1.;
- if (pstar > WR[4]) {
- qR = sqrtf(1. +
- 0.5 * (hydro_gamma + 1.) / hydro_gamma * (pstar / WR[4] - 1.));
+ float qR = 1.f;
+ if (pstar > WR[4] && WR[4] > 0.f) {
+ qR = sqrtf(1.f +
+ 0.5f * hydro_gamma_plus_one * hydro_one_over_gamma *
+ (pstar / WR[4] - 1.f));
}
- SL = uL - aL * qL;
- SR = uR + aR * qR;
- Sstar = (WR[4] - WL[4] + WL[0] * uL * (SL - uL) - WR[0] * uR * (SR - uR)) /
- (WL[0] * (SL - uL) - WR[0] * (SR - uR));
+ const float SL = uL - aL * qL;
+ const float SR = uR + aR * qR;
+ const float Sstar =
+ (WR[4] - WL[4] + WL[0] * uL * (SL - uL) - WR[0] * uR * (SR - uR)) /
+ (WL[0] * (SL - uL) - WR[0] * (SR - uR));
/* STEP 3: HLLC flux in a frame moving with the interface velocity */
- if (Sstar >= 0.) {
+ if (Sstar >= 0.f) {
/* flux FL */
totflux[0] = WL[0] * uL;
/* these are the actual correct fluxes in the boosted lab frame
@@ -99,15 +104,18 @@ __attribute__((always_inline)) INLINE static void riemann_solve_for_flux(
totflux[1] = WL[0] * WL[1] * uL + WL[4] * n[0];
totflux[2] = WL[0] * WL[2] * uL + WL[4] * n[1];
totflux[3] = WL[0] * WL[3] * uL + WL[4] * n[2];
- v2 = WL[1] * WL[1] + WL[2] * WL[2] + WL[3] * WL[3];
- eL = WL[4] / hydro_gamma_minus_one / WL[0] + 0.5 * v2;
+ const float v2 = WL[1] * WL[1] + WL[2] * WL[2] + WL[3] * WL[3];
+ const float eL = WL[4] * hydro_one_over_gamma_minus_one / WL[0] + 0.5f * v2;
totflux[4] = WL[0] * eL * uL + WL[4] * uL;
- if (SL < 0.) {
+ if (SL < 0.f) {
+
+ float UstarL[5];
+
/* add flux FstarL */
- UstarL[0] = 1.;
+ UstarL[0] = 1.f;
/* we need UstarL in the lab frame:
- subtract the velocity in the interface frame from the lab frame
- velocity and then add Sstar in interface frame */
+ * subtract the velocity in the interface frame from the lab frame
+ * velocity and then add Sstar in interface frame */
UstarL[1] = WL[1] + (Sstar - uL) * n[0];
UstarL[2] = WL[2] + (Sstar - uL) * n[1];
UstarL[3] = WL[3] + (Sstar - uL) * n[2];
@@ -129,12 +137,19 @@ __attribute__((always_inline)) INLINE static void riemann_solve_for_flux(
totflux[1] = WR[0] * WR[1] * uR + WR[4] * n[0];
totflux[2] = WR[0] * WR[2] * uR + WR[4] * n[1];
totflux[3] = WR[0] * WR[3] * uR + WR[4] * n[2];
- v2 = WR[1] * WR[1] + WR[2] * WR[2] + WR[3] * WR[3];
- eR = WR[4] / hydro_gamma_minus_one / WR[0] + 0.5 * v2;
+ const float v2 = WR[1] * WR[1] + WR[2] * WR[2] + WR[3] * WR[3];
+ const float eR = WR[4] * hydro_one_over_gamma_minus_one / WR[0] + 0.5f * v2;
totflux[4] = WR[0] * eR * uR + WR[4] * uR;
- if (SR > 0.) {
+ if (SR > 0.f) {
+
+ float UstarR[5];
+
/* add flux FstarR */
- UstarR[0] = 1.;
+ UstarR[0] = 1.f;
+
+ /* we need UstarR in the lab frame:
+ * subtract the velocity in the interface frame from the lab frame
+ * velocity and then add Sstar in interface frame */
UstarR[1] = WR[1] + (Sstar - uR) * n[0];
UstarR[2] = WR[2] + (Sstar - uR) * n[1];
UstarR[3] = WR[3] + (Sstar - uR) * n[2];
@@ -152,19 +167,23 @@ __attribute__((always_inline)) INLINE static void riemann_solve_for_flux(
}
}
- /* deboost to lab frame
- we add the flux contribution due to the movement of the interface
- the density flux is unchanged
+ /* deboost to lab frame we add the flux contribution due to the
+ movement of the interface the density flux is unchanged
we add the extra velocity flux due to the absolute motion of the fluid
similarly, we need to add the energy fluxes due to the absolute motion */
- v2 = vij[0] * vij[0] + vij[1] * vij[1] + vij[2] * vij[2];
- // order is important: we first use the momentum fluxes to update the energy
- // flux and then de-boost the momentum fluxes!
+ const float v2 = vij[0] * vij[0] + vij[1] * vij[1] + vij[2] * vij[2];
+
+ /* order is important: we first use the momentum fluxes to update the energy
+ flux and then de-boost the momentum fluxes! */
totflux[4] += vij[0] * totflux[1] + vij[1] * totflux[2] +
- vij[2] * totflux[3] + 0.5 * v2 * totflux[0];
+ vij[2] * totflux[3] + 0.5f * v2 * totflux[0];
totflux[1] += vij[0] * totflux[0];
totflux[2] += vij[1] * totflux[0];
totflux[3] += vij[2] * totflux[0];
+
+#ifdef SWIFT_DEBUG_CHECKS
+ riemann_check_output(WL, WR, n, vij, totflux);
+#endif
}
#endif /* SWIFT_RIEMANN_HLLC_H */
diff --git a/src/riemann/riemann_trrs.h b/src/riemann/riemann_trrs.h
index 44fbe6ffb8fb1f557045b6e466c87e7d8b591aab..597f6f8edd9690dbdd34a3c9885ae228c2af734c 100644
--- a/src/riemann/riemann_trrs.h
+++ b/src/riemann/riemann_trrs.h
@@ -1,6 +1,7 @@
/*******************************************************************************
* This file is part of SWIFT.
- * Coypright (c) 2015 Bert Vandenbroucke (bert.vandenbroucke@ugent.be)
+ * Copyright (c) 2015 Bert Vandenbroucke (bert.vandenbroucke@ugent.be)
+ * 2018 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as published
@@ -29,6 +30,7 @@
#include "adiabatic_index.h"
#include "error.h"
#include "minmax.h"
+#include "riemann_checks.h"
#include "riemann_vacuum.h"
#ifndef EOS_IDEAL_GAS
@@ -52,7 +54,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 +162,12 @@ __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) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+ riemann_check_input(Wi, Wj, n_unit, vij);
+#endif
float Whalf[5];
float flux[5][3];
@@ -206,6 +213,10 @@ __attribute__((always_inline)) INLINE static void riemann_solve_for_flux(
flux[3][0] * n_unit[0] + flux[3][1] * n_unit[1] + flux[3][2] * n_unit[2];
totflux[4] =
flux[4][0] * n_unit[0] + flux[4][1] * n_unit[1] + flux[4][2] * n_unit[2];
+
+#ifdef SWIFT_DEBUG_CHECKS
+ riemann_check_output(Wi, Wj, n_unit, vij, totflux);
+#endif
}
#endif /* SWIFT_RIEMANN_TRRS_H */
diff --git a/src/riemann/riemann_vacuum.h b/src/riemann/riemann_vacuum.h
index 743abb910193380793ccdf3d7eddbcedc4968691..c1144a2557889df20ece7dea6918af804ee03232 100644
--- a/src/riemann/riemann_vacuum.h
+++ b/src/riemann/riemann_vacuum.h
@@ -24,20 +24,20 @@
* @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;
- }
+ if (!WL[0] || !WR[0]) return 1;
+
/* vacuum generation */
- if (2.0f * aL / hydro_gamma_minus_one + 2.0f * aR / hydro_gamma_minus_one <=
- vR - vL) {
+ else if (hydro_two_over_gamma_minus_one * aL +
+ hydro_two_over_gamma_minus_one * aR <=
+ vR - vL)
return 1;
- }
/* no vacuum */
- return 0;
+ else
+ return 0;
}
/**
@@ -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];
diff --git a/src/sign.h b/src/sign.h
new file mode 100644
index 0000000000000000000000000000000000000000..04ecdb3a99c5baf35e40d6f8286d25102d22e5fd
--- /dev/null
+++ b/src/sign.h
@@ -0,0 +1,46 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2018 Matthieu Schaller (matthieu.schaller@durham.ac.uk).
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_SIGN_H
+#define SWIFT_SIGN_H
+
+/**
+ * @brief Return the sign of a floating point number
+ *
+ * @param x The number of interest.
+ * @return >0 if positive, 0 if negative.
+ */
+__attribute__((always_inline)) INLINE static int signf(float x) {
+#ifdef __GNUC__
+ return !signbit(x);
+#else
+ return (0.f < val) - (val < 0.f);
+#endif
+}
+
+/**
+ * @brief Return 1 if two numbers have the same sign, 0 otherwise
+ *
+ * @param x The first number
+ * @param y The second number
+ */
+__attribute__((always_inline)) INLINE static int same_signf(float x, float y) {
+ return signf(x) == signf(y);
+}
+
+#endif