Started preparing the MPI version for benchmarks and added a solve-independent version of the verification functions to the utils library

src/mpitest/CMakeLists.txt

@@ -9,4 +9,5 @@ find_package(MPI REQUIRED)
 
 add_executable(mpitest main.c)
 target_include_directories(mpitest PRIVATE ${CMAKE_SOURCE_DIR}/src/standalone ${MPI_C_INCLUDE_PATH})
-target_link_libraries(mpitest astaroth_core astaroth_standalone ${MPI_C_LIBRARIES})
+target_link_libraries(mpitest astaroth_core astaroth_utils ${MPI_C_LIBRARIES})
+target_compile_definitions(mpitest PRIVATE -DAC_DEFAULT_CONFIG="${CMAKE_SOURCE_DIR}/config/astaroth.conf")

src/mpitest/main.c

@@ -26,158 +26,35 @@
 #include <string.h>
 
 #include "astaroth.h"
-#include "autotest.h"
 
 #include <mpi.h>
 
-// From Astaroth Standalone
+// From Astaroth Utils
-#include "config_loader.h"
+#include "src/utils/config_loader.h"
-#include "model/host_memory.h"
+#include "src/utils/memory.h"
-
+#include "src/utils/verification.h"
-static void
-distribute_mesh(const AcMesh* src, AcMesh* dst)
-{
-    MPI_Datatype datatype = MPI_FLOAT;
-    if (sizeof(AcReal) == 8)
-        datatype = MPI_DOUBLE;
-
-    int process_id, num_processes;
-    MPI_Comm_rank(MPI_COMM_WORLD, &process_id);
-    MPI_Comm_size(MPI_COMM_WORLD, &num_processes);
-
-    const size_t count = acVertexBufferSize(dst->info);
-    for (int i = 0; i < NUM_VTXBUF_HANDLES; ++i) {
-
-        // Communicate to self
-        if (process_id == 0) {
-            assert(src);
-            assert(dst);
-            memcpy(&dst->vertex_buffer[i][0], //
-                   &src->vertex_buffer[i][0], //
-                   count * sizeof(src->vertex_buffer[i][0]));
-        }
-        // Communicate to others
-        for (int j = 1; j < num_processes; ++j) {
-            if (process_id == 0) {
-                assert(src);
-
-                // Send
-                // TODO RECHECK THESE j INDICES
-                const size_t src_idx = j * dst->info.int_params[AC_mx] *
-                                       dst->info.int_params[AC_my] * src->info.int_params[AC_nz] /
-                                       num_processes;
-
-                MPI_Send(&src->vertex_buffer[i][src_idx], count, datatype, j, 0, MPI_COMM_WORLD);
-            }
-            else {
-                assert(dst);
-
-                // Recv
-                const size_t dst_idx = 0;
-                MPI_Status status;
-                MPI_Recv(&dst->vertex_buffer[i][dst_idx], count, datatype, 0, 0, MPI_COMM_WORLD,
-                         &status);
-            }
-        }
-    }
-}
-
-static void
-gather_mesh(const AcMesh* src, AcMesh* dst)
-{
-    MPI_Datatype datatype = MPI_FLOAT;
-    if (sizeof(AcReal) == 8)
-        datatype = MPI_DOUBLE;
-
-    int process_id, num_processes;
-    MPI_Comm_rank(MPI_COMM_WORLD, &process_id);
-    MPI_Comm_size(MPI_COMM_WORLD, &num_processes);
-
-    size_t count = acVertexBufferSize(src->info);
-
-    for (int i = 0; i < NUM_VTXBUF_HANDLES; ++i) {
-        // Communicate to self
-        if (process_id == 0) {
-            assert(src);
-            assert(dst);
-            memcpy(&dst->vertex_buffer[i][0], //
-                   &src->vertex_buffer[i][0], //
-                   count * sizeof(AcReal));
-        }
-
-        // Communicate to others
-        for (int j = 1; j < num_processes; ++j) {
-            if (process_id == 0) {
-                // Recv
-                // const size_t dst_idx = j * acVertexBufferCompdomainSize(dst->info);
-                const size_t dst_idx = j * dst->info.int_params[AC_mx] *
-                                       dst->info.int_params[AC_my] * dst->info.int_params[AC_nz] /
-                                       num_processes;
-
-                assert(dst_idx + count <= acVertexBufferSize(dst->info));
-                MPI_Status status;
-                MPI_Recv(&dst->vertex_buffer[i][dst_idx], count, datatype, j, 0, MPI_COMM_WORLD,
-                         &status);
-            }
-            else {
-                // Send
-                const size_t src_idx = 0;
-
-                assert(src_idx + count <= acVertexBufferSize(src->info));
-                MPI_Send(&src->vertex_buffer[i][src_idx], count, datatype, 0, 0, MPI_COMM_WORLD);
-            }
-        }
-    }
-}
 
 int
 main(void)
 {
     MPI_Init(NULL, NULL);
 
-    int num_processes, process_id;
+    AcMeshInfo info;
-    MPI_Comm_size(MPI_COMM_WORLD, &num_processes);
+    acLoadConfig(AC_DEFAULT_CONFIG, &info);
-    MPI_Comm_rank(MPI_COMM_WORLD, &process_id);
 
-    char processor_name[MPI_MAX_PROCESSOR_NAME];
+    AcMesh model;
-    int name_len;
+    acMeshCreate(info, &model);
-    MPI_Get_processor_name(processor_name, &name_len);
-    printf("Processor %s. Process %d of %d.\n", processor_name, process_id, num_processes);
 
-    AcMeshInfo mesh_info;
+    AcMesh candidate;
-    load_config(&mesh_info);
+    acMeshCreate(info, &candidate);
-    update_config(&mesh_info);
 
-    AcMesh* main_mesh     = NULL;
+    acMeshSet(1, &model);
-    ModelMesh* model_mesh = NULL;
+    acMeshSet(1.00000005, &candidate);
-    if (process_id == 0) {
-        main_mesh = acmesh_create(mesh_info);
-        acmesh_init_to(INIT_TYPE_RANDOM, main_mesh);
-        model_mesh = modelmesh_create(mesh_info);
-        acmesh_to_modelmesh(*main_mesh, model_mesh);
-    }
 
-    AcMeshInfo submesh_info = mesh_info;
+    acVerifyMesh(model, candidate);
-    submesh_info.int_params[AC_nz] /= num_processes;
-    update_config(&submesh_info);
-
-    AcMesh* submesh = acmesh_create(submesh_info);
-
-    /////////////////////
-    distribute_mesh(main_mesh, submesh);
-    gather_mesh(submesh, main_mesh);
-    /////////////////////////
-    // Autotest
-    bool is_acceptable = verify_meshes(*model_mesh, *main_mesh);
-    /////
-
-    acmesh_destroy(submesh);
-
-    if (process_id == 0) {
-        modelmesh_destroy(model_mesh);
-        acmesh_destroy(main_mesh);
-    }
 
+    acMeshDestroy(&model);
+    acMeshDestroy(&candidate);
     MPI_Finalize();
     return EXIT_SUCCESS;
 }

src/mpitest/main0.c

@@ -0,0 +1,181 @@
+/*
+    Copyright (C) 2014-2019, Johannes Pekkilae, Miikka Vaeisalae.
+
+    This file is part of Astaroth.
+
+    Astaroth is free software: you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    Astaroth is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with Astaroth.  If not, see <http://www.gnu.org/licenses/>.
+*/
+/**
+    Running: mpirun -np <num processes> <executable>
+*/
+#undef NDEBUG // Assert always
+#include <assert.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include "astaroth.h"
+
+#include <mpi.h>
+
+// From Astaroth Utils
+#include "src/utils/config_loader.h"
+#include "src/utils/memory.h"
+
+static void
+distribute_mesh(const AcMesh* src, AcMesh* dst)
+{
+    MPI_Datatype datatype = MPI_FLOAT;
+    if (sizeof(AcReal) == 8)
+        datatype = MPI_DOUBLE;
+
+    int process_id, num_processes;
+    MPI_Comm_rank(MPI_COMM_WORLD, &process_id);
+    MPI_Comm_size(MPI_COMM_WORLD, &num_processes);
+
+    const size_t count = acVertexBufferSize(dst->info);
+    for (int i = 0; i < NUM_VTXBUF_HANDLES; ++i) {
+
+        // Communicate to self
+        if (process_id == 0) {
+            assert(src);
+            assert(dst);
+            memcpy(&dst->vertex_buffer[i][0], //
+                   &src->vertex_buffer[i][0], //
+                   count * sizeof(src->vertex_buffer[i][0]));
+        }
+        // Communicate to others
+        for (int j = 1; j < num_processes; ++j) {
+            if (process_id == 0) {
+                assert(src);
+
+                // Send
+                // TODO RECHECK THESE j INDICES
+                const size_t src_idx = j * dst->info.int_params[AC_mx] *
+                                       dst->info.int_params[AC_my] * src->info.int_params[AC_nz] /
+                                       num_processes;
+
+                MPI_Send(&src->vertex_buffer[i][src_idx], count, datatype, j, 0, MPI_COMM_WORLD);
+            }
+            else {
+                assert(dst);
+
+                // Recv
+                const size_t dst_idx = 0;
+                MPI_Status status;
+                MPI_Recv(&dst->vertex_buffer[i][dst_idx], count, datatype, 0, 0, MPI_COMM_WORLD,
+                         &status);
+            }
+        }
+    }
+}
+
+static void
+gather_mesh(const AcMesh* src, AcMesh* dst)
+{
+    MPI_Datatype datatype = MPI_FLOAT;
+    if (sizeof(AcReal) == 8)
+        datatype = MPI_DOUBLE;
+
+    int process_id, num_processes;
+    MPI_Comm_rank(MPI_COMM_WORLD, &process_id);
+    MPI_Comm_size(MPI_COMM_WORLD, &num_processes);
+
+    size_t count = acVertexBufferSize(src->info);
+
+    for (int i = 0; i < NUM_VTXBUF_HANDLES; ++i) {
+        // Communicate to self
+        if (process_id == 0) {
+            assert(src);
+            assert(dst);
+            memcpy(&dst->vertex_buffer[i][0], //
+                   &src->vertex_buffer[i][0], //
+                   count * sizeof(AcReal));
+        }
+
+        // Communicate to others
+        for (int j = 1; j < num_processes; ++j) {
+            if (process_id == 0) {
+                // Recv
+                // const size_t dst_idx = j * acVertexBufferCompdomainSize(dst->info);
+                const size_t dst_idx = j * dst->info.int_params[AC_mx] *
+                                       dst->info.int_params[AC_my] * dst->info.int_params[AC_nz] /
+                                       num_processes;
+
+                assert(dst_idx + count <= acVertexBufferSize(dst->info));
+                MPI_Status status;
+                MPI_Recv(&dst->vertex_buffer[i][dst_idx], count, datatype, j, 0, MPI_COMM_WORLD,
+                         &status);
+            }
+            else {
+                // Send
+                const size_t src_idx = 0;
+
+                assert(src_idx + count <= acVertexBufferSize(src->info));
+                MPI_Send(&src->vertex_buffer[i][src_idx], count, datatype, 0, 0, MPI_COMM_WORLD);
+            }
+        }
+    }
+}
+
+int
+main(void)
+{
+    MPI_Init(NULL, NULL);
+
+    int num_processes, process_id;
+    MPI_Comm_size(MPI_COMM_WORLD, &num_processes);
+    MPI_Comm_rank(MPI_COMM_WORLD, &process_id);
+
+    char processor_name[MPI_MAX_PROCESSOR_NAME];
+    int name_len;
+    MPI_Get_processor_name(processor_name, &name_len);
+    printf("Processor %s. Process %d of %d.\n", processor_name, process_id, num_processes);
+
+    AcMeshInfo mesh_info;
+    acLoadConfig(AC_DEFAULT_CONFIG, &mesh_info);
+
+    AcMesh* main_mesh     = NULL;
+    ModelMesh* model_mesh = NULL;
+    if (process_id == 0) {
+        main_mesh = acmesh_create(mesh_info);
+        acmesh_init_to(INIT_TYPE_RANDOM, main_mesh);
+        model_mesh = modelmesh_create(mesh_info);
+        acmesh_to_modelmesh(*main_mesh, model_mesh);
+    }
+
+    AcMeshInfo submesh_info = mesh_info;
+    submesh_info.int_params[AC_nz] /= num_processes;
+    update_config(&submesh_info);
+
+    AcMesh* submesh = acmesh_create(submesh_info);
+
+    /////////////////////
+    distribute_mesh(main_mesh, submesh);
+    gather_mesh(submesh, main_mesh);
+    /////////////////////////
+    // Autotest
+    bool is_acceptable = verify_meshes(*model_mesh, *main_mesh);
+    /////
+
+    acmesh_destroy(submesh);
+
+    if (process_id == 0) {
+        modelmesh_destroy(model_mesh);
+        acmesh_destroy(main_mesh);
+    }
+
+    MPI_Finalize();
+    return EXIT_SUCCESS;
+}

src/utils/CMakeLists.txt

@@ -6,5 +6,5 @@ set(CMAKE_C_STANDARD_REQUIRED ON)
 add_compile_options(-Wall -Wextra -Werror -Wdouble-promotion -Wfloat-conversion -Wshadow)
 
 ## Compile
-add_library(astaroth_utils STATIC config_loader.c memory.c)
+add_library(astaroth_utils STATIC config_loader.c memory.c verification.c)
-target_link_libraries(astaroth_utils PRIVATE astaroth_core)
+target_link_libraries(astaroth_utils PRIVATE astaroth_core m)

src/utils/verification.c

@@ -0,0 +1,157 @@
+#include "verification.h"
+
+#include <math.h>
+#include <stdio.h>
+
+#include "astaroth.h"
+
+#define max(a, b) ((a) > (b) ? (a) : (b))
+#define min(a, b) ((a) < (b) ? (a) : (b))
+
+#define fabs(x) ((_Generic((x), float : fabsf, double : fabs, long double : fabsl))(x))
+
+// Defines for colored output
+#define RED "\x1B[31m"
+#define GRN "\x1B[32m"
+#define YEL "\x1B[33m"
+#define BLU "\x1B[34m"
+#define MAG "\x1B[35m"
+#define CYN "\x1B[36m"
+#define WHT "\x1B[37m"
+#define RESET "\x1B[0m"
+
+typedef struct {
+    VertexBufferHandle handle;
+    AcReal model;
+    AcReal candidate;
+    long double abs_error;
+    long double ulp_error;
+    long double rel_error;
+    AcReal maximum_magnitude;
+    AcReal minimum_magnitude;
+} Error;
+
+static inline bool
+is_valid(const AcReal a)
+{
+    return !isnan(a) && !isinf(a);
+}
+
+static Error
+get_error(AcReal model, AcReal candidate)
+{
+    Error error;
+    error.abs_error = 0;
+
+    error.model     = model;
+    error.candidate = candidate;
+
+    if (error.model == error.candidate || fabsl(model - candidate) == 0) { // If exact
+        error.abs_error = 0;
+        error.rel_error = 0;
+        error.ulp_error = 0;
+    }
+    else if (!is_valid(error.model) || !is_valid(error.candidate)) {
+        error.abs_error = INFINITY;
+        error.rel_error = INFINITY;
+        error.ulp_error = INFINITY;
+    }
+    else {
+        const int base = 2;
+        const int p    = sizeof(AcReal) == 4 ? 24 : 53; // Bits in the significant
+
+        const long double e = floorl(logl(fabsl(error.model)) / logl(2));
+
+        const long double ulp             = powl(base, e - (p - 1));
+        const long double machine_epsilon = 0.5 * powl(base, -(p - 1));
+        error.abs_error                   = fabsl(model - candidate);
+        error.ulp_error                   = error.abs_error / ulp;
+        error.rel_error                   = fabsl(1.0l - candidate / model) / machine_epsilon;
+    }
+
+    return error;
+}
+
+static AcReal
+get_maximum_magnitude(const AcReal* field, const AcMeshInfo info)
+{
+    AcReal maximum = -INFINITY;
+
+    for (size_t i = 0; i < acVertexBufferSize(info); ++i)
+        maximum = max(maximum, fabs(field[i]));
+
+    return maximum;
+}
+
+static AcReal
+get_minimum_magnitude(const AcReal* field, const AcMeshInfo info)
+{
+    AcReal minimum = INFINITY;
+
+    for (size_t i = 0; i < acVertexBufferSize(info); ++i)
+        minimum = min(minimum, fabs(field[i]));
+
+    return minimum;
+}
+
+static Error
+get_max_abs_error(const VertexBufferHandle vtxbuf_handle, const AcMesh model_mesh,
+                  const AcMesh candidate_mesh)
+{
+    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 = get_error(model_vtxbuf[i], candidate_vtxbuf[i]);
+
+        if (curr_error.abs_error > error.abs_error)
+            error = curr_error;
+    }
+
+    error.handle            = vtxbuf_handle;
+    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 void
+print_error_to_screen(const Error error)
+{
+    const bool is_acceptable = true;
+
+    printf("\t%-15s... ", vtxbuf_names[error.handle]);
+    if (is_acceptable) {
+        printf(GRN "OK! " RESET);
+    }
+    else {
+        printf(RED "FAIL! " RESET);
+    }
+
+    fprintf(stdout, "| %.2Lg (abs), %.2Lg (ulps), %.2Lg (rel), [%.2g, %.2g] (range)\n", //
+            error.abs_error, error.ulp_error, error.rel_error,                          //
+            (double)error.minimum_magnitude, (double)error.maximum_magnitude);
+}
+
+bool
+acVerifyMesh(const AcMesh model, const AcMesh candidate)
+{
+    for (int i = 0; i < NUM_VTXBUF_HANDLES; ++i) {
+        Error field_error = get_max_abs_error(i, model, candidate);
+        print_error_to_screen(field_error);
+    }
+    return true;
+}

src/utils/verification.h

@@ -0,0 +1,6 @@
+#pragma once
+#include <stdbool.h>
+
+#include "memory.h"
+
+bool acVerifyMesh(const AcMesh model, const AcMesh candidate);