Streamlined verification

include/astaroth_utils.h

@@ -31,11 +31,7 @@ extern "C" {
 
 #include <stdbool.h>
 
-#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);
 
-/** */
+Error acGetError(const AcReal model, const AcReal candidate);
-AcResult acVerifyMesh(const AcMesh model, const AcMesh candidate);
 
-/** */
+bool acEvalError(const char* label, const Error error);
-AcResult acVerifyScalReductions(const AcMesh model, const AcScalReductionTestCase* testCases,
-                                const size_t numCases);
 
-/** */
+AcResult acVerifyMesh(const char* label, const AcMesh model, const AcMesh candidate);
-AcResult acVerifyVecReductions(const AcMesh model, const AcVecReductionTestCase* testCases,
-                               const size_t numCases);
 
 #ifdef __cplusplus
 } // extern "C"

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... ");

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... ");

samples/mpitest/main.cc

@@ -21,12 +21,15 @@
 */
 #include "astaroth.h"
 #include "astaroth_utils.h"
+#include "errchk.h"
 
 #if AC_MPI_ENABLED
 
 #include <mpi.h>
 #include <vector>
 
+#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<AcScalReductionTestCase> 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<AcVecReductionTestCase> 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());
+        printf("---Test: Scalar reductions---\n");
-        acVerifyVecReductions(model, vectorReductionTests.data(), vectorReductionTests.size());
+        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);

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

src/utils/verification.c

@@ -4,6 +4,8 @@
 #include <stdbool.h>
 #include <string.h>
 
+#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,
+get_max_abs_error(const AcReal* model, const AcReal* candidate, const AcMeshInfo info)
-                  const AcMesh candidate_mesh)
 {
-    AcReal* model_vtxbuf     = model_mesh.vertex_buffer[vtxbuf_handle];
+    Error error = {.abs_error = -1};
-    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]);
 
+    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.maximum_magnitude = get_maximum_magnitude(model, info);
-    error.handle = vtxbuf_handle;
+    error.minimum_magnitude = get_minimum_magnitude(model, info);
-    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);
 
     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);
+        const Error error = get_max_abs_error(model.vertex_buffer[i], candidate.vertex_buffer[i],
-        errors_found |= printErrorToScreen(field_error);
+                                              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"
+    if (errors_found > 0)
-                                : "Success. No errors found.");
+        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;
-}