Commit ff49f83e authored by Matthieu Schaller's avatar Matthieu Schaller Committed by Matthieu Schaller
Browse files

Make the Riemann solver unit test more comprehensive by testing many more configurations.

parent 047ce843
...@@ -16,15 +16,20 @@ ...@@ -16,15 +16,20 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>. * 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 <string.h>
#include "error.h"
/* Local headers. */
#include "riemann/riemann_exact.h" #include "riemann/riemann_exact.h"
#include "tools.h" #include "swift.h"
int opposite(float a, float b) { int opposite(float a, float b) {
if ((a - b)) { if ((a - b)) {
return fabs((a + b) / (a - b)) < 1.e-4; return fabsf((a + b) / (a - b)) < 1.e-4f;
} else { } else {
return a == 0.0f; return a == 0.0f;
} }
...@@ -32,7 +37,7 @@ int opposite(float a, float b) { ...@@ -32,7 +37,7 @@ int opposite(float a, float b) {
int equal(float a, float b) { int equal(float a, float b) {
if ((a + b)) { if ((a + b)) {
return fabs((a - b) / (a + b)) < 1.e-4; return fabsf((a - b) / (a + b)) < 1.e-4f;
} else { } else {
return a == 0.0f; return a == 0.0f;
} }
...@@ -46,7 +51,8 @@ int equal(float a, float b) { ...@@ -46,7 +51,8 @@ int equal(float a, float b) {
* @param s String used to identify this check in messages * @param s String used to identify this check in messages
*/ */
void check_value(float a, float b, const char* s) { 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) > 1.e-5f &&
fabsf(a - b) > 1.e-5f) {
error("Values are inconsistent: %g %g (%s)!", a, b, s); error("Values are inconsistent: %g %g (%s)!", a, b, s);
} else { } else {
message("Values are consistent: %g %g (%s).", a, b, s); message("Values are consistent: %g %g (%s).", a, b, s);
...@@ -327,12 +333,24 @@ void check_riemann_symmetry() { ...@@ -327,12 +333,24 @@ void check_riemann_symmetry() {
*/ */
int main() { 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 the exact Riemann solver */
check_riemann_exact(); check_riemann_exact();
/* symmetry test */ /* symmetry test */
int i; for (int i = 0; i < 10000; ++i) check_riemann_symmetry();
for (i = 0; i < 100; ++i) check_riemann_symmetry();
return 0; return 0
} \ No newline at end of file
...@@ -16,13 +16,18 @@ ...@@ -16,13 +16,18 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>. * 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 <string.h>
#include "error.h"
/* Local headers. */
#include "riemann/riemann_hllc.h" #include "riemann/riemann_hllc.h"
#include "tools.h" #include "swift.h"
int consistent_with_zero(float val) { return fabs(val) < 1.e-4; } int consistent_with_zero(float val) { return fabsf(val) < 1.e-4f; }
/** /**
* @brief Check the symmetry of the TRRS Riemann solver * @brief Check the symmetry of the TRRS Riemann solver
...@@ -68,6 +73,12 @@ void check_riemann_symmetry() { ...@@ -68,6 +73,12 @@ void check_riemann_symmetry() {
!consistent_with_zero(totflux1[2] + totflux2[2]) || !consistent_with_zero(totflux1[2] + totflux2[2]) ||
!consistent_with_zero(totflux1[3] + totflux2[3]) || !consistent_with_zero(totflux1[3] + totflux2[3]) ||
!consistent_with_zero(totflux1[4] + totflux2[4])) { !consistent_with_zero(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( message(
"Flux solver asymmetric: [%.3e,%.3e,%.3e,%.3e,%.3e] == " "Flux solver asymmetric: [%.3e,%.3e,%.3e,%.3e,%.3e] == "
"[%.3e,%.3e,%.3e,%.3e,%.3e]\n", "[%.3e,%.3e,%.3e,%.3e,%.3e]\n",
...@@ -88,9 +99,24 @@ void check_riemann_symmetry() { ...@@ -88,9 +99,24 @@ void check_riemann_symmetry() {
*/ */
int main() { 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 */ /* symmetry test */
for (i = 0; i < 100; i++) check_riemann_symmetry(); for (int i = 0; i < 10000; i++) {
check_riemann_symmetry();
}
return 0; 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