diff --git a/src/riemann/riemann_exact.h b/src/riemann/riemann_exact.h
index caff43fb485f271867c27d204829e09a4d9f6b8e..9763d9f0d12da32d9142c24481105b9f139be588 100644
--- a/src/riemann/riemann_exact.h
+++ b/src/riemann/riemann_exact.h
@@ -349,164 +349,152 @@ __attribute__((always_inline)) INLINE static void riemann_solver_solve(
aL = sqrtf(hydro_gamma * WL[4] / WL[0]);
aR = sqrtf(hydro_gamma * WR[4] / WR[0]);
- if (!WL[0] || !WR[0]) {
- /* vacuum: we need a vacuum riemann solver */
+ /* check vacuum (generation) condition */
+ if (riemann_is_vacuum(WL, WR, vL, vR, aL, aR)) {
riemann_solve_vacuum(WL, WR, vL, vR, aL, aR, Whalf, n_unit);
return;
}
- /* check vacuum generation condition */
- if (2.0f * aL / hydro_gamma_minus_one + 2.0f * aR / hydro_gamma_minus_one <
- fabs(vL - vR)) {
- /* vacuum generation: need a vacuum riemann solver */
- riemann_solve_vacuum(WL, WR, vL, vR, aL, aR, Whalf, n_unit);
- return;
- } else {
- /* values are ok: let's find pstar (riemann_f(pstar) = 0)! */
- /* We normally use a Newton-Raphson iteration to find the zeropoint
- of riemann_f(p), but if pstar is close to 0, we risk negative p values.
- Since riemann_f(p) is undefined for negative pressures, we don't
- want this to happen.
- We therefore use Brent's method if riemann_f(0) is larger than some
- value. -5 makes the iteration fail safe while almost never invoking
- the expensive Brent solver. */
- p = 0.;
- /* obtain a first guess for p */
- pguess = riemann_guess_p(WL, WR, vL, vR, aL, aR);
- fp = riemann_f(p, WL, WR, vL, vR, aL, aR);
- fpguess = riemann_f(pguess, WL, WR, vL, vR, aL, aR);
- /* ok, pstar is close to 0, better use Brent's method... */
- /* we use Newton-Raphson until we find a suitable interval */
- if (fp * fpguess >= 0.0f) {
- /* Newton-Raphson until convergence or until suitable interval is found
- to use Brent's method */
- unsigned int counter = 0;
- while (fabs(p - pguess) > 1.e-6f * 0.5f * (p + pguess) &&
- fpguess < 0.0f) {
- p = pguess;
- pguess = pguess - fpguess / riemann_fprime(pguess, WL, WR, aL, aR);
- fpguess = riemann_f(pguess, WL, WR, vL, vR, aL, aR);
- counter++;
- if (counter > 1000) {
- error("Stuck in Newton-Raphson!\n");
- }
- }
- }
- /* As soon as there is a suitable interval: use Brent's method */
- if (1.e6 * fabs(p - pguess) > 0.5f * (p + pguess) && fpguess > 0.0f) {
- p = 0.0f;
- fp = riemann_f(p, WL, WR, vL, vR, aL, aR);
- /* use Brent's method to find the zeropoint */
- p = riemann_solve_brent(p, pguess, fp, fpguess, 1.e-6, WL, WR, vL, vR, aL,
- aR);
- } else {
+ /* values are ok: let's find pstar (riemann_f(pstar) = 0)! */
+ /* We normally use a Newton-Raphson iteration to find the zeropoint
+ of riemann_f(p), but if pstar is close to 0, we risk negative p values.
+ Since riemann_f(p) is undefined for negative pressures, we don't
+ want this to happen.
+ We therefore use Brent's method if riemann_f(0) is larger than some
+ value. -5 makes the iteration fail safe while almost never invoking
+ the expensive Brent solver. */
+ p = 0.;
+ /* obtain a first guess for p */
+ pguess = riemann_guess_p(WL, WR, vL, vR, aL, aR);
+ fp = riemann_f(p, WL, WR, vL, vR, aL, aR);
+ fpguess = riemann_f(pguess, WL, WR, vL, vR, aL, aR);
+ /* ok, pstar is close to 0, better use Brent's method... */
+ /* we use Newton-Raphson until we find a suitable interval */
+ if (fp * fpguess >= 0.0f) {
+ /* Newton-Raphson until convergence or until suitable interval is found
+ to use Brent's method */
+ unsigned int counter = 0;
+ while (fabs(p - pguess) > 1.e-6f * 0.5f * (p + pguess) && fpguess < 0.0f) {
p = pguess;
+ pguess = pguess - fpguess / riemann_fprime(pguess, WL, WR, aL, aR);
+ fpguess = riemann_f(pguess, WL, WR, vL, vR, aL, aR);
+ counter++;
+ if (counter > 1000) {
+ error("Stuck in Newton-Raphson!\n");
+ }
}
+ }
+ /* As soon as there is a suitable interval: use Brent's method */
+ if (1.e6 * fabs(p - pguess) > 0.5f * (p + pguess) && fpguess > 0.0f) {
+ p = 0.0f;
+ fp = riemann_f(p, WL, WR, vL, vR, aL, aR);
+ /* use Brent's method to find the zeropoint */
+ p = riemann_solve_brent(p, pguess, fp, fpguess, 1.e-6, WL, WR, vL, vR, aL,
+ aR);
+ } else {
+ p = pguess;
+ }
- /* calculate the velocity in the intermediate state */
- u = 0.5f * (vL + vR) +
- 0.5f * (riemann_fb(p, WR, aR) - riemann_fb(p, WL, aL));
-
- /* sample the solution */
- /* This corresponds to the flow chart in Fig. 4.14 in Toro */
- if (u < 0.0f) {
- /* advect velocity components */
- Whalf[1] = WR[1];
- Whalf[2] = WR[2];
- Whalf[3] = WR[3];
- pdpR = p / WR[4];
- if (p > WR[4]) {
- /* shockwave */
- SR = vR +
- aR * sqrtf(hydro_gamma_plus_one_over_two_gamma * pdpR +
- hydro_gamma_minus_one_over_two_gamma);
- if (SR > 0.0f) {
- Whalf[0] = WR[0] *
- (pdpR + hydro_gamma_minus_one_over_gamma_plus_one) /
- (hydro_gamma_minus_one_over_gamma_plus_one * pdpR + 1.0f);
+ /* calculate the velocity in the intermediate state */
+ u = 0.5f * (vL + vR) + 0.5f * (riemann_fb(p, WR, aR) - riemann_fb(p, WL, aL));
+
+ /* sample the solution */
+ /* This corresponds to the flow chart in Fig. 4.14 in Toro */
+ if (u < 0.0f) {
+ /* advect velocity components */
+ Whalf[1] = WR[1];
+ Whalf[2] = WR[2];
+ Whalf[3] = WR[3];
+ pdpR = p / WR[4];
+ if (p > WR[4]) {
+ /* shockwave */
+ SR = vR +
+ aR * sqrtf(hydro_gamma_plus_one_over_two_gamma * pdpR +
+ hydro_gamma_minus_one_over_two_gamma);
+ if (SR > 0.0f) {
+ Whalf[0] = WR[0] * (pdpR + hydro_gamma_minus_one_over_gamma_plus_one) /
+ (hydro_gamma_minus_one_over_gamma_plus_one * pdpR + 1.0f);
+ vhalf = u - vR;
+ Whalf[4] = p;
+ } else {
+ Whalf[0] = WR[0];
+ vhalf = 0.0f;
+ Whalf[4] = WR[4];
+ }
+ } else {
+ /* rarefaction wave */
+ SHR = vR + aR;
+ if (SHR > 0.0f) {
+ STR = u + aR * pow_gamma_minus_one_over_two_gamma(pdpR);
+ if (STR <= 0.0f) {
+ Whalf[0] =
+ WR[0] * pow_two_over_gamma_minus_one(
+ hydro_two_over_gamma_plus_one -
+ hydro_gamma_minus_one_over_gamma_plus_one / aR * vR);
+ vhalf = hydro_two_over_gamma_plus_one *
+ (-aR + hydro_gamma_minus_one_over_two * vR) -
+ vR;
+ Whalf[4] =
+ WR[4] * pow_two_gamma_over_gamma_minus_one(
+ hydro_two_over_gamma_plus_one -
+ hydro_gamma_minus_one_over_gamma_plus_one / aR * vR);
+ } else {
+ Whalf[0] = WR[0] * pow_one_over_gamma(pdpR);
vhalf = u - vR;
Whalf[4] = p;
- } else {
- Whalf[0] = WR[0];
- vhalf = 0.0f;
- Whalf[4] = WR[4];
}
} else {
- /* rarefaction wave */
- SHR = vR + aR;
- if (SHR > 0.0f) {
- STR = u + aR * pow_gamma_minus_one_over_two_gamma(pdpR);
- if (STR <= 0.0f) {
- Whalf[0] = WR[0] *
- pow_two_over_gamma_minus_one(
- hydro_two_over_gamma_plus_one -
- hydro_gamma_minus_one_over_gamma_plus_one / aR * vR);
- vhalf = hydro_two_over_gamma_plus_one *
- (-aR + hydro_gamma_minus_one_over_two * vR) -
- vR;
- Whalf[4] = WR[4] *
- pow_two_gamma_over_gamma_minus_one(
- hydro_two_over_gamma_plus_one -
- hydro_gamma_minus_one_over_gamma_plus_one / aR * vR);
- } else {
- Whalf[0] = WR[0] * pow_one_over_gamma(pdpR);
- vhalf = u - vR;
- Whalf[4] = p;
- }
- } else {
- Whalf[0] = WR[0];
- vhalf = 0.0f;
- Whalf[4] = WR[4];
- }
+ Whalf[0] = WR[0];
+ vhalf = 0.0f;
+ Whalf[4] = WR[4];
+ }
+ }
+ } else {
+ Whalf[1] = WL[1];
+ Whalf[2] = WL[2];
+ Whalf[3] = WL[3];
+ pdpL = p / WL[4];
+ if (p > WL[4]) {
+ /* shockwave */
+ SL = vL -
+ aL * sqrtf(hydro_gamma_plus_one_over_two_gamma * pdpL +
+ hydro_gamma_minus_one_over_two_gamma);
+ if (SL < 0.0f) {
+ Whalf[0] = WL[0] * (pdpL + hydro_gamma_minus_one_over_gamma_plus_one) /
+ (hydro_gamma_minus_one_over_gamma_plus_one * pdpL + 1.0f);
+ vhalf = u - vL;
+ Whalf[4] = p;
+ } else {
+ Whalf[0] = WL[0];
+ vhalf = 0.0f;
+ Whalf[4] = WL[4];
}
} else {
- Whalf[1] = WL[1];
- Whalf[2] = WL[2];
- Whalf[3] = WL[3];
- pdpL = p / WL[4];
- if (p > WL[4]) {
- /* shockwave */
- SL = vL -
- aL * sqrtf(hydro_gamma_plus_one_over_two_gamma * pdpL +
- hydro_gamma_minus_one_over_two_gamma);
- if (SL < 0.0f) {
- Whalf[0] = WL[0] *
- (pdpL + hydro_gamma_minus_one_over_gamma_plus_one) /
- (hydro_gamma_minus_one_over_gamma_plus_one * pdpL + 1.0f);
+ /* rarefaction wave */
+ SHL = vL - aL;
+ if (SHL < 0.0f) {
+ STL = u - aL * pow_gamma_minus_one_over_two_gamma(pdpL);
+ if (STL > 0.0f) {
+ Whalf[0] =
+ WL[0] * pow_two_over_gamma_minus_one(
+ hydro_two_over_gamma_plus_one +
+ hydro_gamma_minus_one_over_gamma_plus_one / aL * vL);
+ vhalf = hydro_two_over_gamma_plus_one *
+ (aL + hydro_gamma_minus_one_over_two * vL) -
+ vL;
+ Whalf[4] =
+ WL[4] * pow_two_gamma_over_gamma_minus_one(
+ hydro_two_over_gamma_plus_one +
+ hydro_gamma_minus_one_over_gamma_plus_one / aL * vL);
+ } else {
+ Whalf[0] = WL[0] * pow_one_over_gamma(pdpL);
vhalf = u - vL;
Whalf[4] = p;
- } else {
- Whalf[0] = WL[0];
- vhalf = 0.0f;
- Whalf[4] = WL[4];
}
} else {
- /* rarefaction wave */
- SHL = vL - aL;
- if (SHL < 0.0f) {
- STL = u - aL * pow_gamma_minus_one_over_two_gamma(pdpL);
- if (STL > 0.0f) {
- Whalf[0] = WL[0] *
- pow_two_over_gamma_minus_one(
- hydro_two_over_gamma_plus_one +
- hydro_gamma_minus_one_over_gamma_plus_one / aL * vL);
- vhalf = hydro_two_over_gamma_plus_one *
- (aL + hydro_gamma_minus_one_over_two * vL) -
- vL;
- Whalf[4] = WL[4] *
- pow_two_gamma_over_gamma_minus_one(
- hydro_two_over_gamma_plus_one +
- hydro_gamma_minus_one_over_gamma_plus_one / aL * vL);
- } else {
- Whalf[0] = WL[0] * pow_one_over_gamma(pdpL);
- vhalf = u - vL;
- Whalf[4] = p;
- }
- } else {
- Whalf[0] = WL[0];
- vhalf = 0.0f;
- Whalf[4] = WL[4];
- }
+ Whalf[0] = WL[0];
+ vhalf = 0.0f;
+ Whalf[4] = WL[4];
}
}
}
diff --git a/src/riemann/riemann_hllc.h b/src/riemann/riemann_hllc.h
index feb144589d52d9fb5202e115e7c25028454906e7..fdc22ce05b8d63bdba66e530d1a5a968801a9f10 100644
--- a/src/riemann/riemann_hllc.h
+++ b/src/riemann/riemann_hllc.h
@@ -21,6 +21,7 @@
#define SWIFT_RIEMANN_HLLC_H
#include "adiabatic_index.h"
+#include "riemann_vacuum.h"
__attribute__((always_inline)) INLINE static void riemann_solve_for_flux(
float *WL, float *WR, float *n, float *vij, float *totflux) {
@@ -47,12 +48,9 @@ __attribute__((always_inline)) INLINE static void riemann_solve_for_flux(
aR = sqrtf(hydro_gamma * WR[4] / WR[0]);
/* Handle vacuum: vacuum does not require iteration and is always exact */
- if (!WL[0] || !WR[0]) {
- error("Vacuum not yet supported");
- }
- if (2. * aL / hydro_gamma_minus_one + 2. * aR / hydro_gamma_minus_one <
- fabs(uL - uR)) {
- error("Vacuum not yet supported");
+ if (riemann_is_vacuum(WL, WR, uL, uR, aL, aR)) {
+ riemann_solve_vacuum_flux(WL, WR, uL, uR, aL, aR, n, vij, totflux);
+ return;
}
/* STEP 1: pressure estimate */
diff --git a/src/riemann/riemann_trrs.h b/src/riemann/riemann_trrs.h
index 7a4ac0233f7b7c72850ad441d98ffbebc7687177..b13a76b4c57af548497780e974e5c9ee3a721fac 100644
--- a/src/riemann/riemann_trrs.h
+++ b/src/riemann/riemann_trrs.h
@@ -56,8 +56,9 @@ __attribute__((always_inline)) INLINE static void riemann_solver_solve(
aL = sqrtf(hydro_gamma * WL[4] / WL[0]);
aR = sqrtf(hydro_gamma * WR[4] / WR[0]);
- if (riemann_is_vacuum(WL, WR, vL, vR, aL, aR, Whalf, n_unit)) {
+ if (riemann_is_vacuum(WL, WR, vL, vR, aL, aR)) {
riemann_solve_vacuum(WL, WR, vL, vR, aL, aR, Whalf, n_unit);
+ return;
}
/* calculate the velocity and pressure in the intermediate state */
diff --git a/src/riemann/riemann_vacuum.h b/src/riemann/riemann_vacuum.h
index 9630f312a621f4b54d3dcda8e11382b09e20fa82..8ceac536a3c6160835759c157841ac51b9fff1a1 100644
--- a/src/riemann/riemann_vacuum.h
+++ b/src/riemann/riemann_vacuum.h
@@ -24,16 +24,15 @@
* @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, float* Whalf,
- float* n_unit) {
+ float* WL, float* WR, float vL, float vR, float aL, float aR) {
/* vacuum */
if (!WL[0] || !WR[0]) {
return 1;
}
/* vacuum generation */
- if (2.0f * aL / hydro_gamma_minus_one + 2.0f * aR / hydro_gamma_minus_one <
- fabs(vL - vR)) {
+ if (2.0f * aL / hydro_gamma_minus_one + 2.0f * aR / hydro_gamma_minus_one <=
+ vR - vL) {
return 1;
}
@@ -60,6 +59,12 @@ __attribute__((always_inline)) INLINE static void riemann_solve_vacuum(
float SL, SR;
float vhalf;
+ message(
+ "WL: [%g, %g, %g, %g, %g]\nWR: [%g, %g, %g, %g, %g]\nvL: %g, vR: %g"
+ ", aL: %g, aR: %g",
+ WL[0], WL[1], WL[2], WL[3], WL[4], WR[0], WR[1], WR[2], WR[3], WR[4], vL,
+ vR, aL, aR);
+
if (!WR[0] && !WL[0]) {
/* if both states are vacuum, the solution is also vacuum */
Whalf[0] = 0.0f;
@@ -135,10 +140,12 @@ __attribute__((always_inline)) INLINE static void riemann_solve_vacuum(
Whalf[4] = WR[4];
}
} else {
+ message("Vacuum generation");
/* vacuum generation */
SR = vR - hydro_two_over_gamma_minus_one * aR;
SL = vL + hydro_two_over_gamma_minus_one * aL;
if (SR > 0.0f && SL < 0.0f) {
+ message("Vacuum");
Whalf[0] = 0.0f;
Whalf[1] = 0.0f;
Whalf[2] = 0.0f;
@@ -151,6 +158,7 @@ __attribute__((always_inline)) INLINE static void riemann_solve_vacuum(
Whalf[2] = WL[2];
Whalf[3] = WL[3];
if (aL > vL) {
+ message("Left fan");
Whalf[0] = WL[0] *
pow_two_over_gamma_minus_one(
hydro_two_over_gamma_plus_one +
@@ -163,6 +171,7 @@ __attribute__((always_inline)) INLINE static void riemann_solve_vacuum(
hydro_two_over_gamma_plus_one +
hydro_gamma_minus_one_over_gamma_plus_one / aL * vL);
} else {
+ message("Left state");
Whalf[0] = WL[0];
vhalf = 0.0f;
Whalf[4] = WL[4];
@@ -172,6 +181,7 @@ __attribute__((always_inline)) INLINE static void riemann_solve_vacuum(
Whalf[2] = WR[2];
Whalf[3] = WR[3];
if (-aR < vR) {
+ message("Right fan");
Whalf[0] = WR[0] *
pow_two_over_gamma_minus_one(
hydro_two_over_gamma_plus_one -
@@ -184,6 +194,7 @@ __attribute__((always_inline)) INLINE static void riemann_solve_vacuum(
hydro_two_over_gamma_plus_one -
hydro_gamma_minus_one_over_gamma_plus_one / aR * vR);
} else {
+ message("Right state");
Whalf[0] = WR[0];
vhalf = 0.0f;
Whalf[4] = WR[4];
@@ -199,4 +210,57 @@ __attribute__((always_inline)) INLINE static void riemann_solve_vacuum(
Whalf[3] += vhalf * n_unit[2];
}
+/**
+ * @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) {
+
+ float Whalf[5];
+ float flux[5][3];
+ float vtot[3];
+ float rhoe;
+
+ riemann_solve_vacuum(WL, WR, vL, vR, aL, aR, 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_VACUUM_H */
diff --git a/tests/testRiemannExact.c b/tests/testRiemannExact.c
index 24ed3d08a4ba17768fa57bff4d295fbd27768e81..0c8fba8424f212a827f987ea989a109f1ec752c5 100644
--- a/tests/testRiemannExact.c
+++ b/tests/testRiemannExact.c
@@ -22,6 +22,22 @@
#include "riemann/riemann_exact.h"
#include "tools.h"
+int opposite(float a, float b) {
+ if ((a - b)) {
+ return fabs((a + b) / (a - b)) < 1.e-4;
+ } else {
+ return a == 0.0f;
+ }
+}
+
+int equal(float a, float b) {
+ if ((a + b)) {
+ return fabs((a - b) / (a + b)) < 1.e-4;
+ } else {
+ return a == 0.0f;
+ }
+}
+
/**
* @brief Check that a and b are consistent (up to some error)
*
@@ -248,7 +264,9 @@ void check_riemann_symmetry() {
riemann_solver_solve(WL, WR, Whalf1, n_unit1);
riemann_solver_solve(WR, WL, Whalf2, n_unit2);
- if (memcmp(Whalf1, Whalf2, 5 * sizeof(float))) {
+ if (!equal(Whalf1[0], Whalf2[0]) || !equal(Whalf1[1], Whalf2[1]) ||
+ !equal(Whalf1[2], Whalf2[2]) || !equal(Whalf1[3], Whalf2[3]) ||
+ !equal(Whalf1[4], Whalf2[4])) {
message(
"Solver asymmetric: [%.3e,%.3e,%.3e,%.3e,%.3e] == "
"[%.3e,%.3e,%.3e,%.3e,%.3e]\n",
@@ -277,14 +295,11 @@ void check_riemann_symmetry() {
riemann_solve_for_flux(WL, WR, n_unit1, vij, totflux1);
riemann_solve_for_flux(WR, WL, n_unit2, vij, totflux2);
- /* we expect the fluxes to have a different sign, so we reverse one of them */
- totflux2[0] = -totflux2[0];
- totflux2[1] = -totflux2[1];
- totflux2[2] = -totflux2[2];
- totflux2[3] = -totflux2[3];
- totflux2[4] = -totflux2[4];
-
- if (memcmp(totflux1, totflux2, 5 * sizeof(float))) {
+ if (!opposite(totflux1[0], totflux2[0]) ||
+ !opposite(totflux1[1], totflux2[1]) ||
+ !opposite(totflux1[2], totflux2[2]) ||
+ !opposite(totflux1[3], totflux2[3]) ||
+ !opposite(totflux1[4], totflux2[4])) {
message(
"Flux solver asymmetric: [%.3e,%.3e,%.3e,%.3e,%.3e] == "
"[%.3e,%.3e,%.3e,%.3e,%.3e]\n",
diff --git a/tests/testRiemannHLLC.c b/tests/testRiemannHLLC.c
index bf95aa4b96a46d8d772d1d0f5d34b738cedf2b38..4cf883b68efbcfd795d0b7894adb9e7265b14d14 100644
--- a/tests/testRiemannHLLC.c
+++ b/tests/testRiemannHLLC.c
@@ -22,6 +22,8 @@
#include "riemann/riemann_hllc.h"
#include "tools.h"
+int consistent_with_zero(float val) { return fabs(val) < 1.e-4; }
+
/**
* @brief Check the symmetry of the TRRS Riemann solver
*/
@@ -61,14 +63,11 @@ void check_riemann_symmetry() {
riemann_solve_for_flux(WL, WR, n_unit1, vij, totflux1);
riemann_solve_for_flux(WR, WL, n_unit2, vij, totflux2);
- /* we expect the fluxes to have a different sign, so we reverse one of them */
- totflux2[0] = -totflux2[0];
- totflux2[1] = -totflux2[1];
- totflux2[2] = -totflux2[2];
- totflux2[3] = -totflux2[3];
- totflux2[4] = -totflux2[4];
-
- if (memcmp(totflux1, totflux2, 5 * sizeof(float))) {
+ if (!consistent_with_zero(totflux1[0] + totflux2[0]) ||
+ !consistent_with_zero(totflux1[1] + totflux2[1]) ||
+ !consistent_with_zero(totflux1[2] + totflux2[2]) ||
+ !consistent_with_zero(totflux1[3] + totflux2[3]) ||
+ !consistent_with_zero(totflux1[4] + totflux2[4])) {
message(
"Flux solver asymmetric: [%.3e,%.3e,%.3e,%.3e,%.3e] == "
"[%.3e,%.3e,%.3e,%.3e,%.3e]\n",
diff --git a/tests/testRiemannTRRS.c b/tests/testRiemannTRRS.c
index f284585fe0d093f94209c829e719d3742e24d97e..18ecbdce9173f43674a63b21231322cb01620d29 100644
--- a/tests/testRiemannTRRS.c
+++ b/tests/testRiemannTRRS.c
@@ -22,6 +22,22 @@
#include "riemann/riemann_trrs.h"
#include "tools.h"
+int opposite(float a, float b) {
+ if ((a - b)) {
+ return fabs((a + b) / (a - b)) < 1.e-4;
+ } else {
+ return a == 0.0f;
+ }
+}
+
+int equal(float a, float b) {
+ if ((a + b)) {
+ return fabs((a - b) / (a + b)) < 1.e-4;
+ } else {
+ return a == 0.0f;
+ }
+}
+
/**
* @brief Check that a and b are consistent (up to some error)
*
@@ -248,7 +264,9 @@ void check_riemann_symmetry() {
riemann_solver_solve(WL, WR, Whalf1, n_unit1);
riemann_solver_solve(WR, WL, Whalf2, n_unit2);
- if (memcmp(Whalf1, Whalf2, 5 * sizeof(float))) {
+ if (!equal(Whalf1[0], Whalf2[0]) || !equal(Whalf1[1], Whalf2[1]) ||
+ !equal(Whalf1[2], Whalf2[2]) || !equal(Whalf1[3], Whalf2[3]) ||
+ !equal(Whalf1[4], Whalf2[4])) {
message(
"Solver asymmetric: [%.3e,%.3e,%.3e,%.3e,%.3e] == "
"[%.3e,%.3e,%.3e,%.3e,%.3e]\n",
@@ -270,14 +288,11 @@ void check_riemann_symmetry() {
riemann_solve_for_flux(WL, WR, n_unit1, vij, totflux1);
riemann_solve_for_flux(WR, WL, n_unit2, vij, totflux2);
- /* we expect the fluxes to have a different sign, so we reverse one of them */
- totflux2[0] = -totflux2[0];
- totflux2[1] = -totflux2[1];
- totflux2[2] = -totflux2[2];
- totflux2[3] = -totflux2[3];
- totflux2[4] = -totflux2[4];
-
- if (memcmp(totflux1, totflux2, 5 * sizeof(float))) {
+ if (!opposite(totflux1[0], totflux2[0]) ||
+ !opposite(totflux1[1], totflux2[1]) ||
+ !opposite(totflux1[2], totflux2[2]) ||
+ !opposite(totflux1[3], totflux2[3]) ||
+ !opposite(totflux1[4], totflux2[4])) {
message(
"Solver asymmetric: [%.3e,%.3e,%.3e,%.3e,%.3e] == "
"[%.3e,%.3e,%.3e,%.3e,%.3e]\n",
diff --git a/tests/testSymmetry.c b/tests/testSymmetry.c
index f03bc68eac9d05e1e2c3b07b221bd0b5172e6838..9216602bc166f8d0f7d449f1a8c118878e37e76a 100644
--- a/tests/testSymmetry.c
+++ b/tests/testSymmetry.c
@@ -54,18 +54,11 @@ int main(int argc, char *argv[]) {
pi.primitives.v[1] = random_uniform(-10.0f, 10.0f);
pi.primitives.v[2] = random_uniform(-10.0f, 10.0f);
pi.primitives.P = random_uniform(0.1f, 1.0f);
- /* pj.primitives.rho = random_uniform(0.1f, 1.0f);
- pj.primitives.v[0] = random_uniform(-10.0f, 10.0f);
- pj.primitives.v[1] = random_uniform(-10.0f, 10.0f);
- pj.primitives.v[2] = random_uniform(-10.0f, 10.0f);
- pj.primitives.P = random_uniform(0.1f, 1.0f);*/
- /* make the values for pj the same, since otherwise we suffer from the swap of
- left and right states in the Riemann solver */
- pj.primitives.rho = pi.primitives.rho;
- pj.primitives.v[0] = pi.primitives.v[0];
- pj.primitives.v[1] = pi.primitives.v[1];
- pj.primitives.v[2] = pi.primitives.v[2];
- pj.primitives.P = pi.primitives.P;
+ pj.primitives.rho = random_uniform(0.1f, 1.0f);
+ pj.primitives.v[0] = random_uniform(-10.0f, 10.0f);
+ pj.primitives.v[1] = random_uniform(-10.0f, 10.0f);
+ pj.primitives.v[2] = random_uniform(-10.0f, 10.0f);
+ pj.primitives.P = random_uniform(0.1f, 1.0f);
/* make gradients zero */
pi.primitives.gradients.rho[0] = 0.0f;
pi.primitives.gradients.rho[1] = 0.0f;
@@ -153,9 +146,43 @@ int main(int argc, char *argv[]) {
dx[2] = -dx[2];
runner_iact_nonsym_force(r2, dx, pj2.h, pi2.h, &pj2, &pi2);
- /* Check that the particles are the same */
+/* Check that the particles are the same */
+#if defined(GIZMO_SPH)
+ i_ok = 0;
+ j_ok = 0;
+ for (size_t i = 0; i < sizeof(struct part) / sizeof(float); ++i) {
+ float a = *(((float *)&pi) + i);
+ float b = *(((float *)&pi2) + i);
+ float c = *(((float *)&pj) + i);
+ float d = *(((float *)&pj2) + i);
+
+ int a_is_b;
+ if ((a + b)) {
+ a_is_b = (fabs((a - b) / (a + b)) > 1.e-4);
+ } else {
+ a_is_b = !(a == 0.0f);
+ }
+ int c_is_d;
+ if ((c + d)) {
+ c_is_d = (fabs((c - d) / (c + d)) > 1.e-4);
+ } else {
+ c_is_d = !(c == 0.0f);
+ }
+
+ if (a_is_b) {
+ message("%.8e, %.8e, %lu", a, b, i);
+ }
+ if (c_is_d) {
+ message("%.8e, %.8e, %lu", c, d, i);
+ }
+
+ i_ok |= a_is_b;
+ j_ok |= c_is_d;
+ }
+#else
i_ok = memcmp(&pi, &pi2, sizeof(struct part));
j_ok = memcmp(&pj, &pj2, sizeof(struct part));
+#endif
if (i_ok) {
printParticle_single(&pi, &xpi);