diff --git a/.gitignore b/.gitignore index 7d833a7a..ac2c7228 100644 --- a/.gitignore +++ b/.gitignore @@ -7,6 +7,7 @@ /src/mole_C++/*.a /src/cpp/*.o /src/cpp/*.a +/cpp/build # CMake generated files CMakeCache.txt diff --git a/CMakeLists.txt b/CMakeLists.txt index 2a900e57..5cca56e4 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -2,7 +2,7 @@ cmake_minimum_required(VERSION 3.10) project(TopLevelProject VERSION 1.0 LANGUAGES CXX) -set(CMAKE_CXX_STANDARD 14) +set(CMAKE_CXX_STANDARD 17) set(CMAKE_CXX_STANDARD_REQUIRED True) # Display the detected C++ compiler ID @@ -147,10 +147,9 @@ set(LINK_LIBS ${ARMADILLO_LIBRARIES} ${LAPACK_LIBRARY}) # Add subdirectories -add_subdirectory(src/cpp) -add_subdirectory(tests/cpp) -add_subdirectory(tests/matlab_octave) -add_subdirectory(examples/cpp) +add_subdirectory(src) +add_subdirectory(tests) +add_subdirectory(examples) # Custom target to build everything add_custom_target(all_build DEPENDS mole_C++ tests_C++ examples_C++ tests_matlab_octave) diff --git a/cpp/CMakeLists.txt b/cpp/CMakeLists.txt new file mode 100644 index 00000000..4e02188f --- /dev/null +++ b/cpp/CMakeLists.txt @@ -0,0 +1,170 @@ +cmake_minimum_required(VERSION 3.10) +project(MOLE + VERSION 2.0.0 + DESCRIPTION "Mimetic Operators Library Enhanced (MOLE) - C++ implementation" + LANGUAGES CXX +) + +set(CMAKE_CXX_STANDARD 17) +set(CMAKE_CXX_STANDARD_REQUIRED True) +set(CMAKE_CXX_EXTENSIONS OFF) + +# Display the detected C++ compiler ID +message(STATUS "Detected CXX Compiler ID: ${CMAKE_CXX_COMPILER_ID}") + +# Compiler-specific CXX_FLAGS and linker flags +if (CMAKE_CXX_COMPILER_ID STREQUAL "AppleClang") + set(CMAKE_CXX_FLAGS "-O3 -Xclang -fopenmp -DARMA_DONT_USE_WRAPPER -DARMA_USE_SUPERLU") + set(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} -L/usr/local/opt/libomp/lib -L/opt/homebrew/opt/libomp/lib -lomp") + message(STATUS "Using AppleClang-specific flags.") + include_directories("/usr/local/opt/libomp/include" "/opt/homebrew/opt/libomp/include") + message(STATUS "Adding Eigen3 include directory for AppleClang.") + set(EIGEN3_INCLUDE_DIR "/opt/homebrew/include/eigen3") + include_directories(${EIGEN3_INCLUDE_DIR}) + +elseif (CMAKE_CXX_COMPILER_ID STREQUAL "IntelLLVM") + set(CMAKE_CXX_FLAGS "-O3 -qopenmp -DARMA_DONT_USE_WRAPPER -DARMA_USE_SUPERLU -diag-disable=10430") + set(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS}") + message(STATUS "Using non-Clang compiler flags.") + # Get MKLROOT from environment or fallback to default + if(DEFINED ENV{MKLROOT}) + set(MKLROOT $ENV{MKLROOT}) + else() + set(MKLROOT "/opt/intel/oneapi/mkl/latest") + endif() + # Automatically set Armadillo to link against MKL + set(BLAS_LIBRARIES "${MKLROOT}/lib/intel64/libmkl_rt.so" CACHE STRING "BLAS library path for MKL") + set(LAPACK_LIBRARIES "${MKLROOT}/lib/intel64/libmkl_rt.so" CACHE STRING "LAPACK library path for MKL") + set(ARMA_USE_WRAPPER OFF CACHE BOOL "Disable Armadillo wrapper to directly use MKL") + + message(STATUS "Using MKL from: ${MKLROOT}") + +else() + set(CMAKE_CXX_FLAGS "-O3 -fopenmp -DARMA_DONT_USE_WRAPPER -DARMA_USE_SUPERLU") + set(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS}") + message(STATUS "Using non-Clang compiler flags.") +endif() + +if(UNIX AND NOT APPLE AND NOT MSVC) + message(STATUS "Checking for Fortran compiler (gfortran)...") + enable_language(Fortran OPTIONAL) + if(NOT CMAKE_Fortran_COMPILER) + message(FATAL_ERROR "gfortran is required but was not found. Please install it using 'sudo apt install gfortran'.") + endif() +endif() + + +find_package(Eigen3 3.3.7 REQUIRED) +find_library(OpenBLAS_LIBRARIES NAMES openblas blas PATHS "/usr/lib/x86_64-linux-gnu" "/usr/local/opt/" "usr/local/lib" REQUIRED) +find_library(LAPACK_LIBRARY lapack REQUIRED PATHS "/usr/lib" "/usr/lib/x86_64-linux-gnu" "/usr/local/lib" "/usr/local/opt/") + +if(POLICY CMP0135) + cmake_policy(SET CMP0135 NEW) +endif() + +# Paths for downloading and building libraries +set(ARMADILLO_VERSION "15.2.2") +set(SUPERLU_VERSION "7.0.1") +set(BUILD_DIR "${CMAKE_BINARY_DIR}/third_party_build") +set(INSTALL_DIR "${CMAKE_BINARY_DIR}/third_party_install") + +# SuperLU configuration +set(SUPERLU_TARBALL_URL "https://github.com/xiaoyeli/superlu/archive/refs/tags/v${SUPERLU_VERSION}.tar.gz") +set(SUPERLU_SRC_DIR "${BUILD_DIR}/superlu-${SUPERLU_VERSION}") +set(SUPERLU_BUILD_DIR "${SUPERLU_SRC_DIR}/build") +set(SUPERLU_INSTALL_DIR "${INSTALL_DIR}/superlu-${SUPERLU_VERSION}") + +file(DOWNLOAD ${SUPERLU_TARBALL_URL} ${BUILD_DIR}/superlu-${SUPERLU_VERSION}.tar.gz) +execute_process(COMMAND ${CMAKE_COMMAND} -E tar xzf superlu-${SUPERLU_VERSION}.tar.gz WORKING_DIRECTORY ${BUILD_DIR}) +# Detect platform and apply the correct sed command +if(APPLE) + execute_process( + COMMAND sed -i "" "s/cmake_minimum_required(VERSION 2.8.12)/cmake_minimum_required(VERSION 3.10)/" CMakeLists.txt + WORKING_DIRECTORY ${SUPERLU_SRC_DIR} + ) +else() + execute_process( + COMMAND sed -i "s/cmake_minimum_required(VERSION 2.8.12)/cmake_minimum_required(VERSION 3.10)/" CMakeLists.txt + WORKING_DIRECTORY ${SUPERLU_SRC_DIR} + ) +endif() + +# Confirm the patch +execute_process( + COMMAND grep "cmake_minimum_required" ${SUPERLU_SRC_DIR}/CMakeLists.txt + OUTPUT_VARIABLE cmake_version_line + OUTPUT_STRIP_TRAILING_WHITESPACE +) +message(STATUS "SuperLU CMakeLists.txt after patch: ${cmake_version_line}") + +execute_process(COMMAND ${CMAKE_COMMAND} -E make_directory ${SUPERLU_BUILD_DIR}) + +execute_process(COMMAND ${CMAKE_COMMAND} .. + -DCMAKE_INSTALL_PREFIX=${SUPERLU_INSTALL_DIR} + -Denable_internal_blaslib=NO + -DTPL_BLAS_LIBRARIES=${OpenBLAS_LIBRARIES} + -DCMAKE_POSITION_INDEPENDENT_CODE=ON + -DCMAKE_C_FLAGS="-fPIC" + -DCMAKE_CXX_FLAGS="-fPIC" + WORKING_DIRECTORY ${SUPERLU_BUILD_DIR}) +execute_process(COMMAND make -j4 WORKING_DIRECTORY ${SUPERLU_BUILD_DIR}) +execute_process(COMMAND make install WORKING_DIRECTORY ${SUPERLU_BUILD_DIR}) + +# Armadillo configuration +set(ARMADILLO_TARBALL_URL "https://sourceforge.net/projects/arma/files/armadillo-${ARMADILLO_VERSION}.tar.xz") +set(ARMADILLO_BUILD_DIR "${BUILD_DIR}/armadillo-${ARMADILLO_VERSION}") +set(ARMADILLO_INSTALL_DIR "${INSTALL_DIR}/armadillo-${ARMADILLO_VERSION}") + +file(DOWNLOAD ${ARMADILLO_TARBALL_URL} ${BUILD_DIR}/armadillo-${ARMADILLO_VERSION}.tar.xz) +execute_process(COMMAND ${CMAKE_COMMAND} -E tar xJf armadillo-${ARMADILLO_VERSION}.tar.xz WORKING_DIRECTORY ${BUILD_DIR}) +execute_process(COMMAND ${CMAKE_COMMAND} + -DSuperLU_INCLUDE_DIR=${SUPERLU_INSTALL_DIR}/include + -DSuperLU_LIBRARY=${SUPERLU_INSTALL_DIR}/lib/libsuperlu.a + -DCMAKE_INSTALL_PREFIX=${ARMADILLO_INSTALL_DIR} . + WORKING_DIRECTORY ${ARMADILLO_BUILD_DIR}) +execute_process(COMMAND make -j4 WORKING_DIRECTORY ${ARMADILLO_BUILD_DIR}) +execute_process(COMMAND make install WORKING_DIRECTORY ${ARMADILLO_BUILD_DIR}) + +# Update paths for Armadillo and SuperLU +include_directories(${SUPERLU_INSTALL_DIR}/include ${ARMADILLO_INSTALL_DIR}/include) +link_directories(${SUPERLU_INSTALL_DIR}/lib ${ARMADILLO_INSTALL_DIR}/lib) + +# Find required libraries +set(CMAKE_PREFIX_PATH ${ARMADILLO_INSTALL_DIR} ${CMAKE_PREFIX_PATH}) +set(Armadillo_DIR ${ARMADILLO_INSTALL_DIR}/share/Armadillo) + +find_package(Armadillo REQUIRED) +if(NOT Armadillo_FOUND) + message(FATAL_ERROR "Custom Armadillo not found in ${ARMADILLO_INSTALL_DIR}") +else() + message(STATUS "Using Armadillo from ${ARMADILLO_INSTALL_DIR}") +endif() + +find_package(Eigen3 3.3.7 REQUIRED) +find_library(OpenBLAS_LIBRARIES NAMES openblas blas PATHS "/usr/lib/x86_64-linux-gnu" "/usr/local/opt/" REQUIRED) +find_library(LAPACK_LIBRARY lapack REQUIRED PATHS "/usr/lib" "/usr/lib/x86_64-linux-gnu" "/usr/local/lib" "/usr/local/opt/") + +# Required libraries and link settings +set(LINK_LIBS ${ARMADILLO_LIBRARIES} + ${OpenBLAS_LIBRARIES} + ${SUPERLU_INSTALL_DIR}/lib/libsuperlu.a + ${LAPACK_LIBRARY}) + +# Add subdirectories +add_subdirectory(src) +add_subdirectory(tests) +add_subdirectory(examples/grids) + +# Custom target to build everything +add_custom_target(all_build DEPENDS mole_C++ tests_C++ examples_C++ tests_matlab_octave) +# --------------------------------------------------------------- +# Summary printed at configure time +# --------------------------------------------------------------- +message(STATUS "") +message(STATUS "MOLE ${PROJECT_VERSION} configuration summary:") +message(STATUS " Build type: ${CMAKE_BUILD_TYPE}") +message(STATUS " Armadillo version: ${ARMADILLO_VERSION_STRING}") +message(STATUS " Build examples: ${MOLE_BUILD_EXAMPLES}") +message(STATUS " Build tests: ${MOLE_BUILD_TESTS}") +message(STATUS "") + diff --git a/cpp/README.md b/cpp/README.md new file mode 100644 index 00000000..40b9b6f5 --- /dev/null +++ b/cpp/README.md @@ -0,0 +1,52 @@ + + +# Top subdirectory for MOLE 2.0 C++ files + +Subdirectory and Pathname: **mole/cpp** + +## Purpose + +Top subdirectory for the MOLE 2.0 C++ implementation. +It contains the top CMakefile.list. + +## MOLE C++ 2.0 Directory structure + +```text +mole/ +├── cpp/ +| │── cmake/ +| │── doc/ +| │── examples +| |── src/ +| │ ├── boundaries +| │ └── grids +| │ └── include +| │ └── operators +| │ └── sys +| │ └── utils +| |── tests + +: +: __Other MOLE 2.0 and v1.2.0 files, including other language__ +: __implementations__ +: +├── .gitignore +└── README.md +``` + +## Other MOLE 2.0 C++ Files in this directory + ++ CMakeLists.txt which builds a MOLE 2.0 C++ version of the library +and all its functional modules ++ README.md (this file) diff --git a/cpp/cmake/MOLEConfig.cmake.in b/cpp/cmake/MOLEConfig.cmake.in new file mode 100644 index 00000000..ba8cfd63 --- /dev/null +++ b/cpp/cmake/MOLEConfig.cmake.in @@ -0,0 +1,5 @@ +@PACKAGE_INIT@ + +include("${CMAKE_CURRENT_LIST_DIR}/MOLETargets.cmake") + +check_required_components(MOLE) diff --git a/cpp/examples/README.md b/cpp/examples/README.md new file mode 100644 index 00000000..9d71c966 --- /dev/null +++ b/cpp/examples/README.md @@ -0,0 +1,45 @@ + + +# Subdirectory of Use Cases/Examples for the MOLE 2.0 Implementation + +Subdirectory and Pathname: **mole/cpp/examples** + +## Purpose + +This is the top subdirectory for the MOLE 2.0 C++ example +implementations. Examples are organized by the different modules +contained in the MOLE Library + +## Structure of the MOLE C++ 2.0 Subdirectory of Examples + +```text +mole/ +├── cpp/ +| │── examples +| │ ├── boundaries +| │ └── grids +| │ └── operators +| │ └── sys +| | └── time_integrators +| │ └── utils +``` + +## MOLE 2.0 C++ Files and Subdirectories in this subdirectory + ++ subdirectory **grids**: contains examples on how to declare and + construct MOLE grids ++ subdirectory **sys**: contains examples on how to use MOLE v2.0 + error handling mechanisms ++ file **README.md** this file diff --git a/cpp/examples/grids/CMakeLists.txt b/cpp/examples/grids/CMakeLists.txt new file mode 100644 index 00000000..0cdaabb7 --- /dev/null +++ b/cpp/examples/grids/CMakeLists.txt @@ -0,0 +1,35 @@ +cmake_minimum_required(VERSION 3.14) + +# --------------------------------------------------------------- +# The files in the CMakeLists.txt file are all examples of the +# the use of the MOLE Library +# --------------------------------------------------------------- +set(MOLE_EXAMPLE_SOURCES + grid1D_basic.cpp # basic 1D grid generation + grid2D_basic.cpp # basic 2D grid generation + err_grid2D_basic.cpp # basic 2D grid generation + gridBuilder1D_basic.cpp # 1D grid through gridBuilder + gridBuilder2D_basic.cpp # 2D grid with periodicity + gridBuilder_arg_order.cpp # attribute order independence + gridBuilder_error_handling.cpp # parse-time error reporting +) + +# Builds all the examples listed in MOLE_EXAMPLE_SOURCES +add_custom_target(mole_examples) + +foreach(_mole_example_src IN LISTS MOLE_EXAMPLE_SOURCES) + if(NOT EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/${_mole_example_src}") + message(FATAL_ERROR + "MOLE example source listed in examples/CMakeLists.txt " + "is missing: ${_mole_example_src}") + endif() + get_filename_component(_mole_example_name + ${_mole_example_src} NAME_WE) + add_executable(${_mole_example_name} ${_mole_example_src}) + target_link_libraries(${_mole_example_name} PRIVATE MOLE::mole) + target_compile_options(${_mole_example_name} PRIVATE -Wall -Wextra) + # Group all example binaries under a single "examples" target so + # `cmake --build . --target examples` builds just these, without + # also building the (potentially much larger) test suite. + add_dependencies(mole_examples ${_mole_example_name}) +endforeach() diff --git a/cpp/examples/grids/README.md b/cpp/examples/grids/README.md new file mode 100644 index 00000000..08b6d935 --- /dev/null +++ b/cpp/examples/grids/README.md @@ -0,0 +1,34 @@ + + +# Subdirectory for MOLE 2.0 C++ grid examples + +Subdirectory and Pathname: **mole/cpp/examples/grids/** + +## Purpose + +Subdirectory containing examples of MOLE grid declarations and +generation using MOLE grid classes and the MOLE gridBuilder + +## List of Files in This Subdirectory + ++ **CMakeLists.txt**: support the building of MOLE 2.0 C++ examples of +MOLE grids ++ **err_grid2D_basic.cpp**: a C++ example of how the MOLE error log +mechanisms works in the context of grids ++ **grid1D_basic.cpp**: a C++ example of a 1D grid declaration and +instantiation (checking for errors) ++ **grid2D_basic.cpp**:a C++ example of a 1D grid declaration and +instantiation (checking for errors) ++ **README.md**: (this file) diff --git a/cpp/examples/grids/err_grid2D_basic.cpp b/cpp/examples/grids/err_grid2D_basic.cpp new file mode 100644 index 00000000..e7c5d13a --- /dev/null +++ b/cpp/examples/grids/err_grid2D_basic.cpp @@ -0,0 +1,65 @@ +// This example generates some grid issues +#include "MOLE_grids.h" +#include + +int main() { + gridParams2D p; + p.topology = 'x'; + p.m = 3; p.n = 3; + p.dx = 1.0; p.dy = 1.0; + + cout << "===================================================== " + << endl; + cout << "Attempt to create GRID #1 (fails on invalid topology) " + << endl; + cout << "===================================================== " + << endl; + // This grid has an invalid topology + grid2D g(p); + if (!g.validGrid()) { + g.print_ErrorLog(); + } else { + std::cout << "grid2D built OK, nodes_X is " + << g.grid.nodes_X.data_.n_rows << " x " + << g.grid.nodes_X.data_.n_cols << "\n"; + } + + cout << "===================================================== " + << endl; + cout << "Attempt to create GRID #2 (fails on invalid spacing) " + << endl; + cout << "===================================================== " + << endl; + + // This second grid has an invalid dx + p.topology = 'u'; + p.dx = 0; + grid2D g1(p); + if (!g1.validGrid()) { + g1.print_ErrorLog(); + } else { + std::cout << "grid2D built OK, nodes_X is " + << g1.grid.nodes_X.data_.n_rows << " x " + << g1.grid.nodes_X.data_.n_cols << "\n"; + } + + cout << "======================================================= " + << endl; + cout << "Attempt to create GRID #3 (success - grid instantiated!) " + << endl; + cout << "======================================================= " + << endl; + // This third attempt the grid has valid parameters + p.topology = 'u'; + p.dx = 1.0; + grid2D g2(p); + if (!g2.validGrid()) { + g2.print_ErrorLog(); + } else { + std::cout << "grid2D built OK, nodes_X is " + << g2.grid.nodes_X.data_.n_rows << " x " + << g2.grid.nodes_X.data_.n_cols << "\n"; + } + + return 0; +} diff --git a/cpp/examples/grids/grid1D_basic.cpp b/cpp/examples/grids/grid1D_basic.cpp new file mode 100644 index 00000000..9965ca66 --- /dev/null +++ b/cpp/examples/grids/grid1D_basic.cpp @@ -0,0 +1,18 @@ +#include "MOLE_grids.h" +#include + +int main() { + gridParams1D p; + p.topology = 'u'; + p.m = 4; + p.dx = 0.5; + + grid1D g(p); + if (!g.validGrid()) { + g.print_ErrorLog(); + return 1; + } + std::cout << "grid built OK, " << g.grid.nodes_X.data_.n_elem + << " nodal points\n"; + return 0; +} diff --git a/cpp/examples/grids/grid2D_basic.cpp b/cpp/examples/grids/grid2D_basic.cpp new file mode 100644 index 00000000..c5c81f06 --- /dev/null +++ b/cpp/examples/grids/grid2D_basic.cpp @@ -0,0 +1,19 @@ +#include "MOLE_grids.h" +#include + +int main() { + gridParams2D p; + p.topology = 'u'; + p.m = 3; p.n = 3; + p.dx = 1.0; p.dy = 1.0; + + grid2D g(p); + if (!g.validGrid()) { + g.print_ErrorLog(); + return 1; + } + std::cout << "grid2D built OK, nodes_X is " + << g.grid.nodes_X.data_.n_rows << " x " + << g.grid.nodes_X.data_.n_cols << "\n"; + return 0; +} diff --git a/cpp/examples/grids/gridBuilder1D_basic.cpp b/cpp/examples/grids/gridBuilder1D_basic.cpp new file mode 100644 index 00000000..3eb8a55c --- /dev/null +++ b/cpp/examples/grids/gridBuilder1D_basic.cpp @@ -0,0 +1,23 @@ +// Builds a uniform 1D grid through the gridBuilder utility. +// +// gridBuilder takes pairs in any order and +// returns a gridVar. Compare with grid1D_basic.cpp, which fills a +// gridParams1D and calls the grid1D constructor directly. +#include "grid_builder.h" +#include + +int main() { + gridVar g = gridBuilder("dim", 1, "m", 4, "dx", 0.5, + "topology", 'u'); + + if (!isValidGrid(g)) { + std::visit([](auto&& grid) { grid.print_ErrorLog(); }, g); + return 1; + } + + // dim was 1, so the variant holds a grid1D. + grid1D& g1 = std::get(g); + std::cout << "grid built OK, " << g1.grid.nodes_X.data_.n_elem + << " nodal points\n"; + return 0; +} diff --git a/cpp/examples/grids/gridBuilder2D_basic.cpp b/cpp/examples/grids/gridBuilder2D_basic.cpp new file mode 100644 index 00000000..9c045ec7 --- /dev/null +++ b/cpp/examples/grids/gridBuilder2D_basic.cpp @@ -0,0 +1,35 @@ +// Builds a uniform 2D grid through gridBuilder, with periodicity +// declared on the x-axis only. +// +// isPeriodic is passed as the address of a vector holding one flag +// per dimension. The vector carries its own size, so gridBuilder +// checks it against dim; the pair may appear before "dim" in the +// argument list because both the check and the copy happen after +// parsing. The vector must outlive the gridBuilder call. +#include "grid_builder.h" +#include +#include + +int main() { + const std::vector periodic = {true, false}; + + gridVar g = gridBuilder("dim", 2, + "m", 3, "n", 3, + "dx", 1.0, "dy", 1.0, + "topology", 'u', + "isPeriodic", &periodic); + + if (!isValidGrid(g)) { + std::visit([](auto&& grid) { grid.print_ErrorLog(); }, g); + return 1; + } + + grid2D& g2 = std::get(g); + std::cout << "grid2D built OK, nodes_X is " + << g2.grid.nodes_X.data_.n_rows << " x " + << g2.grid.nodes_X.data_.n_cols << "\n"; + std::cout << "periodic in x: " << g2.grid.bc_isPeriodic[0] + << ", periodic in y: " << g2.grid.bc_isPeriodic[1] + << "\n"; + return 0; +} diff --git a/cpp/examples/grids/gridBuilder_arg_order.cpp b/cpp/examples/grids/gridBuilder_arg_order.cpp new file mode 100644 index 00000000..6135d1be --- /dev/null +++ b/cpp/examples/grids/gridBuilder_arg_order.cpp @@ -0,0 +1,70 @@ +// Attribute pairs may be given in any order. +// +// gridBuilder reads pairs until the nullptr sentinel, +// dispatching on the name, so position carries no meaning. The three +// calls below describe the same 3x3 uniform grid and produce grids +// that compare equal field by field. +// +// The one attribute with an ordering subtlety is isPeriodic: it is +// a pointer to the caller's vector, and dim may not have been read +// yet when the pair is seen. The parser stashes the pointer, then +// checks its size against dim and copies the flags after parsing, +// so "isPeriodic" may still precede "dim". +#include "grid_builder.h" +#include +#include + +// sameGrid compares the parameters two 2D grids were built with. +static bool sameGrid(const grid2D& a, const grid2D& b) { + return a.grid.topology == b.grid.topology + && a.grid.m == b.grid.m + && a.grid.n == b.grid.n + && a.grid.dx == b.grid.dx + && a.grid.dy == b.grid.dy + && a.grid.nodes_X == b.grid.nodes_X + && a.grid.nodes_Y == b.grid.nodes_Y + && a.grid.centers_X == b.grid.centers_X + && a.grid.centers_Y == b.grid.centers_Y + && a.grid.bc_isPeriodic[0] == b.grid.bc_isPeriodic[0] + && a.grid.bc_isPeriodic[1] == b.grid.bc_isPeriodic[1]; +} + +int main() { + const std::vector periodic = {true, false}; + + // Order A: dimensionality, counts, spacings, topology. + gridVar a = gridBuilder("dim", 2, + "m", 3, "n", 3, + "dx", 1.0, "dy", 1.0, + "topology", 'u', + "isPeriodic", &periodic); + + // Order B: the same pairs, reversed. + gridVar b = gridBuilder("isPeriodic", &periodic, + "topology", 'u', + "dy", 1.0, "dx", 1.0, + "n", 3, "m", 3, + "dim", 2); + + // Order C: pairs interleaved by axis rather than by kind. + gridVar c = gridBuilder("topology", 'u', + "m", 3, "dx", 1.0, + "n", 3, "dy", 1.0, + "isPeriodic", &periodic, + "dim", 2); + + if (!isValidGrid(a) || !isValidGrid(b) || !isValidGrid(c)) { + std::cout << "at least one grid failed to build\n"; + return 1; + } + + grid2D& ga = std::get(a); + grid2D& gb = std::get(b); + grid2D& gc = std::get(c); + + std::cout << "order B matches order A: " + << sameGrid(ga, gb) << "\n"; + std::cout << "order C matches order A: " + << sameGrid(ga, gc) << "\n"; + return 0; +} diff --git a/cpp/examples/grids/gridBuilder_error_handling.cpp b/cpp/examples/grids/gridBuilder_error_handling.cpp new file mode 100644 index 00000000..3896b2b2 --- /dev/null +++ b/cpp/examples/grids/gridBuilder_error_handling.cpp @@ -0,0 +1,92 @@ +// This example exercises gridBuilder's parse-time failures. +// +// gridBuilder collects every error it can find rather than stopping +// at the first, pushes MOLE_ERR_INVALID_GRID_ARGS on top of the +// stack, and returns a gridNull carrying the whole log. The caller +// checks isValidGrid and prints the log. +#include "grid_builder.h" +#include +#include + +// report prints the outcome of one gridBuilder call. +static void report(gridVar& g) { + if (isValidGrid(g)) { + std::visit([](auto&& grid) { + std::cout << "grid built OK, dim = " << grid.dim << "\n"; + }, g); + return; + } + std::cout << "holds gridNull: " + << std::holds_alternative(g) << "\n"; + std::visit([](auto&& grid) { grid.print_ErrorLog(); }, g); +} + +int main() { + cout << "====================================================" + << endl; + cout << "GRID #1 (fails on an unknown attribute name) " << endl; + cout << "====================================================" + << endl; + // Parsing stops at the first key that is not a grid attribute, + // because the value type after an unknown key is unknown too. + gridVar g1 = gridBuilder("dim", 1, "m", 5, "spacing", 0.2, + "topology", 'u'); + report(g1); + + cout << "====================================================" + << endl; + cout << "GRID #2 (fails on a missing dim attribute) " << endl; + cout << "====================================================" + << endl; + // dim is required and is never inferred from m, n or o. + gridVar g2 = gridBuilder("m", 5, "dx", 0.2, "topology", 'u'); + report(g2); + + cout << "====================================================" + << endl; + cout << "GRID #3 (accumulates four separate errors) " << endl; + cout << "====================================================" + << endl; + // A 3D grid with no cell counts and a bad topology: three + // missing counts plus the topology, all reported at once. + gridVar g3 = gridBuilder("dim", 3, "topology", 'x'); + report(g3); + + cout << "====================================================" + << endl; + cout << "GRID #4 (fails on a count the dimension cannot use) " + << endl; + cout << "====================================================" + << endl; + // o belongs to a 3D grid only, so supplying it for a 2D grid + // is a cell-count inconsistency. + gridVar g4 = gridBuilder("dim", 2, "m", 5, "n", 5, "o", 5, + "dx", 0.2, "dy", 0.2, "topology", 'u'); + report(g4); + + cout << "====================================================" + << endl; + cout << "GRID #5 (fails on an isPeriodic size mismatch) " + << endl; + cout << "====================================================" + << endl; + // isPeriodic holds one flag per dimension. The vector reports + // its own size, so a 3D grid given two flags is caught rather + // than reading past the end of the vector. + const std::vector tooFew = {true, false}; + gridVar g5 = gridBuilder("dim", 3, "m", 5, "n", 5, "o", 5, + "dx", 0.2, "dy", 0.2, "dz", 0.2, + "topology", 'u', + "isPeriodic", &tooFew); + report(g5); + + cout << "====================================================" + << endl; + cout << "GRID #6 (success - grid instantiated!) " << endl; + cout << "====================================================" + << endl; + gridVar g6 = gridBuilder("dim", 2, "m", 5, "n", 5, + "dx", 0.2, "dy", 0.2, "topology", 'u'); + report(g6); + return 0; +} diff --git a/cpp/src/CMakeLists.txt b/cpp/src/CMakeLists.txt new file mode 100644 index 00000000..d9a8117c --- /dev/null +++ b/cpp/src/CMakeLists.txt @@ -0,0 +1,42 @@ +cmake_minimum_required(VERSION 3.14) + +# --------------------------------------------------------------- +# The MOLE library itself: grids, arrays, error handling, utils. +# -------------------------------------------------------------- +add_library(mole + grids/MOLE_arrays.cpp + grids/MOLE_grids.cpp + grids/grid_builder.cpp + sys/MOLE_Errors.cpp + utils/utils.cpp +) +add_library(MOLE::mole ALIAS mole) + +target_include_directories(mole + PUBLIC + $ + $ +) + +# target_link_libraries(mole PUBLIC ${ARMADILLO_LIBRARIES}) +# target_include_directories(mole PUBLIC ${ARMADILLO_INCLUDE_DIRS}) + +if(MOLE_ENABLE_WARNINGS) + target_compile_options(mole PRIVATE -Wall -Wextra) +endif() + +set_target_properties(mole PROPERTIES + VERSION ${PROJECT_VERSION} + SOVERSION ${PROJECT_VERSION_MAJOR} +) + +install(TARGETS mole + EXPORT MOLETargets + LIBRARY DESTINATION lib + ARCHIVE DESTINATION lib + RUNTIME DESTINATION bin +) + +install(DIRECTORY include/ DESTINATION include + FILES_MATCHING PATTERN "*.h" +) diff --git a/cpp/src/README.md b/cpp/src/README.md new file mode 100644 index 00000000..57b18de8 --- /dev/null +++ b/cpp/src/README.md @@ -0,0 +1,48 @@ + + +# Subdirectory of Source Files for the MOLE 2.0 C++ Implementation + +Subdirectory and Pathname: **mole/cpp/src** + +## Purpose + +This is the top subdirectory for the MOLE 2.0 C++ files containing +the source code implementations of MOLE's core functionalities. The +declarations of public functions implemented in this subdirectory are +in a header files inside the mole/cpp/src/include subdirectory. + +## Structure of the MOLE C++ 2.0 Source Implementations' Subdirectory + +```text +mole/ +├── cpp/ +| │── src +| │ ├── boundaries +| │ └── grids +| │ └── include +| │ └── operators +| | └── sys +| │ └── utils +``` + +## MOLE 2.0 C++ Files and Subdirectories in this subdirectory + ++ **CMakeLists.txt**: CMake file for MOLE sources ++ subdirectory **grids**: contains functional implementations of the + MOLE grid classes and required data structures. ++ subdirectory **operators**: contains functional implementations of + the MOLE operators, including grid interpolators. ++ subdirectory **sys**: contains functional implementations of MOLE + computational support like the error handling mechanisms. ++ file **README.md**: this file diff --git a/cpp/src/grids/MOLE_arrays.cpp b/cpp/src/grids/MOLE_arrays.cpp new file mode 100644 index 00000000..5335e174 --- /dev/null +++ b/cpp/src/grids/MOLE_arrays.cpp @@ -0,0 +1,448 @@ +/* + * SPDX-License-Identifier: GPL-3.0-or-later + * Copyright (c) 2008-2024 San Diego State University Research + * Foundation (SDSURF). + * See LICENSE file or https://www.gnu.org/licenses/gpl-3.0.html + * for details. + */ + +/* + * @file MOLE_arrays.cpp + * + * @brief MOLE Arrays Classes with member function implementations + * + * @date 2026/07/14 + * + */ +#include "MOLE_arrays.h" + +// ----------- +// +// MOLE Arrays Constructors +// +// ----------- +// + + +// array1D::array1D constructor(numelem, fillVal) creates a 1D array +// of Reals of size numelem and fills with fillVal. +// It also records errors during memory allocation or other vector +// operations. +// +array1D::array1D(size_t numelem, Real fillVal) { + try { + data_.set_size(numelem); + data_.fill(fillVal); + } catch (std::bad_alloc& e) { + string wparams = "numelem = " + to_string(numelem); + logArrayError(MOLE_ERR_FAILED_ARRAY_ALLOC, + "array1D Construction", wparams); + } +} + +// +// array2D::array2D constructor(rows, cols, fillVal) creates a 2D +// array of Reals of size rows x cols and fills with fillVal. +// It also records errors during memory allocation or other vector +// operations. +// +array2D::array2D(size_t rows, size_t cols, Real fillVal) { + try { + data_.set_size(rows, cols); + data_.fill(fillVal); + } catch (std::bad_alloc& e) { + string wparams = "rows = " + to_string(rows); + wparams += ", cols = " + to_string(cols); + logArrayError(MOLE_ERR_FAILED_ARRAY_ALLOC, + "array2D Construction", wparams); + } +} + +// +// array3D::array3D constructor(dim1, dim2, dim3, fillVal) creates a 3D array +// of Reals of size rows x cols x depth and fills with fillVal. +// It also records errors during memory allocation or other vector +// operations. +// +array3D::array3D(size_t rows, size_t cols, size_t slices, Real fillVal) { + try { + data_.set_size(rows, cols, slices); + data_.fill(fillVal); + } catch (std::bad_alloc& e) { + string wparams = "rows = " + to_string(rows); + wparams += ", cols = " + to_string(cols); + wparams += ", slices = " + to_string(slices); + logArrayError(MOLE_ERR_FAILED_ARRAY_ALLOC, + "array3D Construction", wparams); + } +} + +// ----------- +// +// Error handling methods for MOLE Arrays +// +// ----------- + +// +// array1D::logArrayErr records errors in the error stack for 1D +//arrays. +// +void array1D::logArrayError(size_t errCode, string errLoc, + string errParm) const { + MOLEerr_log(a_errs, errCode, errLoc, errParm); +} + +// +// array1D::hasArrayErrors check whether are issues or +// errors with the flat2Darray +// +bool array1D::hasArrayErrors() const { + return MOLEerr_haserrors(a_errs); +} + +// +// array1D::print_ErrorLog prints out errors with the array1D +// to standard output +// +void array1D::print_ErrorLog() const { + MOLEerr_print(a_errs); +} + +// +// array1D::write_ErrorLog prints out errors with the array1D +// to an output file with name starting with MOLEArrayErrors - the +// full name of the file also includes a timestamp +// +void array1D::write_ErrorLog() const { + MOLEerr_dumpErrLog(a_errs, "MOLEArrayErrors"); +} + +// +// array2D::logArrayErr records errors in the error stack for 1D +//arrays. +// +void array2D::logArrayError(size_t errCode, string errLoc, + string errParm) const { + MOLEerr_log(a_errs, errCode, errLoc, errParm); +} + +// +// array2D::hasArrayErrors check whether are issues or +// errors with the flat2Darray +// +bool array2D::hasArrayErrors() const { + return MOLEerr_haserrors(a_errs); +} + +// +// array2D::print_ErrorLog prints out errors with the array2D +// to standard output +// +void array2D::print_ErrorLog() const { + MOLEerr_print(a_errs); +} + +// +// array2D::write_ErrorLog prints out errors with the array2D +// to an output file with name starting with MOLEArrayErrors - the +// full name of the file also includes a timestamp +// +void array2D::write_ErrorLog() const { + MOLEerr_dumpErrLog(a_errs, "MOLEArrayErrors"); +} + +// +// array3D::logArrayErr records errors in the error stack for 1D +//arrays. +// +void array3D::logArrayError(size_t errCode, string errLoc, + string errParm) const { + MOLEerr_log(a_errs, errCode, errLoc, errParm); +} + +// +// array3D::hasArrayErrors check whether are issues or +// errors with the flat3Darray +// +bool array3D::hasArrayErrors() const { + return MOLEerr_haserrors(a_errs); +} + +// +// array3D::print_ErrorLog prints out errors with the array3D +// to standard output +// +void array3D::print_ErrorLog() const { + MOLEerr_print(a_errs); +} + +// +// array3D::write_ErrorLog prints out errors with the array3D +// to an output file with name starting with MOLEArrayErrors - the +// full name of the file also includes a timestamp +// +void array3D::write_ErrorLog() const { + MOLEerr_dumpErrLog(a_errs, "MOLEArrayErrors"); +} + +// +// read_ErrorLog methods read the error on top of the stack of an +// array data structure. Its intent is to propagate the error up to a +// user facing class like grids or operators. Users can also print the +// errors to standard output or write them to a file. +// +void array1D::read_ErrorLog(int ErrorCode, std::string &location, + std::string &arrayName) { + if (!a_errs.empty()) { + MOLE_Errors topError = a_errs.top(); + ErrorCode = topError.errCode; + location = topError.errLocation; + arrayName = topError.paramError; + a_errs.pop(); // Remove the top error after reading + } +} + +void array2D::read_ErrorLog(int ErrorCode, std::string &location, + std::string &arrayName) { + if (!a_errs.empty()) { + MOLE_Errors topError = a_errs.top(); + ErrorCode = topError.errCode; + location = topError.errLocation; + arrayName = topError.paramError; + a_errs.pop(); // Remove the top error after reading + } +} + +void array3D::read_ErrorLog(int ErrorCode, std::string &location, + std::string &arrayName) { + if (!a_errs.empty()) { + MOLE_Errors topError = a_errs.top(); + ErrorCode = topError.errCode; + location = topError.errLocation; + arrayName = topError.paramError; + a_errs.pop(); // Remove the top error after reading + } +} +// ------------- +// +// General MOLE Array functions (for flat dense arrays) +// When wrapping other libraries, these functions will need to be +// redefined to work with the other library's data structures. +// 1) validate that an array index is valid (i.e., within bounds) +// 2) array equality, two arrays are == if their dimensions and data +// are the same +// 3) array inequality, two arrays are != if they are stored at +// different memory locations +// 4) array resize, resize an array to a new size, discards old +// memory and allocates new memory (it is not an array reshape). +// +// ------------- + +// +// 1D: valid_index checks whether the i index is valid +// +bool array1D::valid_index(size_t i) const { + return i < data_.n_elem; +} + +// +// 2D: valid_indeces checks whether the i, j indeces are valid +// +bool array2D::valid_indeces(size_t i, size_t j) const { + return i < data_.n_rows && j < data_.n_cols; +} + +// +// 3D: valid_indeces checks whether the i, j, k indeces are valid +// +bool array3D::valid_indeces(size_t i, size_t j, size_t k) const { + return i < data_.n_rows && j < data_.n_cols && k < data_.n_slices; +} + +// +// Array1D::operator == Array equality comparison (dimensions +// and data content) Not using Armadillo's operator == to provide a +// general template for other numerical libraries. +// +bool array1D::operator==(const array1D& other) const { + bool are_equal = true; + if (data_.n_elem != other.data_.n_elem) { // compare sizes + are_equal = false; + } + else if (data_.memptr() != other.data_.memptr()) { // same ptr? + const double* __restrict c = data_.memptr(); // fast compare + const double* __restrict u = other.data_.memptr(); + for (size_t i = 0; i < data_.n_elem; ++i) { + if (c[i] != u[i]) { + are_equal = false; + break; + } + } + } + return are_equal; +} + +// +// Array1D::operator != comparing if these are different arrays (i.e, +// are these arrrays two different memory locations) +// +bool array1D::operator!=(const array1D& other) const { + return (data_.memptr() != other.data_.memptr()); +} + +// +// Array2D::operator == Array equality comparison (dimensions +// and data content) Not using Armadillo's operator == to provide a +// general template for other numerical libraries. +// +bool array2D::operator==(const array2D& other) const { + bool are_equal = true; + if (data_.n_rows != other.data_.n_rows || + data_.n_cols != other.data_.n_cols) { // compare sizes + are_equal = false; + } + else if (data_.memptr() != other.data_.memptr()) { // same ptr? + const double* __restrict c = data_.memptr(); // fast compare + const double* __restrict u = other.data_.memptr(); + for (size_t i = 0; i < data_.n_rows*data_.n_cols; ++i) { + if (c[i] != u[i]) { + are_equal = false; + break; + } + } + } + return are_equal; +} + +// +// Array2D::operator != comparing if these are different arrays (i.e, +// are these arrrays two different memory locations or they don't have +// the same shape) +// +bool array2D::operator!=(const array2D& other) const { + bool are_not_equal = false; + if (data_.n_rows != other.data_.n_rows || + data_.n_cols != other.data_.n_cols) { // compare sizes + are_not_equal = true; + } + else if (data_.memptr() != other.data_.memptr()) { // same ptr? + are_not_equal = true; + } + return are_not_equal; +} + +// +// Array3D::operator == Array equality comparison (dimensions +// and data content) Not using Armadillo's operator == to provide a +// general template for other numerical libraries. +// +bool array3D::operator==(const array3D& other) const { + bool are_equal = true; + if (data_.n_rows != other.data_.n_rows || + data_.n_cols != other.data_.n_cols || + data_.n_slices != other.data_.n_slices) { // compare sizes + are_equal = false; + } + else if (data_.memptr() != other.data_.memptr()) { // same ptr? + const double* __restrict c = data_.memptr(); // fast compare + const double* __restrict u = other.data_.memptr(); + for (size_t i = 0; + i < data_.n_rows*data_.n_cols*data_.n_slices; ++i) { + if (c[i] != u[i]) { + are_equal = false; + break; + } + } + } + return are_equal; +} + +// +// Array3D::operator != comparing if these are different arrays (i.e, +// are these arrrays two different memory locations or they don't have +// the same shape) +// +bool array3D::operator!=(const array3D& other) const { + bool are_not_equal = false; + if (data_.n_rows != other.data_.n_rows || + data_.n_cols != other.data_.n_cols || + data_.n_slices != other.data_.n_slices) { // compare sizes + are_not_equal = true; + } + else if (data_.memptr() != other.data_.memptr()) { // same ptr? + are_not_equal = true; + } + return are_not_equal; +} + +// +// array1D::resize class method to resize a array1D. This method +// guards this operation against an overflow. If an overflow occurs, +// an error is logged + the array size is set to 0 (the resize is +// simply rejected) rather than returning the wrong object. +// +void array1D::resize(size_t numelem, Real fillVal) { + size_t c_size = data_.n_elem; + + try { + data_.set_size(numelem); + data_.fill(fillVal); + } catch (std::bad_alloc& e) { + string wparams = "original numelem = " + to_string(c_size); + wparams += ", requested resize = " + to_string(numelem); + logArrayError(MOLE_ERR_FAILED_ARRAY_RESIZE, + "array1D Resize", wparams); + data_.set_size(0); + } +} + +// +// array2D::resize class method to resize a array2D. This method +// guards this operation against an overflow. If an overflow occurs, +// an error is logged + the array size is set to 0 (the resize is +// simply rejected) rather than returning the wrong object. +// +void array2D::resize(size_t rows, size_t cols, Real fillVal) { + size_t c_rows = data_.n_rows; + size_t c_cols = data_.n_cols; + + try { + data_.set_size(rows, cols); + data_.fill(fillVal); + } catch (std::bad_alloc& e) { + string wparams = "original rows = " + to_string(c_rows); + wparams += ", original cols = " + to_string(c_cols); + wparams += ", requested rows = " + to_string(rows); + wparams += ", requested cols = " + to_string(cols); + logArrayError(MOLE_ERR_FAILED_ARRAY_RESIZE, + "array2D Resize", wparams); + data_.set_size(0, 0); + } +} + +// +// array3D::resize class method to resize a array3D. This method +// guards this operation against an overflow. If an overflow occurs, +// an error is logged + the array size is set to 0 (the resize is +// simply rejected) rather than returning the wrong object. +// +void array3D::resize(size_t rows, size_t cols, size_t slices, + Real fillVal) { + size_t c_rows = data_.n_rows; + size_t c_cols = data_.n_cols; + size_t c_slices = data_.n_slices; + try { + data_.set_size(rows, cols, slices); + data_.fill(fillVal); + } catch (std::bad_alloc& e) { + string wparams = "original rows = " + to_string(c_rows); + wparams += ", original cols = " + to_string(c_cols); + wparams += ", original slices = " + to_string(c_slices); + wparams += ", requested rows = " + to_string(rows); + wparams += ", requested cols = " + to_string(cols); + wparams += ", requested slices = " + to_string(slices); + logArrayError(MOLE_ERR_FAILED_ARRAY_RESIZE, + "array3D Resize", wparams); + data_.set_size(0, 0, 0); + } +} diff --git a/cpp/src/grids/MOLE_grids.cpp b/cpp/src/grids/MOLE_grids.cpp new file mode 100644 index 00000000..f08caf78 --- /dev/null +++ b/cpp/src/grids/MOLE_grids.cpp @@ -0,0 +1,1147 @@ +/* +* SPDX-License-Identifier: GPL-3.0-or-later +* © 2008-2024 San Diego State University Research Foundation (SDSURF). +* See LICENSE file or https://www.gnu.org/licenses/gpl-3.0.html for details. +*/ + +/* + * @file MOLE_grids.cpp + * + * @brief MOLE Grid Implementations + * + * @date 2026/06/24 + * + */ + +#include "MOLE_grids.h" + +#include // abort() for DEBUG_AND_ABORT_MD + +// ------------------------------------------------------------------ +// gridBase Errors Implementation +// Wrappers to the MOLE error handling mechanism that can access +// the gridBase class private error stack +// ------------------------------------------------------------------ +// +// gridBase::logGridErr logs errors for all Grid classes +// +void gridBase::logGridErr(size_t errCode, string errLoc, + string errParm){ + MOLEerr_log(errs, errCode, errLoc, errParm); +} + +// +// gridBase::reportErrors prints out all generated errors +// +bool gridBase::hasGridErrors(){ + return MOLEerr_haserrors(errs); +} + +// +// gridBase::isValidatedGrid checks if a grid has been validated. +// +bool gridBase::isValidatedGrid(){ + // Only checks whether validGrid() has been previously called + return !MOLEerr_contains(errs, MOLE_ERR_GRID_UNCHECKED); +} + +// +// gridBase::setGridValidated is called afer a grid has been +// throroughly validated and sets internal flags with the results. +// +void gridBase::setGridValidated(){ + MOLEerr_remove(errs, MOLE_ERR_GRID_UNCHECKED); +} + +// +// gridBase::print_ErrorLog outputs the contents of a grid's +// errorlog stack to the standard output +// +void gridBase::print_ErrorLog(){ + MOLEerr_print(errs); +} + +// +// gridBase::write_ErrorLog outputs the contents of a grid's +// errorlog stack to a logfile. The name of the output file starts +// with "MOLEGridErrors" and is follow by the timestamp +// +void gridBase::write_ErrorLog(){ + MOLEerr_dumpErrLog(errs, "MOLEGridErrors"); +} + +// +// gridBase::mergeErrors propagates previously-logged errors into +// a grid's error log stack. +// +void gridBase::mergeErrors(const stack& inerrs) { + stack tmp_stk = inerrs; + while (!tmp_stk.empty()) { + logGridErr(tmp_stk.top().errCode, tmp_stk.top().errLocation, + tmp_stk.top().paramError); + tmp_stk.pop(); + } +} + +// +// gridBase::applyDebugMode applies a MOLE debug mode to a grid that +// failed validation. A grid that validated is left untouched. The +// modes are declared in MOLE_errors.h. An unrecognized mode falls +// back to reporting, which is the behaviour that loses the least +// information without ending the user's program. +// +void gridBase::applyDebugMode(size_t debug_mode){ + if (isValidatedGrid()) return; + + switch (debug_mode) { + case DEBUG_DEFAULT_MD: + return; + case DEBUG_REPORTS_STDOUT_MD: + print_ErrorLog(); + return; + case DEBUG_AND_ABORT_MD: + print_ErrorLog(); + abort(); + default: + cout << "Unrecognized MOLE debug mode [" << debug_mode + << "], reporting to standard output." << endl; + print_ErrorLog(); + return; + } +} + +// ---------- +// Set of Auxiliar functions that check for consistency +// in the grid parameters. These are not particular to a grid +// dimensionality nor topology +// ---------- + +// +// validSpacing checks for a valid cell spacing is valid +// +bool validSpacing(Real dh){ + if ( isnan(dh) ) return false; + if ( dh <= 0.0 ) return false; + if ( !isfinite(dh) ) return false; + return true; +} + +// +// generateNodalPts generate arrays used in the generation and +// validation of nodal or normal faces coordinates +// +void generateNodalPts(size_t npts, Real delta, array1D& out_array){ + for (size_t i = 0; i <= npts; ++i) { + out_array.data_[i] = (Real)i * delta; + } +} +// +// generateCenterPts generates cell centered grid coordinates used in +// generation and validation of grid center coordinates +// +void generateCenterPts(size_t npts, Real delta, array1D& out_array){ + // center_X[0] = 0.0, centers_X[1:npts] = (i-0.5)*delta + for (size_t i = 1; i <= npts; ++i) { + out_array.data_[i] = ((Real)i - 0.5) * delta; + } + out_array.data_[npts+1] = (Real)npts * delta; // center_X[npts+1]=npnts*dx +} + +// +// Returns a string with the grid topology name +// +string string_topology (char topology){ + if (topology == 'u') return("Uniform Grid: "); + else if (topology == 'c') return("Curvilinear Grid: "); + else if (topology == 'n') return("Nonuniform Grid"); + else return ("Invalid Topolgy"); +} + +// ---------------------------------------------------------------- +// +// IMPlEMENTATION OF MOLE GRID CLASSES AND RELATED FUNCTIONS +// +// ---------------------------------------------------------------- +// +// gridBase Class +// Default gridBase Constructor (creates only a i grid shell) and +// initializes the error_log for grid validation with and error +// handling. The grid is marked as unvalidated until an explicit +// required call to validGrid() is made, preceding MOLE operations. +gridBase::gridBase(size_t idim){ + MOLEerr_init(errs); // initialized the stack of errors + MOLE_Errors err; + err.errCode = MOLE_ERR_GRID_UNCHECKED; // push error to grid stack + err.errLocation = ""; + err.paramError = ""; + logGridErr(err.errCode, err.errLocation, err.paramError); + if (idim >= 1 && idim <= 3){ // validate grid dimensionality + dim = idim; + } else { + dim = 0; + logGridErr(MOLE_ERR_INVALID_GRID_DIM, + "gridBase declaration", to_string(idim)); + } +} + +// +// valid1DCoordinates validates a user-supplied coordinate array +// against an expected valid grid coordinate one. If the user did not +// provide an array of coordinates, and the grid is uniform, the +// expected value is assigned. +// +bool gridBase::valid1DCoordinates(array1D& userInput, + const array1D& expected, + Real dx, size_t m, + int sizeMismatchErr, + int badCoordsErr) { + if (userInput.data_.is_empty()) { + userInput = expected; // auto-generate + return true; + } + if (userInput.data_.n_elem != expected.data_.n_elem) { + logGridErr(sizeMismatchErr, "grid1D[construct]", + to_string(userInput.data_.n_elem)); + return false; + } + if (!numEqualArray(expected, userInput, 4.0)) {//eps*4.0 precision + string sparams = "dx = " + to_string(dx) + ", m = "; + sparams += to_string(m); + logGridErr(badCoordsErr, "grid1D[construct]", sparams); + return false; + } + return true; +} + +// +// valid2DCoordinates validates a user-supplied coordinate array +// against an expected valid grid coordinate one. If users did not +// provide an array of coordinates, and the grid is uniform, the +// grid is auto-generated. +// +bool gridBase::valid2DCoordinates(array2D& userInput, + const array2D& expected, + Real dx, Real dy, size_t m, size_t n, + int sizeMismatchErr, + int badCoordsErr) { + if (userInput.data_.is_empty()) { + userInput = expected; // auto-generate + return true; + } + if (userInput.data_.n_rows != expected.data_.n_rows || + userInput.data_.n_cols != expected.data_.n_cols) { + string sparams = "m = " + to_string(m); + sparams += ", n = " + to_string(n) + ", dx = "; + sparams += to_string(dx) + ", dy = " + to_string(dy); + logGridErr(sizeMismatchErr, "grid2D[construct]", sparams); + return false; + } + if (!numEqualArray(expected, userInput, 4.0)) {//eps*4.0 precision + string sparams = "m = " + to_string(m); + sparams += ", n = " + to_string(n) + ", dx = "; + sparams += to_string(dx) + ", dy = " + to_string(dy); + logGridErr(badCoordsErr, "grid2D[construct]", sparams); + return false; + } + return true; +} + +// +// buildOrCheck2DCoords implementation (see MOLE_grids.h for details). +// +bool gridBase::buildOrCheck2DCoords(array2D& outX, array2D& outY, + const array1D& xcoord, + const array1D& ycoord, Real dx, + Real dy, size_t m, size_t n, + int sizeMismatchErr, + int badCoordsErr) { + array2D X(xcoord.data_.n_elem, ycoord.data_.n_elem, 0.0); + array2D Y(xcoord.data_.n_elem, ycoord.data_.n_elem, 0.0); + + if (X.data_.n_rows != xcoord.data_.n_elem || + X.data_.n_cols != ycoord.data_.n_elem || + Y.data_.n_rows != xcoord.data_.n_elem || + Y.data_.n_cols != ycoord.data_.n_elem) { + string errmsg = "dim1 = " + to_string(xcoord.data_.n_elem); + errmsg += " X dim2 = " + to_string(ycoord.data_.n_elem); + logGridErr(MOLE_ERR_FAILED_ARRAY_ALLOC, + "gridBase[buildOrCheck2DCoords]", errmsg); + return false; + } + + nd2DGrid(xcoord, ycoord, X, Y); + if (X.hasArrayErrors() || Y.hasArrayErrors()) { + drainArrayErrors(X); + drainArrayErrors(Y); + return false; + } + + bool okX = valid2DCoordinates(outX, X, dx, dy, m, n, + sizeMismatchErr, badCoordsErr); + bool okY = valid2DCoordinates(outY, Y, dx, dy, m, n, + sizeMismatchErr, badCoordsErr); + return okX && okY; +} + +// +// valid3DCoordinates validates a user-supplied coordinate array +// against an expected valid grid coordinate one. If the user did not +// provide an array of coordinates, and the grid is uniform, the +// expected value is assigned. +// +bool gridBase::valid3DCoordinates(array3D& userInput, + const array3D& expected, + Real dx, Real dy, Real dz, + size_t m, size_t n, size_t o, + int sizeMismatchErr, + int badCoordsErr) { + if (userInput.data_.is_empty()) { + userInput = expected; // auto-generate + return true; + } + if (userInput.data_.n_rows != expected.data_.n_rows || + userInput.data_.n_cols != expected.data_.n_cols || + userInput.data_.n_slices != expected.data_.n_slices) { + string sparams = "m = " + to_string(m) + ", n = "; + sparams += to_string(n) + ", o = " + to_string(o) + ", dx = "; + sparams += to_string(dx) + ", dy = " + to_string(dy); + sparams += ", dz = " + to_string(dz); + logGridErr(sizeMismatchErr, "grid3D[construct]", sparams); + return false; + } + if (!numEqualArray(expected, userInput, 4.0)) {//eps*4.0 precision + string sparams = "m = " + to_string(m); + sparams += ", n = " + to_string(n) + ", o = "; + sparams += to_string(o) + ", dx = " + to_string(dx); + sparams += ", dy = " + to_string(dy) + ", dz = "; + sparams += to_string(dz); + logGridErr(badCoordsErr, "grid3D[construct]", sparams); + return false; + } + return true; +} + +// +// buildOrCheck3DCoords implementation (see MOLE_grids.h for details). +// +bool gridBase::buildOrCheck3DCoords(array3D& outX, array3D& outY, + array3D& outZ, + const array1D& xcoord, + const array1D& ycoord, + const array1D& zcoord, + Real dx, Real dy, Real dz, + size_t m, size_t n, size_t o, + int sizeMismatchErr, + int badCoordsErr) { + size_t nx = xcoord.data_.n_elem, ny = ycoord.data_.n_elem, + nz = zcoord.data_.n_elem; + array3D X(nx, ny, nz, 0.0), Y(nx, ny, nz, 0.0), + Z(nx, ny, nz, 0.0); + + if (X.data_.n_rows != nx || X.data_.n_cols != ny || + X.data_.n_slices != nz || Y.data_.n_rows != nx || + Y.data_.n_cols != ny || Y.data_.n_slices != nz || + Z.data_.n_rows != nx || Z.data_.n_cols != ny || + Z.data_.n_slices != nz) { + string errmsg = "dim1 = " + to_string(nx) + " X dim2 = "; + errmsg += to_string(ny) + " X dim3 = " + to_string(nz); + logGridErr(MOLE_ERR_FAILED_ARRAY_ALLOC, + "gridBase[buildOrValidateLayer3D]", errmsg); + return false; + } + + nd3DGrid(xcoord, ycoord, zcoord, X, Y, Z); + if (X.hasArrayErrors() || Y.hasArrayErrors() || + Z.hasArrayErrors()) { + drainArrayErrors(X); + drainArrayErrors(Y); + drainArrayErrors(Z); + return false; + } + + bool okX = valid3DCoordinates(outX, X, dx, dy, dz, m, n, o, + sizeMismatchErr, badCoordsErr); + bool okY = valid3DCoordinates(outY, Y, dx, dy, dz, m, n, o, + sizeMismatchErr, badCoordsErr); + bool okZ = valid3DCoordinates(outZ, Z, dx, dy, dz, m, n, o, + sizeMismatchErr, badCoordsErr); + return okX && okY && okZ; +} + +// ------------------------------------------------------------------ +// +// MOLE 1D Grid Class methods (declarations in MOLE_grids.h) +// +// ------------------------------------------------------------------ + +// +// describeGrid1D generates a diagnostic string used when a grid +// construction fails. +// +static string describeGrid1D(const gridParams1D& grid) { + string errmsg = "1D "; + errmsg += string_topology(grid.topology); + errmsg += "ncells = " + to_string(grid.m); + errmsg += ", dx = " + to_string(grid.dx) + ", Periodic = "; + errmsg += grid.bc_isPeriodic ? "YES." : "NO."; + return errmsg; +} + +// +// Checks whether the member 1D grid is valid or not and reports all +// errors found with the grid in its error stack. When the grid +// topology is uniform, this function also generates coordinate +// arrays not provided by the user. +// +bool grid1D::validGrid() { + bool isValid = true; + + if (grid.m <= 0) { + logGridErr(MOLE_ERR_INVALID_GRID_SIZE, "grid1D[construct]", + to_string(grid.m)); + isValid = false; + } + + // u = uniform, c = curvilinear, n = 'non-uniform + switch (grid.topology) { + case 'u': { + if (!validSpacing(grid.dx)) { + logGridErr(MOLE_ERR_INVALID_GRID_SPACING, + "grid1D[construct]", to_string(grid.dx)); + isValid = false; + break; + } + + array1D xn(grid.m + 1), xc(grid.m + 2); + // check for proper allocation of arrays (safety) + if (xn.data_.n_elem == grid.m+1 && + xc.data_.n_elem == grid.m+2){ + generateNodalPts(grid.m, grid.dx, xn); + generateCenterPts(grid.m, grid.dx, xc); + if (!valid1DCoordinates(grid.nodes_X, xn, grid.dx, + grid.m, MOLE_ERR_GRID_NODAL_SZ_MISMATCH, + MOLE_ERR_INVALID_NODAL_COORDINATES)) + isValid = false; + + if (!valid1DCoordinates(grid.centers_X, xc, grid.dx, + grid.m, MOLE_ERR_GRID_CENTERS_SZ_MISMATCH, + MOLE_ERR_INVALID_CENTER_COORDINATES)) + isValid = false; + } + else{ // problems allocating one or both arrays (xn and cn) + string errmsg = "Either xn_dim = "; + errmsg += to_string(grid.m+1) + ", and xc_dim = "; + errmsg += to_string(grid.m+2) + "could not be allocated"; + logGridErr(MOLE_ERR_FAILED_ARRAY_ALLOC, + "grid1D[construct] generating coordinates", errmsg); + isValid = false; + } + break; + } + case 'c': // 1D curvilinear grids are fundamentally undefined + logGridErr(MOLE_ERR_INVALID_1D_CURVILINEAR, + "grid1D[construct]", ""); + isValid = false; + break; + case 'n': // nonuniform grids require user-supplied nodes_X + if (grid.nodes_X.data_.is_empty()) { + logGridErr(MOLE_ERR_INVALID_NONUNIFORM_GRID, + "grid1D[construct]", ""); + isValid = false; + } + break; + default: // invalid topology + string errmsg = string_topology(grid.topology); + logGridErr(MOLE_ERR_INVALID_GRID_TOPOLOGY, + "grid1D[construct]", errmsg); + isValid = false; + break; + } + if (isValid) setGridValidated(); + return isValid; +} + +// +// This grid1D constructor creates and validates user supplied grid. +// Users use a gridParams1D struct to create a valid MOLE grid with +// this constructore. The minimum grid required attributes from a +// user are: m, dx, and topology. Optionally. users can also provide +// the grid coordinate arrays (nodes_X & centers_X), and whether +// the grid has periodic boundary conditions (default = non-periodic) +// +grid1D::grid1D(const gridParams1D p1): gridBase(1) { + grid = p1; // copy the gridParams1D struct into the grid1D member + + if (!validGrid()){ + logGridErr(MOLE_ERR_GRID_CONSTRUCTION_FAILED, + "grid1D[grid1D constructor]", describeGrid1D(grid)); + } +} + +// +// Like the grid1D constructor above, and applies a MOLE debug mode +// to the result. The delegated-to constructor has already run +// validGrid() by the time the body executes. +// +grid1D::grid1D(const gridParams1D p1, size_t debug_mode) + : grid1D(p1) { + applyDebugMode(debug_mode); +} + +// +// Like the grid1D constructor, this constructor also creates and +// validates a user supplied grid. The similar minimum requirements +// apply to a gridParams1D. The inners contains previously logged +// errors that need to be accumulated with the MOLE grid. +// +grid1D::grid1D(const gridParams1D p1, + const stack& inerrs): gridBase(1) { + grid = p1; + mergeErrors(inerrs); + + if (!validGrid()){ + logGridErr(MOLE_ERR_GRID_CONSTRUCTION_FAILED, + "grid1D[grid1D constructor]", describeGrid1D(grid)); + } +} + +// +// Like the grid1D constructor above, and applies a MOLE debug mode +// to the result. +// +grid1D::grid1D(const gridParams1D p1, + const stack& inerrs, + size_t debug_mode) : grid1D(p1, inerrs) { + applyDebugMode(debug_mode); +} + +// ------------------------------------------------------------------ +// +// MOLE 2D Grid Class methods (declarations in MOLE_grids.h) +// +// ------------------------------------------------------------------ + +// +// Functions that generate 2D Grid Coordinates +// + +// +// describeGrid2D generates a diagnostic string used when a grid +// construction fails. +// +static string describeGrid2D(const gridParams2D& grid) { + string errmsg = "2D "; + errmsg += string_topology(grid.topology); + errmsg += ", m cells = " + to_string(grid.m); + errmsg += ", n cells = " + to_string(grid.n); + errmsg += ", dx = " + to_string(grid.dx); + errmsg += ", dy = " + to_string(grid.dy) + ", x-Periodic = "; + errmsg += grid.bc_isPeriodic[0] ? "YES" : "NO"; + errmsg += ", y-Periodic = "; + errmsg += grid.bc_isPeriodic[1] ? "YES." : "NO."; + return errmsg; +} + +// +// n2DGrid creates a 2D rectangular grid from 2 input coordinate +// vectors x, y, and it outputs the corresponding 2D grid in two +// 2D-arrays, X and Y, with the rectangular grid coordinates +// [X, Y] = n2DGrid(x, y) +// +void nd2DGrid(const array1D& x, const array1D& y, + array2D& OutX, array2D& OutY) { + size_t rows = x.data_.n_elem; + size_t cols = y.data_.n_elem; + + //check that the inputs have the correct dimensions + if (x.data_.n_elem == 0 || y.data_.n_elem == 0) { + OutX.data_.set_size(0, 0); + OutY.data_.set_size(0, 0); + string errmsg = "x size = " + to_string(x.data_.n_elem); + errmsg += ", and y size = " + to_string(y.data_.n_elem); + OutX.logArrayError(MOLE_ERR_INVALID_ARRAY_SIZE, "nd2DGrid", + errmsg); + OutY.logArrayError(MOLE_ERR_INVALID_ARRAY_SIZE, "nd2DGrid", + errmsg); + return; + } + if (OutX.data_.n_rows != rows || OutX.data_.n_cols != cols) { + OutX.data_.set_size(0, 0); + OutY.data_.set_size(0, 0); + string errmsg = "X rows = " + to_string(OutX.data_.n_rows); + errmsg += " and cols = " + to_string(OutX.data_.n_cols); + errmsg += ", but needs rows = " + to_string(rows); + errmsg += ", cols = " + to_string(cols); + OutX.logArrayError(MOLE_ERR_INVALID_ARRAY_SIZE, "nd2DGrid", + errmsg); + return; + } + if (OutY.data_.n_rows != rows || OutY.data_.n_cols != cols) { + OutX.data_.set_size(0, 0); + OutY.data_.set_size(0, 0); + string errmsg = "Y rows = " + to_string(OutY.data_.n_rows); + errmsg += " and cols = " + to_string(OutY.data_.n_cols); + errmsg += ", but needs rows = " + to_string(rows); + errmsg += ", cols = " + to_string(cols); + OutY.logArrayError(MOLE_ERR_INVALID_ARRAY_SIZE, "nd2DGrid", + errmsg); + return; + } + + // Generate the 2D grid coodinates following the Octave + // ndgrid() function + for (size_t r = 0; r < rows; ++r) { + for (size_t c = 0; c < cols; ++c) { + OutX.data_(r,c) = x.data_(r); + OutY.data_(r,c) = y.data_(c); + } + } + return; +} + +// +// Checks whether the member 2D grid is valid or not and reports all +// errors found with the grid in its error stack. When the grid +// topology is uniform, this function also generates coordinate +// arrays not provided by the user. +// +bool grid2D::validGrid() { + bool isValid = true; + + if (grid.m <= 0 || grid.n <= 0) { + logGridErr(MOLE_ERR_INVALID_GRID_SIZE, "grid2D[construct]", + to_string(grid.m)); + isValid = false; + } + + // u = uniform, c = curvilinear, n = 'non-uniform + switch (grid.topology) { + case 'u': { + if (!validSpacing(grid.dx) || !validSpacing(grid.dy)) { + string errmsg = "dx = "; + errmsg += to_string(grid.dx) + ", dy = "; + errmsg += to_string(grid.dy); + logGridErr(MOLE_ERR_INVALID_GRID_SPACING, + "grid2D[construct]", errmsg); + isValid = false; + break; + } + // Generate the four 1D coordinate arrays used for computing + // nodal coordinates: xn = (0:m)*dx, yn = (0:n)*dy, center: + // xc = [0, (0.5:m-0.5)*dx, m*dx], yc likewise, and normal + // faces xv = xc[1:m], yu = yc[1:n] (interior points only). + array1D xn(grid.m + 1), yn(grid.n + 1); + array1D xc(grid.m + 2), yc(grid.n + 2); + array1D xv(grid.m), yu(grid.n); + + if (xn.data_.n_elem == grid.m+1 && + yn.data_.n_elem == grid.n+1 && + xc.data_.n_elem == grid.m+2 && + yc.data_.n_elem == grid.n+2 && + xv.data_.n_elem == grid.m && + yu.data_.n_elem == grid.n) { + generateNodalPts(grid.m, grid.dx, xn); + generateNodalPts(grid.n, grid.dy, yn); + generateCenterPts(grid.m, grid.dx, xc); + generateCenterPts(grid.n, grid.dy, yc); + xv.data_ = xc.data_.subvec(1, grid.m); + yu.data_ = yc.data_.subvec(1, grid.n); + } else { + string errmsg = "one or more 1D coordinate arrays"; + errmsg += " for m = " + to_string(grid.m) + ", n = "; + errmsg += to_string(grid.n) + " could not be allocated"; + logGridErr(MOLE_ERR_FAILED_ARRAY_ALLOC, + "grid2D[construct] 1D coordinates", errmsg); + isValid = false; + break; + } + + // Build/validate nodal coordinates + if (!buildOrCheck2DCoords(grid.nodes_X, grid.nodes_Y, xn, yn, + grid.dx, grid.dy, grid.m+1, grid.n+1, + MOLE_ERR_GRID_NODAL_SZ_MISMATCH, + MOLE_ERR_INVALID_NODAL_COORDINATES)) + isValid = false; + + // Build/validate cell center coordinates + if (!buildOrCheck2DCoords(grid.centers_X, grid.centers_Y, xc, + yc, grid.dx, grid.dy, grid.m+2, grid.n+2, + MOLE_ERR_GRID_CENTERS_SZ_MISMATCH, + MOLE_ERR_INVALID_CENTER_COORDINATES)) + isValid = false; + + // Build/validate normal face coordinates + if (!buildOrCheck2DCoords(grid.faces_u_X, grid.faces_u_Y, xn, + yu, grid.dx, grid.dy, grid.m+1, grid.n, + MOLE_ERR_GRID_FACES_SZ_MISMATCH, + MOLE_ERR_INVALID_NORMAL_FACE_COORDS)) + isValid = false; + + if (!buildOrCheck2DCoords(grid.faces_v_X, grid.faces_v_Y, xv, + yn, grid.dx, grid.dy, grid.m, grid.n+1, + MOLE_ERR_GRID_FACES_SZ_MISMATCH, + MOLE_ERR_INVALID_NORMAL_FACE_COORDS)) + isValid = false; + + break; + } + case 'c': // User must provide at least nodal coordinates + if (grid.nodes_X.data_.is_empty() || + grid.nodes_Y.data_.is_empty()){ + logGridErr(MOLE_ERR_INVALID_CURVILINEAR_GRID, + "grid2D[construct]", ""); + isValid = false; + } + break; + case 'n': // nonuniform grids require user-supplied nodes_X + if (grid.nodes_X.data_.is_empty() || + grid.nodes_Y.data_.is_empty()) { + logGridErr(MOLE_ERR_INVALID_NONUNIFORM_GRID, + "grid2D[construct]", ""); + isValid = false; + } + break; + default: // invalid grid topology + string errmsg = string_topology(grid.topology); + logGridErr(MOLE_ERR_INVALID_GRID_TOPOLOGY, + "grid2D[construct]", errmsg); + isValid = false; + break; + } + + if (isValid) setGridValidated(); + return isValid; +} + +// +// This grid2D constructor creates and validates a 2D grid. +// User needs to input at least m, n, dx, dy and topology in +// a gridParam2D structure, and optionally they can also provide the +// grid coordinate arrays (nodes_X, Nodes_Y, centers_X, centers_Y, +// and faces). This function also initializes the GridBase error_log +// For curvilinear and nonuniform grids, users need to providal +// nodal grid information +// +grid2D::grid2D(gridParams2D p2): gridBase(2) { + grid = p2; // memberwise assignement of the gridParams2D struct + + if (!validGrid()){ + logGridErr(MOLE_ERR_GRID_CONSTRUCTION_FAILED, + "grid2D[grid2D constructor]", describeGrid2D(grid)); + } +} + +// +// Like the grid2D constructor above, and applies a MOLE debug mode +// to the result. +// +grid2D::grid2D(gridParams2D p2, size_t debug_mode) : grid2D(p2) { + applyDebugMode(debug_mode); +} + +// +// Like the grid2D constructor, this constructor also creates and +// validates a user supplied grid. The similar minimum requirements +// apply to a gridParams2D. The inners contains previously logged +// errors that need to be accumulated with the MOLE grid. +// +grid2D::grid2D(gridParams2D p2, + const stack& inerrs): gridBase(2) { + grid = p2; //memberwise assignement of the gridParams2D struct + + mergeErrors(inerrs); + + if (!validGrid()){ + logGridErr(MOLE_ERR_GRID_CONSTRUCTION_FAILED, + "grid2D[grid2D constructor]", describeGrid2D(grid)); + } +} + +// +// Like the grid2D constructor above, and applies a MOLE debug mode +// to the result. +// +grid2D::grid2D(gridParams2D p2, + const stack& inerrs, + size_t debug_mode) : grid2D(p2, inerrs) { + applyDebugMode(debug_mode); +} + +// ------------------------------------------------------------------ +// +// MOLE 3D Grid Class methods (declarations in MOLE_grids.h) +// +// ------------------------------------------------------------------ + +// +// Functions that generate 3D Grid Coordinates +// + +// +// describeGrid3D generates a diagnostic string used when a grid +// construction fails. +// +static string describeGrid3D(const gridParams3D& grid) { + string errmsg = "3D "; + errmsg += string_topology(grid.topology); + errmsg += ", m cells = " + to_string(grid.m); + errmsg += ", n cells = " + to_string(grid.n); + errmsg += ", o cells = " + to_string(grid.o); + errmsg += ", dx = " + to_string(grid.dx); + errmsg += ", dy = " + to_string(grid.dy); + errmsg += ", dz = " + to_string(grid.dz) + ", x-Periodic = "; + errmsg += grid.bc_isPeriodic[0] ? "YES" : "NO"; + errmsg += ", y-Periodic = "; + errmsg += grid.bc_isPeriodic[1] ? "YES" : "NO"; + errmsg += ", z-Periodic = "; + errmsg += grid.bc_isPeriodic[2] ? "YES." : "NO."; + return errmsg; +} + +// +// n3DGrid creates a 3D grid from 3 input coordinate vectors x, y, +// and z. The function outputs the corresponding 3D grid in three +// 3D-arrays, X, Y and Z, with the rectangular grid coordinates +// [X, Y, Z] = n3DGrid(x, y, z) +// +void nd3DGrid(const array1D& x, const array1D& y, const array1D& z, + array3D& OutX, array3D& OutY, array3D& OutZ) { + size_t rows = x.data_.n_elem; + size_t cols = y.data_.n_elem; + size_t slices = z.data_.n_elem; + + //check that the inputs have the correct dimensions + if (x.data_.n_elem == 0 || y.data_.n_elem == 0 || + z.data_.n_elem == 0) { + OutX.data_.set_size(0, 0, 0); + OutY.data_.set_size(0, 0, 0); + OutZ.data_.set_size(0, 0, 0); + string errmsg = "x size = " + to_string(x.data_.n_elem); + errmsg += ", and y size = " + to_string(y.data_.n_elem); + errmsg += ", and z size = " + to_string(z.data_.n_elem); + OutX.logArrayError(MOLE_ERR_INVALID_ARRAY_SIZE, "nd3DGrid", + errmsg); + OutY.logArrayError(MOLE_ERR_INVALID_ARRAY_SIZE, "nd3DGrid", + errmsg); + OutZ.logArrayError(MOLE_ERR_INVALID_ARRAY_SIZE, "nd3DGrid", + errmsg); + return; + } + if (OutX.data_.n_rows != rows || OutX.data_.n_cols != cols || + OutX.data_.n_slices != slices) { + OutX.data_.set_size(0, 0, 0); + OutY.data_.set_size(0, 0, 0); + OutZ.data_.set_size(0, 0, 0); + string errmsg = "X rows = " + to_string(OutX.data_.n_rows); + errmsg += " and cols = " + to_string(OutX.data_.n_cols); + errmsg += " and slices = " + to_string(OutX.data_.n_slices); + errmsg += ", but needs rows = " + to_string(rows); + errmsg += ", cols = " + to_string(cols) + ", slices = "; + errmsg += to_string(slices); + OutX.logArrayError(MOLE_ERR_INVALID_ARRAY_SIZE, "nd3DGrid", + errmsg); + return; + } + if (OutY.data_.n_rows != rows || OutY.data_.n_cols != cols || + OutY.data_.n_slices != slices) { + OutX.data_.set_size(0, 0, 0); + OutY.data_.set_size(0, 0, 0); + OutZ.data_.set_size(0, 0, 0); + string errmsg = "Y rows = " + to_string(OutY.data_.n_rows); + errmsg += " and cols = " + to_string(OutY.data_.n_cols); + errmsg += " and slices = " + to_string(OutY.data_.n_slices); + errmsg += ", but needs rows = " + to_string(rows); + errmsg += ", cols = " + to_string(cols) + ", slices = "; + errmsg += to_string(slices); + OutY.logArrayError(MOLE_ERR_INVALID_ARRAY_SIZE, "nd3DGrid", + errmsg); + return; + } + if (OutZ.data_.n_rows != rows || OutZ.data_.n_cols != cols || + OutZ.data_.n_slices != slices) { + OutX.data_.set_size(0, 0, 0); + OutY.data_.set_size(0, 0, 0); + OutZ.data_.set_size(0, 0, 0); + string errmsg = "Z rows = " + to_string(OutZ.data_.n_rows); + errmsg += " and cols = " + to_string(OutZ.data_.n_cols); + errmsg += " and slices = " + to_string(OutZ.data_.n_slices); + errmsg += ", but needs rows = " + to_string(rows); + errmsg += ", cols = " + to_string(cols) + ", slices = "; + errmsg += to_string(slices); + OutZ.logArrayError(MOLE_ERR_INVALID_ARRAY_SIZE, "nd3DGrid", + errmsg); + return; + } + // Allocate the 3D vectors + + // Populate the 3D workspace + for (size_t r = 0; r < rows; ++r) { + for (size_t c = 0; c < cols; ++c) { + for (size_t s = 0; s < slices; ++s) { + OutX.data_(r, c, s) = x.data_(r); + OutY.data_(r, c, s) = y.data_(c); + OutZ.data_(r, c, s) = z.data_(s); + } + } + } +} + + + +// +// Checks whether the member 3D grid is valid or not, and reports all +// errors found with the grid in its error stack. When the grid +// topology is uniform, this function also generates coordinate +// arrays not provided by the user. +// +bool grid3D::validGrid() { + bool isValid = true; + + if (grid.m <= 0 || grid.n <= 0 || grid.o <=0) { + logGridErr(MOLE_ERR_INVALID_GRID_SIZE, "grid3D[construct]", + to_string(grid.m)); + isValid = false; + } + + // u = uniform, c = curvilinear, n = 'non-uniform + switch (grid.topology) { + case 'u': { + if (!validSpacing(grid.dx) || !validSpacing(grid.dy) || + !validSpacing(grid.dz)) { + string errmsg = "dx = " + to_string(grid.dx); + errmsg += ", dy = " + to_string(grid.dy); + errmsg += ", dz = " + to_string(grid.dz); + logGridErr(MOLE_ERR_INVALID_GRID_SPACING, + "grid3D[construct]", errmsg); + isValid = false; + break; + } + + // Generate the 1D coordinate arrays used for computing + // nodal coordinates (xn, yn, zn), cell centers (xc, yc, zc) + // and normal faces (xv, yu, zu). The latter 3 arrays are + // the interior points of the center arrays. + array1D xn(grid.m + 1), yn(grid.n + 1), zn(grid.o + 1); + array1D xc(grid.m + 2), yc(grid.n + 2), zc(grid.o + 2); + array1D xv(grid.m), yu(grid.n), zu(grid.o); + + if (xn.data_.n_elem == grid.m+1 && + yn.data_.n_elem == grid.n+1 && + zn.data_.n_elem == grid.o+1 && + xc.data_.n_elem == grid.m+2 && + yc.data_.n_elem == grid.n+2 && + zc.data_.n_elem == grid.o+2 && + xv.data_.n_elem == grid.m && + yu.data_.n_elem == grid.n && + zu.data_.n_elem == grid.o) { + generateNodalPts(grid.m, grid.dx, xn); + generateNodalPts(grid.n, grid.dy, yn); + generateNodalPts(grid.o, grid.dz, zn); + generateCenterPts(grid.m, grid.dx, xc); + generateCenterPts(grid.n, grid.dy, yc); + generateCenterPts(grid.o, grid.dz, zc); + xv.data_ = xc.data_.subvec(1, grid.m); + yu.data_ = yc.data_.subvec(1, grid.n); + zu.data_ = zc.data_.subvec(1, grid.o); + } else { + string errmsg = "one or more 1D coordinate arrays for"; + errmsg += " m = " + to_string(grid.m) + ", n = "; + errmsg += to_string(grid.n) + ", o = "; + errmsg += to_string(grid.o) + " could not be allocated"; + logGridErr(MOLE_ERR_FAILED_ARRAY_ALLOC, + "grid3D[construct] 1D coordinates", errmsg); + isValid = false; + break; + } + + // Build/validate nodal coordinates + if (!buildOrCheck3DCoords(grid.nodes_X, grid.nodes_Y, grid.nodes_Z, + xn, yn, zn, grid.dx, grid.dy, grid.dz, + grid.m+1, grid.n+1, grid.o+1, + MOLE_ERR_GRID_NODAL_SZ_MISMATCH, + MOLE_ERR_INVALID_NODAL_COORDINATES)) + isValid = false; + + // Build/validate cell center coordinates + if (!buildOrCheck3DCoords(grid.centers_X, grid.centers_Y, grid.centers_Z, + xc, yc, zc, grid.dx, grid.dy, grid.dz, + grid.m+2, grid.n+2, grid.o+2, + MOLE_ERR_GRID_CENTERS_SZ_MISMATCH, + MOLE_ERR_INVALID_CENTER_COORDINATES)) + isValid = false; + + // Build/validate normal face coordinates + if (!buildOrCheck3DCoords(grid.faces_u_X, grid.faces_u_Y, grid.faces_u_Z, + xn, yu, zu, grid.dx, grid.dy, grid.dz, + grid.m+1, grid.n, grid.o, + MOLE_ERR_GRID_FACES_SZ_MISMATCH, + MOLE_ERR_INVALID_NORMAL_FACE_COORDS)) + isValid = false; + + if (!buildOrCheck3DCoords(grid.faces_v_X, grid.faces_v_Y, grid.faces_v_Z, + xv, yn, zu, grid.dx, grid.dy, grid.dz, + grid.m, grid.n+1, grid.o, + MOLE_ERR_GRID_FACES_SZ_MISMATCH, + MOLE_ERR_INVALID_NORMAL_FACE_COORDS)) + isValid = false; + + if (!buildOrCheck3DCoords(grid.faces_w_X, grid.faces_w_Y, grid.faces_w_Z, + xv, yu, zn, grid.dx, grid.dy, grid.dz, + grid.m, grid.n, grid.o+1, + MOLE_ERR_GRID_FACES_SZ_MISMATCH, + MOLE_ERR_INVALID_NORMAL_FACE_COORDS)) + isValid = false; + + break; + } + case 'c': // User must provide at least nodal coordinates + if (grid.nodes_X.data_.is_empty() || grid.nodes_Y.data_.is_empty() || + grid.nodes_Z.data_.is_empty()){ + logGridErr(MOLE_ERR_INVALID_CURVILINEAR_GRID, + "grid3D[construct]", ""); + isValid = false; + } + break; + + + case 'n': // nonuniform grids require user-supplied nodes_X + if (grid.nodes_X.data_.is_empty() || grid.nodes_Y.data_.is_empty() || + grid.nodes_Z.data_.is_empty()) { + logGridErr(MOLE_ERR_INVALID_NONUNIFORM_GRID, + "grid3D[construct]", ""); + isValid = false; + } + break; + default: // invalid grid topology + string errmsg = string_topology(grid.topology); + logGridErr(MOLE_ERR_INVALID_GRID_TOPOLOGY, + "grid3D[construct]", errmsg); + isValid = false; + break; + } + + if (isValid) setGridValidated(); + return isValid; +} + +// +// This grid3D constructor creates and validates a 3D grid. +// User needs to input at least m, n, o, dx, dy, dz and topology in +// a gridParam3D structure, and optionally they can also provide the +// grid coordinate arrays (nodes_X, Nodes_Y, Nodes_Z, centers_X, +// centers_Y, centers_Z, and faces). This function also initializes +// the GridBase error_log. For curvilinear and nonuniform grids, +// users need to providal nodal grid information +// +grid3D::grid3D(gridParams3D p3): gridBase(3) { + grid = p3; // memberwise assignement of the gridParams2D struct + + if (!validGrid()){ + logGridErr(MOLE_ERR_GRID_CONSTRUCTION_FAILED, + "grid3D[grid3D constructor]", describeGrid3D(grid)); + } +} + +// +// Like the grid3D constructor above, and applies a MOLE debug mode +// to the result. +// +grid3D::grid3D(gridParams3D p3, size_t debug_mode) : grid3D(p3) { + applyDebugMode(debug_mode); +} + +// +// Like the grid3D constructor, this constructor also creates and +// validates a user supplied grid. The similar minimum requirements +// apply to a gridParams3D. The inners contains previously logged +// errors that need to be accumulated with the MOLE grid. +// +grid3D::grid3D(gridParams3D p3, + const stack& inerrs): gridBase(3) { + grid = p3; // memberwise assignement of the gridParams2D struct + mergeErrors(inerrs); + + if (!validGrid()){ + logGridErr(MOLE_ERR_GRID_CONSTRUCTION_FAILED, + "grid3D[grid3D constructor]", describeGrid3D(grid)); + } +} + +// +// Like the grid3D constructor above, and applies a MOLE debug mode +// to the result. +// +grid3D::grid3D(gridParams3D p3, + const stack& inerrs, + size_t debug_mode) : grid3D(p3, inerrs) { + applyDebugMode(debug_mode); +} + +// ------------------------------------------------------------------ +// +// MOLE gridNull Class methods (declarations in MOLE_grids.h) +// +// ------------------------------------------------------------------ + +// +// gridNull sole constructor requires a paramsNull struct and errors. +// It is private; makeGridNull is the only caller. +// +gridNull::gridNull(const paramsNull in_p, + const stack& inerrs): gridBase(0){ + mergeErrors(inerrs); + ErrData.num_errs = in_p.num_errs; + if(ErrData.num_errs == 0){ + ErrData.num_errs = inerrs.size(); + } else { + ErrData.num_errs += inerrs.size(); + } + ErrData.type_errs.push("MOLE Grid"); +} + +// +// makeGridNull is the only entry point to the gridNull constructor. +// It exists so that makeGrid can produce the failure case without +// opening the constructor to users. +// +gridNull makeGridNull(const paramsNull in_p, + const stack& inerrs){ + return gridNull(in_p, inerrs); +} + +// ---------------------------------------------------------------- +// +// gridVar makeGrid is a factory function that works for any of the 3 +// grid dimensionalities, intended for cases when the users need to +// define the dimensionality at runtime. The function takes a variant +// grid structure paramVars and produces a variant grid class gridVar +// +// ---------------------------------------------------------------- +gridVar makeGrid(paramVars params, const stack& errs){ + // you need to include the errors that are already in the + // stack to the + // the corresponding grid structure + return std::visit([&errs](auto&& p) -> gridVar { + using T = std::decay_t; + if constexpr (std::is_same_v) { + return grid1D(p, errs); + } else if constexpr (std::is_same_v) { + return grid2D(p, errs); + } else if constexpr (std::is_same_v){ + return grid3D(p, errs); + } else if constexpr (std::is_same_v){ + return makeGridNull(p, errs); + } + + }, params); +} + +// --------------------------------------------------------- +// Dispatch function which takes a gridVar (generic grid) and +// validates it using the instantiated MOLE grid class (i.e., +// grid1D, grid2D or grid3D) for the actual grid validation +// --------------------------------------------------------- +bool isValidGrid(gridVar& g) { + return std::visit([](auto&& gridObj) + { return gridObj.validGrid(); }, g); +} diff --git a/cpp/src/grids/README.md b/cpp/src/grids/README.md new file mode 100644 index 00000000..7867387b --- /dev/null +++ b/cpp/src/grids/README.md @@ -0,0 +1,34 @@ + + +# Subdirectory for MOLE 2.0 C++ grid classes and their functionality + +Subdirectory and Pathname: **mole/cpp/src/grids/** + +## Purpose + +Subdirectory containing source implementations of MOLE grid classes +and member functions, also flat 2D and 3D arrays to avoid the +construction of multidimensional arrays using the C++ vector class. + +## List of Files in This Subdirectory (in alphabetical order) + ++ **MOLE_arrays.cpp**: C++ implementation of MOLE 1D, 2D, and 3D +array classes, which are wrappers to data classes in other numerical +libraries. It also contains functionality that support the classes, +and custom wrappers to the MOLE error handling mechanisms. ++ **MOLE_grids.cpp**: C++ implementation of MOLE 1D, 2D, and 3D grid +classes. It also contains data structures that support the classes, +and custom wrappers to the MOLE error handling mechanisms. ++ **README.md**: (this file) diff --git a/cpp/src/grids/grid_builder.cpp b/cpp/src/grids/grid_builder.cpp new file mode 100644 index 00000000..f21496cb --- /dev/null +++ b/cpp/src/grids/grid_builder.cpp @@ -0,0 +1,340 @@ +/* + * SPDX-License-Identifier: GPL-3.0-or-later + * Copyright (c) 2008-2024 San Diego State University Research Foundation + * (SDSURF). + * See LICENSE file or https://www.gnu.org/licenses/gpl-3.0.html for details. + */ + +/* + * @file grid_builder.cpp + * + * @brief Parses a variable length user input containing grid + * attributes in a pair , the + * functions in this file, parses the user input and creates + * the corresponding grid classes (actual grid validation is + * part of the MOLE grid constructors) + * + * @date 2026/07/27 + * + */ +#include "grid_builder.h" + +#include +#include +#include +#include +#include + +namespace { + +const std::string kLoc = "gridBuilder"; // error-stack location tag + +// parse_va_args parses an input va_args containing pairs +// of attributes until a nullptr sentinel is found. The key needs to +// match an MOLE grid's attribute name, else an error; +// MAKE_GRID_UNKNOWN_ATTRIBUTE, is generated and the parsing process +// is stopped. NOTE: The actual context of values is check in +// validGrid() +gridRaw parse_va_args(std::stack& errs, + const char* firstName, va_list& ap) { + gridRaw g; + for (const char* name = firstName; name != nullptr; + name = va_arg(ap, const char*)) { + if (std::strcmp(name, "m") == 0) + g.m = va_arg(ap, int); + else if (std::strcmp(name, "n") == 0) + g.n = va_arg(ap, int); + else if (std::strcmp(name, "o") == 0) + g.o = va_arg(ap, int); + else if (std::strcmp(name, "dim") == 0) + g.dim = va_arg(ap, int); + else if (std::strcmp(name, "dx") == 0) + g.dx = va_arg(ap, double); + else if (std::strcmp(name, "dy") == 0) + g.dy = va_arg(ap, double); + else if (std::strcmp(name, "dz") == 0) + g.dz = va_arg(ap, double); + + // The MOLE debug modes are integer macros. An unrecognized + // value is not rejected here; applyDebugMode reports it. + else if (std::strcmp(name, "debug") == 0) + g.debug = va_arg(ap, int); + + // char promotes to int through varargs; read int, store char. + else if (std::strcmp(name, "topology") == 0) + g.topology = static_cast(va_arg(ap, int)); + + else if (std::strcmp(name, "nodes.X") == 0) + g.nodesX = va_arg(ap, const void*); + else if (std::strcmp(name, "nodes.Y") == 0) + g.nodesY = va_arg(ap, const void*); + else if (std::strcmp(name, "nodes.Z") == 0) + g.nodesZ = va_arg(ap, const void*); + else if (std::strcmp(name, "centers.X") == 0) + g.centersX = va_arg(ap, const void*); + else if (std::strcmp(name, "centers.Y") == 0) + g.centersY = va_arg(ap, const void*); + else if (std::strcmp(name, "centers.Z") == 0) + g.centersZ = va_arg(ap, const void*); + else if (std::strcmp(name, "faces.u.X") == 0) + g.facesuX = va_arg(ap, const void*); + else if (std::strcmp(name, "faces.u.Y") == 0) + g.facesuY = va_arg(ap, const void*); + else if (std::strcmp(name, "faces.u.Z") == 0) + g.facesuZ = va_arg(ap, const void*); + else if (std::strcmp(name, "faces.v.X") == 0) + g.facesvX = va_arg(ap, const void*); + else if (std::strcmp(name, "faces.v.Y") == 0) + g.facesvY = va_arg(ap, const void*); + else if (std::strcmp(name, "faces.v.Z") == 0) + g.facesvZ = va_arg(ap, const void*); + else if (std::strcmp(name, "faces.w.X") == 0) + g.faceswX = va_arg(ap, const void*); + else if (std::strcmp(name, "faces.w.Y") == 0) + g.faceswY = va_arg(ap, const void*); + else if (std::strcmp(name, "faces.w.Z") == 0) + g.faceswZ = va_arg(ap, const void*); + + // isPeriodic arrives as the address of the caller's vector, + // which carries its own size. dim may not be parsed yet + // given the variable order of the user input pairs, so the + // pointer is stashed and both the size check and the copy + // happen once dim is known. + else if (std::strcmp(name, "isPeriodic") == 0) + g.isPeriodicSrc = va_arg(ap, const std::vector*); + + else { + MOLEerr_log(errs, MAKE_GRID_UNKNOWN_ATTRIBUTE, kLoc, name); + break; + } + } + return g; +} + +} // namespace + + +// runChecks: checks require grid attributes and their consistency +// for MOLE grid specifications. In case of issues with the grid +// specifications, it records all possible errors with the grid +// attributes before returning. +int runChecks(std::stack& errs, const gridRaw& g) { + // The grid dimensionality is a required parameter and is not + // derived from m, n, o. Without it the remaining checks have + // nothing to compare against, so return early. + if (g.dim == -1) { + MOLEerr_log(errs, MOLE_ERR_INVALID_GRID_DIM, kLoc); + return 0; + } + if (g.dim < 1 || g.dim > 3) { + MOLEerr_log(errs, MOLE_ERR_INVALID_GRID_DIM, kLoc, + std::to_string(g.dim)); + return 0; + } + + // The cell counts supplied must agree with the dimensionality, + // in both directions: a count the dimension needs but did not + // get, and a count the dimension has no use for. + if (g.dim >= 1 && g.m < 0) { + std::string errparam = "m = " + std::to_string(g.m); + MOLEerr_log(errs, MOLE_ERR_INVALID_CELL_COUNT, kLoc, + errparam); + } + if (g.dim >= 2 && g.n < 0) { + std::string errparam = "n = " + std::to_string(g.n); + MOLEerr_log(errs, MOLE_ERR_INVALID_CELL_COUNT, kLoc, + errparam); + } + if (g.dim >= 3 && g.o < 0) { + std::string errparam = "o = " + std::to_string(g.o); + MOLEerr_log(errs, MOLE_ERR_INVALID_CELL_COUNT, kLoc, + errparam); + } + if (g.dim == 1 && g.n > 0) { + std::string errparam = + "1D grid with n = " + std::to_string(g.n); + MOLEerr_log(errs, MOLE_ERR_INVALID_CELL_COUNT, kLoc, + errparam); + } + if (g.dim < 3 && g.o > 0) { + std::string errparam = + "Not a 3D grid but o = " + std::to_string(g.o); + MOLEerr_log(errs, MOLE_ERR_INVALID_CELL_COUNT, kLoc, + errparam); + } + + // isPeriodic holds one flag per dimension. The vector reports + // its own size, so the count is checked here rather than + // trusted at the point the flags are read. + if (g.isPeriodicSrc && + g.isPeriodicSrc->size() != static_cast(g.dim)) { + std::string errparam = + "isPeriodic size = " + + std::to_string(g.isPeriodicSrc->size()) + + ", dim = " + std::to_string(g.dim); + MOLEerr_log(errs, MOLE_ERR_INVALID_ISPERIODIC_DIM, kLoc, + errparam); + } + + // Topology is compulsory and has no default value. + if (g.topology == '\0') + MOLEerr_log(errs, MOLE_ERR_INVALID_GRID_TOPOLOGY, kLoc); + else if (g.topology != 'u' && g.topology != 'c' && + g.topology != 'n') + MOLEerr_log(errs, MOLE_ERR_INVALID_GRID_TOPOLOGY, kLoc, + std::string(1, g.topology)); + + return g.dim; +} + +// --- instantiating MOLE grid structures ---------------------------- +// +// Builds the dimension-specific struct from an parsed gridRaw, +// copying across only the fields that dimension has. Anything a user +// did not supply keeps the struct's own default, so the MOLE library +// sees it as unset during validation +// + +namespace { + +// Casts a type-erased coordinate pointer to the concrete array +// class for this dimension and copies the whole object, which +// carries the array's own error stack along with its data. A null +// pointer means the attribute was not supplied, so dst keeps its +// default and MOLE reads it as unset. +template +void copyArrayAs(const void* src, T& dst) { + if (src) dst = *static_cast(src); +} + +// isPeriodicSrc points at a caller-owned vector. It is null when +// the attribute was omitted, so it is checked before any element +// is read. runChecks has already confirmed the size matches dim. +gridParams1D narrow1D(const gridRaw& g) { + gridParams1D p; + p.topology = g.topology; + p.m = static_cast(g.m); + p.dx = g.dx; + copyArrayAs(g.nodesX, p.nodes_X); + copyArrayAs(g.centersX, p.centers_X); + if (g.isPeriodicSrc) + p.bc_isPeriodic = (*g.isPeriodicSrc)[0]; + return p; +} + +gridParams2D narrow2D(const gridRaw& g) { + gridParams2D p; + p.topology = g.topology; + p.m = static_cast(g.m); + p.n = static_cast(g.n); + p.dx = g.dx; + p.dy = g.dy; + copyArrayAs(g.nodesX, p.nodes_X); + copyArrayAs(g.nodesY, p.nodes_Y); + copyArrayAs(g.centersX, p.centers_X); + copyArrayAs(g.centersY, p.centers_Y); + copyArrayAs(g.facesuX, p.faces_u_X); + copyArrayAs(g.facesuY, p.faces_u_Y); + copyArrayAs(g.facesvX, p.faces_v_X); + copyArrayAs(g.facesvY, p.faces_v_Y); + if (g.isPeriodicSrc) { + p.bc_isPeriodic[0] = (*g.isPeriodicSrc)[0]; + p.bc_isPeriodic[1] = (*g.isPeriodicSrc)[1]; + } + return p; +} + +gridParams3D narrow3D(const gridRaw& g) { + gridParams3D p; + p.topology = g.topology; + p.m = static_cast(g.m); + p.n = static_cast(g.n); + p.o = static_cast(g.o); + p.dx = g.dx; + p.dy = g.dy; + p.dz = g.dz; + copyArrayAs(g.nodesX, p.nodes_X); + copyArrayAs(g.nodesY, p.nodes_Y); + copyArrayAs(g.nodesZ, p.nodes_Z); + copyArrayAs(g.centersX, p.centers_X); + copyArrayAs(g.centersY, p.centers_Y); + copyArrayAs(g.centersZ, p.centers_Z); + copyArrayAs(g.facesuX, p.faces_u_X); + copyArrayAs(g.facesuY, p.faces_u_Y); + copyArrayAs(g.facesuZ, p.faces_u_Z); + copyArrayAs(g.facesvX, p.faces_v_X); + copyArrayAs(g.facesvY, p.faces_v_Y); + copyArrayAs(g.facesvZ, p.faces_v_Z); + copyArrayAs(g.faceswX, p.faces_w_X); + copyArrayAs(g.faceswY, p.faces_w_Y); + copyArrayAs(g.faceswZ, p.faces_w_Z); + if (g.isPeriodicSrc) { + p.bc_isPeriodic[0] = (*g.isPeriodicSrc)[0]; + p.bc_isPeriodic[1] = (*g.isPeriodicSrc)[1]; + p.bc_isPeriodic[2] = (*g.isPeriodicSrc)[2]; + } + return p; +} + +// buildGrid runs the dispatch that turns a parsed gridRaw into a +// gridVar. It is split out of gridBuilder_impl so that every path, +// including the failure paths, funnels through a single point where +// the debug mode is applied. +gridVar buildGrid(std::stack& errs, const gridRaw& g, + int dim) { + // A parse-time failure means the grid cannot be built. the error + // MOLE_ERR_INVALID_GRID_ARGS sits on top of the error log; the + // errors beneath it will provide the details. A paramsNull routes + // makeGrid to gridNull, carrying the whole error stack with it. + if (MOLEerr_haserrors(errs)) { + MOLEerr_log(errs, MOLE_ERR_INVALID_GRID_ARGS, kLoc, ""); + return makeGrid(paramsNull{}, errs); + } + + // makeGrid takes the params variant and the error stack (by const + // reference) and runs the grid1D/2D/3D constructor internally, so + // nothing is heap-allocated here. + switch (dim) { + case 1: + return makeGrid(narrow1D(g), errs); + case 2: + return makeGrid(narrow2D(g), errs); + case 3: + return makeGrid(narrow3D(g), errs); + default: + // invalid grid dimension found + MOLEerr_log(errs, MOLE_ERR_INVALID_GRID_DIM, kLoc, + std::to_string(dim)); + return makeGrid(paramsNull{}, errs); + } +} + +} // namespace + +// gridBuilder_impl is called through the gridBuilder macro, which +// appends the nullptr sentinel. The MOLE error structure is created +// and initialized here for the duration of the build. The narrowed +// this function returns a MOLE grid to the user of type gridVar +// which is a C++ variant for Grid1D, Grid2D or Grid3D +gridVar gridBuilder_impl(const char* firstName, ...) { + std::stack errs; + MOLEerr_init(errs); + + va_list ap; + va_start(ap, firstName); + gridRaw g = parse_va_args(errs, firstName, ap); + va_end(ap); + + const int dim = runChecks(errs, g); + + gridVar out = buildGrid(errs, g, dim); + + // The debug mode governs what gridBuilder reports, not what the + // grid holds: the error log is left intact in every mode, so a + // caller can still print it or write it to a file afterwards. + // A grid that validated ignores the mode. + const size_t dbg = static_cast(g.debug); + std::visit([dbg](auto&& grid) { grid.applyDebugMode(dbg); }, out); + + return out; +} diff --git a/cpp/src/include/MOLE_arrays.h b/cpp/src/include/MOLE_arrays.h new file mode 100644 index 00000000..909e1abf --- /dev/null +++ b/cpp/src/include/MOLE_arrays.h @@ -0,0 +1,167 @@ +/* + * SPDX-License-Identifier: GPL-3.0-or-later + * Copyright (c) 2008-2024 San Diego State University Research Foundation + * (SDSURF). + * See LICENSE file or https://www.gnu.org/licenses/gpl-3.0.html for details. + */ + +/* + * @file MOLE_arrays.h + * + * @brief These classes define flat arrays for the MOLE library. + * These flat C++ arrays/vectors to optimize memory access, instead + * of using nested c++ vectors (e.g., std::vector>) + * To preserve compatibility with exising MOLE code, the flat arrays + * are designed to wrapped around the Armadillo library, which is + * used for some sparse and dense linear algebra operations. These + * interfaces should be extended in the future to work with other + * numerical libraries, such as PETSc (PETSc vectors and matrices), + * and libraries in the Trilinos project (Tpetra, Epetra, etc.). + * + * @date 2026/07/14 + * + */ + +#ifndef MOLE_ARRAYS_H +#define MOLE_ARRAYS_H + +#include +#include +#include +#include +#include "MOLE_errors.h" + +using Real = double; + +// array1D is a class designed to be compatible with other MOLE flat +// arrays and previous versions of MOLE that use Armadillo's vectors. +// The class also records errors during memory allocation or other +// vector operations. +class array1D { +protected: + mutable stack a_errs; // err detection+backtracking +public: + arma::vec data_; // flat array of doubles + array1D() = default; // default constructor (empty array) + // constructor for syntax array1D A(numelem, fillval) + array1D(size_t numelem, Real fillVal = 0.0); + // Array equality comparison (checks data + dimensions) + bool operator==(const array1D& other) const; + // Array inequality (stored at different memory locations) + bool operator!=(const array1D& other) const; + // Checks whether an index is valid for this array + bool valid_index(size_t i) const; + // resize a 1D array + check successful mem reallocation + void resize(size_t numelem, Real fillVal = 0.0); + // error handling methods + void logArrayError(size_t errCode, string errLoc, + string errParm) const; + bool hasArrayErrors() const; + void print_ErrorLog() const; + void write_ErrorLog() const; + void read_ErrorLog(int ErrorCode, std::string &location, + std::string &arrayName); +}; + +// array2D is a class designed to perform better than the 2D nested +// vectors in C++, which stores everything by rows of elements. These +// flat array implementations uses a memory heap to store all the +// elements of the two dimensional array. It also records errors +// during memory allocation. It also records errors during memory +// allocation or other vector operations. It is compatible with +// Armadillo's 2D matrix class (mat). + +class array2D { +protected: + mutable stack a_errs; // err detection+backtracking +public: + arma::mat data_; // flat array of doubles + array2D() = default; // default constructor (empty array) + // constructor for syntax array2D A(rows, cols, fillval) + array2D(size_t rows, size_t cols, Real fillVal = 0.0); + // Array equality comparison (checks data + dimensions) + bool operator==(const array2D& other) const; + // Array inequality (stored at different memory locations) + bool operator!=(const array2D& other) const; + // Checks whether a pair of indeces are valid for this array + bool valid_indeces(size_t i, size_t j) const; + // resize a 2D array + check successful mem reallocation + void resize(size_t rows, size_t cols, Real fillVal = 0.0); + void logArrayError(size_t errCode, string errLoc, + string errParm) const; + bool hasArrayErrors() const; + void print_ErrorLog() const; + void write_ErrorLog() const; + void read_ErrorLog(int ErrorCode, std::string &location, + std::string &arrayName); +}; + +// array3D is a class designed to perform better than the 3D nested +// vectors in C++, which stores everything by rows of elements. These +// flat array implementations uses a memory heap to store all the +// elements of the three dimensional array. It also records errors +// during memory allocation. It also records errors during memory +// allocation or other vector operations. It is compatible with +// Armadillo's 3D matrix class (cube). +class array3D { +protected: + mutable stack a_errs; // err detection+backtracking +public: + arma::cube data_; // flat array of doubles + array3D() = default; // default constructor (empty array) + // constructor for syntax array3D A(dim1, dim2, dim3, fillval) + array3D(size_t dim1, size_t dim2, size_t dim3, Real fillVal = 0.0); + bool operator==(const array3D& other) const; + // Array inequality (stored at different memory locations) + bool operator!=(const array3D& other) const; + // Checks whether a triplet of indeces are valid for this array + bool valid_indeces(size_t i, size_t j, size_t k) const; + // resize a 3D array + check successful mem reallocation + void resize(size_t dim1, size_t dim2, size_t dim3, + Real fillVal = 0.0); + void logArrayError(size_t errCode, string errLoc, + string errParm) const; + bool hasArrayErrors() const; + void print_ErrorLog() const; + void write_ErrorLog() const; + void read_ErrorLog(int ErrorCode, std::string &location, + std::string &arrayName); +}; + +// +// This function checks whether two arrays are numerically equal +// within a specified numerical tolerance. It is a template function +// that works for array1D, array2D, and array3D types. The function +// compares the number of elements and then checks if the absolute +// difference between corresponding elements is less than a tolFactor +// times the machine epsilon for double precision. +// This is useful when validating numerical results that may have +// small floating-point differences due to computation. +// +template +bool numEqualArray(const ArrayType& a1, const ArrayType& a2, + double tolFactor) { + // the array type is checked at compile time only + static_assert(std::is_same_v || + std::is_same_v || + std::is_same_v, + "numEqualArray only supports array1D, array2D, or array3D"); + + // checking elements and abs|a1(i)-a2(i)| < eps*tolFactor + if (a1.data_.n_elem != a2.data_.n_elem) return false; + + const double eps = std::numeric_limits::epsilon(); + const double* __restrict c = a1.data_.memptr(); + const double* __restrict u = a2.data_.memptr(); + + for (arma::uword k = 0; k < a1.data_.n_elem; ++k) { + double diff = std::fabs(c[k] - u[k]); + double mag = std::fabs(c[k]) > std::fabs(u[k]) ? + std::fabs(c[k]) : std::fabs(u[k]); + mag = mag > 1.0 ? mag : 1.0; + if (diff > tolFactor * eps * mag) return false; + } + return true; +} + +#endif // MOLE_ARRAYS_H \ No newline at end of file diff --git a/cpp/src/include/MOLE_errors.h b/cpp/src/include/MOLE_errors.h new file mode 100644 index 00000000..e0c776ae --- /dev/null +++ b/cpp/src/include/MOLE_errors.h @@ -0,0 +1,280 @@ +/* + * SPDX-License-Identifier: GPL-3.0-or-later + * Copyright (c) 2008-2024 San Diego State University Research + * Foundation (SDSURF). + * See LICENSE file or https://www.gnu.org/licenses/gpl-3.0.html + * for details. + */ + +/* + * @file MOLE_errors.h + * + * @brief Error handling for the MOLE library. + * + * @date 2026/06/24 + * + */ +#ifndef MOLE_ERRORS_H +#define MOLE_ERRORS_H + +#include +#include +#include + +using namespace std; +// MOLE errors are stored in a stack to allow tracking and +// backtracking of errors. Every error entry contains the error code, +// the location (MOLE function name) where the error occurred, and +// an optional input value that caused the error (this can be a +// concatenation of parameters). +#include // Using std stack class + +// +// MOLE DEBUGGING MODES +// +// We will define the following MOLE library debug modes. +// The MOLE library has been designed in a way that it does not +// break the execution of a user code, instead it creates and +// maintains an error logging mechanism that allows users to +// check for possible errors after a MOLE operation. In addition, +// the logged errors can be reported to std output or to a file. + +// 1. DEBUG_DEFAULT_MD (default behaviour - returns a grid +// that has errors) +// .hasError (e.g., hasGridErrors()) +// Users and application builders are responsible for checking +// for errors and implementing their error corrections. +// What support does a user get from the MOLE library at this point? +// A user can call MOLE functions to either: +// A: print the errors to standard output +// B: print errors to a file + +// 2. DEBUG_REPORTS_STDOUT_MD (has to be passed to the MOLE API - +// returns an object that has errors ) +// The MOLE library will report errors to standard output and return +// control to the users. What support does a user get from the MOLE +// library at this point? +// A user can call MOLE functions to either: +// A: print the errors to standard output +// B: print errors to a file + +// 3. DEBUG_AND_ABORT_MD (Report, then abort so a debugger stops at +// that particular failure point, the difference is that this one +// aborts execution). The MOLE library will report errors to standard +// output and abort the execution (e.g., code exits). + +// A debug mode governs only what happens to a MOLE object that +// failed validation. An object that validated is unaffected by the +// mode it was built with. The mode is not a property of the object; +// it is applied once, at construction, and is not stored. +#define DEBUG_DEFAULT_MD 0 +#define DEBUG_REPORTS_STDOUT_MD 1 +#define DEBUG_AND_ABORT_MD 2 + +// Predefined error codes for the MOLE library. These codes are used +// to identify specific errors that may occur during the execution +// of MOLE functions. These error symbols are used throughout the +// MOLE library to provide consistent error handling and reporting. +// +// Initialization Errors. 10 - 99 gridBuilder errors +// +#define MOLE_ERR_GRID_UNCHECKED 11 // Grid has not been validated +#define MOLE_ERR_INVALID_GRID_ARGS 12 // Grid is invalid +#define MOLE_ERR_GRID_CONSTRUCTION_FAILED 13 +#define MAKE_GRID_INVALID_INPUT_ARGS 14 +#define MAKE_GRID_MISSING_ARGS 15 +#define MAKE_GRID_UNKNOWN_ATTRIBUTE 16 +#define MAKE_GRID_DUPLICATE_ATTRIBUTES 17 +// +// Grid definition Errors. Error codes 100-199 +// +#define MOLE_ERR_INVALID_GRID_DIM 100 // Invalid grid dimension +#define MOLE_ERR_INVALID_GRID_TOPOLOGY 101 // Invalid grid topology +#define MOLE_ERR_INVALID_GRID_SPACING 102 // Invalid grid spacing +#define MOLE_ERR_INVALID_GRID_SIZE 103 // Invalid grid size +#define MOLE_ERR_GRID_NODAL_SZ_MISMATCH 104 // nodal size mismatch +#define MOLE_ERR_GRID_CENTERS_SZ_MISMATCH 105 // center size mismatch +#define MOLE_ERR_GRID_FACES_SZ_MISMATCH 106 // faces size mismatch +#define MOLE_ERR_INVALID_INPUT_TYPE 107 // Invalid input type +#define MOLE_ERR_ARRAY_HAS_NULL_POINTER 108 // Array has null pointer +#define MOLE_ERR_INVALID_CELL_COUNT 109 // Invalid grid cell count +#define MOLE_ERR_INVALID_ISPERIODIC_DIM 110 // invalid array dim +#define MOLE_ERR_ISPERIODIC_TYPE 111 // invalid array type +#define MOLE_ERR_INVALID_CURVILINEAR_GRID 112 // Need nodal coordines +#define MOLE_ERR_INVALID_1D_CURVILINEAR 113 // 1D Curvilinear invalid +#define MOLE_ERR_INVALID_NONUNIFORM_GRID 114 // Need nodal coordines +#define MOLE_ERR_INVALID_ARRAY_INDEX 115 // invalid array indexing +#define MOLE_ERR_INVALID_NODAL_COORDINATES 116 // invalid user coords +#define MOLE_ERR_INVALID_CENTER_COORDINATES 117 // invalid user coords +#define MOLE_ERR_INVALID_NORMAL_FACE_COORDS 118 // invalid coords + +// +// Flat array defition Errors. Error codes 200-299 +// +#define MOLE_ERR_INVALID_ARRAY_SIZE 200 //invalid array sizes +#define MOLE_ERR_ARRAY_SIZE_OVERFLOW 201 // array allocation overflow +#define MOLE_ERR_ARRAY_INDEX_OUTBOUNDS 202 //array index out of bounds +#define MOLE_ERR_FAILED_ARRAY_ALLOC 203 // array alloc failed +#define MOLE_ERR_FAILED_ARRAY_RESIZE 204 // array resize failed + +// +// MD Operators and Numerical Errors codes 300 - 600 +// +#define MOLE_ERR_DIVISION_BY_ZERO 300 +#define MOLE_ERR_INF_VALUE 301 +#define MOLE_ERR_NAN_VALUE 302 + +// +// Dictionary of error codes and their corresponding messages +// for printing error messages in the MOLE library. +// +static unordered_map MOLE_errors_messages = { + // 011 + {MOLE_ERR_GRID_UNCHECKED, + "Grid has not been validated, call validateGrid() first"}, + // 012 + {MOLE_ERR_INVALID_GRID_ARGS, + "Error(s) in input parameters, resulting MOLE grid is invalid."}, + // 013 + {MOLE_ERR_GRID_CONSTRUCTION_FAILED, + "Grid construction failed, review the list of errors"}, + // 014 + {MAKE_GRID_INVALID_INPUT_ARGS, + "Grid attribute value has the wrong type for that attribute"}, + // 015 + {MAKE_GRID_MISSING_ARGS, + "A required grid attribute was not supplied"}, + // 016 + {MAKE_GRID_UNKNOWN_ATTRIBUTE, + "Not a MOLE grid attribute name. Parsing stops here, because " + "the type of the value following an unknown name is unknown"}, + // 017 + {MAKE_GRID_DUPLICATE_ATTRIBUTES, + "Grid attribute supplied more than once"}, + // 100 + {MOLE_ERR_INVALID_GRID_DIM, + "Invalid grid dimension entered. Valid values are 1, 2, or 3"}, + // 101 + {MOLE_ERR_INVALID_GRID_TOPOLOGY, + "Invalid grid topology entered. Valid values are 'u'=uniform " + "or 'c'=curvilinear) or 'n'=non-uniform"}, + // 102 + {MOLE_ERR_INVALID_GRID_SPACING, + "Grid spacing must be > 0.0 and a valid, finite real number"}, + // 103 + {MOLE_ERR_INVALID_GRID_SIZE, + "Grid size must be a natural number > 0"}, + // 104 + {MOLE_ERR_GRID_NODAL_SZ_MISMATCH, + "Nodal grid coordinates array size mismatch"}, + // 105 + {MOLE_ERR_GRID_CENTERS_SZ_MISMATCH, + "Center grid coordinates array size mismatch"}, + // 106 + {MOLE_ERR_GRID_FACES_SZ_MISMATCH, + "Normal faces coordinates array size mismatch"}, + // 107 + {MOLE_ERR_INVALID_INPUT_TYPE, + "Invalid input type for input name"}, + // 108 + {MOLE_ERR_ARRAY_HAS_NULL_POINTER, + "Array has a null pointer. Array allocation may have failed"}, + // 109 + {MOLE_ERR_INVALID_CELL_COUNT, "Non-positive cell count"}, + // 110 + {MOLE_ERR_INVALID_ISPERIODIC_DIM, + "Wrong size for isPeriodic 1D:(1x1), 2D:(2x1), 3D:(3x1)"}, + // 111 + {MOLE_ERR_ISPERIODIC_TYPE, + "isPeriodic requires a boolean array"}, + // 112 + {MOLE_ERR_INVALID_CURVILINEAR_GRID, + "Curvilinear grids need user-provided nodal coordinates"}, + // 113 + {MOLE_ERR_INVALID_1D_CURVILINEAR, + "Curvilinear grids fundamentally cannot be one dimensional"}, + // 114 + {MOLE_ERR_INVALID_NONUNIFORM_GRID, + "Non-uniform grids need user-provided nodal coordinates"}, + // 115 + {MOLE_ERR_INVALID_ARRAY_INDEX, + "One or more indices to the array are out of bound," + "check array dimensions"}, + // 116 + {MOLE_ERR_INVALID_NODAL_COORDINATES, + "User-provided nodal coordinates do not agree with other " + "uniform grid parameters passed"}, + // 117 + {MOLE_ERR_INVALID_CENTER_COORDINATES, + "User-provided center coordinates do not agree with other " + "uniform grid parameters passed"}, + // 118 + {MOLE_ERR_INVALID_NORMAL_FACE_COORDS, + "User-provided normal face coordinates do not agree with " + "other uniform grid parameters passed"}, + // 200 + {MOLE_ERR_INVALID_ARRAY_SIZE, + "Array dimensions need to be natural numbers >= 1"}, + // 201 + {MOLE_ERR_ARRAY_SIZE_OVERFLOW, + "Array dimensions too large to fit in memory"}, + // 202 + {MOLE_ERR_ARRAY_INDEX_OUTBOUNDS, + "Array index is out of bounds." }, + // 203 + {MOLE_ERR_FAILED_ARRAY_ALLOC, "Array allocation failed."}, + // 204 + {MOLE_ERR_FAILED_ARRAY_RESIZE, "Array resize operation failed."}, + // 300 + {MOLE_ERR_DIVISION_BY_ZERO, "Division by zero."}, + // 301 + {MOLE_ERR_INF_VALUE, "Infinity value detected."}, + // 302 + {MOLE_ERR_NAN_VALUE, "NaN value detected."}, +}; + +struct MOLE_Errors { + int errCode; // Logged error codes + string errLocation; // Report location where the error occurred + string paramError; // Holds additional information on the error +}; + +// MOLE_Error General Stack Operations +// 1. Initializes the error stack +void MOLEerr_init(stack& errorStack); + +// 2. Pushes an error onto the error stack +void MOLEerr_log(stack& errorStack, int errCode, + const string& location, const string& inputParam = ""); + +// 3. Checks whether the error stack contains a specific error code +bool MOLEerr_contains(const stack& errorStack, int targetCode); + +// 4. Checks whether there are any errors in the stack. A grid's +// stack carries MOLE_ERR_GRID_UNCHECKED from construction until +// validation removes it, so this reports true on a grid that has +// not been validated yet, and on a valid grid that merged errors +// from an upstream MOLE object. Use isValidatedGrid() to ask +// whether the object itself is usable. +bool MOLEerr_haserrors(const stack& errorStack); + +// 5. Removes a specific error from the stack +void MOLEerr_remove(stack& errorStack, int targetCode); + +// 6. These functions are used for printing the error stack to +// standard output +void MOLEerr_print(const stack& errorStack); + +// 7. Writes error output to a log file (not implemented yet) +void MOLEerr_dumpErrLog(stack& errorStack, string logType); + +// ParamsNull is a struct used for reporting errors with any +// MOLE objects (classes). The structure can accumulate errors from +// other MOLE classes (i.e, grids, boundaries, etc)in type_errs +struct paramsNull { + size_t num_errs = 0; // Number of errors + stack type_errs; // Types of MOLE errors +}; + + +#endif // MOLE_ERRORS_H diff --git a/cpp/src/include/MOLE_grids.h b/cpp/src/include/MOLE_grids.h new file mode 100644 index 00000000..f87a66d3 --- /dev/null +++ b/cpp/src/include/MOLE_grids.h @@ -0,0 +1,376 @@ +/* + * SPDX-License-Identifier: GPL-3.0-or-later + * Copyright (c) 2008-2024 San Diego State University Research Foundation + * (SDSURF). + * See LICENSE file or https://www.gnu.org/licenses/gpl-3.0.html for details. + */ + +/* + * @file MOLE_grids.h + * + * @brief Generic grid class and data structures for 1D, 2D, and 3D grids. + * + * @date 2026/06/24 + * + */ + +#ifndef MOLE_GRIDS_H +#define MOLE_GRIDS_H + +#include +#include +#include // To handle C++ variant classes and structs +#include //For polymorphism in visit's lambda functions +#include +#include "MOLE_errors.h" +#include "MOLE_arrays.h" +#include "utils.h" + + +using namespace std; + +// ----------------------------------------------------------------// +// +// MOLE GRID DATA STRUCTURES +// +// The following 3 structs are used inside the Grid subclasses +// (i.e., grid1D, grid2D, grid3D). +// +// ----------------------------------------------------------------// +// +// 1D grid specific parameters +// - dx is a Real for uniform grids >= 0.0 +// - for other topologies, the user must provide nodes_X +// +struct gridParams1D { + char topology; // 'u'=uniform, 'c'=curvilinear, 'n'=non-uniform + size_t m = 0; // Number of cells in x-direction + Real dx = 0.0; // Uniform cell spacing in x-direction + array1D nodes_X; // nodes_X = faces_X in 1D + array1D centers_X; // Staggered grid in x-dir for MOLE Ops + // Faces_X is always identical to nodes_X in 1D. The attribute + // is provided for consistency with 2D and 3D grids. Make sure + // to call it as a function in 1D cases. + array1D& Faces_X(){return nodes_X;} // In 1D same as nodes_X + const array1D& Faces_X() const { return nodes_X; } + bool bc_isPeriodic = false; // whether the grid has periodic bcs +}; + +// +// 2D grid specific parameters +// - dx and dy are Reals for uniform grids but for nonuniform +// there are arrays of Reals +// - For curvilinear and nonuniform grids users need to provide +// the nodes_X, and nodes_Y arrays +struct gridParams2D { + char topology; // 'u'=uniform, 'c'=curvilinear, 'n'=non-uniform + size_t m = 0; // Number of cells in x-direction + size_t n = 0; // Number of cells in y-direction + Real dx = 0.0; // Uniform cell spacing in x-direction + Real dy = 0.0; // Uniform cell spacing in y-direction + array2D nodes_X; // Nodal coordinates in x-dir + array2D nodes_Y; // Nodal coordinates in y-dir + array2D centers_X; // Staggered grid in x-dir for MOLE Ops + array2D centers_Y; // Staggered grid in y-dir for MOLE Ops + array2D faces_u_X; // x-normal faces + array2D faces_u_Y; // y-normal faces + array2D faces_v_X; // x-normal faces + array2D faces_v_Y; // y-normal faces + bool bc_isPeriodic[2] = {false, false}; //periodic bcs? +}; + +// +// 3D grid specific parameters +// - dx, dy and dz are Reals for uniform grids but arrays of Reals +// for nonuniform grids +// - For curvilinear and nonuniform grids users need to provide +// the nodes_X, nodes_Y, and nodes_Z arrays +struct gridParams3D { + char topology; // 'u'=uniform, 'c'=curvilinear, 'n'=non-uniform + size_t m = 0; // Number of cells in x-direction + size_t n = 0; // Number of cells in y-direction + size_t o = 0; // Number of cells in z-direction + Real dx = 0.0; // Uniform cell spacing in x-direction + Real dy = 0.0; // Uniform cell spacing in y-direction + Real dz = 0.0; // Uniform cell spacing in z-direction + array3D nodes_X; // Nodal coordinates in x-dir + array3D nodes_Y; // Nodal coordinates in y-dir + array3D nodes_Z; // Nodal coordinates in z-dir + array3D centers_X; // Staggered grid in x-dir for MOLE Ops + array3D centers_Y; // Staggered grid in y-dir for MOLE Ops + array3D centers_Z; // Staggered grid in z-dir for MOLE Ops + array3D faces_u_X; // x-normal faces + array3D faces_u_Y; // y-normal faces + array3D faces_u_Z; // z-normal faces + array3D faces_v_X; // x-normal faces + array3D faces_v_Y; // y-normal faces + array3D faces_v_Z; // z-normal faces + array3D faces_w_X; // x-normal faces + array3D faces_w_Y; // y-normal faces + array3D faces_w_Z; // z-normal faces + bool bc_isPeriodic[3] = {false, false, false}; +}; + +// -----------------------------------------------------------------// +// +// gridBase is the grid superclass that handles common interfaces for +// all grid dimensions (e.g., error detection and propagation). This +// class can also be used for errors that occur prematurally when a +// grid is created. +// +// gridBase is internal to MOLE. Its constructor is protected, so it +// can only be reached through grid1D, grid2D, grid3D, or gridNull. +// +// ----------------------------------------------------------------// +// The error stack is member of this superclass and protected (i.e +// can only be access by member functions in a grid class) +class gridBase{ + protected: + stack errs; // for error detection+backtracking + gridBase(size_t idim); + public: + // gridBase shell members for all grids + size_t dim; // explict declaration of dimensionality in grid + virtual ~gridBase() = default; + // validates user-provided coordinates for uniform grids only + bool valid1DCoordinates(array1D& userInput, + const array1D& expected, + Real dx, size_t m, + int sizeMismatchErr, + int badCoordsErr); + bool valid2DCoordinates(array2D& userInput, + const array2D& expected, + Real dx, Real dy, size_t m, size_t n, + int sizeMismatchErr, + int badCoordErr); + bool valid3DCoordinates(array3D& userInput, + const array3D& expected, + Real dx, Real dy, Real dz, + size_t m, size_t n, size_t o, + int sizeMismatchErr, + int badCoordErr); + // Error handling member functions/methods + void logGridErr(size_t errCode, string errLoc, + string errParm); + bool hasGridErrors(); + bool isValidatedGrid(); + void setGridValidated(); + void print_ErrorLog(); + void write_ErrorLog(); + + // + // applyDebugMode implements the MOLE debug modes for grid + // objects. A grid that passed validation is left untouched; + // the mode governs only what happens to a grid that did not. + // The mode is applied once, by the constructor that received + // it, and is not stored on the grid. + // + // The test is the validation state rather than the presence + // of errors: a grid's stack carries MOLE_ERR_GRID_UNCHECKED + // until validation clears it, and mergeErrors folds upstream + // errors into the same stack, so a valid grid can hold + // errors it did not cause. + // + void applyDebugMode(size_t debug_mode); + + // ---------------------------------------------------------- + // The following member function are used as helpers for the + // MOLE grid classes (grid1D, grid2D, grid3D) to avoid code + // duplication. These are used for grid validatation and + // automatic generation, where the different coordinates need + // to "allocate coordinate array -> build mesh -> propagate + // errors -> validate". + // ---------------------------------------------------------- + + // + // 1) drainArrayErrors propagates errors logged onto an array + // or a different MOLE object up to the corresponding grid's + // error stack. + // + template + void drainArrayErrors(ArrayT& a) { + while (a.hasArrayErrors()) { + int errCode; + string location, msgparam; + a.read_ErrorLog(errCode, location, msgparam); + logGridErr(errCode, location, msgparam); + } + } + + // + // 2) mergeErrors combines a previously-logged error stack + // (inerrs) into a grid's error stack. + // + void mergeErrors(const stack& inerrs); + + // + // 3) buildOrCheck2DCoords builds or validates a user-supplied + // 2D coordinate (nodal, cell-centers, or normal faces). When + // building the coordinates, outX and outY will have the + // output coordinate arrays. However, when validating, the + // user-supplied coordinates are passed in outX and outY + // and validated against the expected coordinates (using + // a tolerance value for floating point comparisons). + // + bool buildOrCheck2DCoords(array2D& outX, array2D& outY, + const array1D& xcoord, + const array1D& ycoord, Real dx, + Real dy, size_t m, size_t n, + int sizeMismatchErr, + int badCoordsErr); + + // + // 4) buildOrCheck3DCoords builds or validates a user-supplied + // 3D coordinate (nodal, cell-centers, or normal faces). When + // building the coordinates, outX, outY, and outZ will have the + // output coordinate arrays. However, when validating, the + // user-supplied coordinates are passed in outX, outY, and outZ + // and validated against the expected coordinates (using + // a tolerance value for floating point comparisons). + // + bool buildOrCheck3DCoords(array3D& outX, array3D& outY, + array3D& outZ, + const array1D& xcoord, + const array1D& ycoord, + const array1D& zcoord, + Real dx, Real dy, Real dz, + size_t m, size_t n, size_t o, + int sizeMismatchErr, + int badCoordsErr); +}; + +// ---------------------------------------------------------------- +// Each grid class offers four constructors. The two that take a +// debug_mode delegate to the matching constructor without one, then +// apply the mode to the result. Valid modes are declared in +// MOLE_errors.h: DEBUG_DEFAULT_MD, DEBUG_REPORTS_STDOUT_MD, and +// DEBUG_AND_ABORT_MD. +// ---------------------------------------------------------------- + +// 1D grid class +class grid1D : public gridBase{ + public: + gridParams1D grid; // structure for all 1D grid parameters + ~grid1D() = default; + // grid1D() constructor, users pass a structure with + // the 1D parameters (gridParams1D). + grid1D(const gridParams1D p1); + // like grid1D but applies a MOLE debug mode to the result + grid1D(const gridParams1D p1, size_t debug_mode); + // like grid1D but also takes an MOLE error object containing + // errors that may have occur prior to the grid construction + grid1D(const gridParams1D p1, + const stack& inerrs); + // like the constructor above but applies a MOLE debug mode + grid1D(const gridParams1D p1, + const stack& inerrs, + size_t debug_mode); + bool validGrid(); +}; + +// 2D Grid class +class grid2D : public gridBase{ + public: + gridParams2D grid; // 2D grid structure + ~grid2D() = default; + // grid2D() constructor, users pass a structure with + // the 2D parameters (gridParams2D). + grid2D(const gridParams2D p2D); + // like grid2D but applies a MOLE debug mode to the result + grid2D(const gridParams2D p2D, size_t debug_mode); + // like grid2D but also takes an MOLE error object containing + // errors that may have occur prior to the grid construction + grid2D(const gridParams2D p2, + const stack& inerrs); + // like the constructor above but applies a MOLE debug mode + grid2D(const gridParams2D p2, + const stack& inerrs, + size_t debug_mode); + // checks if the grid params are valid + bool validGrid(); +}; + +// 3D Grid class +class grid3D : public gridBase{ + public: + gridParams3D grid; // 3D grid structure + ~grid3D() = default; + // grid3D() constructor, users pass a structure with + // the 3D parameters (gridParams2D). + grid3D(const gridParams3D p3D); + // like grid3D but applies a MOLE debug mode to the result + grid3D(const gridParams3D p3D, size_t debug_mode); + // like grid3D but also takes an MOLE error object containing + // errors that may have occur prior to the grid construction + grid3D(const gridParams3D p3, + const stack& inerrs); + // like the constructor above but applies a MOLE debug mode + grid3D(const gridParams3D p3, + const stack& inerrs, + size_t debug_mode); + // checks if the grid params are valid + bool validGrid(); +}; + +// Grid Null Class - This grid object is return when errors are found +// with the grid specifications and the grid could not be built. +// +// gridNull is internal to MOLE. Its constructor is private and only +// makeGridNull can reach it, so a user cannot fabricate one. Copy +// and move stay public because gridVar holds a gridNull alternative +// and has to be copyable. +class gridNull: public gridBase{ + private: + gridNull(const paramsNull in_p, + const stack& inerrs); + friend gridNull makeGridNull(const paramsNull in_p, + const stack& inerrs); + public: + paramsNull ErrData; + ~gridNull() = default; + bool validGrid(){return false;} +}; + +// makeGridNull is the only way to build a gridNull. It exists so +// that makeGrid can produce the failure case without opening the +// gridNull constructor to users. +gridNull makeGridNull(const paramsNull in_p, + const stack& inerrs); + +// ---------------------------------------------------------------- +// VARIANT STRUCTURES AND CLASSES +// variant for data structures holding MOLE's grid information +// or ParamsNull = for reporting error when a grid cannot be +// generated (ParamsNull is declared in MOLE_Errors.h) +using paramVars = std::variant; +// Handling MOLE variant classes and data structures. These are +// variant over the concrete derived types — used purely for +// type-safe dispatch for MOLE grid classes and generic grid +// handling. gridNull is an empty shell created whenever errors +// occur and it is unsafe to use the grid. +// ---------------------------------------------------------------- +using gridVar = std::variant; + +// Auxiliar functions used by some MOLE grids + +bool isValidGrid(gridVar& g); // validates a generic grid +bool validSpacing(Real dh); // checks for valid dx, dy, or dz +void generateNodalPts(size_t npts, Real delta, array1D& out_array); +void generateCenterPts(size_t npts, Real delta, array1D& out_array); + +// +// nd2DGrid/nd3DGrid: GNU Octave ndgrid analogs that build 2D +// and 3D coordinates. +// +void nd2DGrid(const array1D& x, const array1D& y, + array2D& OutX, array2D& OutY); +void nd3DGrid(const array1D& x, const array1D& y, const array1D& z, + array3D& OutX, array3D& OutY, array3D& OutZ); + +// gridVar makeGrid is a factory function that works for any of the 3 +// grid dimensionalities, intended for cases when the users need to +// define the dimensionality at runtime. +gridVar makeGrid(paramVars params, const stack& errs); + +#endif // MOLE_GRIDS_H diff --git a/cpp/src/include/README.md b/cpp/src/include/README.md new file mode 100644 index 00000000..0cd0fc4a --- /dev/null +++ b/cpp/src/include/README.md @@ -0,0 +1,41 @@ + + +# Subdirectory for MOLE 2.0 C++ Header Files + +Subdirectory and Pathname: **mole/cpp/src/include** + +## Purpose + +All header files (.h) for the different MOLE 2.0 C++ functionality +are contain in this directory. MOLE examples and tests have their own +include subdirectories. + +Header files contain public declarations of the MOLE C++ API, which +is the MOLE collection of functionalities available to MOLE users. + +## MOLE Modules and Header Files (alphabetical order) + ++ **MOLE_arrays.h**: MOLE arrays that interface with data structures +from other numerical libraries. The current implementation works with +Armadillo's matrices, vectors and cubes. ++ **MOLE_errors.h**: a MOLE error handling, tracking and reporting +mechanism for the library. Users have full control on how to handle +exceptions, the MOLE library only reports them and does not cause an +application's execution to stop. ++ **MOLE_grids.h**: MOLE grid classes implementation, including error +handling and reporting ++ **README.md**: (this file) ++ **utilis.h**: header file for all MOLE utility functions, which +include implementations that support the MOLE library operations. diff --git a/cpp/src/include/grid_builder.h b/cpp/src/include/grid_builder.h new file mode 100644 index 00000000..15dd4b61 --- /dev/null +++ b/cpp/src/include/grid_builder.h @@ -0,0 +1,142 @@ +#ifndef MOLE_GRID_BUILDER_H +#define MOLE_GRID_BUILDER_H + +// grid_builder.h +// +// MOLE Grid Utility: parses a user's input into +// an internal gridRaw, then narrows gridRaw into a dimensionality +// specific gridParams[1D][2D][3D] used inside MOLE 2.0 to construct +// grids. Grid builder will instantiate a temporaty gridRaw structure +// which will be used in the construction of a Grid1D, Grid2D or +// Grid3D object. MOLE grid objects are created and validated during +// construction. +// The wrapper performs only these parse-time checks: +// 1. invalid key -> MAKE_GRID_UNKNOWN_ATTRIBUTE +// 2. dimension consistency -> MOLE_ERR_INVALID_CELL_COUNT +// 3. invalid dimensionality -> MOLE_ERR_INVALID_GRID_DIM +// 4. invalid topology -> MOLE_ERR_INVALID_GRID_TOPOLOGY +// 5. isPeriodic size -> MOLE_ERR_INVALID_ISPERIODIC_DIM +// +// Users need to input the grid's dimensionaly and topology. Thus, +// these are not inferred nor defaulted. This utility uses MOLE's +// error reporting and tracking mechanism. +// +// SYNTAX: +// ------------------------------------------------------------------ +// gridVar g = gridBuilder(va_arg); +// where va_arg is a list of pairs of the form: +// +// Example: +// gridVar g=gridBuilder("dim",1, "m",20, "dx",0.2, "topology",'u'); +// +// isPeriodic takes the address of a vector holding one flag +// per dimension. The vector must outlive the gridBuilder call: +// std::vector per = {true, false}; +// gridVar g = gridBuilder("dim", 2, ..., "isPeriodic", &per); +// ------------------------------------------------------------------ +// Note I: callers who know the dimension at compile time can call +// the MOLE grid constructors directly by using the corresponding +// gridParams1D/2D/3D directly, and invoking the corresponding grid +// constructor, Grid1D/2D/3D, respectively. +// +// Note II: there are not type-checking for va_arg at compile time, +// users need to be aware of the correct data types: +// ATTRIBUTE // +// TYPE: Name: Valid C++ Input Type // +// -----------------------:------------------:---------------------// +// counters m, n, o, dim const int () +// cell spacing dx, dy, dz const double +// grid topology topology const char +// grid coordinates Nodal, Centers, const vector (doubles) +// Normal faces const &flatNDArray +// grid periodicity isPeriodic const vector* +// (size must equal dim) +// debugging mode debug const int +// (MOLE debug mode) +// +// Note III: the debug attribute takes one of the MOLE debug modes +// declared in MOLE_errors.h (DEBUG_DEFAULT_MD, +// DEBUG_REPORTS_STDOUT_MD, DEBUG_AND_ABORT_MD) and governs what +// gridBuilder does with a grid that fails to build. Parsing stops +// at an unrecognized attribute name, because the type of the value +// after an unknown name is unknown too. A debug pair placed after +// one is therefore never read, on exactly the calls that need it. +// Pass debug as the first pair: +// gridVar g = gridBuilder("debug", DEBUG_REPORTS_STDOUT_MD, +// "dim", 1, "m", 20, "dx", 0.2, +// "topology", 'u'); +// +// ------------------------------------------------------------------------- + +#include +#include +#include "MOLE_grids.h" + +// The variadic list, va_arg, needs to end with a nullptr. A va_arg +// without an nullptr at end fails. Therefore, gridBuilder is +// designed as a macro that always appends a nullpointer sentinel, +#define gridBuilder(...) \ + ::gridBuilder_impl(__VA_ARGS__, static_cast(nullptr)) + +// gridRaw is the gridBuilder's generic structure. +struct gridRaw { + int dim = -1; // required; -1 means the user omitted it + char topology = '\0'; // 'u'|'c'|'n'; required, no default + + int m = -1; + int n = -1; + int o = -1; + + Real dx = 0.0; + Real dy = 0.0; + Real dz = 0.0; + + // MOLE debug mode applied to a grid that fails to build. See + // Note III above for why the debug pair has to come first. + int debug = DEBUG_DEFAULT_MD; + + // points at the caller's vector; null when not supplied. The + // vector carries its own size, which runChecks compares with + // dim regardless of the order the attributes arrived in. + const std::vector* isPeriodicSrc = nullptr; + + const void* nodesX = nullptr; + const void* nodesY = nullptr; + const void* nodesZ = nullptr; + const void* centersX = nullptr; + const void* centersY = nullptr; + const void* centersZ = nullptr; + const void* facesuX = nullptr; + const void* facesuY = nullptr; + const void* facesuZ = nullptr; + const void* facesvX = nullptr; + const void* facesvY = nullptr; + const void* facesvZ = nullptr; + const void* faceswX = nullptr; + const void* faceswY = nullptr; + const void* faceswZ = nullptr; +}; + +// runChecks runs the five parse-time validations, and logs failures +// to errs. It also returns the grid dimensionality, or a 0 whenever +// the dim attribute is missing or <= 0 or > 3. +int runChecks(std::stack& errs, const gridRaw& g); + +// gridBuilder_impl is the entry point behind the gridBuilder macro: +// it parses, checks for errors, and instantiates the gridParams. It +// calls MOLE's makeGrid and returns an instant of a MOLE gridVar a +// c++ variant for the MOLE classes Grid[1D][2D][3D]. When a grid +// cannot be constructed it returns a gridNull, built via +// makeGrid(paramsNull, errs) (user checks errors). In MOLE_grids.h: +// using gridVar = std::variant; +// Before returning, the debug mode is applied to the resulting +// grid. The mode governs only what gridBuilder reports; the grid +// keeps its full error log in every mode. +gridVar gridBuilder_impl(const char* firstName, ...); + +// makeGrid (the factory that turns a paramVars into a gridVar) is +// declared in MOLE_grids.h; gridBuilder_impl calls it for both the +// success path (gridParams1D/2D/3D) and the failure path +// (paramsNull -> gridNull). + +#endif // MOLE_GRID_BUILDER_H diff --git a/src/cpp/utils.h b/cpp/src/include/utils.h similarity index 79% rename from src/cpp/utils.h rename to cpp/src/include/utils.h index 1318bed0..bac9ecd3 100644 --- a/src/cpp/utils.h +++ b/cpp/src/include/utils.h @@ -8,6 +8,7 @@ * @file utils.h * @brief Helpers for sparse operations and MATLAB/Octave analogs * @date 2024/10/15 + * New Implementation for MOLE 2.0 Last Modified 2026/07/22 * */ @@ -17,6 +18,7 @@ #define UTILS_H #include +#include "MOLE_errors.h" using Real = double; using namespace arma; @@ -26,6 +28,8 @@ using namespace arma; * */ class Utils { +protected: + std::stack errs; public: /** @@ -84,7 +88,7 @@ class Utils { * @param Y a sparse matrix, will be filled by the function * */ - void meshgrid(const vec &x, const vec &y, mat &X, mat &Y); + void mesh2Dgrid(const vec &x, const vec &y, mat &X, mat &Y); /** * @brief An analog to the MATLAB/Octave 3D meshgrid operation @@ -101,7 +105,7 @@ class Utils { * @param Z a sparse matrix, will be filled by the function * */ - void meshgrid(const vec &x, const vec &y, const vec &z, cube &X, cube &Y, + void mesh3Dgrid(const vec &x, const vec &y, const vec &z, cube &X, cube &Y, cube &Z); /** * @brief Implements the trapezoidal rule for 1D numerical integration @@ -113,26 +117,13 @@ class Utils { * @param y Vector of y-values at corresponding x * @return Estimated area under the curve */ - static double trapz(const vec &x, const vec &y); -}; - -namespace mole { - -/** - * @brief Validate a cell-spacing argument. - * - * Throws std::invalid_argument if @p h is zero, negative, NaN, or Inf. - * Called at the top of every operator constructor and AddScalarBC free - * function that takes a dx / dy / dz argument, so validation is active - * in both Debug and Release builds (unlike assert(), which vanishes - * under NDEBUG). - * - * @param h Spacing value to check. - * @param name Parameter name for the error message (e.g. "dx"). - * @throws std::invalid_argument if h is not a positive finite number. - */ -void check_spacing(Real h, const char* name); + double trapz(const vec &x, const vec &y); -} // namespace mole + bool hasErrors(); + void print_ErrorLog(); + void write_ErrorLog(); + void read_ErrorLog(int ErrorCode, std::string &location, + std::string &arrayName); +}; #endif // UTILS_H diff --git a/cpp/src/sys/MOLE_Errors.cpp b/cpp/src/sys/MOLE_Errors.cpp new file mode 100644 index 00000000..c1374761 --- /dev/null +++ b/cpp/src/sys/MOLE_Errors.cpp @@ -0,0 +1,215 @@ +/* +* SPDX-License-Identifier: GPL-3.0-or-later +* © 2008-2024 San Diego State University Research Foundation (SDSURF). +* See LICENSE file or https://www.gnu.org/licenses/gpl-3.0.html for details. +*/ + +/* + * @file mole_errors.cpp + * + * @brief MOLE Error Handling and Tracking + * + * @date 2026/06/24 + * + */ + + #include "MOLE_errors.h" + #include + #include + + // This function initializes the error stack by + // clearing any existing errors + void MOLEerr_init(stack& errorStack){ + errorStack = stack(); // Clear the stack +} + +// This function pushes an error onto the error stack +// (default for inputParam is already given in the MOLE_Errors.h +// declaration; C++ forbids repeating it here) +void MOLEerr_log(stack& errorStack, int errCode, + const string& location, const string& inputParam){ + MOLE_Errors err; + err.errCode = errCode; + err.errLocation = location; + err.paramError = inputParam; + errorStack.push(err); +} + +// This is an auxiliary function that checks whether a specific error +// code exists in the error stack. It returns true if the error code +// is found, otherwise false. +bool MOLEerr_contains(const stack& errorStack, int targetCode) { + stack tempStack = errorStack; + while (!tempStack.empty()) { + if (tempStack.top().errCode == targetCode) return true; + tempStack.pop(); + } + return false; +} + +// This function checks whether there are any errors logged in +// the error stack. It returns true if there are errors. +bool MOLEerr_haserrors(const stack& errorStack) +{ + return !errorStack.empty(); +} + +// This function removes a specific error from the error stack. +// errorStack is passed by reference. The targetCode is the error +// code to be removed. +void MOLEerr_remove(stack& errorStack, int targetCode) { + stack tempStack; + + // Pop everything, keep the elements we want, reverse order + // into tempStack + while (!errorStack.empty()) { + MOLE_Errors err = errorStack.top(); + errorStack.pop(); + if (err.errCode != targetCode) { + tempStack.push(err); + } + } + + // tempStack now has kept elements in reverse order — + // push back to restore the backtracing original order + while (!tempStack.empty()) { + errorStack.push(tempStack.top()); + tempStack.pop(); + } +} + +// This auxiliaryfunction prints out any error arguments logged +// in the error stack. It is used inside MOLEerr_print. +void MOLEerr_print_args(const string& strargx) { + if (!strargx.empty()) { + cout << "with arg value(s): " << strargx << endl; + } else { + cout << endl; + } +} + +// writeErrtoStOut writes an error to the log_error file +void writeErrtoStdOut(int errNum, int errCode, string errLocation, + string errParams, const std::string& errMsg) { + cout << "Error #" << errNum << ": MOLE Error code [" << + errCode << "] - "<< errMsg << endl; + cout << "occurred inside:" << errLocation; + MOLEerr_print_args(errParams); +} + +// Produces a timestamp string to create a unique filename to output +// log_errors messages and backtracing information +std::string getDateTimeString() { + std::time_t now = std::time(nullptr); + char buf[100]; + std::strftime(buf, sizeof(buf), "%Y%m%d_%H%M%S", + std::localtime(&now)); + return std::string(buf); +} + +// MOLEerr_write_args is similar to MOLEerr_print_args but it writes +// to a log_error file instead of standard output +void MOLEerr_write_args(std::ofstream& ofile, + const string& strargx) { + if (!strargx.empty()) { + ofile << "with arg value(s): " << strargx << endl; + } else { + ofile << endl; + } +} + +// writeErrtoFile writes an error to the log_error file +void writeErrtoFile(std::ofstream& ofile, int errNum, int errCode, + string errLocation, string errParams, + const std::string& errMsg) { + ofile << "Error #" << errNum << ": MOLE Error code [" << + errCode << "] - "<< errMsg << endl; + ofile << "occurred inside:" << errLocation; + MOLEerr_write_args(ofile, errParams); +} + +// This function prints out to std output any errors logged in a +// stack, preserving the error stack for further enabled debugging. +void MOLEerr_print(const stack& errorStack){ + // Create a copy to preserve the original stack + stack tempStack = errorStack; + if (tempStack.empty()) { + cout << "No errors logged." << endl; + return; + } + cout << "========================================"<< endl; + cout << "Backtracing all logged MOLE Errors :" << endl; + cout << "========================================"<< endl; + + int i = 1; + while (!tempStack.empty()) { + MOLE_Errors err = tempStack.top(); + tempStack.pop(); + if (err.errCode == MOLE_ERR_GRID_UNCHECKED) { + writeErrtoStdOut(i, err.errCode, err.errLocation, + err.paramError, MOLE_errors_messages[err.errCode]); + i++; + continue; + } + auto errMsg = MOLE_errors_messages.find(err.errCode); + if (errMsg != MOLE_errors_messages.end()) { + writeErrtoStdOut(i, err.errCode, err.errLocation, + err.paramError, errMsg->second); + } else { + cout << "Error #" << i << ": Invalid MOLE Error code [" + << err.errCode << "] - "<< "occurred inside:" + << err.errLocation; + MOLEerr_print_args(err.paramError); + } + i++; + } +} + +// This function produces a file with all errors logged in a stack. +// The file name is a string form by the concatanetion of the MOLE +// class type associated with the erros (i.e., grid, operator, or +// bcs) and a unique timestamp +void MOLEerr_dumpErrLog(stack& errorStack, + string logType){ + // create a string for the filename + string filename = logType + getDateTimeString(); + std::ofstream outFile(filename); + stack tempStack = errorStack; + + // Dump the logged errors to an output file + if (tempStack.empty()) { + outFile << "No errors logged." << endl; + return; + } + else { + int i = 1; + outFile << "========================================"<< endl; + outFile << "Backtracing all logged MOLE Errors :" << endl; + outFile << "========================================"<< endl; + while (!tempStack.empty()) { + MOLE_Errors err = tempStack.top(); + tempStack.pop(); + // if the grid has not been validated + if (err.errCode == MOLE_ERR_GRID_UNCHECKED) { + writeErrtoFile(outFile, i, MOLE_ERR_GRID_UNCHECKED, + err.errLocation, err.paramError, + MOLE_errors_messages[err.errCode]); + i++; + continue; + } + // For all other MOLE errors recorded, check for valid + // error code, in case the error log has been corrupted + auto errMsg = MOLE_errors_messages.find(err.errCode); + if (errMsg != MOLE_errors_messages.end()) { + writeErrtoFile(outFile, i, err.errCode, + err.errLocation, err.paramError, errMsg->second); + } + else { + writeErrtoFile(outFile, i, 999, err.errLocation, + to_string(err.errCode), "Unknown MOLE error code "); + } + i++; + } + } + outFile.close(); +} \ No newline at end of file diff --git a/cpp/src/sys/README.md b/cpp/src/sys/README.md new file mode 100644 index 00000000..31d1f0ae --- /dev/null +++ b/cpp/src/sys/README.md @@ -0,0 +1,28 @@ + + +# Subdirectory for MOLE 2.0 C++ grid examples + +Subdirectory and Pathname: **mole/cpp/src/sys/** + +## Purpose + +This subdirectory contains functional implementations of MOLE's +computational support like the error handling mechanisms. + +## List of Files in This Subdirectory (in alphabetical order) + ++ **MOLE_Erros.cpp**: a C++ implementation of MOLE's error handling +and tracking mechanisms ++ **README.md**: (this file) diff --git a/cpp/src/utils/README.md b/cpp/src/utils/README.md new file mode 100644 index 00000000..4b42cfb1 --- /dev/null +++ b/cpp/src/utils/README.md @@ -0,0 +1,29 @@ + + +# Subdirectory for MOLE 2.0 C++ Utilities and Wrappers + +Subdirectory and Pathname: **mole/cpp/src/utils** + +## Purpose + +This subdirectory contains functionality that either extends the MOLE +library's functionality through a third-party software library or +extends existing functionality from other software libraries and +customizes it to the MOLE library. + +## List of files in this subdirectory (alphabetical order) + ++ **utils.cpp**: a class implementation with MOLE basic utilities for +sparse matrices and meshes. diff --git a/cpp/src/utils/utils.cpp b/cpp/src/utils/utils.cpp new file mode 100644 index 00000000..5be03d95 --- /dev/null +++ b/cpp/src/utils/utils.cpp @@ -0,0 +1,340 @@ +/* +* SPDX-License-Identifier: GPL-3.0-or-later +* © 2008-2024 San Diego State University Research Foundation (SDSURF). +* See LICENSE file or https://www.gnu.org/licenses/gpl-3.0.html for details. +*/ + + +/* + * @file utils.cpp + * @brief Helpers for sparse operations and MATLAB/Octave analogs + * @date 2024/10/15 + * New Implementation for MOLE 2.0 Last Modified 2026/07/22 + * + * Sparse operations that repeatedly are needed, but not + * necessarily part of the Armadillo library. Some other MATLAB/Octave + * type functions are also here, like meshgrid. + */ + +#include "utils.h" +#include +#include + +#ifdef EIGEN +#include + +vec Utils::spsolve_eigen(const sp_mat &A, const vec &b) { + Eigen::SparseMatrix eigen_A(A.n_rows, A.n_cols); + std::vector> triplets; + Eigen::SparseLU, Eigen::COLAMDOrdering> solver; + + Eigen::VectorXd eigen_x(A.n_rows); + triplets.reserve(5 * A.n_rows); + + auto it = A.begin(); + while (it != A.end()) { + triplets.push_back(Eigen::Triplet(it.row(), it.col(), *it)); + ++it; + } + + eigen_A.setFromTriplets(triplets.begin(), triplets.end()); + triplets.clear(); + + auto b_ = conv_to>::from(b); + Eigen::Map eigen_b(b_.data(), b_.size()); + + solver.analyzePattern(eigen_A); + solver.factorize(eigen_A); + eigen_x = solver.solve(eigen_b); + + return vec(eigen_x.data(), eigen_x.size()); +} +#endif + +// ------------------------------------------------------------------ +// +// Error handling methods for MOLE Utils +// +// ------------------------------------------------------------------ + +// +// hasErrors checks if there are any errors logged in the Utils class +bool Utils::hasErrors() { + return !errs.empty(); +} + +// +// print_ErrorLog prints the error log to stdout +// +void Utils::print_ErrorLog() { + MOLEerr_print(errs); +} + +// +// write_ErrorLog writes the error log to a file with name starting +// with UtilsErrors, the full file name includes a timestamp. +// +void Utils::write_ErrorLog() { + MOLEerr_dumpErrLog(errs, "UtilsErrors"); // writes the error log to a file +} + +// +// read_ErrorLog reads the error on top of the stack +// +void Utils::read_ErrorLog(int ErrorCode, std::string &location, + std::string &arrayName) { + if (!errs.empty()) { + MOLE_Errors topError = errs.top(); + ErrorCode = topError.errCode; + location = topError.errLocation; + arrayName = topError.paramError; + errs.pop(); // Remove the top error after reading + } +} + +// +// spkron computes the sparse Kronecker tensor product of two +//sparse matrices A and B. This version uses Armadillo's sparse +// matrices and it is based on Octave's kron(A, B) operation. +// +sp_mat Utils::spkron(const sp_mat &A, const sp_mat &B) { + sp_mat::const_iterator itA = A.begin(); + sp_mat::const_iterator endA = A.end(); + sp_mat::const_iterator itB = B.begin(); + sp_mat::const_iterator endB = B.end(); + u32 j = 0; + + vec a = nonzeros(A); + vec b = nonzeros(B); + + umat locations(2, a.n_elem * b.n_elem); + vec values(a.n_elem * b.n_elem); + + while (itA != endA) { + while (itB != endB) { + locations(0, j) = itA.row() * B.n_rows + itB.row(); + locations(1, j) = itA.col() * B.n_cols + itB.col(); + values(j) = (*itA) * (*itB); + ++j; + ++itB; + } + + ++itA; + itB = B.begin(); + } + + sp_mat result(locations, values, A.n_rows * B.n_rows, A.n_cols * B.n_cols, + true); + + return result; +} + + +// +// spjoin_rows joins two sparse matrices by rows, returning a new +// sparse matrix. This version uses Armadillo's sparse matrices and +// it is based on Octave's [A B] operation, where A and B are sparse +// matrices. +// +sp_mat Utils::spjoin_rows(const sp_mat &A, const sp_mat &B) { + sp_mat::const_iterator itA = A.begin(); + sp_mat::const_iterator endA = A.end(); + sp_mat::const_iterator itB = B.begin(); + sp_mat::const_iterator endB = B.end(); + u32 j = 0; + + vec a = nonzeros(A); + vec b = nonzeros(B); + + umat locations(2, a.n_elem + b.n_elem); + vec values(a.n_elem + b.n_elem); + + while (itA != endA) { + locations(0, j) = itA.row(); + locations(1, j) = itA.col(); + values(j) = (*itA); + ++itA; + ++j; + } + + while (itB != endB) { + locations(0, j) = itB.row(); + locations(1, j) = itB.col() + A.n_cols; + values(j) = (*itB); + ++itB; + ++j; + } + + sp_mat result(locations, values, A.n_rows, A.n_cols + B.n_cols, true); + + return result; +} + +// +// spjoin_cols joins two sparse matrices by columns, returning a new +// sparse matrix. This version uses Armadillo's sparse matrices and +// it is based on Octave's [A B]^T operation, where A and B are +// sparse matrices. +// +sp_mat Utils::spjoin_cols(const sp_mat &A, const sp_mat &B) { + sp_mat::const_iterator itA = A.begin(); + sp_mat::const_iterator endA = A.end(); + sp_mat::const_iterator itB = B.begin(); + sp_mat::const_iterator endB = B.end(); + u32 j = 0; + + vec a = nonzeros(A); + vec b = nonzeros(B); + + umat locations(2, a.n_elem + b.n_elem); + vec values(a.n_elem + b.n_elem); + + while (itA != endA) { + locations(0, j) = itA.row(); + locations(1, j) = itA.col(); + values(j) = (*itA); + ++itA; + ++j; + } + + while (itB != endB) { + locations(0, j) = itB.row() + A.n_rows; + locations(1, j) = itB.col(); + values(j) = (*itB); + ++itB; + ++j; + } + + sp_mat result(locations, values, A.n_rows + B.n_rows, A.n_cols, true); + + return result; +} + +// +// meshgrid generates a 2D mesh from two 1D arrays representing the x +// and y coordinates (matches Octave meshgrid implemenation, which is +// incolumn-major order). This version uses Armadillo vecs and mats. +// It also checks for valid input sizes and reports errors. +// +void Utils::mesh2Dgrid(const vec &x, const vec &y, mat &X, mat &Y) { + size_t m = x.n_elem; + size_t n = y.n_elem; + bool valid = true; + + if (m <= 0){ // assert(m > 0) + std::string errmsg = "m = " + std::to_string(m); + MOLEerr_log(errs, MOLE_ERR_INVALID_ARRAY_SIZE, "Utils::meshgrid", + errmsg); + valid = false; + } + if (n <= 0) {// assert(n > 0); + std::string errmsg = "n = " + std::to_string(n); + MOLEerr_log(errs, MOLE_ERR_INVALID_ARRAY_SIZE, "Utils::meshgrid", + errmsg); + valid = false; + } + + if (valid){ + // Build X + vec t(n, fill::ones); + + X.zeros(n, m); + Y.zeros(n, m); + + for (size_t ii = 0; ii < m; ++ii) { + X.col(ii) = x(ii) * t; + t.ones(); + } + + // Build Y + for (size_t ii = 0; ii < m; ++ii) + Y.col(ii) = y; + } +} + +// +// meshgrid generates a 3D mesh from three 1D arrays representing the +// x, y, and z coordinates (matches Octave meshgrid implemenation, +// which is in column-major order). This version uses Armadillo vecs +// and mats. It also checks for valid input sizes and reports errors. +// +void Utils::mesh3Dgrid(const vec &x, const vec &y, const vec &z, + cube &X, cube &Y, cube &Z) { + size_t m = x.n_elem; + size_t n = y.n_elem; + size_t o = z.n_elem; + + bool valid = true; + if (m <= 0){ // assert(m > 0) + std::string errmsg = "m = " + std::to_string(m); + MOLEerr_log(errs, MOLE_ERR_INVALID_ARRAY_SIZE, "Utils::meshgrid", + errmsg); + valid = false; + } + if (n <= 0) {// assert(n > 0); + std::string errmsg = "n = " + std::to_string(n); + MOLEerr_log(errs, MOLE_ERR_INVALID_ARRAY_SIZE, "Utils::meshgrid", + errmsg); + valid = false; + } + if (o <= 0) {// assert(o > 0); + std::string errmsg = "o = " + std::to_string(o); + MOLEerr_log(errs, MOLE_ERR_INVALID_ARRAY_SIZE, "Utils::meshgrid", + errmsg); + valid = false; + } + + if (valid){ + // Temporary Holder of sheet of cube + mat sheet(m, n, fill::zeros); + + // Build X + vec t(n, fill::ones); + + X.zeros(m, n, o); + Y.zeros(m, n, o); + Z.zeros(m, n, o); + + // Sheet that repeats each slice + for (size_t ii = 0; ii < m; ++ii) { + sheet.row(ii) = x(ii) * t.t(); + t.ones(); + } + + for (size_t kk = 0; kk < o; ++kk) + X.slice(kk) = sheet; + + // Y Cube, repeats same sheet as well + for (size_t ii = 0; ii < m; ++ii) + sheet.row(ii) = y.t(); + + for (size_t kk = 0; kk < o; ++kk) + Y.slice(kk) = sheet; + + // Z cube goes by slices each with same value + for (size_t kk = 0; kk < o; ++kk) + Z.slice(kk).fill(z(kk)); + } +} + +// +// Trapezoidal rule (trapz) for 1D integration. This version uses +// Armadillo vecs and mats. It losgs and error when the two input +// vectors are not the same size. +// +double Utils::trapz(const vec &x, const vec &y) { + + if (x.n_elem == y.n_elem){ + double sum = 0.0; + for (uword i = 0; i < x.n_elem - 1; ++i) { + sum += (x(i+1) - x(i)) * (y(i) + y(i+1)); + } + return 0.5 * sum; + } else{ + std::string errmsg = "x._n_elem = " + std::to_string(x.n_elem); + errmsg += ", and y.n_elem = " + std::to_string(y.n_elem); + MOLEerr_log(errs, MOLE_ERR_INVALID_ARRAY_SIZE, "Utils::trapz", errmsg); + return std::nan(""); + } +} + diff --git a/cpp/tests/CMakeLists.txt b/cpp/tests/CMakeLists.txt new file mode 100644 index 00000000..08d204e1 --- /dev/null +++ b/cpp/tests/CMakeLists.txt @@ -0,0 +1,44 @@ +cmake_minimum_required(VERSION 3.14) + +# --------------------------------------------------------------- +# MOLE regression test suite. +# +# Each entry below becomes its own executable AND its own CTest +# test, so `ctest` reports pass/fail per source file (and a failure +# in one doesn't stop the others from running). All test binaries +# link against the same MOLE::mole target the examples use, so tests +# are exercising the exact same build the examples/users see. +# --------------------------------------------------------------- +set(MOLE_TEST_SOURCES + arrays/test_arrays.cpp + errors/test_errors.cpp + utils/test_utils.cpp + grids/test_grid1D.cpp + grids/test_grid2D.cpp + grids/test_grid3D.cpp + grids/test_makeGrid.cpp + grids/test_grid_builder.cpp + grids/test_debug_modes.cpp +) + +add_custom_target(mole_tests) + +foreach(_mole_test_src IN LISTS MOLE_TEST_SOURCES) + if(NOT EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/${_mole_test_src}") + message(FATAL_ERROR + "MOLE test source listed in tests/CMakeLists.txt is " + "missing: ${_mole_test_src}") + endif() + + get_filename_component(_mole_test_name ${_mole_test_src} NAME_WE) + set(_mole_test_bin "mole_${_mole_test_name}") + + add_executable(${_mole_test_bin} ${_mole_test_src}) + target_link_libraries(${_mole_test_bin} PRIVATE MOLE::mole) + target_include_directories(${_mole_test_bin} PRIVATE support) + target_compile_options(${_mole_test_bin} PRIVATE -Wall -Wextra) + + add_test(NAME ${_mole_test_name} COMMAND ${_mole_test_bin}) + + add_dependencies(mole_tests ${_mole_test_bin}) +endforeach() diff --git a/cpp/tests/README.md b/cpp/tests/README.md new file mode 100644 index 00000000..23ad0a50 --- /dev/null +++ b/cpp/tests/README.md @@ -0,0 +1,43 @@ + + +# Top subdirectory for MOLE 2.0 C++ files + +Subdirectory and Pathname: **mole/cpp/tests** + +## Purpose + +MOLE 2.0 C++ Top subdirectory for the test suite (regression testing) + +## MOLE C++ 2.0 Directory structure + +```text +mole/ +├── cpp/ +| │── cmake/ +| │── doc/ +| │── examples +| |── src/ +| |── tests +| | ├── arrays +| │ └── errors +| │ └── grids +| │ └── support +| │ └── utils +``` + +## Other MOLE 2.0 C++ Files in this directory + ++ CMakeLists.txt which builds tests for MOLE 2.0 C++ ++ README.md (this file) diff --git a/cpp/tests/arrays/README.md b/cpp/tests/arrays/README.md new file mode 100644 index 00000000..1922652b --- /dev/null +++ b/cpp/tests/arrays/README.md @@ -0,0 +1,26 @@ + + +# Subdirectory for MOLE 2.0 C++ Tests for the array classes + +Subdirectory and Pathname: **mole/cpp/tests/arrays/** + +## Purpose + +Subdirectory containing tests for the MOLE array classes + +## List of Files in This Subdirectory + ++ **test_arrays.cpp**: C++ tests of MOLE array classes ++ **README.md**: (this file) diff --git a/cpp/tests/arrays/test_arrays.cpp b/cpp/tests/arrays/test_arrays.cpp new file mode 100644 index 00000000..e1117c95 --- /dev/null +++ b/cpp/tests/arrays/test_arrays.cpp @@ -0,0 +1,241 @@ +// Regression tests for array1D, array2D, array3D, and numEqualArray. +#include "MOLE_arrays.h" +#include "mole_test.h" + +// --------------------------------------------------------------- +// Default construction / empty-object semantics +// --------------------------------------------------------------- + +TEST_CASE("array1D default construction is empty") { + array1D a; + CHECK(a.data_.is_empty()); + CHECK(a.data_.n_elem == 0); +} + +TEST_CASE("array2D default construction is empty") { + array2D a; + CHECK(a.data_.is_empty()); + CHECK(a.data_.n_rows == 0); + CHECK(a.data_.n_cols == 0); +} + +TEST_CASE("array3D default construction is empty") { + array3D a; + CHECK(a.data_.is_empty()); + CHECK(a.data_.n_rows == 0); + CHECK(a.data_.n_cols == 0); + CHECK(a.data_.n_slices == 0); +} + +// --------------------------------------------------------------- +// Parameterized constructors +// --------------------------------------------------------------- + +TEST_CASE("array1D(n, fill) allocates and fills") { + array1D a(5, 3.0); + REQUIRE(a.data_.n_elem == 5); + for (arma::uword i = 0; i < a.data_.n_elem; ++i) { + CHECK(a.data_(i) == 3.0); + } +} + +TEST_CASE("array1D default fill value is 0.0") { + array1D a(4); + for (arma::uword i = 0; i < a.data_.n_elem; ++i) { + CHECK(a.data_(i) == 0.0); + } +} + +TEST_CASE("array2D(rows, cols, fill) allocates and fills") { + array2D a(3, 4, 7.0); + REQUIRE(a.data_.n_rows == 3); + REQUIRE(a.data_.n_cols == 4); + CHECK(a.data_(0, 0) == 7.0); + CHECK(a.data_(2, 3) == 7.0); +} + +TEST_CASE("array3D(d1, d2, d3, fill) allocates and fills") { + array3D a(2, 3, 4, 9.0); + REQUIRE(a.data_.n_rows == 2); + REQUIRE(a.data_.n_cols == 3); + REQUIRE(a.data_.n_slices == 4); + CHECK(a.data_(1, 2, 3) == 9.0); +} + +// --------------------------------------------------------------- +// valid_index / valid_indeces +// --------------------------------------------------------------- + +TEST_CASE("array1D::valid_index respects bounds") { + array1D a(3, 0.0); + CHECK(a.valid_index(0)); + CHECK(a.valid_index(2)); + CHECK(!a.valid_index(3)); // one-past-the-end is invalid +} + +TEST_CASE("array2D::valid_indeces respects bounds") { + array2D a(2, 3, 0.0); + CHECK(a.valid_indeces(0, 0)); + CHECK(a.valid_indeces(1, 2)); + CHECK(!a.valid_indeces(2, 0)); + CHECK(!a.valid_indeces(0, 3)); +} + +TEST_CASE("array3D::valid_indeces respects bounds") { + array3D a(2, 2, 2, 0.0); + CHECK(a.valid_indeces(1, 1, 1)); + CHECK(!a.valid_indeces(2, 0, 0)); + CHECK(!a.valid_indeces(0, 2, 0)); + CHECK(!a.valid_indeces(0, 0, 2)); +} + +// --------------------------------------------------------------- +// resize: preserves-vs-overwrites semantics as actually implemented +// (resize reallocates + fills with fillVal; it does NOT preserve +// prior contents the way arma::mat::resize does. This test locks in +// the class's actual documented behavior rather than assuming +// Armadillo's semantics carry over.) +// --------------------------------------------------------------- + +TEST_CASE("array1D::resize reallocates and fills with fillVal") { + array1D a(3, 1.0); + a.resize(5, 2.0); + REQUIRE(a.data_.n_elem == 5); + for (arma::uword i = 0; i < a.data_.n_elem; ++i) { + CHECK(a.data_(i) == 2.0); + } + CHECK(!a.hasArrayErrors()); +} + +TEST_CASE("array2D::resize reallocates and fills with fillVal") { + array2D a(2, 2, 1.0); + a.resize(3, 3, 5.0); + REQUIRE(a.data_.n_rows == 3); + REQUIRE(a.data_.n_cols == 3); + CHECK(a.data_(2, 2) == 5.0); + CHECK(!a.hasArrayErrors()); +} + +TEST_CASE("array3D::resize reallocates and fills with fillVal") { + array3D a(2, 2, 2, 1.0); + a.resize(3, 3, 3, 6.0); + REQUIRE(a.data_.n_rows == 3); + REQUIRE(a.data_.n_cols == 3); + REQUIRE(a.data_.n_slices == 3); + CHECK(a.data_(2, 2, 2) == 6.0); + CHECK(!a.hasArrayErrors()); +} + +TEST_CASE("array1D::resize to 0 produces an empty array, not an error") { + array1D a(4, 1.0); + a.resize(0); + CHECK(a.data_.is_empty()); + CHECK(a.data_.n_elem == 0); + CHECK(!a.hasArrayErrors()); +} + +// --------------------------------------------------------------- +// operator== +// +// operator== is a proper value comparison (size + contents), and +// correctly reports equal for two independently-constructed arrays +// that happen to hold the same values. +// --------------------------------------------------------------- + +TEST_CASE("array1D::operator== compares by value, not identity") { + array1D a(3, 2.5), b(3, 2.5); + CHECK(a.data_.memptr() != b.data_.memptr()); // genuinely distinct + CHECK(a == b); // but equal by value +} + +TEST_CASE("array1D::operator== detects differing size") { + array1D a(3, 1.0), b(4, 1.0); + CHECK(!(a == b)); +} + +TEST_CASE("array1D::operator== detects differing content") { + array1D a(3, 1.0), b(3, 2.0); + CHECK(!(a == b)); +} + +TEST_CASE("array2D::operator== compares by value") { + array2D a(2, 2, 4.0), b(2, 2, 4.0); + CHECK(a == b); + array2D c(2, 2, 4.0), d(2, 3, 4.0); + CHECK(!(c == d)); +} + +TEST_CASE("array3D::operator== compares by value") { + array3D a(2, 2, 2, 4.0), b(2, 2, 2, 4.0); + CHECK(a == b); + array3D c(2, 2, 2, 4.0), d(2, 2, 3, 4.0); + CHECK(!(c == d)); +} + +// --------------------------------------------------------------- +// operator!= -- KNOWN DIVERGENCE FROM operator== +// +// As currently implemented, operator!= is an identity check (does +// this array live at a different memory address / have a different +// shape?), NOT the logical negation of operator==. This means two +// independently-constructed, value-equal arrays report BOTH +// (a == b) == true AND (a != b) == true simultaneously, which is +// a genuine logical inconsistency between the two operators. +// +// This test intentionally documents that CURRENT behavior (so a +// future accidental "fix" that changes it doesn't silently pass +// unnoticed) while flagging it clearly as inconsistent with ==. +// See the accompanying regression-test report for a suggested fix +// (operator!= should just be `return !(*this == other);`). +// --------------------------------------------------------------- + +TEST_CASE("KNOWN ISSUE: array1D::operator!= is pointer-identity, " + "not the logical negation of operator==") { + array1D a(3, 2.5), b(3, 2.5); + REQUIRE(a == b); // equal by value ... + CHECK_MSG(a != b, // ... yet also reported "not equal" + "operator!= currently returns true here because it only " + "compares memptr()/shape, not values. If operator!= is ever " + "fixed to be !(a==b), this CHECK will (correctly) start " + "failing and should be updated to CHECK(!(a != b))."); +} + +// --------------------------------------------------------------- +// numEqualArray +// --------------------------------------------------------------- + +TEST_CASE("numEqualArray: identical array1D compares equal") { + array1D a(4, 1.0), b(4, 1.0); + CHECK(numEqualArray(a, b, 4.0)); +} + +TEST_CASE("numEqualArray: differing sizes compare unequal") { + array1D a(4, 1.0), b(5, 1.0); + CHECK(!numEqualArray(a, b, 4.0)); +} + +TEST_CASE("numEqualArray: tiny floating point noise within tolerance") { + array1D a(3, 1.0); + // epsilon for double is ~2.22e-16, so 4*epsilon ~8.88e-16. + // 5e-17 is safely below that; 1e-15 (an earlier version of this + // test used that) is actually ABOVE 4*epsilon and correctly + // fails -- don't widen this without re-checking the tolerance. + array1D b(3, 1.0 + 5e-17); + CHECK(numEqualArray(a, b, 4.0)); +} + +TEST_CASE("numEqualArray: difference beyond tolerance is rejected") { + array1D a(3, 1.0); + array1D b(3, 1.01); // far beyond 4*epsilon + CHECK(!numEqualArray(a, b, 4.0)); +} + +TEST_CASE("numEqualArray works for array2D and array3D") { + array2D a2(2, 2, 3.0), b2(2, 2, 3.0); + CHECK(numEqualArray(a2, b2, 4.0)); + + array3D a3(2, 2, 2, 3.0), b3(2, 2, 2, 3.0); + CHECK(numEqualArray(a3, b3, 4.0)); +} + +MOLE_TEST_MAIN() diff --git a/cpp/tests/errors/README.md b/cpp/tests/errors/README.md new file mode 100644 index 00000000..0fde6ace --- /dev/null +++ b/cpp/tests/errors/README.md @@ -0,0 +1,28 @@ + + +# Subdirectory for MOLE 2.0 C++ Tests for the array classes + +Subdirectory and Pathname: **mole/cpp/tests/errors/** + +## Purpose + +Subdirectory containing tests for the MOLE error handling and +reporting services + +## List of Files in This Subdirectory + ++ **test_errors.cpp**: C++ tests of the MOLE error logging +and reporting mechanisms works ++ **README.md**: (this file) diff --git a/cpp/tests/errors/test_errors.cpp b/cpp/tests/errors/test_errors.cpp new file mode 100644 index 00000000..5107a6f2 --- /dev/null +++ b/cpp/tests/errors/test_errors.cpp @@ -0,0 +1,90 @@ +// Regression tests for MOLE_Errors.cpp (the free-function error +// stack API used by every MOLE class). +#include "MOLE_errors.h" +#include "mole_test.h" +#include + +TEST_CASE("MOLEerr_log pushes an entry, MOLEerr_haserrors sees it") { + std::stack errs; + CHECK(!MOLEerr_haserrors(errs)); + MOLEerr_log(errs, MOLE_ERR_INVALID_GRID_SPACING, "test", "dx=-1"); + CHECK(MOLEerr_haserrors(errs)); + REQUIRE(!errs.empty()); + CHECK(errs.top().errCode == MOLE_ERR_INVALID_GRID_SPACING); + CHECK(errs.top().errLocation == "test"); +} + +TEST_CASE("MOLEerr_contains finds a logged code and ignores others") { + std::stack errs; + MOLEerr_log(errs, MOLE_ERR_INVALID_GRID_SPACING, "a", ""); + MOLEerr_log(errs, MOLE_ERR_GRID_NODAL_SZ_MISMATCH, "b", ""); + CHECK(MOLEerr_contains(errs, MOLE_ERR_INVALID_GRID_SPACING)); + CHECK(MOLEerr_contains(errs, MOLE_ERR_GRID_NODAL_SZ_MISMATCH)); + CHECK(!MOLEerr_contains(errs, MOLE_ERR_INVALID_GRID_TOPOLOGY)); +} + +TEST_CASE("MOLEerr_contains on an empty stack returns false") { + std::stack errs; + CHECK(!MOLEerr_contains(errs, MOLE_ERR_GRID_UNCHECKED)); +} + +TEST_CASE("MOLEerr_remove removes only the targeted code") { + std::stack errs; + MOLEerr_log(errs, MOLE_ERR_INVALID_GRID_SPACING, "a", ""); + MOLEerr_log(errs, MOLE_ERR_GRID_NODAL_SZ_MISMATCH, "b", ""); + MOLEerr_log(errs, MOLE_ERR_INVALID_GRID_SPACING, "c", ""); + + MOLEerr_remove(errs, MOLE_ERR_INVALID_GRID_SPACING); + + CHECK(!MOLEerr_contains(errs, MOLE_ERR_INVALID_GRID_SPACING)); + CHECK(MOLEerr_contains(errs, MOLE_ERR_GRID_NODAL_SZ_MISMATCH)); +} + +TEST_CASE("MOLEerr_print reports an unrecognized code exactly once") { + std::stack errs; + // 99999 is not a key in MOLE_errors_messages. + MOLEerr_log(errs, 99999, "bogus_location", "param"); + + std::ostringstream captured; + std::streambuf* old_cout = std::cout.rdbuf(captured.rdbuf()); + MOLEerr_print(errs); + std::cout.rdbuf(old_cout); + + const std::string out = captured.str(); + size_t count = 0, pos = 0; + while ((pos = out.find("bogus_location", pos)) != std::string::npos) { + ++count; + pos += 1; + } + CHECK_MSG(count == 1, + "expected 'bogus_location' to appear exactly once in " + "MOLEerr_print output, got " << count + << ". Full output:\n" << out); +} + +TEST_CASE("MOLEerr_print reports a recognized code with its message") { + std::stack errs; + MOLEerr_log(errs, MOLE_ERR_INVALID_GRID_SPACING, "test_loc", "dx=0"); + + std::ostringstream captured; + std::streambuf* old_cout = std::cout.rdbuf(captured.rdbuf()); + MOLEerr_print(errs); + std::cout.rdbuf(old_cout); + + const std::string out = captured.str(); + CHECK(out.find("test_loc") != std::string::npos); + CHECK(out.find( + MOLE_errors_messages[MOLE_ERR_INVALID_GRID_SPACING]) + != std::string::npos); +} + +TEST_CASE("MOLEerr_init clears any existing entries") { + std::stack errs; + MOLEerr_log(errs, MOLE_ERR_INVALID_GRID_SPACING, "a", ""); + MOLEerr_log(errs, MOLE_ERR_GRID_NODAL_SZ_MISMATCH, "b", ""); + REQUIRE(MOLEerr_haserrors(errs)); + MOLEerr_init(errs); + CHECK(!MOLEerr_haserrors(errs)); +} + +MOLE_TEST_MAIN() diff --git a/cpp/tests/grids/README.md b/cpp/tests/grids/README.md new file mode 100644 index 00000000..1db6a13b --- /dev/null +++ b/cpp/tests/grids/README.md @@ -0,0 +1,29 @@ + + +# Subdirectory for MOLE 2.0 C++ Tests for the array classes + +Subdirectory and Pathname: **mole/cpp/tests/grids/** + +## Purpose + +Subdirectory containing tests for the MOLE grids + +## List of Files in This Subdirectory + ++ **test_grid1D.cpp**: C++ tests of the 1D Grid Class ++ **test_grid2D.cpp**: C++ tests of the 2D Grid Class ++ **test_grid3D.cpp**: C++ tests of the 3D Grid Class ++ **test_makeGrid.cpp**: C++ tests of the makeGrid functionalities ++ **README.md**: (this file) diff --git a/cpp/tests/grids/test_debug_modes.cpp b/cpp/tests/grids/test_debug_modes.cpp new file mode 100644 index 00000000..907f7222 --- /dev/null +++ b/cpp/tests/grids/test_debug_modes.cpp @@ -0,0 +1,310 @@ +/* + * SPDX-License-Identifier: GPL-3.0-or-later + * Copyright (c) 2008-2024 San Diego State University Research Foundation + * (SDSURF). + * See LICENSE file or https://www.gnu.org/licenses/gpl-3.0.html for details. + */ + +// Regression tests for the MOLE debug modes. +// +// The modes are declared in MOLE_errors.h and applied by +// gridBase::applyDebugMode, which the grid constructors and +// gridBuilder call. Two properties are under test: +// +// 1. a mode changes only what is written to standard output; the +// grid's error log is identical afterwards in every mode +// 2. a grid that passed validation ignores the mode entirely, +// even when it carries errors merged from upstream +// +// DEBUG_AND_ABORT_MD is only exercised on grids that validated, +// since a real abort would end the test binary. If property 2 +// regresses, these cases abort and ctest reports the failure. +#include "MOLE_grids.h" +#include "grid_builder.h" +#include "mole_test.h" + +#include +#include +#include +#include + +// capture redirects standard output for the duration of fn and +// returns whatever was written. The grid error reporting goes to +// cout; mole_test.h writes failures to cerr, so assertion output is +// not swallowed. +template +static std::string capture(F fn) { + std::ostringstream buf; + std::streambuf* old = std::cout.rdbuf(buf.rdbuf()); + fn(); + std::cout.rdbuf(old); + return buf.str(); +} + +// badParams1D returns a gridParams1D that cannot validate: 'z' is +// not a MOLE topology. +static gridParams1D badParams1D() { + gridParams1D p; + p.topology = 'z'; + p.m = 5; + p.dx = 1.0; + return p; +} + +// goodParams1D returns a gridParams1D that validates. +static gridParams1D goodParams1D() { + gridParams1D p; + p.topology = 'u'; + p.m = 5; + p.dx = 0.5; + return p; +} + +// --------------------------------------------------------------- +// What each mode writes +// --------------------------------------------------------------- + +TEST_CASE("DEBUG_DEFAULT_MD writes nothing for an invalid grid") { + std::string out = capture([]{ + grid1D g(badParams1D(), DEBUG_DEFAULT_MD); + CHECK(!g.isValidatedGrid()); + }); + CHECK_MSG(out.empty(), + "expected no output in the default mode, got: " << out); +} + +TEST_CASE("DEBUG_REPORTS_STDOUT_MD writes the log for an invalid " + "grid and returns control") { + bool reached_next_line = false; + std::string out = capture([&]{ + grid1D g(badParams1D(), DEBUG_REPORTS_STDOUT_MD); + reached_next_line = true; + CHECK(!g.isValidatedGrid()); + }); + CHECK(reached_next_line); + CHECK(out.find("MOLE Error code") != std::string::npos); +} + +TEST_CASE("an unrecognized debug mode falls back to reporting") { + std::string out = capture([]{ + grid1D g(badParams1D(), 99); + CHECK(!g.isValidatedGrid()); + }); + CHECK(out.find("Unrecognized MOLE debug mode") + != std::string::npos); + CHECK(out.find("MOLE Error code") != std::string::npos); +} + +// --------------------------------------------------------------- +// A mode must not consume the error log +// +// Reporting has to be non-destructive. The validation flag lives in +// the same stack as the errors (MOLE_ERR_GRID_UNCHECKED), so a mode +// that drained the stack would leave an invalid grid claiming to be +// validated. It would also break the promise that a user can still +// print or write the log after the library has reported it. +// --------------------------------------------------------------- + +TEST_CASE("reporting leaves the same error log behind as the " + "default mode") { + grid1D quiet(badParams1D(), DEBUG_DEFAULT_MD); + + std::string reported = capture([]{ + grid1D loud(badParams1D(), DEBUG_REPORTS_STDOUT_MD); + (void)loud; + }); + + grid1D loud(badParams1D(), DEBUG_DEFAULT_MD); + + std::string quiet_log = capture([&]{ quiet.print_ErrorLog(); }); + std::string loud_log = capture([&]{ loud.print_ErrorLog(); }); + + CHECK_MSG(quiet_log == loud_log, + "the two modes produced different error logs"); + CHECK_MSG(reported == loud_log, + "what the mode reported differs from what the grid kept"); +} + +TEST_CASE("a grid can still be asked for its log after the mode " + "already reported it") { + std::string first; + grid1D g(badParams1D(), DEBUG_DEFAULT_MD); + first = capture([&]{ g.print_ErrorLog(); }); + std::string second = capture([&]{ g.print_ErrorLog(); }); + CHECK(!first.empty()); + CHECK_MSG(first == second, + "print_ErrorLog is not repeatable"); + CHECK(!g.isValidatedGrid()); +} + +// --------------------------------------------------------------- +// A validated grid ignores the mode +// +// These are the regression tests for using isValidatedGrid() rather +// than hasGridErrors() as the trigger. A freshly built grid always +// has MOLE_ERR_GRID_UNCHECKED on its stack until validation clears +// it, and mergeErrors folds upstream errors into the same stack, so +// hasGridErrors() is true for grids that are perfectly usable. +// Under DEBUG_AND_ABORT_MD the wrong trigger ends the process. +// --------------------------------------------------------------- + +TEST_CASE("a valid grid ignores DEBUG_AND_ABORT_MD") { + std::string out = capture([]{ + grid1D g(goodParams1D(), DEBUG_AND_ABORT_MD); + CHECK(g.isValidatedGrid()); + }); + CHECK_MSG(out.empty(), + "a validated grid should produce no output, got: " << out); +} + +TEST_CASE("a valid grid carrying upstream errors ignores " + "DEBUG_AND_ABORT_MD") { + std::stack inerrs; + MOLEerr_init(inerrs); + MOLEerr_log(inerrs, MOLE_ERR_INVALID_INPUT_TYPE, "upstream", ""); + + std::string out = capture([&]{ + grid1D g(goodParams1D(), inerrs, DEBUG_AND_ABORT_MD); + CHECK(g.isValidatedGrid()); + // the upstream error is still on the grid's stack + CHECK(g.hasGridErrors()); + }); + CHECK(out.empty()); +} + +TEST_CASE("the inerrs constructor applies the mode to an invalid " + "grid") { + std::stack inerrs; + MOLEerr_init(inerrs); + MOLEerr_log(inerrs, MOLE_ERR_INVALID_INPUT_TYPE, "upstream", ""); + + std::string out = capture([&]{ + grid1D g(badParams1D(), inerrs, DEBUG_REPORTS_STDOUT_MD); + CHECK(!g.isValidatedGrid()); + }); + CHECK(out.find("MOLE Error code") != std::string::npos); + // the upstream error travelled into the report + CHECK(out.find("upstream") != std::string::npos); +} + +// --------------------------------------------------------------- +// 2D and 3D take the same path +// --------------------------------------------------------------- + +TEST_CASE("grid2D honours the debug modes") { + gridParams2D bad; + bad.topology = 'z'; + bad.m = 3; bad.n = 3; bad.dx = 1.0; bad.dy = 1.0; + + CHECK(capture([&]{ grid2D g(bad, DEBUG_DEFAULT_MD); }).empty()); + CHECK(!capture([&]{ + grid2D g(bad, DEBUG_REPORTS_STDOUT_MD); + }).empty()); + + gridParams2D good; + good.topology = 'u'; + good.m = 3; good.n = 3; good.dx = 1.0; good.dy = 1.0; + CHECK(capture([&]{ + grid2D g(good, DEBUG_AND_ABORT_MD); + CHECK(g.isValidatedGrid()); + }).empty()); +} + +TEST_CASE("grid3D honours the debug modes") { + gridParams3D bad; + bad.topology = 'z'; + bad.m = 2; bad.n = 3; bad.o = 2; + bad.dx = 1.0; bad.dy = 1.0; bad.dz = 1.0; + + CHECK(capture([&]{ grid3D g(bad, DEBUG_DEFAULT_MD); }).empty()); + CHECK(!capture([&]{ + grid3D g(bad, DEBUG_REPORTS_STDOUT_MD); + }).empty()); + + gridParams3D good; + good.topology = 'u'; + good.m = 2; good.n = 3; good.o = 2; + good.dx = 1.0; good.dy = 1.0; good.dz = 1.0; + CHECK(capture([&]{ + grid3D g(good, DEBUG_AND_ABORT_MD); + CHECK(g.isValidatedGrid()); + }).empty()); +} + +// --------------------------------------------------------------- +// The gridBuilder debug attribute +// --------------------------------------------------------------- + +TEST_CASE("gridBuilder without a debug attribute stays quiet on a " + "failed build") { + std::string out = capture([]{ + gridVar g = gridBuilder("dim", 1, "m", 5, "dx", 0.2, + "topology", 'z'); + CHECK(std::holds_alternative(g)); + }); + CHECK(out.empty()); +} + +TEST_CASE("gridBuilder reports a failed build when debug is " + "DEBUG_REPORTS_STDOUT_MD") { + std::string out = capture([]{ + gridVar g = gridBuilder("debug", DEBUG_REPORTS_STDOUT_MD, + "dim", 1, "m", 5, "dx", 0.2, + "topology", 'z'); + CHECK(std::holds_alternative(g)); + }); + CHECK(out.find("MOLE Error code") != std::string::npos); +} + +TEST_CASE("gridBuilder stays quiet when the grid builds, whatever " + "the debug mode") { + std::string out = capture([]{ + gridVar g = gridBuilder("debug", DEBUG_AND_ABORT_MD, + "dim", 1, "m", 5, "dx", 0.2, + "topology", 'u'); + REQUIRE(std::holds_alternative(g)); + CHECK(std::get(g).isValidatedGrid()); + }); + CHECK_MSG(out.empty(), + "a grid that built should produce no output, got: " << out); +} + +TEST_CASE("gridBuilder keeps the full log in every mode") { + gridVar quiet = gridBuilder("dim", 1, "m", 5, "dx", 0.2, + "topology", 'z'); + + std::string reported = capture([]{ + gridVar loud = gridBuilder("debug", DEBUG_REPORTS_STDOUT_MD, + "dim", 1, "m", 5, "dx", 0.2, + "topology", 'z'); + (void)loud; + }); + + std::string quiet_log = capture([&]{ + std::visit([](auto&& g){ g.print_ErrorLog(); }, quiet); + }); + + CHECK(!quiet_log.empty()); + CHECK_MSG(reported == quiet_log, + "gridBuilder reported something other than the log it " + "handed back"); +} + +TEST_CASE("NOTE: a debug pair placed after an unknown attribute is " + "never read") { + // Parsing stops at the first name that is not a grid attribute, + // because the type of the value following it is unknown. The + // debug pair below is never reached, so the build falls back to + // DEBUG_DEFAULT_MD and reports nothing. This is the reason the + // debug pair has to come first. + std::string out = capture([]{ + gridVar g = gridBuilder("dim", 1, "m", 5, "spacing", 0.2, + "debug", DEBUG_REPORTS_STDOUT_MD, + "topology", 'u'); + CHECK(std::holds_alternative(g)); + }); + CHECK_MSG(out.empty(), + "a trailing debug pair was read; the parser changed"); +} + +MOLE_TEST_MAIN() diff --git a/cpp/tests/grids/test_grid1D.cpp b/cpp/tests/grids/test_grid1D.cpp new file mode 100644 index 00000000..1795ac03 --- /dev/null +++ b/cpp/tests/grids/test_grid1D.cpp @@ -0,0 +1,134 @@ +// Regression tests for grid1D. +#include "MOLE_grids.h" +#include "mole_test.h" + +TEST_CASE("grid1D uniform: valid parameters build a valid grid") { + gridParams1D p; + p.topology = 'u'; + p.m = 5; + p.dx = 0.5; + grid1D g(p); + CHECK(g.isValidatedGrid()); + CHECK(g.grid.nodes_X.data_.n_elem == p.m + 1); + CHECK(g.grid.centers_X.data_.n_elem == p.m + 2); +} + +TEST_CASE("grid1D uniform: nodal coordinates have correct values") { + gridParams1D p; + p.topology = 'u'; + p.m = 4; + p.dx = 2.0; + grid1D g(p); + REQUIRE(g.isValidatedGrid()); + for (size_t i = 0; i <= p.m; ++i) { + CHECK_MSG(g.grid.nodes_X.data_(i) == static_cast(i) * p.dx, + "nodes_X(" << i << ") = " << g.grid.nodes_X.data_(i) + << ", expected " << (i * p.dx)); + } +} + +TEST_CASE("grid1D: invalid (non-positive) spacing is rejected") { + gridParams1D p; + p.topology = 'u'; + p.m = 5; + p.dx = -1.0; + grid1D g(p); + CHECK(!g.isValidatedGrid()); +} + +TEST_CASE("grid1D: invalid topology character is rejected") { + gridParams1D p; + p.topology = 'z'; + p.m = 5; + p.dx = 1.0; + grid1D g(p); + CHECK(!g.isValidatedGrid()); +} + +TEST_CASE("grid1D: curvilinear topology is fundamentally invalid in 1D") { + gridParams1D p; + p.topology = 'c'; + p.m = 5; + grid1D g(p); + CHECK(!g.isValidatedGrid()); +} + +TEST_CASE("grid1D: nonuniform topology without user nodes is rejected") { + gridParams1D p; + p.topology = 'n'; + p.m = 5; + grid1D g(p); + CHECK(!g.isValidatedGrid()); +} + +TEST_CASE("grid1D: Faces_X() aliases nodes_X (same object, not a copy)") { + gridParams1D p; + p.topology = 'u'; + p.m = 3; + p.dx = 1.0; + grid1D g(p); + REQUIRE(g.isValidatedGrid()); + CHECK(&g.grid.Faces_X() == &g.grid.nodes_X); + CHECK(g.grid.Faces_X() == g.grid.nodes_X); +} + +TEST_CASE("grid1D: Faces_X() tracks nodes_X correctly across struct " + "copies (does not stay bound to the original)") { + gridParams1D p; + p.topology = 'u'; + p.m = 3; + p.dx = 1.0; + grid1D g(p); + REQUIRE(g.isValidatedGrid()); + + gridParams1D copy = g.grid; // struct copy + copy.nodes_X = array1D(3, 99.0); + + // The copy's Faces_X() must reflect the COPY's own nodes_X, not + // the original grid's. This is the exact hazard a true C++ + // reference member would introduce (see prior discussion) -- + // Faces_X() is implemented as a function precisely to avoid it. + CHECK(copy.Faces_X().data_(0) == 99.0); + CHECK(g.grid.nodes_X.data_(0) != 99.0); +} + +TEST_CASE("grid1D: user-supplied, correctly-valued nodes_X validates") { + size_t m = 4; + double dx = 0.5; + array1D nx(m + 1, 0.0); + for (size_t i = 0; i <= m; ++i) nx.data_(i) = i * dx; + + gridParams1D p; + p.topology = 'u'; + p.m = m; + p.dx = dx; + p.nodes_X = nx; + grid1D g(p); + CHECK(g.isValidatedGrid()); +} + +TEST_CASE("grid1D: user-supplied nodes_X with wrong size is rejected") { + gridParams1D p; + p.topology = 'u'; + p.m = 4; + p.dx = 0.5; + p.nodes_X = array1D(3, 0.0); // wrong size: should be m+1 = 5 + grid1D g(p); + CHECK(!g.isValidatedGrid()); +} + +TEST_CASE("grid1D: constructor with inbound errors merges them in") { + std::stack inerrs; + MOLEerr_log(inerrs, MOLE_ERR_INVALID_INPUT_TYPE, "caller", "note"); + + gridParams1D p; + p.topology = 'u'; + p.m = 4; + p.dx = 1.0; + grid1D g(p, inerrs); + + g.print_ErrorLog(); // not asserted on; just confirm it doesn't crash + CHECK(g.hasGridErrors()); +} + +MOLE_TEST_MAIN() diff --git a/cpp/tests/grids/test_grid2D.cpp b/cpp/tests/grids/test_grid2D.cpp new file mode 100644 index 00000000..25238c07 --- /dev/null +++ b/cpp/tests/grids/test_grid2D.cpp @@ -0,0 +1,145 @@ +// Regression tests for grid2D. +// +// The asymmetric-size tests (m != n) here specifically guard against +// the axis-order/transposition bug found earlier in this codebase's +// history, where Utils::mesh2Dgrid's (len(y), len(x)) convention +// didn't match grid2D's (m, n) convention -- a bug invisible under +// square (m == n) grids. Do not "simplify" these to square grids. +#include "MOLE_grids.h" +#include "mole_test.h" + +TEST_CASE("grid2D uniform: valid asymmetric grid validates") { + gridParams2D p; + p.topology = 'u'; + p.m = 3; p.n = 5; + p.dx = 1.0; p.dy = 1.0; + grid2D g(p); + CHECK(g.isValidatedGrid()); +} + +TEST_CASE("grid2D uniform: layer shapes follow the (m,n) convention, " + "not a transposed one (regression guard)") { + gridParams2D p; + p.topology = 'u'; + p.m = 3; p.n = 5; // deliberately asymmetric + p.dx = 1.0; p.dy = 1.0; + grid2D g(p); + REQUIRE(g.isValidatedGrid()); + + CHECK(g.grid.nodes_X.data_.n_rows == p.m + 1); + CHECK(g.grid.nodes_X.data_.n_cols == p.n + 1); + CHECK(g.grid.nodes_Y.data_.n_rows == p.m + 1); + CHECK(g.grid.nodes_Y.data_.n_cols == p.n + 1); + + CHECK(g.grid.centers_X.data_.n_rows == p.m + 2); + CHECK(g.grid.centers_X.data_.n_cols == p.n + 2); + + CHECK(g.grid.faces_u_X.data_.n_rows == p.m + 1); + CHECK(g.grid.faces_u_X.data_.n_cols == p.n); + + CHECK(g.grid.faces_v_X.data_.n_rows == p.m); + CHECK(g.grid.faces_v_X.data_.n_cols == p.n + 1); +} + +TEST_CASE("grid2D uniform: nodal coordinate values are correct " + "(not just shapes)") { + gridParams2D p; + p.topology = 'u'; + p.m = 3; p.n = 4; + p.dx = 1.5; p.dy = 0.5; + grid2D g(p); + REQUIRE(g.isValidatedGrid()); + + for (size_t i = 0; i <= p.m; ++i) { + for (size_t j = 0; j <= p.n; ++j) { + CHECK(g.grid.nodes_X.data_(i, j) == static_cast(i) * p.dx); + CHECK(g.grid.nodes_Y.data_(i, j) == static_cast(j) * p.dy); + } + } +} + +TEST_CASE("grid2D uniform: user-supplied, correctly-shaped AND " + "correctly-valued nodes round-trip successfully") { + size_t m = 3, n = 5; + double dx = 1.0, dy = 1.0; + array2D nx(m + 1, n + 1, 0.0), ny(m + 1, n + 1, 0.0); + for (size_t i = 0; i <= m; ++i) + for (size_t j = 0; j <= n; ++j) { + nx.data_(i, j) = i * dx; + ny.data_(i, j) = j * dy; + } + + gridParams2D p; + p.topology = 'u'; + p.m = m; p.n = n; p.dx = dx; p.dy = dy; + p.nodes_X = nx; p.nodes_Y = ny; + grid2D g(p); + CHECK(g.isValidatedGrid()); +} + +TEST_CASE("grid2D uniform: user-supplied nodes with wrong shape " + "are rejected (not silently transposed)") { + gridParams2D p; + p.topology = 'u'; + p.m = 3; p.n = 5; + p.dx = 1.0; p.dy = 1.0; + // Deliberately transposed shape: (n+1) x (m+1) instead of (m+1) x (n+1) + p.nodes_X = array2D(p.n + 1, p.m + 1, 0.0); + p.nodes_Y = array2D(p.n + 1, p.m + 1, 0.0); + grid2D g(p); + CHECK(!g.isValidatedGrid()); + CHECK(g.hasGridErrors()); +} + +TEST_CASE("grid2D uniform: user-supplied nodes with correct shape but " + "wrong values are rejected") { + gridParams2D p; + p.topology = 'u'; + p.m = 3; p.n = 5; + p.dx = 1.0; p.dy = 1.0; + p.nodes_X = array2D(p.m + 1, p.n + 1, 0.0); // all zeros: wrong values + p.nodes_Y = array2D(p.m + 1, p.n + 1, 0.0); + grid2D g(p); + CHECK(!g.isValidatedGrid()); +} + +TEST_CASE("grid2D: invalid topology character is rejected") { + gridParams2D p; + p.topology = 'x'; + p.m = 3; p.n = 3; + p.dx = 1.0; p.dy = 1.0; + grid2D g(p); + CHECK(!g.isValidatedGrid()); +} + +TEST_CASE("grid2D: invalid (zero) spacing is rejected") { + gridParams2D p; + p.topology = 'u'; + p.m = 3; p.n = 3; + p.dx = 0.0; p.dy = 1.0; + grid2D g(p); + CHECK(!g.isValidatedGrid()); +} + +TEST_CASE("grid2D: curvilinear without user nodes is rejected") { + gridParams2D p; + p.topology = 'c'; + p.m = 3; p.n = 3; + grid2D g(p); + CHECK(!g.isValidatedGrid()); +} + +TEST_CASE("grid2D: periodic BC flags round-trip through the grid") { + gridParams2D p; + p.topology = 'u'; + p.m = 3; p.n = 3; + p.dx = 1.0; p.dy = 1.0; + p.bc_isPeriodic[0] = true; + p.bc_isPeriodic[1] = false; + grid2D g(p); + REQUIRE(g.isValidatedGrid()); + CHECK(g.grid.bc_isPeriodic[0] == true); + CHECK(g.grid.bc_isPeriodic[1] == false); +} + +MOLE_TEST_MAIN() diff --git a/cpp/tests/grids/test_grid3D.cpp b/cpp/tests/grids/test_grid3D.cpp new file mode 100644 index 00000000..726a8089 --- /dev/null +++ b/cpp/tests/grids/test_grid3D.cpp @@ -0,0 +1,110 @@ +// Regression tests for grid3D. +// +// As with grid2D, sizes here are deliberately distinct (m != n != o) +// to catch axis-order regressions that a symmetric grid would hide. +#include "MOLE_grids.h" +#include "mole_test.h" + +TEST_CASE("grid3D uniform: valid asymmetric grid validates") { + gridParams3D p; + p.topology = 'u'; + p.m = 3; p.n = 4; p.o = 2; + p.dx = 1.0; p.dy = 1.0; p.dz = 1.0; + grid3D g(p); + CHECK(g.isValidatedGrid()); +} + +TEST_CASE("grid3D uniform: layer shapes follow the (m,n,o) convention") { + gridParams3D p; + p.topology = 'u'; + p.m = 3; p.n = 4; p.o = 2; + p.dx = 1.0; p.dy = 1.0; p.dz = 1.0; + grid3D g(p); + REQUIRE(g.isValidatedGrid()); + + CHECK(g.grid.nodes_X.data_.n_rows == p.m + 1); + CHECK(g.grid.nodes_X.data_.n_cols == p.n + 1); + CHECK(g.grid.nodes_X.data_.n_slices == p.o + 1); + + CHECK(g.grid.centers_X.data_.n_rows == p.m + 2); + CHECK(g.grid.centers_X.data_.n_cols == p.n + 2); + CHECK(g.grid.centers_X.data_.n_slices == p.o + 2); + + CHECK(g.grid.faces_u_X.data_.n_rows == p.m + 1); + CHECK(g.grid.faces_u_X.data_.n_cols == p.n); + CHECK(g.grid.faces_u_X.data_.n_slices == p.o); + + CHECK(g.grid.faces_v_X.data_.n_rows == p.m); + CHECK(g.grid.faces_v_X.data_.n_cols == p.n + 1); + CHECK(g.grid.faces_v_X.data_.n_slices == p.o); + + CHECK(g.grid.faces_w_X.data_.n_rows == p.m); + CHECK(g.grid.faces_w_X.data_.n_cols == p.n); + CHECK(g.grid.faces_w_X.data_.n_slices == p.o + 1); +} + +TEST_CASE("grid3D uniform: user-supplied, correctly-valued nodes " + "round-trip successfully") { + size_t m = 3, n = 4, o = 2; + double dx = 1.0, dy = 2.0, dz = 0.5; + array3D nx(m + 1, n + 1, o + 1, 0.0); + array3D ny(m + 1, n + 1, o + 1, 0.0); + array3D nz(m + 1, n + 1, o + 1, 0.0); + for (size_t i = 0; i <= m; ++i) + for (size_t j = 0; j <= n; ++j) + for (size_t k = 0; k <= o; ++k) { + nx.data_(i, j, k) = i * dx; + ny.data_(i, j, k) = j * dy; + nz.data_(i, j, k) = k * dz; + } + + gridParams3D p; + p.topology = 'u'; + p.m = m; p.n = n; p.o = o; + p.dx = dx; p.dy = dy; p.dz = dz; + p.nodes_X = nx; p.nodes_Y = ny; p.nodes_Z = nz; + grid3D g(p); + CHECK(g.isValidatedGrid()); +} + +TEST_CASE("grid3D uniform: user-supplied nodes with wrong shape " + "are rejected") { + gridParams3D p; + p.topology = 'u'; + p.m = 3; p.n = 4; p.o = 2; + p.dx = 1.0; p.dy = 1.0; p.dz = 1.0; + p.nodes_X = array3D(p.n + 1, p.m + 1, p.o + 1, 0.0); // wrong shape + grid3D g(p); + CHECK(!g.isValidatedGrid()); +} + +TEST_CASE("grid3D: invalid spacing is rejected") { + gridParams3D p; + p.topology = 'u'; + p.m = 2; p.n = 2; p.o = 2; + p.dx = 1.0; p.dy = -1.0; p.dz = 1.0; + grid3D g(p); + CHECK(!g.isValidatedGrid()); +} + +TEST_CASE("grid3D: curvilinear without user nodes is rejected") { + gridParams3D p; + p.topology = 'c'; + p.m = 2; p.n = 2; p.o = 2; + grid3D g(p); + CHECK(!g.isValidatedGrid()); +} + +TEST_CASE("grid3D: constructor with inbound errors merges them in") { + std::stack inerrs; + MOLEerr_log(inerrs, MOLE_ERR_INVALID_INPUT_TYPE, "caller", "note"); + + gridParams3D p; + p.topology = 'u'; + p.m = 2; p.n = 2; p.o = 2; + p.dx = 1.0; p.dy = 1.0; p.dz = 1.0; + grid3D g(p, inerrs); + CHECK(g.hasGridErrors()); +} + +MOLE_TEST_MAIN() diff --git a/cpp/tests/grids/test_grid_builder.cpp b/cpp/tests/grids/test_grid_builder.cpp new file mode 100644 index 00000000..42f847ed --- /dev/null +++ b/cpp/tests/grids/test_grid_builder.cpp @@ -0,0 +1,320 @@ +/* + * SPDX-License-Identifier: GPL-3.0-or-later + * Copyright (c) 2008-2024 San Diego State University Research Foundation + * (SDSURF). + * See LICENSE file or https://www.gnu.org/licenses/gpl-3.0.html for details. + */ + +// Regression tests for grid_builder. +// +// runChecks is exercised directly through gridRaw so every parse-time +// validation branch is observable. The gridBuilder macro is exercised +// for dimension dispatch through the real makeGrid. +#include "grid_builder.h" +#include "mole_test.h" + +#include +#include + +// makeRaw builds a gridRaw carrying only the counts a test sets; +// every unset field keeps its header default. +static gridRaw makeRaw(int dim, char topology, int m = -1, + int n = -1, int o = -1) { + gridRaw g; + g.dim = dim; + g.topology = topology; + g.m = m; + g.n = n; + g.o = o; + return g; +} + +// countErrs reports how many entries in the stack carry a given +// error symbol, so a test can assert on the errors it is about +// rather than on the total size of the stack. +static size_t countErrs(const stack& errs, int code) { + stack tmp = errs; + size_t n = 0; + while (!tmp.empty()) { + if (tmp.top().errCode == code) ++n; + tmp.pop(); + } + return n; +} + +// --------------------------------------------------------------- +// runChecks: dimension validity +// --------------------------------------------------------------- + +TEST_CASE("runChecks: a missing dim returns zero and stops early") { + stack errs; + gridRaw g; // dim defaults to -1 + CHECK(runChecks(errs, g) == 0); + CHECK(MOLEerr_contains(errs, MOLE_ERR_INVALID_GRID_DIM)); + // A missing dim short-circuits before the cell-count checks. + CHECK(!MOLEerr_contains(errs, MOLE_ERR_INVALID_CELL_COUNT)); +} + +TEST_CASE("runChecks: an out-of-range dim returns zero") { + for (int dim : {0, 4}) { + stack errs; + gridRaw g = makeRaw(dim, 'u', 5); + CHECK_MSG(runChecks(errs, g) == 0, "dim = " << dim); + CHECK(MOLEerr_contains(errs, MOLE_ERR_INVALID_GRID_DIM)); + CHECK(!MOLEerr_contains(errs, MOLE_ERR_INVALID_CELL_COUNT)); + } +} + +// --------------------------------------------------------------- +// runChecks: cell-count consistency +// --------------------------------------------------------------- + +TEST_CASE("runChecks: a 1D grid without m is rejected") { + stack errs; + gridRaw g = makeRaw(1, 'u'); // m omitted + CHECK(runChecks(errs, g) == 1); + CHECK(MOLEerr_contains(errs, MOLE_ERR_INVALID_CELL_COUNT)); +} + +TEST_CASE("runChecks: a 2D grid without n is rejected") { + stack errs; + gridRaw g = makeRaw(2, 'u', 5); // n omitted + CHECK(runChecks(errs, g) == 2); + CHECK(MOLEerr_contains(errs, MOLE_ERR_INVALID_CELL_COUNT)); +} + +TEST_CASE("runChecks: a 3D grid without o is rejected") { + stack errs; + gridRaw g = makeRaw(3, 'u', 5, 5); // o omitted + CHECK(runChecks(errs, g) == 3); + CHECK(MOLEerr_contains(errs, MOLE_ERR_INVALID_CELL_COUNT)); +} + +TEST_CASE("runChecks: n supplied for a 1D grid is rejected") { + stack errs; + gridRaw g = makeRaw(1, 'u', 5, 3); + CHECK(runChecks(errs, g) == 1); + CHECK(MOLEerr_contains(errs, MOLE_ERR_INVALID_CELL_COUNT)); +} + +TEST_CASE("runChecks: o supplied for a 2D grid is rejected") { + stack errs; + gridRaw g = makeRaw(2, 'u', 5, 5, 3); + CHECK(runChecks(errs, g) == 2); + CHECK(MOLEerr_contains(errs, MOLE_ERR_INVALID_CELL_COUNT)); +} + +// --------------------------------------------------------------- +// runChecks: topology +// --------------------------------------------------------------- + +TEST_CASE("runChecks: a missing topology is rejected") { + stack errs; + gridRaw g = makeRaw(1, '\0', 5); + runChecks(errs, g); + CHECK(MOLEerr_contains(errs, MOLE_ERR_INVALID_GRID_TOPOLOGY)); +} + +TEST_CASE("runChecks: an unknown topology character is rejected") { + stack errs; + gridRaw g = makeRaw(1, 'x', 5); + runChecks(errs, g); + CHECK(MOLEerr_contains(errs, MOLE_ERR_INVALID_GRID_TOPOLOGY)); +} + +// --------------------------------------------------------------- +// runChecks: valid grids clear every check (also covers the three +// valid topology chars 'u', 'c', 'n') +// --------------------------------------------------------------- + +TEST_CASE("runChecks: a valid uniform 1D grid logs nothing") { + stack errs; + gridRaw g = makeRaw(1, 'u', 5); + CHECK(runChecks(errs, g) == 1); + CHECK(!MOLEerr_haserrors(errs)); +} + +TEST_CASE("runChecks: a valid curvilinear 2D grid logs nothing") { + stack errs; + gridRaw g = makeRaw(2, 'c', 5, 5); + CHECK(runChecks(errs, g) == 2); + CHECK(!MOLEerr_haserrors(errs)); +} + +TEST_CASE("runChecks: a valid nonuniform 3D grid logs nothing") { + stack errs; + gridRaw g = makeRaw(3, 'n', 5, 5, 5); + CHECK(runChecks(errs, g) == 3); + CHECK(!MOLEerr_haserrors(errs)); +} + +// --------------------------------------------------------------- +// runChecks: errors accumulate without short-circuit +// --------------------------------------------------------------- + +TEST_CASE("runChecks: every problem is reported, not just the first") { + stack errs; + gridRaw g = makeRaw(3, 'z'); // all counts missing, bad topo + CHECK(runChecks(errs, g) == 3); + // One cell-count error per missing count, and one for the + // unrecognised topology character. + CHECK_MSG(countErrs(errs, MOLE_ERR_INVALID_CELL_COUNT) == 3, + "got " << countErrs(errs, MOLE_ERR_INVALID_CELL_COUNT) + << " cell-count errors, expected 3"); + CHECK(countErrs(errs, MOLE_ERR_INVALID_GRID_TOPOLOGY) == 1); +} + +// --------------------------------------------------------------- +// runChecks: isPeriodic size against dimension +// +// isPeriodic arrives as a pointer to the caller's vector, so its +// size travels with it and can be compared with dim regardless of +// the order the attributes were supplied in. +// --------------------------------------------------------------- + +TEST_CASE("runChecks: an omitted isPeriodic is accepted") { + stack errs; + gridRaw g = makeRaw(3, 'u', 5, 5, 5); // isPeriodicSrc null + CHECK(runChecks(errs, g) == 3); + CHECK(!MOLEerr_contains(errs, MOLE_ERR_INVALID_ISPERIODIC_DIM)); +} + +TEST_CASE("runChecks: one isPeriodic flag suits a 1D grid") { + stack errs; + std::vector periodic = {true}; + gridRaw g = makeRaw(1, 'u', 5); + g.isPeriodicSrc = &periodic; + CHECK(runChecks(errs, g) == 1); + CHECK(!MOLEerr_haserrors(errs)); +} + +TEST_CASE("runChecks: two isPeriodic flags suit a 2D grid") { + stack errs; + std::vector periodic = {true, false}; + gridRaw g = makeRaw(2, 'u', 5, 5); + g.isPeriodicSrc = &periodic; + CHECK(runChecks(errs, g) == 2); + CHECK(!MOLEerr_haserrors(errs)); +} + +TEST_CASE("runChecks: three isPeriodic flags suit a 3D grid") { + stack errs; + std::vector periodic = {true, false, true}; + gridRaw g = makeRaw(3, 'u', 5, 5, 5); + g.isPeriodicSrc = &periodic; + CHECK(runChecks(errs, g) == 3); + CHECK(!MOLEerr_haserrors(errs)); +} + +TEST_CASE("runChecks: too few isPeriodic flags are rejected") { + stack errs; + std::vector periodic = {true}; // 3D grid needs three + gridRaw g = makeRaw(3, 'u', 5, 5, 5); + g.isPeriodicSrc = &periodic; + runChecks(errs, g); + CHECK(MOLEerr_contains(errs, MOLE_ERR_INVALID_ISPERIODIC_DIM)); +} + +TEST_CASE("runChecks: too many isPeriodic flags are rejected") { + stack errs; + std::vector periodic = {true, false}; // 1D needs one + gridRaw g = makeRaw(1, 'u', 5); + g.isPeriodicSrc = &periodic; + runChecks(errs, g); + CHECK(MOLEerr_contains(errs, MOLE_ERR_INVALID_ISPERIODIC_DIM)); +} + +TEST_CASE("runChecks: an empty isPeriodic vector is rejected") { + stack errs; + std::vector periodic; + gridRaw g = makeRaw(2, 'u', 5, 5); + g.isPeriodicSrc = &periodic; + runChecks(errs, g); + CHECK(MOLEerr_contains(errs, MOLE_ERR_INVALID_ISPERIODIC_DIM)); +} + +// A bad dimension stops the checks before isPeriodic is reached, so +// no size error is reported against a dimension that is itself +// invalid. +TEST_CASE("runChecks: a bad dim skips the isPeriodic size check") { + stack errs; + std::vector periodic = {true, false, true}; + gridRaw g = makeRaw(7, 'u', 5); + g.isPeriodicSrc = &periodic; + CHECK(runChecks(errs, g) == 0); + CHECK(MOLEerr_contains(errs, MOLE_ERR_INVALID_GRID_DIM)); + CHECK(!MOLEerr_contains(errs, MOLE_ERR_INVALID_ISPERIODIC_DIM)); +} + +// --------------------------------------------------------------- +// gridBuilder macro: dimension dispatch through makeGrid +// --------------------------------------------------------------- + +TEST_CASE("gridBuilder: dim 1 yields a grid1D") { + gridVar g = gridBuilder("dim", 1, "m", 5, "dx", 0.2, + "topology", 'u'); + CHECK(std::holds_alternative(g)); +} + +TEST_CASE("gridBuilder: dim 2 yields a grid2D") { + gridVar g = gridBuilder("dim", 2, "m", 5, "n", 5, "dx", 0.2, + "dy", 0.2, "topology", 'u'); + CHECK(std::holds_alternative(g)); +} + +TEST_CASE("gridBuilder: dim 3 yields a grid3D") { + gridVar g = gridBuilder("dim", 3, "m", 5, "n", 5, "o", 5, + "dx", 0.2, "dy", 0.2, "dz", 0.2, + "topology", 'u'); + CHECK(std::holds_alternative(g)); +} + +// --------------------------------------------------------------- +// gridBuilder macro: failure paths return gridNull +// --------------------------------------------------------------- + +TEST_CASE("gridBuilder: an unknown attribute yields a gridNull") { + gridVar g = gridBuilder("dim", 1, "m", 5, "dx", 0.2, + "topology", 'u', "bogus", 7); + CHECK(std::holds_alternative(g)); +} + +TEST_CASE("gridBuilder: a missing dim yields a gridNull") { + gridVar g = gridBuilder("m", 5, "dx", 0.2, "topology", 'u'); + CHECK(std::holds_alternative(g)); +} + +// --------------------------------------------------------------- +// gridBuilder macro: isPeriodic through the full parse path +// --------------------------------------------------------------- + +TEST_CASE("gridBuilder: a matching isPeriodic vector builds a grid") { + std::vector periodic = {true, false}; + gridVar g = gridBuilder("dim", 2, "m", 5, "n", 5, "dx", 0.2, + "dy", 0.2, "topology", 'u', + "isPeriodic", &periodic); + REQUIRE(std::holds_alternative(g)); + grid2D& g2 = std::get(g); + CHECK(g2.grid.bc_isPeriodic[0] == true); + CHECK(g2.grid.bc_isPeriodic[1] == false); +} + +TEST_CASE("gridBuilder: a mismatched isPeriodic vector is rejected") { + std::vector periodic = {true}; + gridVar g = gridBuilder("dim", 2, "m", 5, "n", 5, "dx", 0.2, + "dy", 0.2, "topology", 'u', + "isPeriodic", &periodic); + CHECK(std::holds_alternative(g)); +} + +// The size check reads the vector rather than the argument order, +// so isPeriodic supplied before dim validates the same way. +TEST_CASE("gridBuilder: isPeriodic may precede dim") { + std::vector periodic = {true}; + gridVar g = gridBuilder("isPeriodic", &periodic, "dim", 2, + "m", 5, "n", 5, "dx", 0.2, "dy", 0.2, + "topology", 'u'); + CHECK(std::holds_alternative(g)); +} + +MOLE_TEST_MAIN() diff --git a/cpp/tests/grids/test_makeGrid.cpp b/cpp/tests/grids/test_makeGrid.cpp new file mode 100644 index 00000000..2ea371d2 --- /dev/null +++ b/cpp/tests/grids/test_makeGrid.cpp @@ -0,0 +1,113 @@ +// Regression tests for makeGrid(), gridNull, and isValidGrid(). +#include "MOLE_grids.h" +#include "mole_test.h" +#include + +TEST_CASE("makeGrid dispatches gridParams1D to a valid grid1D") { + gridParams1D p; + p.topology = 'u'; p.m = 4; p.dx = 1.0; + std::stack errs; + gridVar g = makeGrid(p, errs); + REQUIRE(std::holds_alternative(g)); + CHECK(std::get(g).isValidatedGrid()); +} + +TEST_CASE("makeGrid dispatches gridParams2D to a valid grid2D") { + gridParams2D p; + p.topology = 'u'; p.m = 3; p.n = 4; p.dx = 1.0; p.dy = 1.0; + std::stack errs; + gridVar g = makeGrid(p, errs); + REQUIRE(std::holds_alternative(g)); + CHECK(std::get(g).isValidatedGrid()); +} + +TEST_CASE("makeGrid dispatches gridParams3D to a valid grid3D") { + gridParams3D p; + p.topology = 'u'; p.m = 2; p.n = 3; p.o = 2; + p.dx = 1.0; p.dy = 1.0; p.dz = 1.0; + std::stack errs; + gridVar g = makeGrid(p, errs); + REQUIRE(std::holds_alternative(g)); + CHECK(std::get(g).isValidatedGrid()); +} + +TEST_CASE("makeGrid dispatches paramsNull to gridNull, which is " + "always invalid") { + paramsNull np; + np.num_errs = 1; + std::stack errs; + MOLEerr_log(errs, MOLE_ERR_GRID_CONSTRUCTION_FAILED, "prior", ""); + gridVar g = makeGrid(np, errs); + REQUIRE(std::holds_alternative(g)); + CHECK(!std::get(g).validGrid()); +} + +TEST_CASE("makeGrid propagates a pre-existing error stack into the " + "resulting grid") { + gridParams1D p; + p.topology = 'u'; p.m = 4; p.dx = 1.0; + std::stack errs; + MOLEerr_log(errs, MOLE_ERR_INVALID_INPUT_TYPE, "upstream", ""); + gridVar g = makeGrid(p, errs); + REQUIRE(std::holds_alternative(g)); + CHECK(std::get(g).hasGridErrors()); +} + +TEST_CASE("isValidGrid dispatches to the underlying grid's validGrid()") { + gridParams2D p; + p.topology = 'u'; p.m = 3; p.n = 3; p.dx = 1.0; p.dy = 1.0; + std::stack errs; + gridVar g = makeGrid(p, errs); + CHECK(isValidGrid(g)); +} + +TEST_CASE("isValidGrid on an invalid gridVar (invalid topology) " + "returns false") { + gridParams2D p; + p.topology = 'q'; p.m = 3; p.n = 3; p.dx = 1.0; p.dy = 1.0; + std::stack errs; + gridVar g = makeGrid(p, errs); + CHECK(!isValidGrid(g)); +} + +TEST_CASE("NOTE: isValidGrid re-runs validGrid() on an already-" + "constructed grid, which re-appends to its error log") { + // makeGrid() already calls validGrid() once inside the grid's + // constructor. Calling isValidGrid() on the result runs + // validGrid() a SECOND time. For a grid that fails validation, + // this means its error stack accumulates two copies of every + // validation error rather than one. This test documents that + // behavior explicitly, using captured print_ErrorLog() output as + // an observable proxy for "how many times was this error logged" + // (the error stack itself is protected and has no public size + // query). If this is ever changed to be idempotent, this test + // should start failing loudly rather than silently. + gridParams2D p; + p.topology = 'q'; // invalid topology -> validGrid() will fail + p.m = 3; p.n = 3; p.dx = 1.0; p.dy = 1.0; + std::stack errs; + gridVar g = makeGrid(p, errs); // 1st validGrid() call, inside ctor + + CHECK(!isValidGrid(g)); // 2nd validGrid() call; still invalid + + std::ostringstream captured; + std::streambuf* old_cout = std::cout.rdbuf(captured.rdbuf()); + std::visit([](auto&& gridObj) { gridObj.print_ErrorLog(); }, g); + std::cout.rdbuf(old_cout); + + const std::string out = captured.str(); + size_t count = 0, pos = 0; + const std::string needle = "grid2D[construct]"; + while ((pos = out.find(needle, pos)) != std::string::npos) { + ++count; + pos += needle.size(); + } + CHECK_MSG(count == 2, + "expected the invalid-topology error to appear twice after " + "two validGrid() calls (once from the constructor, once from " + "isValidGrid()), got " << count << ". If this is now 1, " + "validGrid() has been made idempotent -- update this test " + "to CHECK(count == 1) and remove this NOTE."); +} + +MOLE_TEST_MAIN() diff --git a/cpp/tests/support/mole_test.h b/cpp/tests/support/mole_test.h new file mode 100644 index 00000000..38fbea12 --- /dev/null +++ b/cpp/tests/support/mole_test.h @@ -0,0 +1,158 @@ +/* + * SPDX-License-Identifier: GPL-3.0-or-later + * + * mole_test.h - a tiny, dependency-free unit test framework for the + * MOLE regression suite. + * + * Why not Catch2/GoogleTest? MOLE's build has no other third-party + * dependency besides Armadillo, and the test suite doesn't need + * anything more than "run named test cases, report pass/fail, + * non-zero exit on failure" to work well with CTest. Each test + * executable is self-contained: include this header, write one or + * more TEST_CASE blocks, and end the file with MOLE_TEST_MAIN(). + * + * Usage: + * + * #include "mole_test.h" + * + * TEST_CASE("array1D default construction is empty") { + * array1D a; + * CHECK(a.data_.is_empty()); + * CHECK(a.data_.n_elem == 0); + * } + * + * MOLE_TEST_MAIN() + * + * CHECK(cond) records a failure and keeps running the rest of the + * test case (use this for independent assertions you want full + * visibility into). REQUIRE(cond) aborts the current test case + * immediately on failure (use this when later checks in the same + * case would be meaningless or crash after a failed precondition, + * e.g. dereferencing something that might not exist). + */ +#ifndef MOLE_TEST_H +#define MOLE_TEST_H + +#include +#include +#include +#include +#include + +namespace moletest { + +struct AssertionFailure { + std::string message; +}; + +struct TestCase { + std::string name; + std::function fn; +}; + +inline std::vector& registry() { + static std::vector r; + return r; +} + +struct Registrar { + Registrar(std::string name, std::function fn) { + registry().push_back({std::move(name), std::move(fn)}); + } +}; + +// Per-test-case failure counter. Reset before each case runs. +inline int& failures_in_case() { + static int n = 0; + return n; +} + +inline void report_failure(const char* kind, const char* condText, + const char* file, int line, + const std::string& extra = "") { + ++failures_in_case(); + std::cerr << " [" << kind << " failed] " << condText + << " (" << file << ":" << line << ")"; + if (!extra.empty()) std::cerr << "\n " << extra; + std::cerr << "\n"; +} + +} // namespace moletest + +#define MOLE_TEST_CONCAT_(a, b) a##b +#define MOLE_TEST_CONCAT(a, b) MOLE_TEST_CONCAT_(a, b) + +#define TEST_CASE(name) \ + static void MOLE_TEST_CONCAT(mole_test_fn_, __LINE__)(); \ + static ::moletest::Registrar MOLE_TEST_CONCAT(mole_test_reg_, \ + __LINE__)(name, MOLE_TEST_CONCAT(mole_test_fn_, __LINE__)); \ + static void MOLE_TEST_CONCAT(mole_test_fn_, __LINE__)() + +#define CHECK(cond) \ + do { \ + if (!(cond)) { \ + ::moletest::report_failure("CHECK", #cond, __FILE__, __LINE__);\ + } \ + } while (0) + +// CHECK_MSG lets a failing check carry a dynamically-built message +// (e.g. the actual vs. expected values), which is often the +// difference between "test failed" and "test failed, here's why". +#define CHECK_MSG(cond, msg) \ + do { \ + if (!(cond)) { \ + std::ostringstream oss_; oss_ << msg; \ + ::moletest::report_failure("CHECK", #cond, __FILE__, __LINE__, \ + oss_.str()); \ + } \ + } while (0) + +#define REQUIRE(cond) \ + do { \ + if (!(cond)) { \ + ::moletest::report_failure("REQUIRE", #cond, __FILE__, \ + __LINE__); \ + throw ::moletest::AssertionFailure{ \ + std::string("REQUIRE failed: ") + #cond}; \ + } \ + } while (0) + +inline int mole_run_all_tests() { + using namespace moletest; + int failed_cases = 0; + const auto& tests = registry(); + std::cout << "Running " << tests.size() << " MOLE test case(s)\n"; + std::cout << "----------------------------------------------------\n"; + for (const auto& tc : tests) { + failures_in_case() = 0; + std::cout << "[ RUN ] " << tc.name << "\n"; + try { + tc.fn(); + } catch (const AssertionFailure&) { + // Already reported by REQUIRE; just stop this case. + } catch (const std::exception& e) { + std::cerr << " [uncaught std::exception] " << e.what() + << "\n"; + ++failures_in_case(); + } catch (...) { + std::cerr << " [uncaught unknown exception]\n"; + ++failures_in_case(); + } + if (failures_in_case() == 0) { + std::cout << "[ OK ] " << tc.name << "\n"; + } else { + std::cout << "[ FAILED ] " << tc.name << "\n"; + ++failed_cases; + } + } + std::cout << "----------------------------------------------------\n"; + std::cout << tests.size() << " test case(s) run, " + << (tests.size() - failed_cases) << " passed, " + << failed_cases << " failed.\n"; + return failed_cases == 0 ? 0 : 1; +} + +#define MOLE_TEST_MAIN() \ + int main() { return mole_run_all_tests(); } + +#endif // MOLE_TEST_H diff --git a/cpp/tests/utils/README.md b/cpp/tests/utils/README.md new file mode 100644 index 00000000..5eec0406 --- /dev/null +++ b/cpp/tests/utils/README.md @@ -0,0 +1,26 @@ + + +# Subdirectory for MOLE 2.0 C++ Tests for the array classes + +Subdirectory and Pathname: **mole/cpp/tests/utils/** + +## Purpose + +Subdirectory containing tests for the MOLE utilities + +## List of Files in This Subdirectory + ++ **test_utils.cpp**: C++ tests of the MOLE utilities ++ **README.md**: (this file) diff --git a/cpp/tests/utils/test_utils.cpp b/cpp/tests/utils/test_utils.cpp new file mode 100644 index 00000000..948d431a --- /dev/null +++ b/cpp/tests/utils/test_utils.cpp @@ -0,0 +1,34 @@ +// Regression tests for utils/utils.cpp. +#include "utils.h" +#include "mole_test.h" +#include + +TEST_CASE("Utils::trapz integrates a known linear function exactly") { + // y = x on [0,4], exact area = 0.5*4*4 = 8 + arma::vec x = {0.0, 1.0, 2.0, 3.0, 4.0}; + arma::vec y = {0.0, 1.0, 2.0, 3.0, 4.0}; + Utils u; + double area = u.trapz(x, y); + CHECK_MSG(std::fabs(area - 8.0) < 1e-9, "got area = " << area); + CHECK(!u.hasErrors()); +} + +TEST_CASE("Utils::trapz integrates a constant function") { + // y = 2 on [0,3], exact area = 6 + arma::vec x = {0.0, 1.0, 2.0, 3.0}; + arma::vec y = {2.0, 2.0, 2.0, 2.0}; + Utils u; + double area = u.trapz(x, y); + CHECK_MSG(std::fabs(area - 6.0) < 1e-9, "got area = " << area); +} + +TEST_CASE("Utils::trapz on mismatched sizes returns NaN and logs an error") { + arma::vec x = {0.0, 1.0, 2.0}; + arma::vec y = {0.0, 1.0}; + Utils u; + double result = u.trapz(x, y); + CHECK(std::isnan(result)); + CHECK(u.hasErrors()); +} + +MOLE_TEST_MAIN() diff --git a/src/cpp/utils.cpp b/src/cpp/utils.cpp deleted file mode 100644 index a35369f3..00000000 --- a/src/cpp/utils.cpp +++ /dev/null @@ -1,258 +0,0 @@ -/* -* SPDX-License-Identifier: GPL-3.0-or-later -* © 2008-2024 San Diego State University Research Foundation (SDSURF). -* See LICENSE file or https://www.gnu.org/licenses/gpl-3.0.html for details. -*/ - - -/* - * @file utils.cpp - * @brief Helpers for sparse operations and MATLAB/Octave analogs - * @date 2024/10/15 - * - * Sparse operations that repeatedly are needed, but not - * necessarily part of the Armadillo library. Some other MATLAB/Octave - * type functions are also here, like meshgrid. - */ - -#include "utils.h" -#include -#include -#include -#include - -#ifdef EIGEN -#include - -vec Utils::spsolve_eigen(const sp_mat &A, const vec &b) { - Eigen::SparseMatrix eigen_A(A.n_rows, A.n_cols); - std::vector> triplets; - Eigen::SparseLU, Eigen::COLAMDOrdering> solver; - - Eigen::VectorXd eigen_x(A.n_rows); - triplets.reserve(5 * A.n_rows); - - auto it = A.begin(); - while (it != A.end()) { - triplets.push_back(Eigen::Triplet(it.row(), it.col(), *it)); - ++it; - } - - eigen_A.setFromTriplets(triplets.begin(), triplets.end()); - triplets.clear(); - - auto b_ = conv_to>::from(b); - Eigen::Map eigen_b(b_.data(), b_.size()); - - solver.analyzePattern(eigen_A); - solver.factorize(eigen_A); - eigen_x = solver.solve(eigen_b); - - return vec(eigen_x.data(), eigen_x.size()); -} -#endif - -// Basic implementation of Kronecker product -/* -sp_mat Utils::spkron(const sp_mat &A, const sp_mat &B) -{ - sp_mat result; - - for (u32 i = 0; i < A.n_rows; i++) { - sp_mat BLOCK; - for (u32 j = 0; j < A.n_cols; j++) { - BLOCK = join_rows(BLOCK, A(i, j)*B); - } - result = join_cols(result, BLOCK); - } - - return result; -} -*/ - -sp_mat Utils::spkron(const sp_mat &A, const sp_mat &B) { - sp_mat::const_iterator itA = A.begin(); - sp_mat::const_iterator endA = A.end(); - sp_mat::const_iterator itB = B.begin(); - sp_mat::const_iterator endB = B.end(); - u32 j = 0; - - vec a = nonzeros(A); - vec b = nonzeros(B); - - umat locations(2, a.n_elem * b.n_elem); - vec values(a.n_elem * b.n_elem); - - while (itA != endA) { - while (itB != endB) { - locations(0, j) = itA.row() * B.n_rows + itB.row(); - locations(1, j) = itA.col() * B.n_cols + itB.col(); - values(j) = (*itA) * (*itB); - ++j; - ++itB; - } - - ++itA; - itB = B.begin(); - } - - sp_mat result(locations, values, A.n_rows * B.n_rows, A.n_cols * B.n_cols, - true); - - return result; -} - - -sp_mat Utils::spjoin_rows(const sp_mat &A, const sp_mat &B) { - sp_mat::const_iterator itA = A.begin(); - sp_mat::const_iterator endA = A.end(); - sp_mat::const_iterator itB = B.begin(); - sp_mat::const_iterator endB = B.end(); - u32 j = 0; - - vec a = nonzeros(A); - vec b = nonzeros(B); - - umat locations(2, a.n_elem + b.n_elem); - vec values(a.n_elem + b.n_elem); - - while (itA != endA) { - locations(0, j) = itA.row(); - locations(1, j) = itA.col(); - values(j) = (*itA); - ++itA; - ++j; - } - - while (itB != endB) { - locations(0, j) = itB.row(); - locations(1, j) = itB.col() + A.n_cols; - values(j) = (*itB); - ++itB; - ++j; - } - - sp_mat result(locations, values, A.n_rows, A.n_cols + B.n_cols, true); - - return result; -} - - -sp_mat Utils::spjoin_cols(const sp_mat &A, const sp_mat &B) { - sp_mat::const_iterator itA = A.begin(); - sp_mat::const_iterator endA = A.end(); - sp_mat::const_iterator itB = B.begin(); - sp_mat::const_iterator endB = B.end(); - u32 j = 0; - - vec a = nonzeros(A); - vec b = nonzeros(B); - - umat locations(2, a.n_elem + b.n_elem); - vec values(a.n_elem + b.n_elem); - - while (itA != endA) { - locations(0, j) = itA.row(); - locations(1, j) = itA.col(); - values(j) = (*itA); - ++itA; - ++j; - } - - while (itB != endB) { - locations(0, j) = itB.row() + A.n_rows; - locations(1, j) = itB.col(); - values(j) = (*itB); - ++itB; - ++j; - } - - sp_mat result(locations, values, A.n_rows + B.n_rows, A.n_cols, true); - - return result; -} - - -void Utils::meshgrid(const vec &x, const vec &y, mat &X, mat &Y) { - int m = x.n_elem; - int n = y.n_elem; - - assert(m > 0); - assert(n > 0); - - // Build X - vec t(n, fill::ones); - - X.zeros(n, m); - Y.zeros(n, m); - - for (int ii = 0; ii < m; ++ii) { - X.col(ii) = x(ii) * t; - t.ones(); - } - - for (int ii = 0; ii < m; ++ii) - Y.col(ii) = y; -} - - -void Utils::meshgrid(const vec &x, const vec &y, const vec &z, cube &X, cube &Y, - cube &Z) { - int m = x.n_elem; - int n = y.n_elem; - int o = z.n_elem; - - assert(m > 0); - assert(n > 0); - assert(o > 0); - - // Temporary Holder of sheet of cube - mat sheet(m, n, fill::zeros); - - // Build X - vec t(n, fill::ones); - - X.zeros(m, n, o); - Y.zeros(m, n, o); - Z.zeros(m, n, o); - - // Sheet that repeats each slice - for (int ii = 0; ii < m; ++ii) { - sheet.row(ii) = x(ii) * t.t(); - t.ones(); - } - - for (int kk = 0; kk < o; ++kk) - X.slice(kk) = sheet; - - // Y Cube, repeats same sheet as well - for (int ii = 0; ii < m; ++ii) - sheet.row(ii) = y.t(); - - for (int kk = 0; kk < o; ++kk) - Y.slice(kk) = sheet; - - // Z cube goes by slices each with same value - for (int kk = 0; kk < o; ++kk) - Z.slice(kk).fill(z(kk)); -} - -// Trapezoidal rule (trapz) for 1D integration -double Utils::trapz(const vec &x, const vec &y) { - assert(x.n_elem == y.n_elem); - double sum = 0.0; - for (uword i = 0; i < x.n_elem - 1; ++i) { - sum += (x(i+1) - x(i)) * (y(i) + y(i+1)); - } - return 0.5 * sum; -} - -// Spacing validation shared across every operator entry point. -void mole::check_spacing(Real h, const char* name) { - if (!std::isfinite(h) || h <= 0.0) { - throw std::invalid_argument( - std::string("MOLE: ") + name + - " must be a positive finite number, got " + std::to_string(h)); - } -} -