From 226de326513e2a1eba1137faa4f606f64f749028 Mon Sep 17 00:00:00 2001 From: jpekkila Date: Wed, 3 Jun 2020 13:37:00 +0300 Subject: [PATCH] Added model solution for reductions and functions for automated testing --- include/astaroth.h | 11 +++ include/astaroth_utils.h | 28 ++++++ samples/mpitest/main.cc | 14 +++ src/core/device.cc | 24 +++-- src/utils/CMakeLists.txt | 2 +- src/utils/modelreduce.c | 209 +++++++++++++++++++++++++++++++++++++++ src/utils/verification.c | 23 ++--- 7 files changed, 287 insertions(+), 24 deletions(-) create mode 100644 src/utils/modelreduce.c diff --git a/include/astaroth.h b/include/astaroth.h index b170796..47beb88 100644 --- a/include/astaroth.h +++ b/include/astaroth.h @@ -320,6 +320,17 @@ AcResult acGridIntegrate(const Stream stream, const AcReal dt); /** */ AcResult acGridPeriodicBoundconds(const Stream stream); + +/** TODO */ +AcResult +acGridReduceScal(const Stream stream, const ReductionType rtype, + const VertexBufferHandle vtxbuf_handle, AcReal* result); + +/** TODO */ +AcResult +acGridReduceVec(const Stream stream, const ReductionType rtype, + const VertexBufferHandle vtxbuf0, const VertexBufferHandle vtxbuf1, + const VertexBufferHandle vtxbuf2, AcReal* result); #endif // AC_MPI_ENABLED /* diff --git a/include/astaroth_utils.h b/include/astaroth_utils.h index 740e7f1..0422b7b 100644 --- a/include/astaroth_utils.h +++ b/include/astaroth_utils.h @@ -29,6 +29,25 @@ extern "C" { #endif + #include + +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; + +/** 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); @@ -56,6 +75,15 @@ 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); diff --git a/samples/mpitest/main.cc b/samples/mpitest/main.cc index 8d0a4fc..7b12fe2 100644 --- a/samples/mpitest/main.cc +++ b/samples/mpitest/main.cc @@ -53,6 +53,13 @@ main(void) acGridIntegrate(STREAM_DEFAULT, FLT_EPSILON); acGridPeriodicBoundconds(STREAM_DEFAULT); + // Do reductions + AcReal cand_reduce_res = 0; + VertexBufferHandle vtxbuf = VTXBUF_UUX; + ReductionType rtype = RTYPE_MAX; + acGridReduceScal(STREAM_DEFAULT, rtype, vtxbuf, &cand_reduce_res); // TODO + + acGridStoreMesh(STREAM_DEFAULT, &candidate); acGridQuit(); @@ -62,6 +69,13 @@ main(void) acMeshApplyPeriodicBounds(&model); acVerifyMesh(model, candidate); + + // Check reductions + AcReal model_reduce_res = acModelReduceScal(model, RTYPE_MAX, vtxbuf); + Error error = acGetError(model_reduce_res, cand_reduce_res); + error.handle = vtxbuf; + printErrorToScreen(error); + acMeshDestroy(&model); acMeshDestroy(&candidate); } diff --git a/src/core/device.cc b/src/core/device.cc index 4846a44..c82e80f 100644 --- a/src/core/device.cc +++ b/src/core/device.cc @@ -1641,8 +1641,8 @@ acMPIReduceScal(const AcReal local_result, const ReductionType rtype, AcReal* re MPI_Datatype datatype = MPI_DOUBLE; #else MPI_Datatype datatype = MPI_FLOAT; - #endif - + #endif + /* int rank; MPI_Comm_rank(MPI_COMM_WORLD, &rank); @@ -1664,31 +1664,43 @@ acMPIReduceScal(const AcReal local_result, const ReductionType rtype, AcReal* re } AcResult -acGridReduceScal(const Device device, const Stream stream, const ReductionType rtype, +acGridReduceScal(const Stream stream, const ReductionType rtype, const VertexBufferHandle vtxbuf_handle, AcReal* result) { + ERRCHK(grid.initialized); + // acGridSynchronizeStream(stream); + + const Device device = grid.device; + //const int3 nn = grid.nn; + acGridSynchronizeStream(STREAM_ALL); MPI_Barrier(MPI_COMM_WORLD); AcReal local_result; acDeviceReduceScal(device, stream, rtype, vtxbuf_handle, &local_result); - return acMPIReduceScal(local_result,rtype,result); + return acMPIReduceScal(local_result,rtype,result); } AcResult -acGridReduceVec(const Device device, const Stream stream, const ReductionType rtype, +acGridReduceVec(const Stream stream, const ReductionType rtype, const VertexBufferHandle vtxbuf0, const VertexBufferHandle vtxbuf1, const VertexBufferHandle vtxbuf2, AcReal* result) { + ERRCHK(grid.initialized); + // acGridSynchronizeStream(stream); + + const Device device = grid.device; + //const int3 nn = grid.nn; + acGridSynchronizeStream(STREAM_ALL); MPI_Barrier(MPI_COMM_WORLD); AcReal local_result; acDeviceReduceVec(device, stream, rtype, vtxbuf0, vtxbuf1, vtxbuf2, &local_result); - return acMPIReduceScal(local_result,rtype,result); + return acMPIReduceScal(local_result,rtype,result); } #endif // AC_MPI_ENABLED diff --git a/src/utils/CMakeLists.txt b/src/utils/CMakeLists.txt index 63e918c..47f1116 100644 --- a/src/utils/CMakeLists.txt +++ b/src/utils/CMakeLists.txt @@ -1,3 +1,3 @@ ## Astaroth Utils -add_library(astaroth_utils STATIC config_loader.c memory.c verification.c modelsolver.c) +add_library(astaroth_utils STATIC config_loader.c memory.c verification.c modelsolver.c modelreduce.c) add_dependencies(astaroth_utils dsl_headers) diff --git a/src/utils/modelreduce.c b/src/utils/modelreduce.c new file mode 100644 index 0000000..d95fc8c --- /dev/null +++ b/src/utils/modelreduce.c @@ -0,0 +1,209 @@ +/* + Copyright (C) 2014-2020, Johannes Pekkila, Miikka Vaisala. + + 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 . +*/ + +/** + * @file + * \brief Brief info. + * + * Detailed info. + * + */ + #include "astaroth.h" + +#include + +#include "errchk.h" + +#if AC_DOUBLE_PRECISION == 0 // HACK TODO fix, make cleaner (purkkaratkaisu) +#define fabs fabsf +#define exp expf +#define sqrt sqrtf +#endif + +// Function pointer definitions +typedef AcReal (*ReduceFunc)(const AcReal, const AcReal); +typedef AcReal (*ReduceInitialScalFunc)(const AcReal); +typedef AcReal (*ReduceInitialVecFunc)(const AcReal, const AcReal, + const AcReal); + +// clang-format off +/* Comparison funcs */ +static inline AcReal +max(const AcReal a, const AcReal b) { return a > b ? a : b; } + +static inline AcReal +min(const AcReal a, const AcReal b) { return a < b ? a : b; } + +static inline AcReal +sum(const AcReal a, const AcReal b) { return a + b; } + +/* Function used to determine the values used during reduction */ +static inline AcReal +length_scal(const AcReal a) { return (AcReal)(a); } + +static inline AcReal +length_vec(const AcReal a, const AcReal b, const AcReal c) { return sqrt(a*a + b*b + c*c); } + +static inline AcReal +squared_scal(const AcReal a) { return (AcReal)(a*a); } + +static inline AcReal +squared_vec(const AcReal a, const AcReal b, const AcReal c) { return squared_scal(a) + squared_scal(b) + squared_scal(c); } + +static inline AcReal +exp_squared_scal(const AcReal a) { return exp(a)*exp(a); } + +static inline AcReal +exp_squared_vec(const AcReal a, const AcReal b, const AcReal c) { return exp_squared_scal(a) + exp_squared_scal(b) + exp_squared_scal(c); } +// clang-format on + +AcReal +acModelReduceScal(const AcMesh mesh, const ReductionType rtype, const VertexBufferHandle a) +{ + ReduceInitialScalFunc reduce_initial; + ReduceFunc reduce; + + bool solve_mean = false; + + switch (rtype) { + case RTYPE_MAX: + reduce_initial = length_scal; + reduce = max; + break; + case RTYPE_MIN: + reduce_initial = length_scal; + reduce = min; + break; + case RTYPE_RMS: + reduce_initial = squared_scal; + reduce = sum; + solve_mean = true; + break; + case RTYPE_RMS_EXP: + reduce_initial = exp_squared_scal; + reduce = sum; + solve_mean = true; + break; + case RTYPE_SUM: + reduce_initial = length_scal; + reduce = sum; + break; + default: + ERROR("Unrecognized RTYPE"); + } + + const int initial_idx = acVertexBufferIdx(mesh.info.int_params[AC_nx_min], + mesh.info.int_params[AC_ny_min], + mesh.info.int_params[AC_nz_min], mesh.info); + + AcReal res; + if (rtype == RTYPE_MAX || rtype == RTYPE_MIN) + res = reduce_initial(mesh.vertex_buffer[a][initial_idx]); + else + res = 0; + + for (int k = mesh.info.int_params[AC_nz_min]; k < mesh.info.int_params[AC_nz_max]; ++k) { + for (int j = mesh.info.int_params[AC_ny_min]; j < mesh.info.int_params[AC_ny_max]; ++j) { + for (int i = mesh.info.int_params[AC_nx_min]; i < mesh.info.int_params[AC_nx_max]; + ++i) { + const int idx = acVertexBufferIdx(i, j, k, mesh.info); + const AcReal curr_val = reduce_initial(mesh.vertex_buffer[a][idx]); + res = reduce(res, curr_val); + } + } + } + + if (solve_mean) { + const AcReal inv_n = (AcReal)1.0 / mesh.info.int_params[AC_nxyz]; + return sqrt(inv_n * res); + } + else { + return res; + } +} + +AcReal +acModelReduceVec(const AcMesh mesh, const ReductionType rtype, const VertexBufferHandle a, + const VertexBufferHandle b, const VertexBufferHandle c) +{ + // AcReal (*reduce_initial)(AcReal, AcReal, AcReal); + ReduceInitialVecFunc reduce_initial; + ReduceFunc reduce; + + bool solve_mean = false; + + switch (rtype) { + case RTYPE_MAX: + reduce_initial = length_vec; + reduce = max; + break; + case RTYPE_MIN: + reduce_initial = length_vec; + reduce = min; + break; + case RTYPE_RMS: + reduce_initial = squared_vec; + reduce = sum; + solve_mean = true; + break; + case RTYPE_RMS_EXP: + reduce_initial = exp_squared_vec; + reduce = sum; + solve_mean = true; + break; + case RTYPE_SUM: + reduce_initial = length_vec; + reduce = sum; + break; + default: + ERROR("Unrecognized RTYPE"); + } + + const int initial_idx = acVertexBufferIdx(mesh.info.int_params[AC_nx_min], + mesh.info.int_params[AC_ny_min], + mesh.info.int_params[AC_nz_min], mesh.info); + + AcReal res; + if (rtype == RTYPE_MAX || rtype == RTYPE_MIN) + res = reduce_initial(mesh.vertex_buffer[a][initial_idx], mesh.vertex_buffer[b][initial_idx], + mesh.vertex_buffer[c][initial_idx]); + else + res = 0; + + for (int k = mesh.info.int_params[AC_nz_min]; k < mesh.info.int_params[AC_nz_max]; k++) { + for (int j = mesh.info.int_params[AC_ny_min]; j < mesh.info.int_params[AC_ny_max]; j++) { + for (int i = mesh.info.int_params[AC_nx_min]; i < mesh.info.int_params[AC_nx_max]; + i++) { + const int idx = acVertexBufferIdx(i, j, k, mesh.info); + const AcReal curr_val = reduce_initial(mesh.vertex_buffer[a][idx], + mesh.vertex_buffer[b][idx], + mesh.vertex_buffer[c][idx]); + res = reduce(res, curr_val); + } + } + } + + if (solve_mean) { + const AcReal inv_n = (AcReal)1.0 / mesh.info.int_params[AC_nxyz]; + return sqrt(inv_n * res); + } + else { + return res; + } +} diff --git a/src/utils/verification.c b/src/utils/verification.c index 3277660..a3ddb18 100644 --- a/src/utils/verification.c +++ b/src/utils/verification.c @@ -18,25 +18,14 @@ #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 +acGetError(AcReal model, AcReal candidate) { Error error; error.abs_error = 0; @@ -109,7 +98,7 @@ get_max_abs_error(const VertexBufferHandle vtxbuf_handle, const AcMesh model_mes for (size_t i = 0; i < acVertexBufferSize(model_mesh.info); ++i) { - Error curr_error = get_error(model_vtxbuf[i], candidate_vtxbuf[i]); + Error curr_error = acGetError(model_vtxbuf[i], candidate_vtxbuf[i]); if (curr_error.abs_error > error.abs_error) error = curr_error; @@ -147,8 +136,8 @@ is_acceptable(const Error error) return false; } -static bool -print_error_to_screen(const Error error) +bool +printErrorToScreen(const Error error) { bool errors_found = false; @@ -177,7 +166,7 @@ acVerifyMesh(const AcMesh model, const AcMesh candidate) bool errors_found = false; for (int i = 0; i < NUM_VTXBUF_HANDLES; ++i) { Error field_error = get_max_abs_error(i, model, candidate); - errors_found |= print_error_to_screen(field_error); + errors_found |= printErrorToScreen(field_error); } printf("%s\n", errors_found ? "Failure. Found errors in one or more vertex buffers"