diff --git a/CMakeLists.txt b/CMakeLists.txt index 64c8e9a..2eaa7e3 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,28 +1,38 @@ # SPDX-FileCopyrightText: 2015 UCLA -# SPDX-FileCopyrightText: 2025 Bernhard Haas (GFZ) +# SPDX-FileCopyrightText: 2025 GFZ Helmholtz Centre for Geosciences # # SPDX-License-Identifier: BSD-3-Clause -# CMakeLists.txt -# Buildsystem for VERB4D - See README.md for build instructions -# -# Author: Anthony Miyaguchi (acmiyaguchi@gmail.com) +cmake_minimum_required(VERSION 3.14) +project(VERB4D LANGUAGES C CXX) -cmake_minimum_required(VERSION 3.9) -project(VERB4D C CXX) +# --- Project Settings --- set(CMAKE_CXX_STANDARD 17) +set(CMAKE_CXX_STANDARD_REQUIRED ON) +list(APPEND CMAKE_MODULE_PATH "${CMAKE_CURRENT_SOURCE_DIR}/cmake") + +# --- Options --- +option(DATA_ASSIMILATION "Enable data assimilation" OFF) +option(DATA_ASSIMILATION_DEBUG "Enable debug output for data assimilation" OFF) +option(USE_PPFV "Compile positive preserving finite volume method" OFF) +option(FAST_CONVECTION "Use larger time steps for convection (less accurate)" OFF) +option(PYTHON_BINDINGS "Build Python bindings using pybind11" OFF) +option(SAVE_PSD_LOST_CONV "Save lost PSD during convection step" OFF) +option(LU_CACHING "Cache LU decompositions during Lapack calls" OFF) + +# --- Dependencies --- +find_package(OpenMP) +include(GetBlasType) # this defines BLAS_LIBRARIES_TO_LINK and BLAS_INCLUDE_DIRS + +if(USE_PPFV) + message(STATUS "PPFV enabled: Finding or downloading Eigen3 and xtensor...") + include(FindOrDownloadEigen) + include(FindOrDownloadXTensor) +endif() -# Enable/disable data assimilation -set(DATA_ASSIMILATION FALSE CACHE BOOL "Enable data assimilation") - -# Math Constants are not defined in Standard C++. To use them, we must first define _USE_MATH_DEFINES -add_compile_definitions(_USE_MATH_DEFINES) - -# Recursively search for header files to include for project files. This isn't -# necessary to compile, but its nice to have in Visual Studio solutions. -file(GLOB_RECURSE VERB4D_HEADERS *.h) +# --- Source Files --- +file(GLOB_RECURSE VERB4D_HEADERS src/*.h) -# Add new source files here and re-run cmake set(VERB4D_SOURCES src/VERB4D_Solver.cpp src/Convection_1D_ULTIMATE_QUICKEST6.cpp @@ -42,165 +52,95 @@ set(VERB4D_SOURCES src/MonotCubicInterpolator.cpp src/Parameters.cpp src/ReadInitialData.cpp - src/UpdatableMatrix.cpp) + src/UpdatableMatrix.cpp +) -# Add data assimilation source files and macros if(DATA_ASSIMILATION) - add_definitions(-DDATA_ASSIMILATION) - - set(VERB4D_SOURCES - ${VERB4D_SOURCES} + message(STATUS "Data Assimilation enabled") + list(APPEND VERB4D_SOURCES src/MatrixOperations.cpp src/CustomDate.cpp src/DataAssimilationHelper.cpp - src/DataAssimilation.cpp) - - if(MATLAB_FOUND) - set(VERB4D_SOURCES - ${VERB4D_SOURCES} - src/PMF.cpp) - endif(MATLAB_FOUND) - - if(DATA_ASSIMILATION_DEBUG) - message("Saving DA debug output!") - add_definitions(-DDATA_ASSIMILATION_DEBUG) - endif(DATA_ASSIMILATION_DEBUG) + src/DataAssimilation.cpp + ) +endif() -endif(DATA_ASSIMILATION) +if(USE_PPFV) + list(APPEND VERB4D_SOURCES src/Diffusion_2D_PPFV.cpp) +endif() -# Build and link +# --- Main Executable Target --- add_executable(VERB4D_Solver ${VERB4D_SOURCES} ${VERB4D_HEADERS}) -# Apply different compiler options based on the compiler +# --- Compiler Definitions and Options --- +target_compile_definitions(VERB4D_Solver PRIVATE + _USE_MATH_DEFINES + $<$:DATA_ASSIMILATION> + $<$:DATA_ASSIMILATION_DEBUG> + $<$:USE_PPFV> + $<$:FAST_CONVECTION> + $<$:SAVE_PSD_LOST_CONV> + $<$:LU_CACHING> + $<$:RELEASE> +) + if(MSVC) - # For MSVC, use /W4 for warnings and /WX to treat warnings as errors. - # /wd4804 - T should not be bool in some cases. TODO: solve later - # /wd4996 'fpclassify': ambiguous call to overloaded function [C:/Program Files (x86)/Windows Kits/10/Include/10.0.22621.0/ucrt/corecrt_math.h"] - # /wd4849 OpenMP 'collapse' clause ignored in 'parallel for' directive - set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} /W4 /WX /permissive /wd4804 /wd4996 /wd4849") + target_compile_options(VERB4D_Solver PRIVATE /W4) else() - # For GCC/Clang, use -Wall -Wextra -Wpedantic (-Werror removed) - set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -Wextra -Wpedantic") + target_compile_options(VERB4D_Solver PRIVATE -Wall -Wextra -Wpedantic) endif() -# Add directories for include files -target_include_directories(VERB4D_Solver PUBLIC ${PROJECT_SOURCE_DIR}/src) - - -# Check if interprocedural optimization (link time optimization) is supported by the compiler -include(CheckIPOSupported) -check_ipo_supported(RESULT ipo_supported OUTPUT ipo_error) - -if(NOT DEFINED BLAS_TYPE) - message("No BLAS_TYPE specified, choosing LAPACK") - set(BLAS_TYPE LAPACK) -endif() -string(TOUPPER ${BLAS_TYPE} BLAS_TYPE) -if(${BLAS_TYPE} STREQUAL "LAPACK") - if(WIN32) - # The libraries for lapack and BLAS for windows are in the lib folder - set(LAPACK_LIBRARIES ${CMAKE_SOURCE_DIR}/lib/clapack.lib) - set(BLAS_LIBRARIES ${CMAKE_SOURCE_DIR}/lib/blas.lib) - set(F2C_LIBRARIES ${CMAKE_SOURCE_DIR}/lib/libf2c.lib) - list(APPEND DEPLIBS ${F2C_LIBRARIES}) - else(WIN32) - # Installed by distribution's package manager - find_package(LAPACK REQUIRED) - endif(WIN32) - list(APPEND DEPLIBS ${LAPACK_LIBRARIES} ${BLAS_LIBRARIES}) -elseif(${BLAS_TYPE} STREQUAL "OPENBLAS") - add_library(openblas SHARED IMPORTED) - if(WIN32) - if(NOT DEFINED OPENBLAS_LIBPATH) - set(OPENBLAS_LIBPATH $ENV{HOMEPATH}/AppData/Local/OpenBLAS/openblas.lib) - message("No OPENBLAS_LIBPATH specified, choosing ${OPENBLAS_LIBPATH}") - endif() - set_target_properties(openblas PROPERTIES IMPORTED_IMPLIB ${OPENBLAS_LIBPATH}) - else() - if(NOT DEFINED OPENBLAS_LIBPATH) - set(OPENBLAS_LIBPATH $ENV{HOME}/.local/lib/libopenblas.so) - message("No OPENBLAS_LIBPATH specified, choosing ${OPENBLAS_LIBPATH}") - endif() - set_target_properties(openblas PROPERTIES IMPORTED_LOCATION ${OPENBLAS_LIBPATH}) - endif() - list(APPEND DEPLIBS openblas) -elseif(${BLAS_TYPE} STREQUAL "MKL") - set(MKL_THREADING sequential) - set(MKL_INTERFACE lp64) - find_package(MKL CONFIG REQUIRED) - list(APPEND DEPLIBS ${MKL_IMPORTED_TARGETS}) - target_include_directories(VERB4D_Solver PUBLIC ${MKL_INCLUDE}) - add_definitions(-DMKL_FOUND) -else() - message( FATAL_ERROR "Unknown BLAS_TYPE. Choose either LAPACK, OPENBLAS or MKL" ) +# --- Include Directories --- +target_include_directories(VERB4D_Solver PUBLIC src) +if(BLAS_INCLUDE_DIRS) + target_include_directories(VERB4D_Solver PRIVATE ${BLAS_INCLUDE_DIRS}) endif() -# Find OpenMP -find_package(OpenMP REQUIRED) -if(OPENMP_FOUND) - set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${OpenMP_C_FLAGS}") - set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}") - set(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} ${OpenMP_EXE_LINKER_FLAGS}") -endif(OPENMP_FOUND) - -# Find Matlab -find_package(Matlab COMPONENTS ENG_LIBRARY MAT_LIBRARY) - -if(Matlab_FOUND) - include_directories(${Matlab_INCLUDE_DIRS}) - list(APPEND DEPLIBS - ${Matlab_MAT_LIBRARY} - ${Matlab_MX_LIBRARY} - ${Matlab_ENG_LIBRARY}) - # Set the definition for MATLAB_CAPABLE - add_definitions(-DMATLAB_CAPABLE=true) -endif() +# --- Link Libraries --- +target_link_libraries(VERB4D_Solver PRIVATE ${BLAS_LIBRARIES_TO_LINK}) -# use interprocidural optimization if supported -if( ipo_supported ) - message(STATUS "IPO enabled") - set_property(TARGET VERB4D_Solver PROPERTY INTERPROCEDURAL_OPTIMIZATION TRUE) +if(OpenMP_CXX_FOUND) + message(STATUS "OpenMP found! Parallelization enabled.") + target_link_libraries(VERB4D_Solver PRIVATE OpenMP::OpenMP_CXX) else() - message(STATUS "IPO not supported: <${ipo_error}>") + message(WARNING "OpenMP not found. VERB4D will be compiled in SERIAL mode.") endif() -target_link_libraries(VERB4D_Solver ${DEPLIBS}) - -# Flag compile type feedback -if(DEFINED CMAKE_BUILD_TYPE) - if(CMAKE_BUILD_TYPE STREQUAL "") - message(WARNING "\nCMAKE_BUILD_TPYE not specified. Run cmake with -DCMAKE_BUILD_TYPE=Release or -DCMAKE_BUILD_TYPE=Debug") - else(CMAKE_BUILD_TYPE STREQUAL "") - string(TOUPPER ${CMAKE_BUILD_TYPE} CMAKE_BUILD_TYPE_upper) - if(CMAKE_BUILD_TYPE_upper STREQUAL "RELEASE") - add_definitions(-DRELEASE) - endif() - endif() +# --- PPFV --- +if(USE_PPFV) + set(PPFV_SOURCES + src/ppfv/Mesh.cpp + src/ppfv/Solver.cpp + src/ppfv/Cases/VERB_ppfv2d.cpp + ) + add_library(ppfv SHARED ${PPFV_SOURCES} ${VERB4D_HEADERS}) + + target_link_libraries(ppfv PUBLIC Eigen3::Eigen) + target_link_libraries(ppfv PUBLIC xtensor) + + # Link the internal library to the main solver + target_link_libraries(VERB4D_Solver PRIVATE ppfv) endif() -# Flag to decide convection step size in Convection_2D.cpp -if(${FAST_CONVECTION}) - add_definitions(-DFAST_CONVECTION) +# --- Optimization (IPO/LTO) --- +include(CheckIPOSupported) +check_ipo_supported(RESULT ipo_supported OUTPUT ipo_error) +if(ipo_supported) + set_property(TARGET VERB4D_Solver PROPERTY INTERPROCEDURAL_OPTIMIZATION TRUE) + message(STATUS "IPO/LTO enabled for VERB4D_Solver") endif() -unset(FAST_CONVECTION CACHE) -# Python bindings -if (${PYTHON_BINDINGS}) - include_directories(${PROJECT_SOURCE_DIR}/extern/pybind11/include ${PROJECT_SOURCE_DIR}/src/python_bindings ${PROJECT_SOURCE_DIR}/src) +# --- Python Bindings --- +if(PYTHON_BINDINGS) + message(STATUS "Building Python bindings with pybind11...") add_subdirectory(extern/pybind11) - pybind11_add_module(verb4d_solver src/python_bindings/Define_python_modules.cpp ${VERB4D_SOURCES}) - target_link_libraries(verb4d_solver PUBLIC ${DEPLIBS}) -endif() -unset(PYTHON_BINDINGS CACHE) -# save lost PSD during Convection_2D -if(${SAVE_PSD_LOST_CONV}) - add_definitions(-DSAVE_PSD_LOST_CONV) -endif() -unset(SAVE_PSD_LOST CACHE) + pybind11_add_module(verb4d_solver src/python_bindings/Define_python_modules.cpp ${VERB4D_SOURCES}) + target_link_libraries(verb4d_solver PUBLIC OpenMP::OpenMP_CXX ${BLAS_LIBRARIES_TO_LINK}) + include_directories(${PROJECT_SOURCE_DIR}/src/python_bindings ${PROJECT_SOURCE_DIR}/src) -# cache LU decomposition during Lapack call in Diffusion_2D.cpp -if(${LU_CACHING}) - add_definitions(-DLU_CACHING) -endif() -unset(SAVE_PSD_LOST CACHE) + if(USE_PPFV) + target_compile_definitions(verb4d_solver PRIVATE USE_PPFV) + target_link_libraries(verb4d_solver PUBLIC ppfv) + endif() +endif() \ No newline at end of file diff --git a/README.md b/README.md index 0ad76a8..4a24f05 100644 --- a/README.md +++ b/README.md @@ -1,6 +1,6 @@ @@ -10,16 +10,15 @@ Source code for the VERB-4D solver. Need to be used together with Matlab or Pyth create all necessary input. The main file is in `./src/VERB4D_Solver.cpp` ## Requirements -Download cmake 3.7 or above. Note that any new source files should be added to +Download cmake 3.14 or above. Note that any new source files should be added to the VERB4D_SOURCES line in CMakeLists.txt, and the build files must be regenerated with the below commands. The following dependencies are required and will be automatically searched via cmake. Use homebrew or the package manager of your choice to install them: -* OpenMP +* OpenMP (optional) * Lapack, OpenBLAS (added as a submodule) or Intel's MKL -* MatLab (optional) A C++ compiler is necessary, a Fortran compiler (only for OpenBLAS) recommended. If no Fortran compiler is installed, OpenBLAS will compile from C code instead. @@ -32,63 +31,39 @@ If CMAKE cannot find the correct g++ version, you can set the environmental vari `export CXX=/u/local/compilers/gcc/10.2.0/bin/g++` #### Windows -##### Visual Studio -The MSVC C++ compiler can be installed using [Visual Studio](https://visualstudio.microsoft.com/downloads/). -When choosing OpenBLAS, the only option in Visual Studio is to compile OpenBLAS from C code, -which might be less optimized than building from the Fortran source code. - -##### MinGW and Ninja -The recommended way of building OpenBLAS on Windows is by using the [Ninja](https://ninja-build.org/) -build system instead of Visual Studio. Additionally, [MSYS2](https://www.msys2.org/) is necessary -to install MinGW's gcc and gfortran compiler. To install both compilers, open -an MSYS2 MinGW terminal and run -``` -pacman -S mingw-w64-x86_64-gcc -pacman -S mingw-w64-x86_64-gcc-fortran -pacman -Syu -``` -Make sure that the location of `ninja.exe` and the MinGW binary folder, e.g. `C:\msys64\mingw64\bin\`, are on your PATH. -#### OSX -You will need to install gcc or g++ in order to support compilation -with OpenMP. The version of gcc on OSX 10.10 is actually linked to clang. You -will need to install a version of gcc from homebrew or ports. To test that it -works, run the command `gcc-5 --version` and verify that it is indeed gcc and -not clang. +*coming soon* ## Build instructions Clone the verb4d_solver repository #### via ssh ``` -git clone git@rbm9.gfz-potsdam.de:verb/verb4d_solver.git +git clone git@git.gfz-potsdam.de:spsw/verb/verb4d/verb4d_solver.git ``` #### via https ``` -git clone https://rbm9.gfz-potsdam.de/verb/verb4d_solver.git +git clone https://git.gfz-potsdam.de/spsw/verb/verb4d/verb4d_solver.git ``` -Add the `--recursive` flag to also clone the OpenBLAS submodule. Enter the verb4d_solver directory with `cd verb4d_solver`. If you cloned verb4d_solver non-recursively, you can update the OpenBLAS submodule by executing `git submodule update --init`. - -### Compile and install the OpenBLAS library (optional) -#### Linux / OSX -Build the shared OpenBLAS library and install it to `~/.local/lib` by calling -``` -./install_openblas.sh -``` -### Windows -For the Visual Studio build, open a PowerShell terminal and execute -``` -install_openblas_vs.ps1 -``` -For the Ninja build, make sure that both gcc and gfortran are installed (e.g. `gfortran --version`). -Call the Ninja install script -``` -install_openblas_ninja.ps1 -``` -Both build types compile OpenBLAS as a dynamic library and copy the library files to `:\Users\\AppData\Local\OpenBLAS\`. -This location is automatically appended to the user's PATH environment variable, so VERB-4D can use it at runtime. +Add the `--recursive` flag to also clone the OpenBLAS and Pybind11 submodule. Enter the verb4d_solver directory with `cd verb4d_solver`. If you cloned verb4d_solver non-recursively, you can update the submodules by executing `git submodule update --init`. ### Compile the VERB-4D solver + +#### List of compilation flags + +| Command line argument | Description | Possible values +| --------------------- | ----------- | -------------- | +| -DCMAKE_BUILD_TYPE | Standard cmake option enabling debugging symbols or optimization flags | Release, Debug | +| -DBLAS_TYPE | Library type of LAPACK and BLAS routines | LAPACK, OpenBLAS, MKL | +| -DOPENBLAS_LIBPATH | Path towards compiled openblas.lib | | +| -DDATA_ASSIMILATION | Bool whether the data assimilation routines will be compiled | True, False | +| -DDATA_ASSIMILATION_DEBUG | Bool whether the data assimilation debug output will be saved (!huge file sizes!) | True, False | +| -DFAST_CONVECTION | Bool to decide convection step size in Convection_2D | True, False | +| -DPYTHON_BINDINGS | Bool whether the python-bindings-library will be compiled | True, False | +| -DSAVE_PSD_LOST_CONV | Bool whether the lost PSD during Convection_2d will be saved | True, False | +| -DLU_CACHING | Bool whether LU caching will be used for local diffusion. It's faster, but needs more memory. | True, False | +| -DUSE_PPFV | Bool wheter the positive-preserving finite volume method should be compiled. | True, False | + #### Linux / OSX ##### CMake Manually compile the VERB-4D solver by creating a build folder @@ -110,118 +85,19 @@ The executable `VERB4D_Solver` can be found in `./build`. #### Windows ##### CMake (PowerShell) -Manually compile the VERB-4D solver by creating a build folder -``` -mkdir build && cd build -``` -and call CMake for a release built and the BLAS library of your choice, either Lapack (default), OpenBLAS or MKL. -``` -cmake .. -DBLAS_TYPE=OpenBLAS -cmake --build . -j --config Release -``` -For the first cmake call, when choosing OpenBLAS, an additional parameter `-DOPENBLAS_LIBPATH=` can be specified
- (default `:\Users\\AppData\Local\OpenBLAS\openblas.lib`) -##### PowerShell script -Automatically compile the VERB-4D solver (using OpenBLAS) by calling +*coming soon* + +### Compile and install the OpenBLAS library (optional) +#### Linux / OSX +Build the shared OpenBLAS library and install it to `~/.local/lib` by calling ``` -install_verb.ps1 +./install_openblas.sh ``` -The executable `VERB4D_Solver.exe` can be found in `.\build\Release`. - -#### List of compilation flags - -| Command line argument | Description | Possible values -| --------------------- | ----------- | -------------- | -| -DCMAKE_BUILD_TYPE | Standard cmake option enabling debugging symbols or optimization flags | Release, Debug | -| -DBLAS_TYPE | Library type of LAPACK and BLAS routines | LAPACK, OpenBLAS, MKL | -| -DOPENBLAS_LIBPATH | Path towards compiled openblas.lib | | -| -DDATA_ASSIMILATION | Bool whether the data assimilation routines will be compiled | True, False | -| -DDATA_ASSIMILATION_DEBUG | Bool whether the data assimilation debug output will be saved (!huge file sizes!) | True, False | -| -DFAST_CONVECTION | Bool to decide convection step size in Convection_2D | True, False | -| -DPYTHON_BINDINGS | Bool whether the python-bindings-library will be compiled | True, False | -| -DSAVE_PSD_LOST_CONV | Bool whether the lost PSD during Convection_2d will be saved | True, False | -| -DLU_CACHING | Bool whether LU caching will be used for local diffusion. It's faster, but needs more memory. | True, False | - -### Use the VERB-4D solver - -Copy the executable to a simulation folder that contains a `parameters.ini` file and run it to start a simulation. - - -## Old Notes (Possibly outdated) -ReadMe - COMPILING ON MAC - -Note: All tabbed items should be executed on the command line - -Xcode does not come with openmp support so a different compiler is necessary. -However, you will still need Xcode for some tools in order to build the new -compiler - -download Xcode with developer tools - located in app store (free) (developer -tools might need to be downloaded separately) - -in terminal write: xcode-select --install - -download gcc-5.1.0 as explained at -https://solarianprogrammer.com/2015/05/01/compiling-gcc-5-mac-os-x/ NOTE: some -of the links seem to not work anymore. Just copy and paste the name of the -desired file to download into google and you should be able to find it. - -Get the verb4d folder from gitlab - -In the terminal go to the src directory which contains VERB4D_Solver.cpp from -verb4d cd ~/…(your path if any).../verb4d/VERB4D_Solver/src - -add paths for both gcc and omp to work export PATH=/usr/gcc-5.1.0/bin:$PATH -Assuming you kept the gcc-5.1.0 file in your downloads like in the above -tutorial export PATH=~/Downloads/gcc-5.1.0/libgomp:$PATH - -Might need to change to in the file Matrix.h - -Create the executable make - -You now have an executable name VERB4D_Solver - -To test that it works: - -Copy the executable “VERB4D_Solver” into the 1.0 1D_Alpha_Example folder cp -VERB4D_Solver ../../Examples/1.0\ 1D_Alpha_Example/VERB4D_Solver - -Open Conv_Dif.m from that file in matlab - -Scroll down to the run portion and change system('VERB4D_Solver.exe’) to -system('./VERB4D_Solver') - -press run (play button) at top of page - -COMPILING not on mac - -Must have OPENMP2.0 or above - comes with visual studio etc. Must have c++ -compiler such as gcc or g++ etc. Change the CC variable in Makefile to -whatever compiler you are using - -MATLAB CAPABILITY - -If you have matlab installed on your machine then you can opt to read/write -files into .mat format which is faster to use and more precise than text files. - -Make sure MATLAB_CAPABLE is set to true in Matrix.h - -In all Conv_Dif.m files there is a variable named use_matlab. Set this to true -as well. - -Also to compile the solver you must link to the matlab/c++ shared libraries: -for mac they can be found at /Applications/MATLAB_R2015a.app/bin/maci64 -assuming you downloaded matlab into applications folder. Look at the Makefile -for an idea of how to link the proper Matlab libraries for compiling. - -For correct path type matlabroot into matlab, and then add on the necessary -extensions such as extern/include and bin/os +#### Windows +*coming soon* -If you do not have matlab capabilities use the MakefileNoMatlab to compile -after setting MATLAB_CAPABLE to false in Matrix.h. +## Contributors and Acknowledgements -NOTES: +We thank all the contributors of the VERB-4D Solver. -For the standard testing of an example include_boundary should be true and -inversion method should be lapack. The ADI inversion methods are not working -properly as of Sept 2015 +We acknowledge Xin Tao and all developers of the [Sayram code](https://github.com/xtaohub/Sayram-2D/) for providing the C++ code for the positive-preserving finite volume scheme. \ No newline at end of file diff --git a/cmake/FindOrDownloadEigen.cmake b/cmake/FindOrDownloadEigen.cmake new file mode 100644 index 0000000..e74ee0b --- /dev/null +++ b/cmake/FindOrDownloadEigen.cmake @@ -0,0 +1,20 @@ +include(FetchContent) + +find_package(Eigen QUIET) + +if(NOT Eigen_FOUND) + FetchContent_Declare( + Eigen + GIT_REPOSITORY https://gitlab.com/libeigen/eigen.git + GIT_TAG 3.4 + GIT_SHALLOW TRUE + GIT_PROGRESS TRUE + ) + + set(BUILD_TESTING OFF) + set(EIGEN_BUILD_TESTING OFF) + set(EIGEN_MPL2_ONLY ON) + set(EIGEN_BUILD_PKGCONFIG OFF) + set(EIGEN_BUILD_DOC OFF) + FetchContent_MakeAvailable(Eigen) +endif() \ No newline at end of file diff --git a/cmake/FindOrDownloadXTensor.cmake b/cmake/FindOrDownloadXTensor.cmake new file mode 100644 index 0000000..30d1f14 --- /dev/null +++ b/cmake/FindOrDownloadXTensor.cmake @@ -0,0 +1,18 @@ +include(FetchContent) + +find_package(xtensor QUIET) + +if(NOT xtensor_FOUND) + FetchContent_Declare(xtl GIT_REPOSITORY https://github.com/xtensor-stack/xtl.git GIT_TAG 0.7.6) + FetchContent_MakeAvailable(xtl) + + FetchContent_Declare( + xtensor + GIT_REPOSITORY https://github.com/xtensor-stack/xtensor.git + GIT_TAG 0.24.7 + GIT_SHALLOW TRUE + GIT_PROGRESS TRUE + ) + FetchContent_MakeAvailable(xtensor) + +endif() \ No newline at end of file diff --git a/cmake/GetBlasType.cmake b/cmake/GetBlasType.cmake new file mode 100644 index 0000000..3c20d46 --- /dev/null +++ b/cmake/GetBlasType.cmake @@ -0,0 +1,60 @@ +if(NOT DEFINED BLAS_TYPE) + message("No BLAS_TYPE specified, choosing LAPACK") + set(BLAS_TYPE LAPACK) +endif() + +string(TOUPPER ${BLAS_TYPE} BLAS_TYPE) + +# Initialize a local list for dependencies +set(BLAS_LIBRARIES_TO_LINK "") + +if(${BLAS_TYPE} STREQUAL "LAPACK") + + if(WIN32) + # The libraries for lapack and BLAS for windows are in the lib folder + set(LAPACK_LIBRARIES ${CMAKE_SOURCE_DIR}/lib/clapack.lib) + set(BLAS_LIBRARIES ${CMAKE_SOURCE_DIR}/lib/blas.lib) + set(F2C_LIBRARIES ${CMAKE_SOURCE_DIR}/lib/libf2c.lib) + list(APPEND BLAS_LIBRARIES_TO_LINK ${LAPACK_LIBRARIES} ${BLAS_LIBRARIES} ${F2C_LIBRARIES}) + else() + # Installed by distribution's package manager + find_package(LAPACK REQUIRED QUIET) + list(APPEND BLAS_LIBRARIES_TO_LINK ${LAPACK_LIBRARIES} ${BLAS_LIBRARIES}) + endif() + +elseif(${BLAS_TYPE} STREQUAL "OPENBLAS") + if(NOT DEFINED OPENBLAS_LIBPATH) + if(WIN32) + set(OPENBLAS_LIBPATH "$ENV{HOMEPATH}/AppData/Local/OpenBLAS/openblas.lib") + else() + set(OPENBLAS_LIBPATH "$ENV{HOME}/.local/lib/libopenblas.so") + endif() + message(STATUS "No OPENBLAS_LIBPATH specified, using: ${OPENBLAS_LIBPATH}") + endif() + + if(NOT TARGET OpenBLAS::OpenBLAS) + add_library(OpenBLAS::OpenBLAS SHARED IMPORTED) + if(WIN32) + set_target_properties(OpenBLAS::OpenBLAS PROPERTIES IMPORTED_IMPLIB "${OPENBLAS_LIBPATH}") + else() + set_target_properties(OpenBLAS::OpenBLAS PROPERTIES IMPORTED_LOCATION "${OPENBLAS_LIBPATH}") + endif() + endif() + list(APPEND BLAS_LIBRARIES_TO_LINK OpenBLAS::OpenBLAS) + +elseif(${BLAS_TYPE} STREQUAL "MKL") + set(MKL_THREADING sequential CACHE STRING "MKL threading mode") + set(MKL_INTERFACE lp64 CACHE STRING "MKL interface mode") + + find_package(MKL CONFIG REQUIRED) + list(APPEND BLAS_LIBRARIES_TO_LINK ${MKL_IMPORTED_TARGETS}) + + # Store MKL include dir in a variable to be used after target creation + set(BLAS_INCLUDE_DIRS ${MKL_INCLUDE}) + +else() + message(FATAL_ERROR "Unknown BLAS_TYPE: ${BLAS_TYPE}. Choose LAPACK, OPENBLAS, or MKL") +endif() + +# Export the list of libraries to the main scope +set(BLAS_LIBRARIES_TO_LINK ${BLAS_LIBRARIES_TO_LINK}) \ No newline at end of file diff --git a/src/Diffusion_2D_PPFV.cpp b/src/Diffusion_2D_PPFV.cpp new file mode 100644 index 0000000..c5d3a53 --- /dev/null +++ b/src/Diffusion_2D_PPFV.cpp @@ -0,0 +1,93 @@ +// SPDX-FileCopyrightText: 2026 GFZ Helmholtz Centre for Geosciences +// +// SPDX-License-Identifier: BSD-3-Clause + +#include "Diffusion_2D_PPFV.h" + +#include + +Eigen::VectorXd convertToEigen(const Matrix1D &mat) +{ + Eigen::VectorXd eigen_vec(mat.size_q1); + for (size_t i = 0; i < mat.size_q1; ++i) + { + eigen_vec(i) = mat.plane_array[i]; + } + return eigen_vec; +} + +xt::xtensor convertToXt(const Matrix2D &mat) +{ + xt::xtensor xt_mat = xt::zeros({mat.size_q1, mat.size_q2}); + for (size_t i = 0; i < mat.size_q1; ++i) + { + for (size_t j = 0; j < mat.size_q2; ++j) + { + xt_mat(i, j) = mat.plane_array[i * mat.size_q2 + j]; + } + } + return xt_mat; +} + +Matrix2D convertToMatrix(const xt::xtensor &xt_mat) +{ + size_t size_q1 = xt_mat.shape(0); + size_t size_q2 = xt_mat.shape(1); + + Matrix2D mat(size_q1, size_q2); + + for (size_t i = 0; i < size_q1; ++i) + { + for (size_t j = 0; j < size_q2; ++j) + { + mat.plane_array[i * size_q2 + j] = xt_mat(i, j); + } + } + + return mat; +} + +bool Diffusion_2D_PPFV( + Matrix2D &psd, + const Matrix1D &x, const Matrix1D &y, + int &x_size, int &y_size, + const Matrix1D &x_LBC, const Matrix1D &x_UBC, + const Matrix1D &y_LBC, const Matrix1D &y_UBC, + BoundaryConditionType x_LBC_type, BoundaryConditionType x_UBC_type, + BoundaryConditionType y_LBC_type, BoundaryConditionType y_UBC_type, + const Matrix2D &Dxx, const Matrix2D &Dyy, const Matrix2D &Dxy, + const Matrix2D &G, const Matrix2D &Sources, const Matrix2D &Losses, double total_time, double dt) +{ + + Eigen::VectorXd x_fvm = convertToEigen(x); + Eigen::VectorXd y_fvm = convertToEigen(y); + Eigen::VectorXd xlbc_fvm = convertToEigen(x_LBC); + Eigen::VectorXd xubc_fvm = convertToEigen(x_UBC); + Eigen::VectorXd ylbc_fvm = convertToEigen(y_LBC); + Eigen::VectorXd yubc_fvm = convertToEigen(y_UBC); + + xt::xtensor PSD0_fvm = convertToXt(psd); + xt::xtensor Dxx_fvm = convertToXt(Dxx); + xt::xtensor Dyy_fvm = convertToXt(Dyy); + xt::xtensor Dxy_fvm = convertToXt(Dxy); + xt::xtensor G_fvm = convertToXt(G); + xt::xtensor Sources_fvm = convertToXt(Sources); + xt::xtensor Losses_fvm = convertToXt(Losses); + + + const int num_substeps = std::max(1, int(total_time / dt)); + dt = total_time / num_substeps; + bool success; + + Mesh mesh(x_fvm, y_fvm, dt); + VERB_ppfv2d eq(mesh, PSD0_fvm, Dxx_fvm, Dxy_fvm, Dyy_fvm, xlbc_fvm, xubc_fvm, ylbc_fvm, yubc_fvm, G_fvm, x_LBC_type, x_UBC_type, y_LBC_type, y_UBC_type, Losses_fvm); + Solver solver(mesh, &eq); + + for (int i = 0; i < num_substeps; ++i) + { + solver.update(); + } + PSD0_fvm = solver.PSD(); + psd = convertToMatrix(PSD0_fvm); + return success; +} \ No newline at end of file diff --git a/src/Diffusion_2D_PPFV.h b/src/Diffusion_2D_PPFV.h new file mode 100644 index 0000000..76c76a7 --- /dev/null +++ b/src/Diffusion_2D_PPFV.h @@ -0,0 +1,29 @@ +// SPDX-FileCopyrightText: 2026 GFZ Helmholtz Centre for Geosciences +// +// SPDX-License-Identifier: BSD-3-Clause + +#pragma once + +#include "Matrix.h" +#include "MatrixSolver.h" +#include "BoundaryConditionType.hpp" +#include +#include +#include "ppfv/Mesh.h" +#include "ppfv/Cases/VERB_ppfv2d.h" +#include "ppfv/Solver.h" + +Eigen::VectorXd convertToEigen(const Matrix1D &mat); +xt::xtensor convertToXt(const Matrix2D &mat); + +bool Diffusion_2D_PPFV( + Matrix2D &psd, + const Matrix1D& x, const Matrix1D& y, + int &x_size, int &y_size, + const Matrix1D& x_LBC, const Matrix1D& x_UBC, + const Matrix1D& y_LBC, const Matrix1D& y_UBC, + BoundaryConditionType x_LBC_type, BoundaryConditionType x_UBC_type, + BoundaryConditionType y_LBC_type, BoundaryConditionType y_UBC_type, + const Matrix2D& Dxx, const Matrix2D& Dyy, const Matrix2D& Dxy, + const Matrix2D& G, const Matrix2D& Sources, const Matrix2D& Losses, double total_time ,double dt + ); \ No newline at end of file diff --git a/src/Matrix.cpp b/src/Matrix.cpp index b146d17..7e38153 100644 --- a/src/Matrix.cpp +++ b/src/Matrix.cpp @@ -237,8 +237,8 @@ Matrix1D::Matrix1D( const Matrix1D &M ) { } template -Matrix1D::Matrix1D(size_t size_q1, const T *data) { - this->size_q1 = size_q1; +Matrix1D::Matrix1D(long int size_q1, const T *data) { + this->size_q1 = static_cast(size_q1); this->AllocateMemory(this->size_q1); // fast values copying as a memory range @@ -653,7 +653,7 @@ void Matrix1D::readFromFile(const std::string& filename) { * This is the same as the readFromFile() function although only compatible with .mat files instead of .plt or other text files */ template -void Matrix1D::readFromMatlabFile(const std::string& file , int columnNumber) +void Matrix1D::readFromMatlabFile([[maybe_unused]] const std::string& file, [[maybe_unused]] int columnNumber) { #if (MATLAB_CAPABLE) @@ -1021,7 +1021,7 @@ void Matrix1D::readFromFile(const std::string& filename, const Matrix1D& g * This is the same as the readFromFile() function although only compatible with .mat files instead of .plt or other text files */ template -void Matrix1D::readFromMatlabFile(const std::string& file , const Matrix1D& grid_x) +void Matrix1D::readFromMatlabFile([[maybe_unused]] const std::string& file, [[maybe_unused]] const Matrix1D& grid_x) { #if (MATLAB_CAPABLE) @@ -1459,7 +1459,7 @@ Matrix2D::Matrix2D( const Matrix2D &M ) { } template -Matrix2D::Matrix2D(const size_t *sizes, const T *data) { +Matrix2D::Matrix2D(const long int *sizes, const T *data) { size_q1 = sizes[0]; size_q2 = sizes[1]; AllocateMemory(size_q1, size_q2); @@ -1913,7 +1913,7 @@ void Matrix2D::readFromFile(const std::string& filename, int read_column) { * This is the same as the readFromFile() function although only compatible with .mat files instead of .plt or other text files */ template -void Matrix2D::readFromMatlabFile(const std::string& file , int columnNumber) +void Matrix2D::readFromMatlabFile([[maybe_unused]] const std::string& file, [[maybe_unused]] int columnNumber) { #if (MATLAB_CAPABLE) @@ -2285,7 +2285,7 @@ void Matrix2D::readFromFile(const std::string& filename, const Matrix2D& g * This is the same as the readFromFile() function although only compatible with .mat files instead of .plt or other text files */ template -void Matrix2D::readFromMatlabFile(const std::string& file , const Matrix2D& grid_x, const Matrix2D& grid_y) +void Matrix2D::readFromMatlabFile([[maybe_unused]] const std::string& file, [[maybe_unused]] const Matrix2D& grid_x, [[maybe_unused]] const Matrix2D& grid_y) { #if (MATLAB_CAPABLE) @@ -3131,7 +3131,7 @@ void Matrix3D::readFromFile(const std::string& filename, int read_column) { * This is the same as the readFromFile() function although only compatible with .mat files instead of .plt or other text files */ template -void Matrix3D::readFromMatlabFile(const std::string& file , int columnNumber) +void Matrix3D::readFromMatlabFile([[maybe_unused]]const std::string& file, [[maybe_unused]] int columnNumber) { #if (MATLAB_CAPABLE) @@ -3515,7 +3515,7 @@ void Matrix3D::readFromFile(const std::string& filename, const Matrix3D& g * This is the same as the readFromFile() function although only compatible with .mat files instead of .plt or other text files */ template -void Matrix3D::readFromMatlabFile(const std::string& file , const Matrix3D& grid_x, const Matrix3D& grid_y, const Matrix3D& grid_z) +void Matrix3D::readFromMatlabFile([[maybe_unused]]const std::string& file, [[maybe_unused]] const Matrix3D& grid_x, [[maybe_unused]] const Matrix3D& grid_y, [[maybe_unused]] const Matrix3D& grid_z) { #if (MATLAB_CAPABLE) @@ -4738,7 +4738,7 @@ mxArray* Matrix4D::createStructMatrix(const std::string& info) const * Struct has 7 fields in including - arr time size, size1, size2, size3, size4 */ template -void Matrix4D::writeToMatlabFile(const std::string& filename, const std::string& info) const { +void Matrix4D::writeToMatlabFile([[maybe_unused]] const std::string& filename, [[maybe_unused]] const std::string& info) const { #if (MATLAB_CAPABLE) @@ -4901,7 +4901,7 @@ void Matrix4D::writeToFile(const std::string& filename, const Matrix4D &gr * Uses the createStructMatrix() function to pack all the grid dimensions into seperate variables and then combines these variables into a single matlab structure to save in .mat */ template -void Matrix4D::writeToMatlabFile(const std::string& file, const Matrix4D &grid_w, const Matrix4D &grid_x, const Matrix4D &grid_y, const Matrix4D &grid_z) const { +void Matrix4D::writeToMatlabFile([[maybe_unused]] const std::string& file, [[maybe_unused]] const Matrix4D &grid_w, [[maybe_unused]] const Matrix4D &grid_x, [[maybe_unused]] const Matrix4D &grid_y, [[maybe_unused]] const Matrix4D &grid_z) const { #if (MATLAB_CAPABLE) @@ -5013,7 +5013,7 @@ void Matrix4D::readFromFile(const std::string& filename, int read_column) { * This is the same as the readFromFile() function although only compatible with .mat files instead of .plt or other text files */ template -void Matrix4D::readFromMatlabFile(const std::string& file , int columnNumber) +void Matrix4D::readFromMatlabFile([[maybe_unused]]const std::string& file, [[maybe_unused]] int columnNumber) { #if (MATLAB_CAPABLE) @@ -5464,7 +5464,7 @@ void Matrix4D::readFromFile(const std::string& filename, const Matrix4D& g * This is the same as the readFromFile() function although only compatible with .mat files instead of .plt or other text files */ template -void Matrix4D::readFromMatlabFile(const std::string& file , const Matrix4D& grid_w, const Matrix4D& grid_x, const Matrix4D& grid_y, const Matrix4D& grid_z) +void Matrix4D::readFromMatlabFile([[maybe_unused]] const std::string& file , [[maybe_unused]] const Matrix4D& grid_w, [[maybe_unused]] const Matrix4D& grid_x, [[maybe_unused]] const Matrix4D& grid_y, [[maybe_unused]] const Matrix4D& grid_z) { #if (MATLAB_CAPABLE) diff --git a/src/Matrix.h b/src/Matrix.h index e7293e1..e0e0730 100644 --- a/src/Matrix.h +++ b/src/Matrix.h @@ -67,7 +67,7 @@ class Matrix1D { Matrix1D() = default; Matrix1D( size_t size_q1); Matrix1D( const Matrix1D &M ); - Matrix1D(size_t size_q1, const T* data); // pybind is using this constructor + Matrix1D(long int size_q1, const T* data); // pybind is using this constructor ~Matrix1D(); virtual void AllocateMemory( size_t size_q1 ); @@ -167,7 +167,7 @@ template class Matrix2D { Matrix2D() = default; Matrix2D( const Matrix2D &M ); Matrix2D( size_t size_q1, size_t size_q2 ); - Matrix2D(const size_t* sizes, const T* data); // pybind is using this constructor + Matrix2D(const long int* sizes, const T* data); // pybind is using this constructor ~Matrix2D(); virtual void AllocateMemory(size_t size_q1, size_t size_q2); diff --git a/src/ppfv/Cases/VERB_ppfv2d.cpp b/src/ppfv/Cases/VERB_ppfv2d.cpp new file mode 100644 index 0000000..2d76e7f --- /dev/null +++ b/src/ppfv/Cases/VERB_ppfv2d.cpp @@ -0,0 +1,92 @@ +// SPDX-FileCopyrightText: 2025 Xin Tao +// SPDX-FileCopyrightText: 2026 GFZ Helmholtz Centre Potsdam +// SPDX-FileContributor: Bernhard Haas +// SPDX-FileContributor: Peng Peng +// SPDX-FileContributor: Xin Tao +// +// SPDX-License-Identifier: BSD-3 +// SPDX-License-Identifier: MIT +// +// This code is adapted from the original version published at https://github.com/xtaohub/Sayram-2D/tree/v1.0.0. + +#include "VERB_ppfv2d.h" + +VERB_ppfv2d::VERB_ppfv2d(const Mesh &m_in, + const Xtensor2d &PSD, + const Xtensor2d &Dxx, + const Xtensor2d &Dxy, + const Xtensor2d &Dyy, + const Eigen::VectorXd &x_LBC, + const Eigen::VectorXd &x_UBC, + const Eigen::VectorXd &y_LBC, + const Eigen::VectorXd &y_UBC, + const Xtensor2d &G, + const BoundaryConditionType &x_LBC_type, + const BoundaryConditionType &x_UBC_type, + const BoundaryConditionType &y_LBC_type, + const BoundaryConditionType &y_UBC_type, + const Xtensor2d &Losses) + : Equation(m_in), m(m_in), + PSD(PSD), + x_LBC_(x_LBC), x_UBC_(x_UBC), + y_LBC_(y_LBC), y_UBC_(y_UBC), + x_LBC_type(x_LBC_type), + x_UBC_type(x_UBC_type), + y_LBC_type(y_LBC_type), + y_UBC_type(y_UBC_type) +{ + + Equation::G_ = G; + Equation::Dxx_ = Dxx; + Equation::Dyy_ = Dyy; + Equation::Dxy_ = Dxy; + Equation::Losses_ = Losses; +} + +double VERB_ppfv2d::init_PSD(const Ind &ind) const +{ + return calculate_init_PSD(ind); +} + +BoundaryConditionType VERB_ppfv2d::get_x_LBC_type() const +{ + return x_LBC_type; +} +BoundaryConditionType VERB_ppfv2d::get_x_UBC_type() const +{ + return x_UBC_type; +} +BoundaryConditionType VERB_ppfv2d::get_y_LBC_type() const +{ + return y_LBC_type; +} +BoundaryConditionType VERB_ppfv2d::get_y_UBC_type() const +{ + return y_UBC_type; +} + +void VERB_ppfv2d::apply_bcs(Xtensor2d *vertex_fp) const +{ + Xtensor2d &vertex_f = *vertex_fp; + + // i == 0 and m.nx() boundary condition case + for (std::size_t j = 0; j < m.ny() + 1; ++j) + { + if (x_LBC_type == BoundaryConditionType::ConstantValue) vertex_f(0, j) = x_LBC_(j); + else if (x_LBC_type == BoundaryConditionType::ConstantDerivative) vertex_f(0, j) = PSD(0, j); + + if (x_UBC_type == BoundaryConditionType::ConstantValue) vertex_f(m.nx(), j) = x_UBC_(j); + else if (x_UBC_type == BoundaryConditionType::ConstantDerivative) vertex_f(m.nx(), j) = PSD(m.nx()-1, j); + + } + // j == 0 and j == m.ny() boundary + for (std::size_t i = 0; i < m.nx() + 1; ++i) + { + if (y_LBC_type == BoundaryConditionType::ConstantValue) vertex_f(i, 0) = y_LBC_(i); + else if (y_LBC_type == BoundaryConditionType::ConstantDerivative) vertex_f(i, 0) = PSD(i, 0); + if (y_UBC_type == BoundaryConditionType::ConstantValue) vertex_f(i, m.ny()) = y_UBC_(i); + else if (y_UBC_type == BoundaryConditionType::ConstantDerivative) vertex_f(i, m.ny()) = PSD(i, m.ny()-1); + } + +} + diff --git a/src/ppfv/Cases/VERB_ppfv2d.h b/src/ppfv/Cases/VERB_ppfv2d.h new file mode 100644 index 0000000..6450d86 --- /dev/null +++ b/src/ppfv/Cases/VERB_ppfv2d.h @@ -0,0 +1,67 @@ +/* + * SPDX-FileCopyrightText: 2025 Xin Tao + * SPDX-FileCopyrightText: 2026 GFZ Helmholtz Centre Potsdam + * SPDX-FileContributor: Bernhard Haas + * SPDX-FileContributor: Peng Peng + * SPDX-FileContributor: Xin Tao + * + * SPDX-License-Identifier: BSD-3 + * SPDX-License-Identifier: MIT + * + * This code is adapted from the original version published at https://github.com/xtaohub/Sayram-2D/tree/v1.0.0. + */ + +#pragma once + +#include "../Equation.h" +#include "../../BoundaryConditionType.hpp" + +class VERB_ppfv2d : public Equation +{ +public: + VERB_ppfv2d(const Mesh &m_in, + const Xtensor2d &PSD, + const Xtensor2d &Dxx, + const Xtensor2d &Dxy, + const Xtensor2d &Dyy, + const Eigen::VectorXd &x_LBC, + const Eigen::VectorXd &x_UBC, + const Eigen::VectorXd &y_LBC, + const Eigen::VectorXd &y_UBC, + const Xtensor2d &G, + const BoundaryConditionType &x_LBC_type, + const BoundaryConditionType &x_UBC_type, + const BoundaryConditionType &y_LBC_type, + const BoundaryConditionType &y_UBC_type, + const Xtensor2d &Losses); + + double init_PSD(const Ind &ind) const; + void apply_bcs(Xtensor2d *vertex_fp) const; + BoundaryConditionType x_LBC_type; + BoundaryConditionType x_UBC_type; + BoundaryConditionType y_LBC_type; + BoundaryConditionType y_UBC_type; + BoundaryConditionType get_x_LBC_type() const; + BoundaryConditionType get_x_UBC_type() const; + BoundaryConditionType get_y_LBC_type() const; + BoundaryConditionType get_y_UBC_type() const; + + private: + const Mesh &m; + const Xtensor2d &PSD; + const Eigen::VectorXd &x_LBC_, &x_UBC_, &y_LBC_, &y_UBC_; + double calculate_init_PSD(const Ind &ind) const + { + // convert from cell indices (ind) to vertex indices + + int idx_left = ind.i; + int idx_right = ind.i + 1; + + int idx_bottom = ind.j; + int idx_top = ind.j+1; + + return PSD(ind.i, ind.j); + } + + void init(); +}; diff --git a/src/ppfv/Edge.h b/src/ppfv/Edge.h new file mode 100644 index 0000000..4af6556 --- /dev/null +++ b/src/ppfv/Edge.h @@ -0,0 +1,67 @@ +/* + * SPDX-FileCopyrightText: 2025 Xin Tao + * SPDX-FileCopyrightText: 2026 GFZ Helmholtz Centre Potsdam + * + * SPDX-FileContributor: Bernhard Haas + * SPDX-FileContributor: Peng Peng + * SPDX-FileContributor: Xin Tao + * + * SPDX-License-Identifier: BSD-3 + * SPDX-License-Identifier: MIT + * + * This code is adapted from the original version published at https://github.com/xtaohub/Sayram-2D/tree/v1.0.0. + */ + +#ifndef EDGE_H_ +#define EDGE_H_ + +#include "common.h" + +enum Direction {XPOS, XNEG, YPOS, YNEG}; + +typedef std::array Edge_Vertices; + +class Edge{ + public: + static const int NT = 2; // number of points for each edge. + const Point& v(int i) const { return vs_[i]; } // vertices + int v_size() const { return vs_.size(); } + double length() const { return length_; } + const Vector2& n() const { return n_; } // normal vector, depending on the edge direction + + void set_vs_dir(const Edge_Vertices& vs, Direction dir) { + vs_ = vs; + length_ = (vs[1] - vs[0]).norm(); + + switch (dir) { + case XPOS: // Edge perpendicular to positive x-axis + n_ = {1, 0}; + break; + + case XNEG: + n_ = {-1, 0}; + break; + + case YPOS: // Edge perpendicular to y-axis + n_ = {0, 1}; + break; + + case YNEG: + n_ = {0, -1}; + break; + + default: // incorrect dir + n_ = {0, 0}; + } + + } + + private: + + Edge_Vertices vs_; + double length_; + Vector2 n_; // the unit normal vector +}; + +#endif /* EDGE_H */ + diff --git a/src/ppfv/Equation.h b/src/ppfv/Equation.h new file mode 100644 index 0000000..0beef22 --- /dev/null +++ b/src/ppfv/Equation.h @@ -0,0 +1,70 @@ +/* + * SPDX-FileCopyrightText: 2025 Xin Tao + * SPDX-FileCopyrightText: 2026 GFZ Helmholtz Centre Potsdam + * + * SPDX-FileContributor: Bernhard Haas + * SPDX-FileContributor: Peng Peng + * SPDX-FileContributor: Xin Tao + * + * SPDX-License-Identifier: BSD-3 + * SPDX-License-Identifier: MIT + * + * This code is adapted from the original version published at https://github.com/xtaohub/Sayram-2D/tree/v1.0.0. + */ + +#ifndef EQUATION_H_ +#define EQUATION_H_ + +#include "common.h" +#include "Mesh.h" +#include "../BoundaryConditionType.hpp" + +// +// Definitions of the Equation are given here +// including D and boundary conditions. +// To handle different cases, we Equation as an abstract class. +// + +class Equation{ + public: + Equation(const Mesh& m) { + std::size_t nx = m.nx(); + std::size_t ny = m.ny(); + + G_.resize({nx,ny}); + Losses_.resize({nx,ny}); + + Dxx_.resize({nx,ny}); + Dyy_.resize({nx,ny}); + Dxy_.resize({nx,ny}); + } + + double Losses(const Ind& ind) const {return Losses_(ind.i, ind.j); } + + double G(const Ind& ind) const { return G_(ind.i, ind.j); } + + double Dxx(const Ind& ind) const { return Dxx_(ind.i, ind.j); } + double Dyy(const Ind& ind) const { return Dyy_(ind.i, ind.j); } + double Dxy(const Ind& ind) const { return Dxy_(ind.i, ind.j); } + + + // pure virtual functions to be defined by the CASE of interest. + virtual double init_PSD(const Ind& ind) const = 0; + virtual void apply_bcs(Xtensor2d* vertex_fp) const = 0; + virtual BoundaryConditionType get_x_LBC_type() const = 0; + virtual BoundaryConditionType get_x_UBC_type() const = 0; + virtual BoundaryConditionType get_y_LBC_type() const = 0; + virtual BoundaryConditionType get_y_UBC_type() const = 0; + + protected: + Xtensor2d G_; + Xtensor2d Losses_; + + Xtensor2d Dxx_; + Xtensor2d Dyy_; + Xtensor2d Dxy_; +}; + + +#endif /* EQUATION_H */ + diff --git a/src/ppfv/Mesh.cpp b/src/ppfv/Mesh.cpp new file mode 100644 index 0000000..8726e24 --- /dev/null +++ b/src/ppfv/Mesh.cpp @@ -0,0 +1,154 @@ +// SPDX-FileCopyrightText: 2025 Xin Tao +// SPDX-FileCopyrightText: 2026 GFZ Helmholtz Centre Potsdam +// +// SPDX-FileContributor: Bernhard Haas +// SPDX-FileContributor: Peng Peng +// SPDX-FileContributor: Xin Tao +// +// SPDX-License-Identifier: BSD-3 +// SPDX-License-Identifier: MIT +// +// This code is adapted from the original version published at https://github.com/xtaohub/Sayram-2D/tree/v1.0.0. + +#include "Mesh.h" + +Mesh::Mesh(const Eigen::VectorXd& x, const Eigen::VectorXd& y, + double dt) +: x_centers_(x), y_centers_(y), dt_(dt) { + + nx_ = x_centers_.size(); + ny_ = y_centers_.size(); + + dt_ = dt; + + dx_.resize(nx_); + dy_.resize(ny_); + + calculate_mesh_spacing(); + + xO_ = x_centers_(0) - dx(0)/2.0; + yO_ = y_centers_(0) - dy(0)/2.0; + + x_vertices_ = Eigen::VectorXd::Zero(nx()+1); + y_vertices_ = Eigen::VectorXd::Zero(ny()+1); + + x_vertices_(0) = xO(); + y_vertices_(0) = yO(); + + for (std::size_t i=1; i x_vertices_(1); + bool is_flipped_y = y_vertices_(0) > y_vertices_(1); + + for (std::size_t i=0; i +#include "common.h" +// #include "Parameters.h" +#include "Edge.h" + +struct Ind{ + std::size_t i; + std::size_t j; +}; + +class Mesh { + public: + Mesh(const Eigen::VectorXd& x, const Eigen::VectorXd& y, + double dt); + + const Eigen::VectorXd& x() const { return x_centers_; } + const Eigen::VectorXd& y() const { return y_centers_; } + + double x(int i) const { return x_centers_(i); } + double y(int j) const { return y_centers_(j); } + + double xO() const { return xO_; } + double yO() const { return yO_; } + + std::size_t nx() const { return nx_; } + std::size_t ny() const { return ny_; } + + double dx(int i) const { return dx_(i); } + double dy(int j) const { return dy_(j); } + double dt() const { return dt_; } + + double cell_area(int i, int j) const { return std::abs(dx(i) * dy(j)); } + double cell_area_dt(int i, int j) const { return cell_area(i, j) / dt();} + double cell_area_dt(Ind ind) const { return cell_area(ind.i, ind.j) / dt();} + + void indO(const Point& A, Ind* indp) const { + // calculate the i,j coordinate relative to the Origin + // Note: not the cell index. + // This function is useful to calculate fA and fB from interpolation + + int i = 0; + while (A(0) != x_vertices_(i)) { + i++; + } + + int j = 0; + while (A(1) != y_vertices_(j)) { + j++; + } + + indp->i = i; + indp->j = j; + + } + + std::size_t flatten_index(const Ind& ind) const { // map 2d indices to 1 + return ind.j*nx() + ind.i; + } + + std::size_t nnbrs() const { return 4; } // each cell has 4 nbrs + + // define the neighbor # of 4 adjacent cells + // im -- (i-1, j, k); jp -- (i, j+1, k) + // ip -- (i+1, j, k); jm -- (i, j-1, k) + int inbr_im() const { return 0; } + int inbr_jp() const { return 1; } + int inbr_ip() const { return 2; } + int inbr_jm() const { return 3; } + + int rinbr(int inbr) const { return rinbr_(inbr); } + + void get_nbr_ind(const Ind& ind, int inbr, Ind* nbr_indp) const { + *nbr_indp = nbr_inds(ind.i,ind.j,inbr); + } + + void get_nbr_edge(const Ind& ind, int inbr, Edge* edgep) const { + *edgep = edges(ind.i,ind.j,inbr); + } + + private: + std::size_t nx_; + std::size_t ny_; + Eigen::VectorXd dx_; // Individual dx for each x-cell + Eigen::VectorXd dy_; // Individual dy for each y-cell + double dt_; + + // coordinate origin: corresponds to i-0.5, j-0.5 + double xO_; + double yO_; + + Eigen::VectorXd x_centers_; + Eigen::VectorXd y_centers_; + + Eigen::VectorXd x_vertices_; + Eigen::VectorXd y_vertices_; + + Eigen::Matrix rinbr_; + + xt::xtensor nbr_inds; + xt::xtensor edges; + + void build_connectivity(); + void calculate_mesh_spacing(); +}; + +#endif /* MESH_H */ + diff --git a/src/ppfv/Solver.cpp b/src/ppfv/Solver.cpp new file mode 100644 index 0000000..452ce48 --- /dev/null +++ b/src/ppfv/Solver.cpp @@ -0,0 +1,329 @@ +// SPDX-FileCopyrightText: 2025 Xin Tao +// SPDX-FileCopyrightText: 2026 GFZ Helmholtz Centre Potsdam +// +// SPDX-FileContributor: Bernhard Haas +// SPDX-FileContributor: Peng Peng +// SPDX-FileContributor: Xin Tao +// +// SPDX-License-Identifier: BSD-3 +// SPDX-License-Identifier: MIT +// +// This code is adapted from the original version published at https://github.com/xtaohub/Sayram-2D/tree/v1.0.0. + +#include "Solver.h" +#include +#include // for setw +#include "../BoundaryConditionType.hpp" + +Solver::Solver(const Mesh &m_in, Equation *eqp) + : m(m_in), eq(*eqp) +{ + + size_t nx = m.nx(); + size_t ny = m.ny(); + + M_.resize(nx * ny, nx * ny); + + PSD_.resize({nx, ny}); + + R_.resize(nx * ny); + ftmp_.resize(nx * ny); + + Lambda_.resize({nx, ny}); + vertex_f_.resize({nx + 1, ny + 1}); + + init(); +} + +void Solver::init() +{ + + // iterSolver.setTolerance(1e-8); + + for (size_t i = 0; i < m.nx(); i++) + { + for (size_t j = 0; j < m.ny(); j++) + { + PSD_(i, j) = eq.init_PSD({i, j}); + } + } + + // apply boundary conditions + + update_Lambda(); + update_vertex_f(); + +} + +void Solver::print_matrix(const Xtensor2d &mat, const string &name, const int iter) const +{ + if (iter > 0) + { + cout << name << iter << endl; + } + else + { + cout << name << endl; + } + + const int width = 8; + + for (size_t i = 0; i < m.nx(); i++) + { + for (size_t j = 0; j < m.ny(); j++) + { + cout << setw(width) << std::fixed << setprecision(5) << mat(i, j) << " "; + } + cout << endl; + } +} + +void Solver::print_vertex_matrix(const Xtensor2d &mat, const string &name) const +{ + cout << name << endl; + const int width = 6; + + for (size_t i = 1; i < m.nx(); i++) + { + for (size_t j = 1; j < m.ny(); j++) + { + cout << setw(width) << mat(i, j) << " "; + } + cout << endl; + } +} + +void Solver::update_Lambda() +{ + + for (size_t i = 0; i < m.nx(); i++) + { + for (size_t j = 0; j < m.ny(); j++) + { + Lambda_(i, j) << eq.Dxx({i, j}) * eq.G({i, j}), eq.Dxy({i, j}) * eq.G({i, j}), + eq.Dxy({i, j}) * eq.G({i, j}), eq.Dyy({i, j}) * eq.G({i, j}); + } + } +} + +void Solver::a_sigma_func(const Ind &ind, int inbr, Vector2 *a_sigma_i_p, double *a_sigmap) +{ + + Edge edge; + m.get_nbr_edge(ind, inbr, &edge); + + Point K = {m.x(ind.i), m.y(ind.j)}; + + Point A = edge.v(0); + Point B = edge.v(1); + + Ind indA, indB; + + Vector2 vkb = B - K; // as Vector KB (from K point to B) + Vector2 vka = A - K; + + Vector2 rvkb = {vkb(1), -vkb(0)}; // (x,y) rotated by 90 clockwise, becoming (y, -x) + Vector2 rvka = {vka(1), -vka(0)}; + + double a_sigma_a = edge.length() * (edge.n().transpose() * Lambda(ind) * rvkb)(0) / (vka.transpose() * rvkb)(0); + double a_sigma_b = edge.length() * (edge.n().transpose() * Lambda(ind) * rvka)(0) / (vkb.transpose() * rvka)(0); + + if (std::isnan(a_sigma_a) || std::isnan(a_sigma_b)) { + exit(-1); + } + + (*a_sigma_i_p)(0) = a_sigma_a; + (*a_sigma_i_p)(1) = a_sigma_b; + + m.indO(A, &indA); + m.indO(B, &indB); + + *a_sigmap = a_sigma_a * vertex_f(indA) + a_sigma_b * vertex_f(indB); +} + +void Solver::update_coeff_inner_pair(const Ind &ind, int inbr) +{ // add coefficient from a inner (no-a-boundary) face + + Ind nbr_ind; + + int rinb; + + double a_sigma_K, a_sigma_L; + double mu_K, mu_L; + + double B_sigma, B_sigma_abs, B_sigma_plus, B_sigma_minus; + double A_K, A_L; + + Vector2 a_sigma_i_K; + Vector2 a_sigma_i_L; + + // calculate A_K + a_sigma_func(ind, inbr, &a_sigma_i_K, &a_sigma_K); + + // calculate A_L + m.get_nbr_ind(ind, inbr, &nbr_ind); + rinb = m.rinbr(inbr); + a_sigma_func(nbr_ind, rinb, &a_sigma_i_L, &a_sigma_L); + + calculate_mu(a_sigma_K, a_sigma_L, &mu_K, &mu_L); + + B_sigma = mu_L * a_sigma_L - mu_K * a_sigma_K; + + B_sigma_abs = abs(B_sigma); + B_sigma_plus = (B_sigma_abs + B_sigma) / 2.0; + B_sigma_minus = (B_sigma_abs - B_sigma) / 2.0; + + A_K = mu_K * (a_sigma_i_K[0] + a_sigma_i_K[1]) + B_sigma_plus / (PSD(ind) + gEPS); + A_L = mu_L * (a_sigma_i_L[0] + a_sigma_i_L[1]) + B_sigma_minus / (PSD(nbr_ind) + gEPS); + + if (std::isnan(A_K) || std::isnan(A_L)) { + exit(-1); + } + + size_t K = m.flatten_index(ind), L = m.flatten_index(nbr_ind); + + M_coeffs_.push_back(T(K, K, A_K)); + M_coeffs_.push_back(T(K, L, -A_L)); + M_coeffs_.push_back(T(L, L, A_L)); + M_coeffs_.push_back(T(L, K, -A_K)); +} + +void Solver::update_coeff_dirbc(const Ind &ind, int inbr) +{ // add coefficient from a Dirichlet face + + double a_sigma_K; + + double B_sigma, B_sigma_abs, B_sigma_plus, B_sigma_minus; + double A_K; + Vector2 a_sigma_i_K; + + // calculate A_K + a_sigma_func(ind, inbr, &a_sigma_i_K, &a_sigma_K); + + B_sigma = -a_sigma_K; + B_sigma_abs = abs(B_sigma); + B_sigma_plus = (B_sigma_abs + B_sigma) / 2.0; + B_sigma_minus = (B_sigma_abs - B_sigma) / 2.0; + + A_K = a_sigma_i_K[0] + a_sigma_i_K[1] + B_sigma_plus / (PSD(ind) + gEPS); + + size_t K = m.flatten_index(ind); + M_coeffs_.push_back(T(K, K, A_K)); + R_(K) += B_sigma_minus; +} + +void Solver::assemble() +{ // obtain M and R + + size_t i, j; + + // i direction + // i==0: dirichlet boundary condition + if (eq.get_x_LBC_type() == BoundaryConditionType::ConstantValue) + { + i = 0; + for (j = 0; j < m.ny(); ++j) + update_coeff_dirbc({i, j}, m.inbr_im()); + } + + for (i = 1; i < m.nx(); ++i) + for (j = 0; j < m.ny(); ++j) + update_coeff_inner_pair({i, j}, m.inbr_im()); + + if (eq.get_x_UBC_type() == BoundaryConditionType::ConstantValue) + { + i = m.nx() - 1; + for (j = 0; j < m.ny(); ++j) + update_coeff_dirbc({i, j}, m.inbr_ip()); + } + + if (eq.get_y_LBC_type() == BoundaryConditionType::ConstantValue) + { + j = 0; + for (i = 0; i < m.nx(); ++i) + update_coeff_dirbc({i, j}, m.inbr_jm()); + } + + for (i = 0; i < m.nx(); ++i) + for (j = 1; j < m.ny(); ++j) + update_coeff_inner_pair({i, j}, m.inbr_jm()); + + if (eq.get_y_UBC_type() == BoundaryConditionType::ConstantValue) + { + j = m.ny() - 1; + for (i = 0; i < m.nx(); ++i) + update_coeff_dirbc({i, j}, m.inbr_jp()); + } + + // no updates for BoundaryConditionType::ConstantDerivative condition + + + long K; + double UKK; + + for (size_t i = 0; i < m.nx(); ++i) + { + for (size_t j = 0; j < m.ny(); ++j) + { + K = m.flatten_index({i, j}); + UKK = eq.G({i, j}) * m.cell_area_dt({i,j}); + M_coeffs_.push_back(T(K, K, UKK)); + + M_coeffs_.push_back(T(K, K, -m.cell_area(i,j) * eq.G({i, j}) * eq.Losses({i, j}))); + + // M_(K, K) += UKK; + UKK = eq.G({i, j}) * m.cell_area_dt(i,j); + R_(K) += UKK * PSD_(i, j); + } + } + + M_.setFromTriplets(M_coeffs_.begin(), M_coeffs_.end()); +} + +void Solver::update() +{ + R_.setZero(); + M_coeffs_.clear(); + + assemble(); + + // iterSolver.compute(M_); + // ftmp_ = iterSolver.solve(R_); + + if (!already_factorized_) { + solver.analyzePattern(M_); + solver.factorize(M_); + already_factorized_ = true; + } + + ftmp_ = solver.solve(R_); + + // cout << "iterSolver.info() = " << solver.info() << endl; + + for (size_t i = 0; i < m.nx(); ++i) + { + for (size_t j = 0; j < m.ny(); ++j) + { + PSD_(i, j) = ftmp_(m.flatten_index({i, j})); + } + } + + update_Lambda(); + update_vertex_f(); + + // print_vertex_matrix(vertex_f_, "vertex_f_"); +} + +void Solver::update_vertex_f() +{ + + for (size_t i = 1; i < m.nx(); ++i) + { + for (size_t j = 1; j < m.ny(); ++j) + { + vertex_f_(i, j) = (PSD_(i - 1, j - 1) + PSD_(i - 1, j) + PSD_(i, j - 1) + PSD_(i, j)) / 4.0; + } + } + + eq.apply_bcs(&vertex_f_); +} diff --git a/src/ppfv/Solver.h b/src/ppfv/Solver.h new file mode 100644 index 0000000..c4328e6 --- /dev/null +++ b/src/ppfv/Solver.h @@ -0,0 +1,86 @@ +/* + * SPDX-FileCopyrightText: 2025 Xin Tao + * SPDX-FileCopyrightText: 2026 GFZ Helmholtz Centre Potsdam + * + * SPDX-FileContributor: Bernhard Haas + * SPDX-FileContributor: Peng Peng + * SPDX-FileContributor: Xin Tao + * + * SPDX-License-Identifier: BSD-3 + * SPDX-License-Identifier: MIT + * + * This code is adapted from the original version published at https://github.com/xtaohub/Sayram-2D/tree/v1.0.0. + */ + +#ifndef SOLVER_H +#define SOLVER_H + +#include "common.h" +#include "Mesh.h" +#include "Equation.h" + +class Solver +{ +public: + Solver(const Mesh &m_in, Equation *eqp); + + void update(); + const Xtensor2d &PSD() const { return PSD_; } + const Xtensor2d &PSD_vertices() const { return vertex_f_; } + double PSD(const Ind &ind) const { return PSD_(ind.i, ind.j); } + double PSD(int i, int j) const { return PSD_(i, j); } + void print_matrix(const Xtensor2d &mat, const std::string &name, const int iter) const; + void print_vertex_matrix(const Xtensor2d &mat, const std::string &name) const; + +private: + const Mesh &m; + + Equation &eq; + bool already_factorized_ = false; + + // M f = R + Eigen::SparseLU> solver; + SpMat M_; + std::vector M_coeffs_; + + Xtensor2d PSD_; + Eigen::VectorXd R_; + Eigen::VectorXd ftmp_; // used to store results from the Eigen Solver. + + xt::xtensor Lambda_; // Lambda_(i,j) would be the Lambda Matrix at cell (i,j) + const Eigen::Matrix2d &Lambda(int i, int j) const { return Lambda_(i, j); } + const Eigen::Matrix2d &Lambda(const Ind &ind) const { return Lambda_(ind.i, ind.j); } + void update_Lambda(); + + // + // f at vertices to build a lookup table + // + Xtensor2d vertex_f_; + + double vertex_f(int i, int j) const { return vertex_f_(i, j); } + double vertex_f(const Ind &ind) const { return vertex_f_(ind.i, ind.j); } + void update_vertex_f(); + + void assemble(); + + // to calculate coefficients alpha_sigma_i and a_sigma_i + void a_sigma_func(const Ind &ind, int inbr, Vector2 *a_sigma_i_p, double *a_sigmap); + + // update coefficients in M for cell Ind, and its neighbor. + // note that coefficients for both cell Ind and its neighbor are updated. + void update_coeff_inner_pair(const Ind &ind, int inbr); + + // Here: the inbr neighbor is a Dirichlet boundary. + void update_coeff_dirbc(const Ind &ind, int inbr); + + void calculate_mu(double a_sigma_K, double a_sigma_L, double *mu_Kp, double *mu_Lp) + { + double denom = std::abs(a_sigma_K) + std::abs(a_sigma_L) + 2 * gEPS; + (*mu_Kp) = (std::abs(a_sigma_L) + gEPS) / denom; + (*mu_Lp) = 1 - (*mu_Kp); + } + + void init(); +}; + +#endif /* SOLVER_H */ diff --git a/src/ppfv/common.h b/src/ppfv/common.h new file mode 100644 index 0000000..8b94a6b --- /dev/null +++ b/src/ppfv/common.h @@ -0,0 +1,59 @@ +/* + * SPDX-FileCopyrightText: 2025 Xin Tao + * SPDX-FileCopyrightText: 2026 GFZ Helmholtz Centre Potsdam + * + * SPDX-FileContributor: Bernhard Haas + * SPDX-FileContributor: Peng Peng + * SPDX-FileContributor: Xin Tao + * + * SPDX-License-Identifier: BSD-3 + * SPDX-License-Identifier: MIT + * + * This code is adapted from the original version published at https://github.com/xtaohub/Sayram-2D/tree/v1.0.0. + */ + +#ifndef COMMON_H_ +#define COMMON_H_ + +#define EIGEN_NO_DEBUG + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + + +typedef Eigen::SparseMatrix SpMat; +typedef Eigen::Triplet T; +typedef Eigen::Vector2d Vector2; +typedef Eigen::Vector2d Point; + +typedef xt::xarray Xarray1d; +typedef xt::xtensor Xtensor2d; + +struct Loc{ + int i0; + int j0; + double wi; + double wj; +}; + +// constants +const double gEPS = std::numeric_limits::epsilon(); +const double gPI = 3.141592653589793238462; +const double gD2R = gPI / 180.0; // convert degree to radian +const double gC = 1; +const double gE0 = 0.511875; // MeV +const double gME = gE0 / (gC * gC); + +using namespace std; //bad practice, used here for convenience + +#endif diff --git a/test/test_builds.py b/test/test_builds.py new file mode 100644 index 0000000..36aaea5 --- /dev/null +++ b/test/test_builds.py @@ -0,0 +1,133 @@ +# SPDX-FileCopyrightText: 2026 GFZ Helmholtz Centre for Geosciences +# +# SPDX-License-Identifier: BSD-3-Clause +# +# /// script +# dependencies = [ +# "rich", +# ] +# /// + +import subprocess +import shutil +from pathlib import Path +from rich.console import Console + +console = Console() + +default_target = {} # use defaults + +open_blas_target = {"BLAS_TYPE": "OpenBLAS"} + +data_assimilation_target = { + "BLAS_TYPE": "OpenBLAS", + "DATA_ASSIMILATION": "ON", +} + +data_assimilation_debug_target = { + "BLAS_TYPE": "OpenBLAS", + "DATA_ASSIMILATION": "ON", + "DATA_ASSIMILATION_DEBUG": "ON", +} + +ppfv_target = {"USE_PPFV": "ON"} + +python_bindings_target = {"PYTHON_BINDINGS": "ON"} + +ppfv_python_bindings_target = { + "USE_PPFV": "ON", + "PYTHON_BINDINGS": "ON", +} + +fast_convection_target = {"FAST_CONVECTION": "ON"} + +save_psd_lost_target = {"SAVE_PSD_LOST_CONV": "ON"} + +lu_caching_target = {"LU_CACHING": "ON"} + +all_enabled_target = { + "BLAS_TYPE": "OpenBLAS", + "DATA_ASSIMILATION": "ON", + "DATA_ASSIMILATION_DEBUG": "ON", + "USE_PPFV": "ON", + "PYTHON_BINDINGS": "ON", + "FAST_CONVECTION": "ON", + "SAVE_PSD_LOST_CONV": "ON", + "LU_CACHING": "ON", +} + +test_targets = [ + default_target, open_blas_target, data_assimilation_target, data_assimilation_debug_target, + ppfv_target, python_bindings_target, ppfv_python_bindings_target, fast_convection_target, + save_psd_lost_target, lu_caching_target, all_enabled_target +] + +# uncomment to limit tests +# test_targets = [ppfv_python_bindings_target] + +def run_build(config, index, base_path: Path): + build_dir = base_path / f"build_test_{index}" + log_dir = base_path / "logs" + + if build_dir.exists(): + shutil.rmtree(build_dir) + build_dir.mkdir(parents=True, exist_ok=True) + log_dir.mkdir(exist_ok=True) + + cmake_flags = [f"-D{k}={v}" for k, v in config.items()] + cmake_flags.append("-DCMAKE_BUILD_TYPE=Release") + + # Run Configure + conf = subprocess.run( + ["cmake", "..", "-B", str(build_dir)] + cmake_flags, + capture_output=True, + text=True, + ) + if conf.returncode != 0: + (log_dir / f"config_fail_{index}.log").write_text(conf.stderr) + return False + + # Run Build + build = subprocess.run( + ["cmake", "--build", str(build_dir), "--parallel", "4"], + capture_output=True, + text=True, + ) + if build.returncode != 0: + (log_dir / f"build_fail_{index}.log").write_text(build.stderr) + return False + + return True + + +root_path = Path(__file__).parent.resolve() + +console.print("[bold red]Cleaning up old build directories and logs...[/bold red]") +# Find all directories matching build_test_* and remove them +for old_dir in root_path.glob("build_test_*"): + if old_dir.is_dir(): + shutil.rmtree(old_dir) + +# Remove old logs +old_logs = root_path / "logs" +if old_logs.exists(): + shutil.rmtree(old_logs) + +console.print( + f"[bold blue]Starting Matrix Build for {len(test_targets)} configurations...[/bold blue]\n" +) + +with console.status("", spinner="bouncingBar") as status: + for i, config in enumerate(test_targets): + # Format the current options for the display + if len(config) == 0: + opt_str = "default" + else: + opt_str = ", ".join([f"[cyan]{k}[/]={v}" for k, v in config.items()]) + status.update(f"[bold yellow]Config {i}:[/bold yellow] {opt_str}") + + success = run_build(config, i, root_path) + + # Log result above the sticky status bar + status_icon = "[green]PASS[/green]" if success else "[red]FAIL[/red]" + console.log(f"{status_icon} Config {i}: {opt_str}")