diff --git a/CMakeLists.txt b/CMakeLists.txt index 93e29a96ff0090554ba3021094a2b61e3a660dc8..6b4f0a54d493750762038fb9ca4230b4293f5655 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -2,6 +2,7 @@ cmake_minimum_required (VERSION 3.10) # CMake utils list(APPEND CMAKE_MODULE_PATH "${CMAKE_CURRENT_SOURCE_DIR}/cmake") +list(APPEND CMAKE_MODULE_PATH "${CMAKE_CURRENT_SOURCE_DIR}/packages/cmake-modules") # Forbids in-source builds include(CheckNotInSources) @@ -113,13 +114,12 @@ if (NOT PUGS_ENABLE_MPI MATCHES "^(AUTO|ON|OFF)$") message(FATAL_ERROR "PUGS_ENABLE_MPI='${PUGS_ENABLE_MPI}'. Must be set to one of AUTO, ON or OFF") endif() -# checks for MPI +# Check for MPI if (PUGS_ENABLE_MPI MATCHES "^(AUTO|ON)$") set(MPI_DETERMINE_LIBRARY_VERSION TRUE) find_package(MPI) endif() - #------------------------------------------------------ # Search for ParMETIS @@ -137,6 +137,16 @@ if(${MPI_FOUND}) endif() endif() +#------------------------------------------------------ +# Check for PETSc +find_package(PETSc) +set(PUGS_HAS_PETSC ${PETSC_FOUND}) + +if (${PETSC_FOUND}) + include_directories(SYSTEM ${PETSC_INCLUDES}) +endif() + + # ----------------------------------------------------- if (${MPI_FOUND}) @@ -414,15 +424,18 @@ add_executable( target_link_libraries( pugs PugsMesh + PugsAlgebra PugsUtils PugsLanguage PugsLanguageAST PugsLanguageModules PugsLanguageAlgorithms PugsMesh + PugsAlgebra PugsUtils PugsLanguageUtils kokkos + ${PETSC_LIBRARIES} ${PARMETIS_LIBRARIES} ${MPI_CXX_LINK_FLAGS} ${MPI_CXX_LIBRARIES} ${KOKKOS_CXX_FLAGS} @@ -450,37 +463,43 @@ endif() message("") message("====== pugs build options ======") -message(STATUS "version: ${PUGS_VERSION}") -message(STATUS "build type: ${CMAKE_BUILD_TYPE}") -message(STATUS "compiler: ${CMAKE_CXX_COMPILER_ID} ${CMAKE_CXX_COMPILER_VERSION}") -message(STATUS "kokkos devices: ${PUGS_BUILD_KOKKOS_DEVICES}") +message(" version: ${PUGS_VERSION}") +message(" build type: ${CMAKE_BUILD_TYPE}") +message(" compiler: ${CMAKE_CXX_COMPILER_ID} ${CMAKE_CXX_COMPILER_VERSION}") +message(" kokkos devices: ${PUGS_BUILD_KOKKOS_DEVICES}") if (MPI_FOUND) - message(STATUS "MPI: ${MPI_CXX_LIBRARY_VERSION_STRING}") + message(" MPI: ${MPI_CXX_LIBRARY_VERSION_STRING}") else() if(NOT PARMETIS_LIBRARIES) - message(STATUS "MPI: deactivated: ParMETIS cannot be found!") + message(" MPI: deactivated: ParMETIS cannot be found!") else() - message(STATUS "MPI: not found!") + message(" MPI: not found!") endif() endif() +if (PETSC_FOUND) + message(" PETSc: ${PETSC_VERSION}") +else() + message(" PETSc: not found!") +endif() + if(CLANG_FORMAT) - message(STATUS "clang-format: ${CLANG_FORMAT}") + message(" clang-format: ${CLANG_FORMAT}") else() - message(STATUS "clang-format: not found!") + message(" clang-format: not found!") endif() if(CLAZY_STANDALONE) - message(STATUS "clazy-standalone: ${CLAZY_STANDALONE}") + message(" clazy-standalone: ${CLAZY_STANDALONE}") else() - message(STATUS "clazy-standalone: no found!") + message(" clazy-standalone: no found!") endif() if (DOXYGEN_FOUND) - message(STATUS "doxygen: ${DOXYGEN_EXECUTABLE}") + message(" doxygen: ${DOXYGEN_EXECUTABLE}") else() - message(STATUS "doxygen: no found!") + message(" doxygen: no found!") endif() message("================================") diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 0d6e8cbb996d4f6c18c9098bb9af7051c11b0dbf..1b21f79a7545c442be5448792132acbdf840947a 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -10,7 +10,7 @@ add_subdirectory(utils) add_subdirectory(language) # Pugs algebra -#add_subdirectory(algebra) +add_subdirectory(algebra) # Pugs mesh add_subdirectory(mesh) diff --git a/src/algebra/CMakeLists.txt b/src/algebra/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..bbd2d01bf417baaacb13c6c7362a825337a180be --- /dev/null +++ b/src/algebra/CMakeLists.txt @@ -0,0 +1,6 @@ +# ------------------- Source files -------------------- + +add_library( + PugsAlgebra + LinearSystemSolver.cpp + PETScWrapper.cpp) diff --git a/src/algebra/LinearSystemSolver.cpp b/src/algebra/LinearSystemSolver.cpp new file mode 100644 index 0000000000000000000000000000000000000000..488d6763d1838cc4d5edc596589310dc2ef66038 --- /dev/null +++ b/src/algebra/LinearSystemSolver.cpp @@ -0,0 +1,331 @@ +#include <algebra/LinearSystemSolver.hpp> +#include <utils/pugs_config.hpp> + +#ifdef PUGS_HAS_PETSC +#include <petsc.h> +#endif // PUGS_HAS_PETSC + +struct LinearSystemSolver::Internals +{ + static bool + hasLibrary(LSLibrary library) + { + switch (library) { + case LSLibrary::builtin: { + return true; + } + case LSLibrary::petsc: { +#ifdef PUGS_HAS_PETSC + return true; +#else + return false; +#endif + } + default: { + throw UnexpectedError("Linear system library (" + name(library) + ") was not set!"); + } + } + } + + static void + checkHasLibrary(LSLibrary library) + { + if (not hasLibrary(library)) { + throw NormalError(name(library) + " is not linked to pugs. Cannot use it!"); + } + } + + static std::string + name(LSLibrary library) + { + switch (library) { + case LSLibrary::builtin: { + return "builtin"; + } + case LSLibrary::petsc: { + return "PETSc"; + } + default: { + throw UnexpectedError("Linear system library name is not defined!"); + } + } + } + + static std::string + name(LSMethod method) + { + switch (method) { + case LSMethod::cg: { + return "conjugate gradient"; + } + case LSMethod::bicgstab: { + return "bicg-stab"; + } + case LSMethod::gmres: { + return "gmres"; + } + default: { + throw UnexpectedError("Linear system method name is not defined!"); + } + } + } + + static std::string + name(LSPrecond precond) + { + switch (precond) { + case LSPrecond::none: { + return "none"; + } + case LSPrecond::diagonal: { + return "Jacobi (diagonal)"; + } + case LSPrecond::incomplete_choleski: { + return "incomplete Choleski"; + } + case LSPrecond::incomplete_LU: { + return "incomplete LU"; + } + default: { + throw UnexpectedError("Linear system preconditioner name is not defined!"); + } + } + } + + static void + checkBuiltinMethod(LSMethod method) + { + switch (method) { + case LSMethod::cg: + case LSMethod::bicgstab: { + break; + } + default: { + throw NormalError(name(method) + " is not a builtin linear solver!"); + } + } + } + + static void + checkPETScMethod(LSMethod method) + { + switch (method) { + case LSMethod::cg: + case LSMethod::bicgstab: + case LSMethod::gmres: { + break; + } + default: { + throw NormalError(name(method) + " is not a builtin linear solver!"); + } + } + } + + static void + checkBuiltinPrecond(LSPrecond precond) + { + switch (precond) { + case LSPrecond::none: { + break; + } + default: { + throw NormalError(name(precond) + " is not a builtin preconditioner!"); + } + } + } + + static void + checkPETScPrecond(LSPrecond precond) + { + switch (precond) { + case LSPrecond::none: + case LSPrecond::diagonal: + case LSPrecond::incomplete_choleski: + case LSPrecond::incomplete_LU: { + break; + } + default: { + throw NormalError(name(precond) + " is not a PETSc preconditioner!"); + } + } + } + + static void + checkOptions(LSLibrary library, LSMethod method, LSPrecond precond) + { + switch (library) { + case LSLibrary::builtin: { + checkBuiltinMethod(method); + checkBuiltinPrecond(precond); + break; + } + case LSLibrary::petsc: { + checkPETScMethod(method); + checkPETScPrecond(precond); + break; + } + default: { + throw UnexpectedError("Undefined options compatibility for this library (" + name(library) + ")!"); + } + } + } + +#ifdef PUGS_HAS_PETSC + + static int + petscMonitor(KSP, int it, double norm, void*) + { + std::cout << "iteration: " << it << " residual norm = " << norm << '\n'; + return 0; + } + + static void + petscSolveLocalSystem(const SparseMatrixDescriptor<double, size_t>& A, + Vector<double>& x, + const Vector<double>& b, + LSMethod method, + LSPrecond precond, + double epsilon, + size_t maximum_iteration) + { + Assert(x.size() == b.size() and x.size() == A.numberOfRows()); + + Vec petscB; + VecCreateMPIWithArray(PETSC_COMM_WORLD, 1, b.size(), b.size(), &b[0], &petscB); + Vec petscX; + VecCreateMPIWithArray(PETSC_COMM_WORLD, 1, x.size(), x.size(), &x[0], &petscX); + + Mat petscMat; + MatCreate(PETSC_COMM_WORLD, &petscMat); + MatSetSizes(petscMat, PETSC_DECIDE, PETSC_DECIDE, x.size(), x.size()); + + MatSetType(petscMat, MATAIJ); + MatSetUp(petscMat); + + for (size_t i = 0; i < A.numberOfRows(); ++i) { + const auto& row = A.row(i); + PetscInt I = i; + for (auto [j, value] : row.columnIdToValueMap()) { + PetscInt J = j; + MatSetValues(petscMat, 1, &I, 1, &J, &value, INSERT_VALUES); + } + } + + MatAssemblyBegin(petscMat, MAT_FINAL_ASSEMBLY); + MatAssemblyEnd(petscMat, MAT_FINAL_ASSEMBLY); + + PetscPrintf(PETSC_COMM_WORLD, "Matrix _before_ MatScaleUserImpl() operation\n"); + MatView(petscMat, PETSC_VIEWER_STDOUT_SELF); + + KSP ksp; + + KSPCreate(PETSC_COMM_WORLD, &ksp); + KSPSetTolerances(ksp, epsilon, 1E-100, 1E5, maximum_iteration); + + KSPSetOperators(ksp, petscMat, petscMat); + + switch (method) { + case LSMethod::bicgstab: { + KSPSetType(ksp, KSPBCGS); + break; + } + case LSMethod::cg: { + KSPSetType(ksp, KSPCG); + break; + } + case LSMethod::gmres: { + KSPSetType(ksp, KSPGMRES); + break; + } + default: { + throw UnexpectedError("unexpected method: " + name(method)); + } + } + + PC pc; + KSPGetPC(ksp, &pc); + + switch (precond) { + case LSPrecond::diagonal: { + PCSetType(pc, PCJACOBI); + break; + } + case LSPrecond::incomplete_LU: { + PCSetType(pc, PCILU); + break; + } + case LSPrecond::incomplete_choleski: { + PCSetType(pc, PCICC); + break; + } + case LSPrecond::none: { + PCSetType(pc, PCNONE); + break; + } + default: { + throw UnexpectedError("unexpected preconditioner: " + name(precond)); + } + } + + KSPMonitorSet(ksp, petscMonitor, 0, 0); + + KSPSolve(ksp, petscB, petscX); + + VecView(petscX, PETSC_VIEWER_STDOUT_SELF); + + MatDestroy(&petscMat); + } + +#else + + static void + petscSolveLocalSystem(const SparseMatrixDescriptor<double, size_t>&, + Vector<double>&, + const Vector<double>&, + LSMethod, + LSPrecond, + double, + size_t) + { + checkHasLibrary(LSLibrary::petsc); + } + +#endif // PUGS_HAS_PETSC +}; + +bool +LinearSystemSolver::hasLibrary(LSLibrary library) const +{ + return Internals::hasLibrary(library); +} + +LinearSystemSolver::LinearSystemSolver(LSLibrary library, + LSMethod method, + LSPrecond precond, + double epsilon, + size_t maximum_iteration) + : m_library{library}, m_method{method}, m_precond{precond}, m_epsilon{epsilon}, m_maximum_iteration{maximum_iteration} +{ + Internals::checkHasLibrary(library); + Internals::checkOptions(library, method, precond); +} + +void +LinearSystemSolver::solveLocalSystem(const SparseMatrixDescriptor<double, size_t>& A, + Vector<double>& x, + const Vector<double>& b) +{ + switch (m_library) { + case LSLibrary::builtin: { + throw NotImplementedError("builtin solvers are not plugged yet!"); + break; + } + case LSLibrary::petsc: { + Internals::petscSolveLocalSystem(A, x, b, m_method, m_precond, m_epsilon, m_maximum_iteration); + break; + } + default: { + throw UnexpectedError(Internals::name(m_library) + " cannot solve local systems for sparse matrices"); + } + } +} diff --git a/src/algebra/LinearSystemSolver.hpp b/src/algebra/LinearSystemSolver.hpp new file mode 100644 index 0000000000000000000000000000000000000000..f689f96ddba12ca8f7f4e4c19bbd78621dd1981f --- /dev/null +++ b/src/algebra/LinearSystemSolver.hpp @@ -0,0 +1,79 @@ +#ifndef LINEAR_SYSTEM_SOLVER_HPP +#define LINEAR_SYSTEM_SOLVER_HPP + +#include <algebra/SparseMatrixDescriptor.hpp> +#include <algebra/TinyMatrix.hpp> +#include <algebra/TinyVector.hpp> +#include <algebra/Vector.hpp> + +#include <utils/Exceptions.hpp> + +enum class LSLibrary +{ + builtin, + petsc +}; + +enum class LSMethod +{ + cg, + bicgstab, + gmres +}; + +enum class LSPrecond +{ + none, + diagonal, + incomplete_choleski, + incomplete_LU +}; + +class LinearSystemSolver +{ + private: + struct Internals; + + const LSLibrary m_library; + const LSMethod m_method; + const LSPrecond m_precond; + + const double m_epsilon; + const size_t m_maximum_iteration; + + bool m_is_verbose = false; + + void _solveLocalDense(size_t N, const double* A, double* x, const double* b); + + public: + bool + isVerbose() const + { + return m_is_verbose; + } + + void + setVerbose(bool verbose) + { + m_is_verbose = verbose; + } + + bool hasLibrary(LSLibrary library) const; + + void solveLocalSystem(const SparseMatrixDescriptor<double, size_t>& A, Vector<double>& x, const Vector<double>& b); + + template <size_t N> + void + solveLocalSystem(const TinyMatrix<N>& A, TinyVector<N>& x, const TinyVector<N>& b) + { + this->_solveLocalDense(N, &A(0, 0), &x[0], &b[0]); + } + + LinearSystemSolver(LSLibrary library, + LSMethod method, + LSPrecond precond, + double epsilon = 1e-6, + size_t maximum_iteration = 100); +}; + +#endif // LINEAR_SYSTEM_SOLVER_HPP diff --git a/src/algebra/PETScWrapper.cpp b/src/algebra/PETScWrapper.cpp new file mode 100644 index 0000000000000000000000000000000000000000..76e56692a67e86bb61f6333c09680a30bc9e07ae --- /dev/null +++ b/src/algebra/PETScWrapper.cpp @@ -0,0 +1,26 @@ +#include <algebra/PETScWrapper.hpp> + +#include <utils/pugs_config.hpp> + +#ifdef PUGS_HAS_PETSC +#include <petsc.h> +#endif // PUGS_HAS_PETSC + +namespace PETScWrapper +{ +void +initialize([[maybe_unused]] int& argc, [[maybe_unused]] char* argv[]) +{ +#ifdef PUGS_HAS_PETSC + PetscInitialize(&argc, &argv, 0, 0); +#endif // PUGS_HAS_PETSC +} + +void +finalize() +{ +#ifdef PUGS_HAS_PETSC + PetscFinalize(); +#endif // PUGS_HAS_PETSC +} +} // namespace PETScWrapper diff --git a/src/algebra/PETScWrapper.hpp b/src/algebra/PETScWrapper.hpp new file mode 100644 index 0000000000000000000000000000000000000000..2f40b37f8c4562b7abeccffeb6e8cfb27688c7c4 --- /dev/null +++ b/src/algebra/PETScWrapper.hpp @@ -0,0 +1,10 @@ +#ifndef PETSC_WRAPPER_HPP +#define PETSC_WRAPPER_HPP + +namespace PETScWrapper +{ +void initialize(int& argc, char* argv[]); +void finalize(); +} // namespace PETScWrapper + +#endif // PETSC_WRAPPER_HPP diff --git a/src/algebra/SparseMatrixDescriptor.hpp b/src/algebra/SparseMatrixDescriptor.hpp index 100000366d3bd0e5a072084864e82cdb4b73ba29..09f13085b4881d2d602d41a42e0d8d0918cca633 100644 --- a/src/algebra/SparseMatrixDescriptor.hpp +++ b/src/algebra/SparseMatrixDescriptor.hpp @@ -22,6 +22,12 @@ class SparseMatrixDescriptor std::map<IndexType, DataType> m_id_value_map; public: + const auto& + columnIdToValueMap() const + { + return m_id_value_map; + } + IndexType numberOfValues() const noexcept { diff --git a/src/utils/BuildInfo.cpp b/src/utils/BuildInfo.cpp index 13b8633c2664add44d91b11b0a6d55afe79d2955..7162aeccf69a7b5e6420e23ae7dca6cea6ced9b1 100644 --- a/src/utils/BuildInfo.cpp +++ b/src/utils/BuildInfo.cpp @@ -8,6 +8,10 @@ #include <mpi.h> #endif // PUGS_HAS_MPI +#ifdef PUGS_HAS_PETSC +#include <petsc.h> +#endif // PUGS_HAS_PETSC + std::string BuildInfo::type() { @@ -42,3 +46,14 @@ BuildInfo::mpiLibrary() return "none"; #endif // PUGS_HAS_MPI } + +std::string +BuildInfo::petscLibrary() +{ +#ifdef PUGS_HAS_PETSC + return std::to_string(PETSC_VERSION_MAJOR) + "." + std::to_string(PETSC_VERSION_MINOR) + "." + + std::to_string(PETSC_VERSION_SUBMINOR); +#else + return "none"; +#endif // PUGS_HAS_PETSC +} diff --git a/src/utils/BuildInfo.hpp b/src/utils/BuildInfo.hpp index 81d2b4960085fbeffff667fb92e8a4069848fb8d..67134a782f595ccf68dc89cd29a6f273c7ce6140 100644 --- a/src/utils/BuildInfo.hpp +++ b/src/utils/BuildInfo.hpp @@ -9,6 +9,7 @@ struct BuildInfo static std::string compiler(); static std::string kokkosDevices(); static std::string mpiLibrary(); + static std::string petscLibrary(); }; #endif // BUILD_INFO_HPP diff --git a/src/utils/PugsUtils.cpp b/src/utils/PugsUtils.cpp index 7ea50b0ed95845d2cee2633477e33080453f0f30..8f388c5b492dcd8065ed77740d322556bdc76db3 100644 --- a/src/utils/PugsUtils.cpp +++ b/src/utils/PugsUtils.cpp @@ -9,6 +9,8 @@ #include <utils/FPEManager.hpp> #include <utils/SignalManager.hpp> +#include <algebra/PETScWrapper.hpp> + #include <rang.hpp> #include <Kokkos_Core.hpp> @@ -41,7 +43,8 @@ initialize(int& argc, char* argv[]) std::cout << "type: " << rang::style::bold << BuildInfo::type() << rang::style::reset << '\n'; std::cout << "compiler: " << rang::style::bold << BuildInfo::compiler() << rang::style::reset << '\n'; std::cout << "kokkos: " << rang::style::bold << BuildInfo::kokkosDevices() << rang::style::reset << '\n'; - std::cout << "mpi: " << rang::style::bold << BuildInfo::mpiLibrary() << rang::style::reset << '\n'; + std::cout << "MPI: " << rang::style::bold << BuildInfo::mpiLibrary() << rang::style::reset << '\n'; + std::cout << "PETSc: " << rang::style::bold << BuildInfo::petscLibrary() << rang::style::reset << '\n'; std::cout << "-------------------------------------------------------\n"; std::string filename; @@ -81,6 +84,8 @@ initialize(int& argc, char* argv[]) SignalManager::init(enable_signals); } + PETScWrapper::initialize(argc, argv); + Kokkos::initialize(argc, argv); std::cout << "-------------------- " << rang::fg::green << "exec info" << rang::fg::reset << " ------------------------" << '\n'; @@ -97,5 +102,6 @@ void finalize() { Kokkos::finalize(); + PETScWrapper::finalize(); parallel::Messenger::destroy(); } diff --git a/src/utils/pugs_config.hpp.in b/src/utils/pugs_config.hpp.in index 45878e64327797cd6549b18f4d6922f8cde8a29d..326884f97f401f0d70688fa4ad7259bdb64ab967 100644 --- a/src/utils/pugs_config.hpp.in +++ b/src/utils/pugs_config.hpp.in @@ -3,6 +3,7 @@ #cmakedefine PUGS_HAS_FENV_H #cmakedefine PUGS_HAS_MPI +#cmakedefine PUGS_HAS_PETSC #cmakedefine SYSTEM_IS_LINUX #cmakedefine SYSTEM_IS_DARWIN @@ -10,4 +11,4 @@ #define PUGS_BUILD_TYPE "@CMAKE_BUILD_TYPE@" -#endif // PUGS_CONFIG_HPP +#endif // PUGS_CONFIG_HPP diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 7cd1ed89023e619af478625b58441002ecae0af2..8b2a7c5516394444debece5e32413dde7804c56d 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -93,6 +93,7 @@ target_link_libraries (unit_tests PugsLanguageAlgorithms PugsLanguageUtils PugsMesh + PugsAlgebra PugsUtils kokkos ${PARMETIS_LIBRARIES} @@ -102,6 +103,7 @@ target_link_libraries (unit_tests target_link_libraries (mpi_unit_tests PugsUtils + PugsAlgebra PugsMesh kokkos ${PARMETIS_LIBRARIES}