Commit 2643736b authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Merge branch 'thorough_riemann_solver_tests' into 'master'

Update to the Riemann solver unit tests

See merge request !529
parents 7db83d29 c89526be
......@@ -16,26 +16,83 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
#include "../config.h"
/* Some standard headers. */
#include <fenv.h>
#include <stdio.h>
#include <string.h>
#include "error.h"
/* Local headers. */
#include "riemann/riemann_exact.h"
#include "tools.h"
#include "swift.h"
int opposite(float a, float b) {
if ((a - b)) {
return fabs((a + b) / (a - b)) < 1.e-4;
} else {
return a == 0.0f;
const float max_abs_error = 1e-3f;
const float max_rel_error = 1e-2f;
const float min_threshold = 1e-2f;
/**
* @brief Checks whether two numbers are opposite of each others.
*/
int are_symmetric(float a, float b) {
/* Check that the signs are different */
if ((a * b) > 0.f) {
message("Identical signs a=%.8e b=%.8e", a, b);
return 0;
}
const float abs_a = fabsf(a);
const float abs_b = fabsf(b);
const float abs_error = fabsf(abs_a - abs_b);
/* Avoid FPEs... */
if (fabsf(abs_a + abs_b) == 0.f) {
return 1;
}
/* Avoid things close to 0 */
if ((abs_a < min_threshold) || (abs_b < min_threshold)) {
return 1;
}
const float rel_error = 0.5f * abs_error / fabsf(abs_a + abs_b);
/* Check that we do not breach the relative error limit */
if (rel_error > max_rel_error) {
message("Relative error too large a=%.8e b=%.8e rel=%.8e", a, b, rel_error);
return 0;
}
/* All good */
return 1;
}
int equal(float a, float b) {
if ((a + b)) {
return fabs((a - b) / (a + b)) < 1.e-4;
} else {
return a == 0.0f;
const float abs_error = fabsf(a - b);
/* Avoid FPEs... */
if (fabsf(a + b) == 0.f) {
return 1;
}
/* Avoid things close to 0 */
if ((fabsf(a) < min_threshold) || (fabsf(b) < min_threshold)) {
return 1;
}
const float rel_error = 0.5f * abs_error / fabsf(a + b);
/* Check that we do not breach the relative error limit */
if (rel_error > max_rel_error) {
message("Relative error too large a=%.8e b=%.8e rel=%.8e", a, b, rel_error);
return 0;
}
/* All good */
return 1;
}
/**
......@@ -46,7 +103,8 @@ int equal(float a, float b) {
* @param s String used to identify this check in messages
*/
void check_value(float a, float b, const char* s) {
if (fabsf(a - b) / fabsf(a + b) > 1.e-5f && fabsf(a - b) > 1.e-5f) {
if (fabsf(a + b) != 0.f && fabsf(a - b) / fabsf(a + b) > max_rel_error &&
fabsf(a - b) > max_abs_error) {
error("Values are inconsistent: %g %g (%s)!", a, b, s);
} else {
message("Values are consistent: %g %g (%s).", a, b, s);
......@@ -295,17 +353,22 @@ void check_riemann_symmetry() {
riemann_solve_for_flux(WL, WR, n_unit1, vij, totflux1);
riemann_solve_for_flux(WR, WL, n_unit2, vij, totflux2);
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])) {
if (!are_symmetric(totflux1[0], totflux2[0]) ||
!are_symmetric(totflux1[1], totflux2[1]) ||
!are_symmetric(totflux1[2], totflux2[2]) ||
!are_symmetric(totflux1[3], totflux2[3]) ||
!are_symmetric(totflux1[4], totflux2[4])) {
message("WL=[%.8e, %.8e, %.8e, %.8e, %.8e]", WL[0], WL[1], WL[2], WL[3],
WL[4]);
message("WR=[%.8e, %.8e, %.8e, %.8e, %.8e]", WR[0], WR[1], WR[2], WR[3],
WR[4]);
message("n_unit1=[%.8e, %.8e, %.8e]", n_unit1[0], n_unit1[1], n_unit1[2]);
message("vij=[%.8e, %.8e, %.8e]\n", vij[0], vij[1], vij[2]);
message(
"Flux solver asymmetric: [%.3e,%.3e,%.3e,%.3e,%.3e] == "
"[%.3e,%.3e,%.3e,%.3e,%.3e]\n",
totflux1[0], totflux1[1], totflux1[2], totflux1[3], totflux1[4],
totflux2[0], totflux2[1], totflux2[2], totflux2[3], totflux2[4]);
message("Asymmetry in flux solution!");
/* This asymmetry is to be expected, since we do an iteration. Are the
results at least consistent? */
check_value(totflux1[0], totflux2[0], "Mass flux");
......@@ -313,6 +376,7 @@ void check_riemann_symmetry() {
check_value(totflux1[2], totflux2[2], "Momentum[1] flux");
check_value(totflux1[3], totflux2[3], "Momentum[2] flux");
check_value(totflux1[4], totflux2[4], "Energy flux");
error("Asymmetry in flux solution!");
} else {
/* message( */
/* "Flux solver symmetric: [%.3e,%.3e,%.3e,%.3e,%.3e] == " */
......@@ -327,12 +391,27 @@ void check_riemann_symmetry() {
*/
int main() {
/* Initialize CPU frequency, this also starts time. */
unsigned long long cpufreq = 0;
clocks_set_cpufreq(cpufreq);
/* Choke on FP-exceptions */
#ifdef HAVE_FE_ENABLE_EXCEPT
feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
#endif
/* Get some randomness going */
const int seed = time(NULL);
message("Seed = %d", seed);
srand(seed);
/* check the exact Riemann solver */
check_riemann_exact();
/* symmetry test */
int i;
for (i = 0; i < 100; ++i) check_riemann_symmetry();
for (int i = 0; i < 100000; ++i) {
check_riemann_symmetry();
}
return 0;
}
......@@ -16,16 +16,68 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
#include "../config.h"
/* Some standard headers. */
#include <fenv.h>
#include <math.h>
#include <stdio.h>
#include <string.h>
#include "error.h"
/* Local headers. */
#include "riemann/riemann_hllc.h"
#include "tools.h"
#include "swift.h"
int consistent_with_zero(float val) { return fabs(val) < 1.e-4; }
const float max_abs_error = 1e-3f;
const float max_rel_error = 1e-3f;
const float min_threshold = 1e-2f;
/**
* @brief Check the symmetry of the TRRS Riemann solver
* @brief Checks whether two numbers are opposite of each others.
*/
int are_symmetric(float a, float b) {
/* Check that the signs are different */
if ((a * b) > 0.f) {
message("Identical signs a=%.8e b=%.8e", a, b);
return 0;
}
const float abs_a = fabsf(a);
const float abs_b = fabsf(b);
const float abs_error = fabsf(abs_a - abs_b);
/* Check that we do not breach the absolute error limit */
if (abs_error > max_abs_error) {
message("Absolute error too large a=%.8e b=%.8e abs=%.8e", a, b, abs_error);
return 0;
}
/* Avoid FPEs... */
if (fabsf(abs_a + abs_b) == 0.f) {
return 1;
}
/* Avoid things close to 0 */
if ((abs_a < min_threshold) || (abs_b < min_threshold)) {
return 1;
}
const float rel_error = 0.5f * abs_error / fabsf(abs_a + abs_b);
/* Check that we do not breach the relative error limit */
if (rel_error > max_rel_error) {
message("Relative error too large a=%.8e b=%.8e rel=%.8e", a, b, rel_error);
return 0;
}
/* All good */
return 1;
}
/**
* @brief Check the symmetry of the HLLC Riemann solver for a random setup
*/
void check_riemann_symmetry() {
float WL[5], WR[5], n_unit1[3], n_unit2[3], n_norm, vij[3], totflux1[5],
......@@ -63,14 +115,20 @@ void check_riemann_symmetry() {
riemann_solve_for_flux(WL, WR, n_unit1, vij, totflux1);
riemann_solve_for_flux(WR, WL, n_unit2, vij, totflux2);
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])) {
if (!are_symmetric(totflux1[0], totflux2[0]) ||
!are_symmetric(totflux1[1], totflux2[1]) ||
!are_symmetric(totflux1[2], totflux2[2]) ||
!are_symmetric(totflux1[3], totflux2[3]) ||
!are_symmetric(totflux1[4], totflux2[4])) {
message("WL=[%.8e, %.8e, %.8e, %.8e, %.8e]", WL[0], WL[1], WL[2], WL[3],
WL[4]);
message("WR=[%.8e, %.8e, %.8e, %.8e, %.8e]", WR[0], WR[1], WR[2], WR[3],
WR[4]);
message("n_unit1=[%.8e, %.8e, %.8e]", n_unit1[0], n_unit1[1], n_unit1[2]);
message("vij=[%.8e, %.8e, %.8e]\n", vij[0], vij[1], vij[2]);
message(
"Flux solver asymmetric: [%.3e,%.3e,%.3e,%.3e,%.3e] == "
"[%.3e,%.3e,%.3e,%.3e,%.3e]\n",
"Flux solver asymmetric: [%.6e,%.6e,%.6e,%.6e,%.6e] == "
"[%.6e,%.6e,%.6e,%.6e,%.6e]\n",
totflux1[0], totflux1[1], totflux1[2], totflux1[3], totflux1[4],
totflux2[0], totflux2[1], totflux2[2], totflux2[3], totflux2[4]);
error("Asymmetry in flux solution!");
......@@ -88,9 +146,24 @@ void check_riemann_symmetry() {
*/
int main() {
int i;
/* Initialize CPU frequency, this also starts time. */
unsigned long long cpufreq = 0;
clocks_set_cpufreq(cpufreq);
/* Choke on FP-exceptions */
#ifdef HAVE_FE_ENABLE_EXCEPT
feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
#endif
/* Get some randomness going */
const int seed = time(NULL);
message("Seed = %d", seed);
srand(seed);
/* symmetry test */
for (i = 0; i < 100; i++) check_riemann_symmetry();
for (int i = 0; i < 100000; i++) {
check_riemann_symmetry();
}
return 0;
}
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment