diff --git a/include/astaroth_utils.h b/include/astaroth_utils.h index ee5ebd4..632d7ce 100644 --- a/include/astaroth_utils.h +++ b/include/astaroth_utils.h @@ -31,11 +31,7 @@ extern "C" { #include -#define ERROR_LABEL_LENGTH 30 - typedef struct { - char label[ERROR_LABEL_LENGTH]; - VertexBufferHandle handle; AcReal model; AcReal candidate; long double abs_error; @@ -45,42 +41,9 @@ typedef struct { AcReal minimum_magnitude; } Error; -typedef struct { - char label[ERROR_LABEL_LENGTH]; - VertexBufferHandle vtxbuf; - ReductionType rtype; - AcReal candidate; -} AcScalReductionTestCase; - -typedef struct { - char label[ERROR_LABEL_LENGTH]; - VertexBufferHandle a; - VertexBufferHandle b; - VertexBufferHandle c; - ReductionType rtype; - AcReal candidate; -} AcVecReductionTestCase; - -/** TODO comment */ -Error acGetError(AcReal model, AcReal candidate); - -/** TODO comment */ -bool printErrorToScreen(const Error error); - /** Loads data from the config file */ AcResult acLoadConfig(const char* config_path, AcMeshInfo* config); -/** */ -AcScalReductionTestCase acCreateScalReductionTestCase(const char* label, - const VertexBufferHandle vtxbuf, - const ReductionType rtype); - -/** */ -AcVecReductionTestCase acCreateVecReductionTestCase(const char* label, const VertexBufferHandle a, - const VertexBufferHandle b, - const VertexBufferHandle c, - const ReductionType rtype); - /** */ AcResult acMeshSet(const AcReal value, AcMesh* mesh); @@ -96,23 +59,18 @@ AcResult acMeshClear(AcMesh* mesh); /** */ AcResult acModelIntegrateStep(AcMesh mesh, const AcReal dt); -/** TODO */ +/** */ AcReal acModelReduceScal(const AcMesh mesh, const ReductionType rtype, const VertexBufferHandle a); -/** TODO */ +/** */ AcReal acModelReduceVec(const AcMesh mesh, const ReductionType rtype, const VertexBufferHandle a, const VertexBufferHandle b, const VertexBufferHandle c); -/** */ -AcResult acVerifyMesh(const AcMesh model, const AcMesh candidate); +Error acGetError(const AcReal model, const AcReal candidate); -/** */ -AcResult acVerifyScalReductions(const AcMesh model, const AcScalReductionTestCase* testCases, - const size_t numCases); +bool acEvalError(const char* label, const Error error); -/** */ -AcResult acVerifyVecReductions(const AcMesh model, const AcVecReductionTestCase* testCases, - const size_t numCases); +AcResult acVerifyMesh(const char* label, const AcMesh model, const AcMesh candidate); #ifdef __cplusplus } // extern "C" diff --git a/samples/cpptest/main.cc b/samples/cpptest/main.cc index 8e8e9b8..3dde831 100644 --- a/samples/cpptest/main.cc +++ b/samples/cpptest/main.cc @@ -42,7 +42,7 @@ main(void) acInit(info); acLoad(model); acStore(&candidate); - acVerifyMesh(model, candidate); + acVerifyMesh("Load/Store", model, candidate); // Attempt to integrate and check max and min printf("Integrating... "); diff --git a/samples/ctest/main.c b/samples/ctest/main.c index 8bf4e2d..b340154 100644 --- a/samples/ctest/main.c +++ b/samples/ctest/main.c @@ -41,7 +41,7 @@ main(void) acInit(info); acLoad(model); acStore(&candidate); - acVerifyMesh(model, candidate); + acVerifyMesh("Load/Store", model, candidate); // Attempt to integrate and check max and min printf("Integrating... "); diff --git a/samples/mpitest/main.cc b/samples/mpitest/main.cc index 492c1df..606ea3b 100644 --- a/samples/mpitest/main.cc +++ b/samples/mpitest/main.cc @@ -21,12 +21,15 @@ */ #include "astaroth.h" #include "astaroth_utils.h" +#include "errchk.h" #if AC_MPI_ENABLED #include #include +#define ARRAY_SIZE(arr) (sizeof(arr) / sizeof(*arr)) + int main(void) { @@ -35,6 +38,9 @@ main(void) MPI_Comm_size(MPI_COMM_WORLD, &nprocs); MPI_Comm_rank(MPI_COMM_WORLD, &pid); + // Set random seed for reproducibility + srand(321654987); + // CPU alloc AcMeshInfo info; acLoadConfig(AC_DEFAULT_CONFIG, &info); @@ -50,57 +56,69 @@ main(void) // GPU alloc & compute acGridInit(info); - // INTEGRATION TESTS START --------------------------------------------------------------------- + // Boundconds + acGridLoadMesh(STREAM_DEFAULT, model); + acGridPeriodicBoundconds(STREAM_DEFAULT); + acGridStoreMesh(STREAM_DEFAULT, &candidate); + if (pid == 0) { + acMeshApplyPeriodicBounds(&model); + const AcResult res = acVerifyMesh("Boundconds", model, candidate); + ERRCHK_ALWAYS(res == AC_SUCCESS); + acMeshRandomize(&model); + } + + // Integration acGridLoadMesh(STREAM_DEFAULT, model); acGridIntegrate(STREAM_DEFAULT, FLT_EPSILON); acGridPeriodicBoundconds(STREAM_DEFAULT); acGridStoreMesh(STREAM_DEFAULT, &candidate); - if (pid == 0) { acModelIntegrateStep(model, FLT_EPSILON); acMeshApplyPeriodicBounds(&model); - acVerifyMesh(model, candidate); + const AcResult res = acVerifyMesh("Integration", model, candidate); + ERRCHK_ALWAYS(res == AC_SUCCESS); + acMeshRandomize(&model); } - // INTEGRATION TESTS END ----------------------------------------------------------------------- - // REDUCTION TESTS START ----------------------------------------------------------------------- + // Scalar reductions acGridLoadMesh(STREAM_DEFAULT, model); - std::vector scalarReductionTests{ - acCreateScalReductionTestCase("Scalar MAX", VTXBUF_UUX, RTYPE_MAX), - acCreateScalReductionTestCase("Scalar MIN", VTXBUF_UUX, RTYPE_MIN), - /* - acCreateScalReductionTestCase("Scalar RMS", VTXBUF_UUX, RTYPE_RMS), - acCreateScalReductionTestCase("Scalar RMS_EXP", VTXBUF_UUX, RTYPE_RMS_EXP), - acCreateScalReductionTestCase("Scalar SUM", VTXBUF_UUX, RTYPE_SUM), - */ - }; - std::vector vectorReductionTests{ - acCreateVecReductionTestCase("Vector MAX", VTXBUF_UUX, VTXBUF_UUY, VTXBUF_UUZ, RTYPE_MAX), - acCreateVecReductionTestCase("Vector MIN", VTXBUF_UUX, VTXBUF_UUY, VTXBUF_UUZ, RTYPE_MIN), - /* - acCreateVecReductionTestCase("Vector RMS", VTXBUF_UUX, VTXBUF_UUY, VTXBUF_UUZ, RTYPE_RMS), - acCreateVecReductionTestCase("Vector RMS_EXP", VTXBUF_UUX, VTXBUF_UUY, VTXBUF_UUZ, - RTYPE_RMS_EXP), - acCreateVecReductionTestCase("Vector SUM", VTXBUF_UUX, VTXBUF_UUY, VTXBUF_UUZ, RTYPE_SUM), - */ - }; - // False positives due to too strict error bounds, skip the tests until we can determine a - // proper error bound - fprintf(stderr, "WARNING: RTYPE_RMS, RTYPE_RMS_EXP, and RTYPE_SUM tests skipped\n"); - - for (auto& testCase : scalarReductionTests) { - acGridReduceScal(STREAM_DEFAULT, testCase.rtype, testCase.vtxbuf, &testCase.candidate); - } - for (auto& testCase : vectorReductionTests) { - acGridReduceVec(STREAM_DEFAULT, testCase.rtype, testCase.a, testCase.b, testCase.c, - &testCase.candidate); - } if (pid == 0) { - acVerifyScalReductions(model, scalarReductionTests.data(), scalarReductionTests.size()); - acVerifyVecReductions(model, vectorReductionTests.data(), vectorReductionTests.size()); + printf("---Test: Scalar reductions---\n"); + printf("Warning: testing only RTYPE_MAX and RTYPE_MIN\n"); + } + for (size_t i = 0; i < 2; ++i) { // NOTE: 2 instead of NUM_RTYPES + const VertexBufferHandle v0 = VTXBUF_UUX; + AcReal candval; + acGridReduceScal(STREAM_DEFAULT, (ReductionType)i, v0, &candval); + if (pid == 0) { + const AcReal modelval = acModelReduceScal(model, (ReductionType)i, v0); + Error error = acGetError(modelval, candval); + error.maximum_magnitude = acModelReduceScal(model, RTYPE_MAX, v0); + error.minimum_magnitude = acModelReduceScal(model, RTYPE_MIN, v0); + ERRCHK_ALWAYS(acEvalError(rtype_names[i], error)); + } + } + + // Vector reductions + if (pid == 0) { + printf("---Test: Vector reductions---\n"); + printf("Warning: testing only RTYPE_MAX and RTYPE_MIN\n"); + } + for (size_t i = 0; i < 2; ++i) { // NOTE: 2 instead of NUM_RTYPES + const VertexBufferHandle v0 = VTXBUF_UUX; + const VertexBufferHandle v1 = VTXBUF_UUY; + const VertexBufferHandle v2 = VTXBUF_UUZ; + AcReal candval; + acGridReduceVec(STREAM_DEFAULT, (ReductionType)i, v0, v1, v2, &candval); + if (pid == 0) { + const AcReal modelval = acModelReduceVec(model, (ReductionType)i, v0, v1, v2); + Error error = acGetError(modelval, candval); + error.maximum_magnitude = acModelReduceVec(model, RTYPE_MAX, v0, v1, v2); + error.minimum_magnitude = acModelReduceVec(model, RTYPE_MIN, v0, v1, v1); + ERRCHK_ALWAYS(acEvalError(rtype_names[i], error)); + } } - // REDUCTION TESTS END ------------------------------------------------------------------------- if (pid == 0) { acMeshDestroy(&model); diff --git a/samples/standalone/autotest.cc b/samples/standalone/autotest.cc index 6362965..14ae784 100644 --- a/samples/standalone/autotest.cc +++ b/samples/standalone/autotest.cc @@ -245,7 +245,7 @@ check_reductions(const AcMeshInfo& config) acLoad(*mesh); - for (int rtype = 0; rtype < NUM_REDUCTION_TYPES; ++rtype) { + for (int rtype = 0; rtype < NUM_RTYPES; ++rtype) { if (rtype == RTYPE_SUM) { // Skip SUM test for now. The failure is either caused by floating-point diff --git a/src/utils/verification.c b/src/utils/verification.c index 393fbc3..5a5cad1 100644 --- a/src/utils/verification.c +++ b/src/utils/verification.c @@ -4,6 +4,8 @@ #include #include +#include "errchk.h" + #define max(a, b) ((a) > (b) ? (a) : (b)) #define min(a, b) ((a) < (b) ? (a) : (b)) @@ -26,7 +28,7 @@ is_valid(const AcReal a) } Error -acGetError(AcReal model, AcReal candidate) +acGetError(const AcReal model, const AcReal candidate) { Error error; error.abs_error = 0; @@ -57,9 +59,48 @@ acGetError(AcReal model, AcReal candidate) error.rel_error = fabsl(1.0l - candidate / model) / machine_epsilon; } + error.maximum_magnitude = error.minimum_magnitude = 0; + return error; } +static inline void +print_error_to_file(const char* label, const Error error, const char* path) +{ + FILE* file = fopen(path, "a"); + ERRCHK_ALWAYS(file); + fprintf(file, "%s, %Lg, %Lg, %Lg, %g, %g\n", label, error.ulp_error, error.abs_error, + error.rel_error, (double)error.maximum_magnitude, (double)error.minimum_magnitude); + fclose(file); +} + +/** Returns true if the error is acceptable, false otherwise. */ +bool +acEvalError(const char* label, const Error error) +{ + // Accept the error if the relative error is < max_ulp_error ulps. + // Also consider the error zero if it is less than the minimum value in the mesh scaled to + // machine epsilon + const long double max_ulp_error = 5; + + bool acceptable; + if (error.ulp_error < max_ulp_error) + acceptable = true; + else if (error.abs_error < error.minimum_magnitude * AC_REAL_EPSILON) + acceptable = true; + else + acceptable = false; + + printf("%-15s... %s ", label, acceptable ? GRN "OK! " RESET : RED "FAIL! " RESET); + + printf("| %.3Lg (abs), %.3Lg (ulps), %.3Lg (rel). Range: [%.3g, %.3g]\n", // + error.abs_error, error.ulp_error, error.rel_error, // + (double)error.minimum_magnitude, (double)error.maximum_magnitude); + print_error_to_file(label, error, "verification.out"); + + return acceptable; +} + static AcReal get_maximum_magnitude(const AcReal* field, const AcMeshInfo info) { @@ -88,168 +129,39 @@ get_minimum_magnitude(const AcReal* field, const AcMeshInfo info) // floating-point precision range and gives huge errors with values that should be considered // zero (f.ex. 1e-19 and 1e-22 give error of around 1e4 ulps) static Error -get_max_abs_error(const VertexBufferHandle vtxbuf_handle, const AcMesh model_mesh, - const AcMesh candidate_mesh) +get_max_abs_error(const AcReal* model, const AcReal* candidate, const AcMeshInfo info) { - AcReal* model_vtxbuf = model_mesh.vertex_buffer[vtxbuf_handle]; - AcReal* candidate_vtxbuf = candidate_mesh.vertex_buffer[vtxbuf_handle]; - - Error error; - error.abs_error = -1; - - for (size_t i = 0; i < acVertexBufferSize(model_mesh.info); ++i) { - - Error curr_error = acGetError(model_vtxbuf[i], candidate_vtxbuf[i]); + Error error = {.abs_error = -1}; + for (size_t i = 0; i < acVertexBufferSize(info); ++i) { + Error curr_error = acGetError(model[i], candidate[i]); if (curr_error.abs_error > error.abs_error) error = curr_error; } - - error.handle = vtxbuf_handle; - strncpy(error.label, vtxbuf_names[vtxbuf_handle], ERROR_LABEL_LENGTH - 1); - error.label[ERROR_LABEL_LENGTH - 1] = '\0'; - error.maximum_magnitude = get_maximum_magnitude(model_vtxbuf, model_mesh.info); - error.minimum_magnitude = get_minimum_magnitude(model_vtxbuf, model_mesh.info); + error.maximum_magnitude = get_maximum_magnitude(model, info); + error.minimum_magnitude = get_minimum_magnitude(model, info); return error; } -static inline void -print_error_to_file(const char* path, const int n, const Error error) -{ - FILE* file = fopen(path, "a"); - fprintf(file, "%d, %Lg, %Lg, %Lg, %g, %g\n", n, error.ulp_error, error.abs_error, - error.rel_error, (double)error.maximum_magnitude, (double)error.minimum_magnitude); - fclose(file); -} - -static bool -is_acceptable(const Error error) -{ - // Accept the error if the relative error is < max_ulp_error ulps. - // Also consider the error zero if it is less than the minimum value in the mesh scaled to - // machine epsilon - const long double max_ulp_error = 5; - - if (error.ulp_error < max_ulp_error) - return true; - else if (error.abs_error < error.minimum_magnitude * AC_REAL_EPSILON) - return true; - else - return false; -} - -bool -printErrorToScreen(const Error error) -{ - bool errors_found = false; - - printf("\t%-15s... ", error.label); - if (is_acceptable(error)) { - printf(GRN "OK! " RESET); - } - else { - printf(RED "FAIL! " RESET); - errors_found = true; - } - - fprintf(stdout, "| %.3Lg (abs), %.3Lg (ulps), %.3Lg (rel). Range: [%.3g, %.3g]\n", // - error.abs_error, error.ulp_error, error.rel_error, // - (double)error.minimum_magnitude, (double)error.maximum_magnitude); - - return errors_found; -} - /** Returns true when successful, false if errors were found. */ AcResult -acVerifyMesh(const AcMesh model, const AcMesh candidate) +acVerifyMesh(const char* label, const AcMesh model, const AcMesh candidate) { + printf("---Test: %s---\n", label); printf("Errors at the point of the maximum absolute error:\n"); - bool errors_found = false; + int errors_found = 0; for (int i = 0; i < NUM_VTXBUF_HANDLES; ++i) { - Error field_error = get_max_abs_error(i, model, candidate); - errors_found |= printErrorToScreen(field_error); + const Error error = get_max_abs_error(model.vertex_buffer[i], candidate.vertex_buffer[i], + model.info); + const bool acceptable = acEvalError(vtxbuf_names[i], error); + if (!acceptable) + ++errors_found; } - printf("%s\n", errors_found ? "Failure. Found errors in one or more vertex buffers" - : "Success. No errors found."); + if (errors_found > 0) + printf("Failure. Found %d errors\n", errors_found); return errors_found ? AC_FAILURE : AC_SUCCESS; } - -/** Verification function for scalar reductions*/ -AcResult -acVerifyScalReductions(const AcMesh model, const AcScalReductionTestCase* testCases, - const size_t numCases) -{ - printf("\nTesting scalar reductions:\n"); - - bool errors_found = false; - for (size_t i = 0; i < numCases; i++) { - AcReal model_reduction = acModelReduceScal(model, testCases[i].rtype, testCases[i].vtxbuf); - Error error = acGetError(model_reduction, testCases[i].candidate); - strncpy(error.label, testCases[i].label, ERROR_LABEL_LENGTH - 1); - error.label[ERROR_LABEL_LENGTH - 1] = '\0'; - errors_found |= printErrorToScreen(error); - } - printf("%s\n", errors_found ? "Failure. Found errors in one or more scalar reductions" - : "Success. No errors found."); - - return errors_found ? AC_FAILURE : AC_SUCCESS; -} - -/** Verification function for vector reductions*/ -AcResult -acVerifyVecReductions(const AcMesh model, const AcVecReductionTestCase* testCases, - const size_t numCases) -{ - printf("\nTesting vector reductions:\n"); - - bool errors_found = false; - for (size_t i = 0; i < numCases; i++) { - AcReal model_reduction = acModelReduceVec(model, testCases[i].rtype, testCases[i].a, - testCases[i].b, testCases[i].c); - Error error = acGetError(model_reduction, testCases[i].candidate); - strncpy(error.label, testCases[i].label, ERROR_LABEL_LENGTH - 1); - error.label[ERROR_LABEL_LENGTH - 1] = '\0'; - errors_found |= printErrorToScreen(error); - } - printf("%s\n", errors_found ? "Failure. Found errors in one or more vector reductions" - : "Success. No errors found."); - - return errors_found ? AC_FAILURE : AC_SUCCESS; -} - -/** Constructor for scalar reduction test case */ -AcScalReductionTestCase -acCreateScalReductionTestCase(const char* label, const VertexBufferHandle vtxbuf, const ReductionType rtype) -{ - AcScalReductionTestCase testCase; - - strncpy(testCase.label,label,ERROR_LABEL_LENGTH - 1); - testCase.label[ERROR_LABEL_LENGTH - 1] = '\0'; - testCase.vtxbuf = vtxbuf; - testCase.rtype = rtype; - testCase.candidate = 0; - - return testCase; -} - -/** Constructor for vector reduction test case */ -AcVecReductionTestCase -acCreateVecReductionTestCase(const char* label, const VertexBufferHandle a, - const VertexBufferHandle b, const VertexBufferHandle c, const ReductionType rtype) -{ - AcVecReductionTestCase testCase; - - strncpy(testCase.label,label,ERROR_LABEL_LENGTH - 1); - testCase.label[ERROR_LABEL_LENGTH - 1] = '\0'; - testCase.a = a; - testCase.b = b; - testCase.c = c; - testCase.rtype = rtype; - testCase.candidate = 0; - - return testCase; -}