diff --git a/CMakeLists.txt b/CMakeLists.txt index 90af1219b..30cfe278f 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -35,6 +35,7 @@ option(USE_DAMARIS_LIB "Use the Damaris library for asynchronous I/O?" OFF) option(USE_GPU_BRIDGE "Enable the GPU bridge (GPU/AMGCL solvers)" ON) option(USE_TRACY_PROFILER "Enable tracy profiling" OFF) option(CONVERT_CUDA_TO_HIP "Convert CUDA code to HIP (to run on AMD cards)" OFF) +option(USE_HYPRE "Use the Hypre library for linear solvers?" OFF) set(OPM_COMPILE_COMPONENTS "2;3;4;5;6;7" CACHE STRING "The components to compile support for") # Ensure OPM_COMPILE_COMPONENTS is a list of at least one element @@ -648,8 +649,13 @@ add_custom_target(extra_test ${CMAKE_CTEST_COMMAND} -C ExtraTests) if(CUDA_FOUND) if (NOT USE_HIP) - target_link_libraries( opmsimulators PUBLIC ${CUDA_cusparse_LIBRARY} ) - target_link_libraries( opmsimulators PUBLIC ${CUDA_cublas_LIBRARY} ) + target_link_libraries(opmsimulators + PUBLIC + ${CUDA_cusparse_LIBRARY} + ${CUDA_cublas_LIBRARY} + ${CUDA_nvptxcompiler_static_LIBRARY} + ) + foreach(tgt test_gpu_safe_call test_cuda_check_last_error test_GpuVector) target_link_libraries(${tgt} CUDA::cudart) endforeach() @@ -715,3 +721,19 @@ endif() install(DIRECTORY doc/man1 DESTINATION ${CMAKE_INSTALL_MANDIR} FILES_MATCHING PATTERN "*.1") + +if(USE_HYPRE) + find_package(HYPRE) + if(HYPRE_FOUND) + set(HAVE_HYPRE 1) + target_link_libraries(opmsimulators PUBLIC HYPRE::HYPRE) + if (HYPRE_USING_CUDA) + set_tests_properties(HyprePreconditionerGPU PROPERTIES LABELS gpu_cuda) + elseif (HYPRE_USING_HIP) + set_tests_properties(HyprePreconditionerGPU PROPERTIES LABELS gpu_hip) + endif() + else() + message(WARNING "Hypre requested but not found. Continuing without Hypre support.") + set(USE_HYPRE OFF) + endif() +endif() diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 353567f16..e67d7667e 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -492,6 +492,13 @@ if(HDF5_FOUND) list(APPEND TEST_SOURCE_FILES tests/test_HDF5Serializer.cpp) endif() +if(HYPRE_FOUND) + list(APPEND TEST_SOURCE_FILES tests/test_HyprePreconditionerCPU.cpp) + if(HYPRE_USING_CUDA OR HYPRE_USING_HIP) + list(APPEND TEST_SOURCE_FILES tests/test_HyprePreconditionerGPU.cpp) + endif() +endif() + list (APPEND TEST_DATA_FILES tests/equil_base.DATA tests/equil_capillary.DATA @@ -1156,3 +1163,9 @@ if(dune-alugrid_FOUND) examples/fracture_discretefracture.cpp ) endif() + +if(HYPRE_FOUND) + list(APPEND PUBLIC_HEADER_FILES + opm/simulators/linalg/HyprePreconditioner.hpp + ) +endif() diff --git a/opm-simulators-prereqs.cmake b/opm-simulators-prereqs.cmake index f2c6fb9be..2af64605f 100644 --- a/opm-simulators-prereqs.cmake +++ b/opm-simulators-prereqs.cmake @@ -23,6 +23,7 @@ set (opm-simulators_CONFIG_VAR HAVE_SUITESPARSE_UMFPACK HAVE_DAMARIS HAVE_HDF5 + HAVE_HYPRE USE_HIP USE_TRACY FLOW_INSTANTIATE_FLOAT @@ -56,6 +57,8 @@ set (opm-simulators_DEPS # packages from ROCm framework "rocblas" "rocsparse" + # HYPRE solver library + "HYPRE" # OPM dependency "opm-common REQUIRED" "opm-grid REQUIRED" diff --git a/opm/simulators/flow/Main.cpp b/opm/simulators/flow/Main.cpp index 1d40593ed..91864fba9 100644 --- a/opm/simulators/flow/Main.cpp +++ b/opm/simulators/flow/Main.cpp @@ -35,6 +35,11 @@ #include #endif +#if HAVE_HYPRE +#include +#include +#endif + namespace Opm { Main::Main(int argc, char** argv, bool ownMPI) @@ -88,6 +93,10 @@ Main::~Main() } #endif // HAVE_MPI +#if HAVE_HYPRE + HYPRE_Finalize(); +#endif + if (ownMPI_) { FlowGenericVanguard::setCommunication(nullptr); } @@ -163,6 +172,14 @@ void Main::initMPI() #endif #endif // HAVE_MPI + +#if HAVE_HYPRE +#if HYPRE_RELEASE_NUMBER >= 22900 + HYPRE_Initialize(); +#else + HYPRE_Init(); +#endif +#endif } void Main::handleVersionCmdLine_(int argc, char** argv, diff --git a/opm/simulators/linalg/HyprePreconditioner.hpp b/opm/simulators/linalg/HyprePreconditioner.hpp new file mode 100644 index 000000000..067fe6e66 --- /dev/null +++ b/opm/simulators/linalg/HyprePreconditioner.hpp @@ -0,0 +1,407 @@ +/* + Copyright 2024 SINTEF AS + Copyright 2024 Equinor ASA + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#ifndef OPM_HYPRE_PRECONDITIONER_HEADER_INCLUDED +#define OPM_HYPRE_PRECONDITIONER_HEADER_INCLUDED + +#include +#include +#include +#include + +#include +#include + +#include +#include +#include +#include <_hypre_utilities.h> + +#include +#include + +namespace Hypre { + +/** + * @brief Wrapper for Hypre's BoomerAMG preconditioner. + * + * This class provides an interface to the BoomerAMG preconditioner from the Hypre library. + * It is designed to work with matrices, update vectors, and defect vectors specified by the template parameters. + * + * @tparam M The matrix type the preconditioner is for. + * @tparam X The type of the update vector. + * @tparam Y The type of the defect vector. + */ +template +class HyprePreconditioner : public Dune::PreconditionerWithUpdate { +public: + + /** + * @brief Constructor for the HyprePreconditioner class. + * + * Initializes the preconditioner with the given matrix and property tree. + * + * @param A The matrix for which the preconditioner is constructed. + * @param prm The property tree containing configuration parameters. + */ + HyprePreconditioner (const M& A, const Opm::PropertyTree prm) + : A_(A) + { + OPM_TIMEBLOCK(prec_construct); + + int size; + MPI_Comm_size(MPI_COMM_WORLD, &size); + if (size > 1) { + OPM_THROW(std::runtime_error, "HyprePreconditioner is currently only implemented for sequential runs"); + } + + use_gpu_ = prm.get("use_gpu", false); + + // Set memory location and execution policy +#if HYPRE_USING_CUDA || HYPRE_USING_HIP + if (use_gpu_) { + HYPRE_SetMemoryLocation(HYPRE_MEMORY_DEVICE); + HYPRE_SetExecutionPolicy(HYPRE_EXEC_DEVICE); + // use hypre's SpGEMM instead of vendor implementation + HYPRE_SetSpGemmUseVendor(false); + // use cuRand for PMIS + HYPRE_SetUseGpuRand(1); + HYPRE_DeviceInitialize(); + HYPRE_PrintDeviceInfo(); + } + else +#endif + { + HYPRE_SetMemoryLocation(HYPRE_MEMORY_HOST); + HYPRE_SetExecutionPolicy(HYPRE_EXEC_HOST); + } + + // Create the solver (BoomerAMG) + HYPRE_BoomerAMGCreate(&solver_); + + // Set parameters from property tree with defaults + HYPRE_BoomerAMGSetPrintLevel(solver_, prm.get("print_level", 0)); + HYPRE_BoomerAMGSetMaxIter(solver_, prm.get("max_iter", 1)); + HYPRE_BoomerAMGSetStrongThreshold(solver_, prm.get("strong_threshold", 0.5)); + HYPRE_BoomerAMGSetAggTruncFactor(solver_, prm.get("agg_trunc_factor", 0.3)); + HYPRE_BoomerAMGSetInterpType(solver_, prm.get("interp_type", 6)); + HYPRE_BoomerAMGSetMaxLevels(solver_, prm.get("max_levels", 15)); + HYPRE_BoomerAMGSetTol(solver_, prm.get("tolerance", 0.0)); + + if (use_gpu_) { + HYPRE_BoomerAMGSetRelaxType(solver_, 16); + HYPRE_BoomerAMGSetCoarsenType(solver_, 8); + HYPRE_BoomerAMGSetAggNumLevels(solver_, 0); + HYPRE_BoomerAMGSetAggInterpType(solver_, 6); + // Keep transpose to avoid SpMTV + HYPRE_BoomerAMGSetKeepTranspose(solver_, true); + } + else { + HYPRE_BoomerAMGSetRelaxType(solver_, prm.get("relax_type", 13)); + HYPRE_BoomerAMGSetCoarsenType(solver_, prm.get("coarsen_type", 10)); + HYPRE_BoomerAMGSetAggNumLevels(solver_, prm.get("agg_num_levels", 1)); + HYPRE_BoomerAMGSetAggInterpType(solver_, prm.get("agg_interp_type", 4)); + } + + // Create Hypre vectors + N_ = A_.N(); + nnz_ = A_.nonzeroes(); + HYPRE_IJVectorCreate(MPI_COMM_SELF, 0, N_-1, &x_hypre_); + HYPRE_IJVectorCreate(MPI_COMM_SELF, 0, N_-1, &b_hypre_); + HYPRE_IJVectorSetObjectType(x_hypre_, HYPRE_PARCSR); + HYPRE_IJVectorSetObjectType(b_hypre_, HYPRE_PARCSR); + HYPRE_IJVectorInitialize(x_hypre_); + HYPRE_IJVectorInitialize(b_hypre_); + // Create indices vector + indices_.resize(N_); + std::iota(indices_.begin(), indices_.end(), 0); + if (use_gpu_) { + indices_device_ = hypre_CTAlloc(HYPRE_BigInt, N_, HYPRE_MEMORY_DEVICE); + hypre_TMemcpy(indices_device_, indices_.data(), HYPRE_BigInt, N_, HYPRE_MEMORY_DEVICE, HYPRE_MEMORY_HOST); + // Allocate device vectors + x_values_device_ = hypre_CTAlloc(HYPRE_Real, N_, HYPRE_MEMORY_DEVICE); + b_values_device_ = hypre_CTAlloc(HYPRE_Real, N_, HYPRE_MEMORY_DEVICE); + } + + // Create Hypre matrix + HYPRE_IJMatrixCreate(MPI_COMM_SELF, 0, N_-1, 0, N_-1, &A_hypre_); + HYPRE_IJMatrixSetObjectType(A_hypre_, HYPRE_PARCSR); + HYPRE_IJMatrixInitialize(A_hypre_); + + setupSparsityPattern(); + update(); + } + + /** + * @brief Destructor for the HyprePreconditioner class. + * + * Cleans up resources allocated by the preconditioner. + */ + ~HyprePreconditioner() { + if (solver_) { + HYPRE_BoomerAMGDestroy(solver_); + } + if (A_hypre_) { + HYPRE_IJMatrixDestroy(A_hypre_); + } + if (x_hypre_) { + HYPRE_IJVectorDestroy(x_hypre_); + } + if (b_hypre_) { + HYPRE_IJVectorDestroy(b_hypre_); + } + if (values_device_) { + hypre_TFree(values_device_, HYPRE_MEMORY_DEVICE); + } + if (x_values_device_) { + hypre_TFree(x_values_device_, HYPRE_MEMORY_DEVICE); + } + if (b_values_device_) { + hypre_TFree(b_values_device_, HYPRE_MEMORY_DEVICE); + } + if (indices_device_) { + hypre_TFree(indices_device_, HYPRE_MEMORY_DEVICE); + } + } + + /** + * @brief Updates the preconditioner with the current matrix values. + * + * This method should be called whenever the matrix values change. + */ + void update() override { + OPM_TIMEBLOCK(prec_update); + copyMatrixToHypre(); + HYPRE_BoomerAMGSetup(solver_, parcsr_A_, par_b_, par_x_); + } + + /** + * @brief Pre-processing step before applying the preconditioner. + * + * This method is currently a no-op. + * + * @param v The update vector. + * @param d The defect vector. + */ + void pre(X& /*v*/, Y& /*d*/) override { + } + + /** + * @brief Applies the preconditioner to a vector. + * + * Performs one AMG V-cycle to solve the system. + * + * @param v The update vector. + * @param d The defect vector. + */ + void apply(X& v, const Y& d) override { + OPM_TIMEBLOCK(prec_apply); + + // Copy vectors to Hypre format + copyVectorsToHypre(v, d); + + // Apply the preconditioner (one AMG V-cycle) + HYPRE_BoomerAMGSolve(solver_, parcsr_A_, par_b_, par_x_); + + // Copy result back + copyVectorFromHypre(v); + } + + /** + * @brief Post-processing step after applying the preconditioner. + * + * This method is currently a no-op. + * + * @param v The update vector. + */ + void post(X& /*v*/) override { + } + + /** + * @brief Returns the solver category. + * + * @return The solver category, which is sequential. + */ + Dune::SolverCategory::Category category() const override { + return Dune::SolverCategory::sequential; + } + + /** + * @brief Checks if the preconditioner has a perfect update. + * + * @return True, indicating that the preconditioner can be perfectly updated. + */ + bool hasPerfectUpdate() const override + { + // The Hypre preconditioner can depend on the values of the matrix so it does not have perfect update. + // However, copying the matrix to Hypre requires to setup the solver again, so this is handled internally. + // So for ISTLSolver, we can return true. + return true; + } + +private: + /** + * @brief Sets up the sparsity pattern for the Hypre matrix. + * + * Allocates and initializes arrays required by Hypre. + */ + void setupSparsityPattern() { + // Allocate arrays required by Hypre + ncols_.resize(N_); + rows_.resize(N_); + cols_.resize(nnz_); + + // Setup arrays and fill column indices + int pos = 0; + for (auto row = A_.begin(); row != A_.end(); ++row) { + const int rowIdx = row.index(); + rows_[rowIdx] = rowIdx; + ncols_[rowIdx] = row->size(); + + for (auto col = row->begin(); col != row->end(); ++col) { + cols_[pos++] = col.index(); + } + } + if (use_gpu_) { + // Allocate device arrays + ncols_device_ = hypre_CTAlloc(HYPRE_Int, N_, HYPRE_MEMORY_DEVICE); + rows_device_ = hypre_CTAlloc(HYPRE_BigInt, N_, HYPRE_MEMORY_DEVICE); + cols_device_ = hypre_CTAlloc(HYPRE_BigInt, nnz_, HYPRE_MEMORY_DEVICE); + values_device_ = hypre_CTAlloc(HYPRE_Real, nnz_, HYPRE_MEMORY_DEVICE); + + // Copy to device + hypre_TMemcpy(ncols_device_, ncols_.data(), HYPRE_Int, N_, HYPRE_MEMORY_DEVICE, HYPRE_MEMORY_HOST); + hypre_TMemcpy(rows_device_, rows_.data(), HYPRE_BigInt, N_, HYPRE_MEMORY_DEVICE, HYPRE_MEMORY_HOST); + hypre_TMemcpy(cols_device_, cols_.data(), HYPRE_BigInt, nnz_, HYPRE_MEMORY_DEVICE, HYPRE_MEMORY_HOST); + } + } + + /** + * @brief Copies the matrix values to the Hypre matrix. + * + * This method transfers the matrix data from the host to the Hypre matrix. + * It assumes that the values of the matrix are stored in a contiguous array. + * If GPU is used, the data is transferred to the device. + */ + void copyMatrixToHypre() { + OPM_TIMEBLOCK(prec_copy_matrix); + // Get pointer to matrix values array + const HYPRE_Real* values = &(A_[0][0][0][0]); + // Indexing explanation: + // A_[0] - First row of the matrix + // [0] - First block in that row + // [0] - First row within the 1x1 block + // [0] - First column within the 1x1 block + + if (use_gpu_) { + hypre_TMemcpy(values_device_, values, HYPRE_Real, nnz_, HYPRE_MEMORY_DEVICE, HYPRE_MEMORY_HOST); + HYPRE_IJMatrixSetValues(A_hypre_, N_, ncols_device_, rows_device_, cols_device_, values_device_); + } + else { + HYPRE_IJMatrixSetValues(A_hypre_, N_, ncols_.data(), rows_.data(), cols_.data(), values); + } + + HYPRE_IJMatrixAssemble(A_hypre_); + HYPRE_IJMatrixGetObject(A_hypre_, (void**)&parcsr_A_); + } + + /** + * @brief Copies vectors to the Hypre format. + * + * Transfers the update and defect vectors to Hypre. + * If GPU is used, the data is transferred from the host to the device. + * + * @param v The update vector. + * @param d The defect vector. + */ + void copyVectorsToHypre(const X& v, const Y& d) { + OPM_TIMEBLOCK(prec_copy_vectors_to_hypre); + const HYPRE_Real* x_vals = &(v[0][0]); + const HYPRE_Real* b_vals = &(d[0][0]); + + if (use_gpu_) { + hypre_TMemcpy(x_values_device_, x_vals, HYPRE_Real, N_, HYPRE_MEMORY_DEVICE, HYPRE_MEMORY_HOST); + hypre_TMemcpy(b_values_device_, b_vals, HYPRE_Real, N_, HYPRE_MEMORY_DEVICE, HYPRE_MEMORY_HOST); + + HYPRE_IJVectorSetValues(x_hypre_, N_, indices_device_, x_values_device_); + HYPRE_IJVectorSetValues(b_hypre_, N_, indices_device_, b_values_device_); + } + else { + HYPRE_IJVectorSetValues(x_hypre_, N_, indices_.data(), x_vals); + HYPRE_IJVectorSetValues(b_hypre_, N_, indices_.data(), b_vals); + } + + HYPRE_IJVectorAssemble(x_hypre_); + HYPRE_IJVectorAssemble(b_hypre_); + HYPRE_IJVectorGetObject(x_hypre_, (void**)&par_x_); + HYPRE_IJVectorGetObject(b_hypre_, (void**)&par_b_); + } + + /** + * @brief Copies the solution vector from Hypre. + * + * Transfers the solution vector from Hypre back to the host. + * If GPU is used, the data is transferred from the device to the host. + * + * @param v The update vector. + */ + void copyVectorFromHypre(X& v) { + OPM_TIMEBLOCK(prec_copy_vector_from_hypre); + HYPRE_Real* values = &(v[0][0]); + if (use_gpu_) { + HYPRE_IJVectorGetValues(x_hypre_, N_, indices_device_, x_values_device_); + hypre_TMemcpy(values, x_values_device_, HYPRE_Real, N_, HYPRE_MEMORY_HOST, HYPRE_MEMORY_DEVICE); + } + else { + HYPRE_IJVectorGetValues(x_hypre_, N_, indices_.data(), values); + } + } + + const M& A_; //!< The matrix for which the preconditioner is constructed. + bool use_gpu_ = false; //!< Flag indicating whether to use GPU acceleration. + + HYPRE_Solver solver_ = nullptr; //!< The Hypre solver object. + HYPRE_IJMatrix A_hypre_ = nullptr; //!< The Hypre matrix object. + HYPRE_ParCSRMatrix parcsr_A_ = nullptr; //!< The parallel CSR matrix object. + HYPRE_IJVector x_hypre_ = nullptr; //!< The Hypre solution vector. + HYPRE_IJVector b_hypre_ = nullptr; //!< The Hypre right-hand side vector. + HYPRE_ParVector par_x_ = nullptr; //!< The parallel solution vector. + HYPRE_ParVector par_b_ = nullptr; //!< The parallel right-hand side vector. + + std::vector ncols_; //!< Number of columns per row. + std::vector rows_; //!< Row indices. + std::vector cols_; //!< Column indices. + HYPRE_Int* ncols_device_ = nullptr; //!< Device array for number of columns per row. + HYPRE_BigInt* rows_device_ = nullptr; //!< Device array for row indices. + HYPRE_BigInt* cols_device_ = nullptr; //!< Device array for column indices. + HYPRE_Real* values_device_ = nullptr; //!< Device array for matrix values. + + std::vector indices_; //!< Indices vector for copying vectors to/from Hypre. + HYPRE_BigInt* indices_device_ = nullptr; //!< Device array for indices. + HYPRE_Int N_ = -1; //!< Number of rows in the matrix. + HYPRE_Int nnz_ = -1; //!< Number of non-zero elements in the matrix. + + HYPRE_Real* x_values_device_ = nullptr; //!< Device array for solution vector values. + HYPRE_Real* b_values_device_ = nullptr; //!< Device array for right-hand side vector values. +}; + +} // namespace Hypre + +#endif diff --git a/opm/simulators/linalg/OwningTwoLevelPreconditioner.hpp b/opm/simulators/linalg/OwningTwoLevelPreconditioner.hpp index 5f32cccc6..98ab1cac1 100644 --- a/opm/simulators/linalg/OwningTwoLevelPreconditioner.hpp +++ b/opm/simulators/linalg/OwningTwoLevelPreconditioner.hpp @@ -174,7 +174,7 @@ public: } virtual bool hasPerfectUpdate() const override { - return false; + return twolevel_method_.hasPerfectUpdate(); } private: diff --git a/opm/simulators/linalg/PreconditionerFactory_impl.hpp b/opm/simulators/linalg/PreconditionerFactory_impl.hpp index c53c069e7..04bedf0e2 100644 --- a/opm/simulators/linalg/PreconditionerFactory_impl.hpp +++ b/opm/simulators/linalg/PreconditionerFactory_impl.hpp @@ -49,6 +49,10 @@ // Include all cuistl/GPU preconditioners inside of this headerfile #include +#if HAVE_HYPRE +#include +#endif + namespace Opm { @@ -543,6 +547,15 @@ struct StandardPreconditioners { return getRebuildOnUpdateWrapper>(op, crit, parms); } }); +#if HAVE_HYPRE + // Only add Hypre for scalar matrices + if constexpr (M::block_type::rows == 1 && M::block_type::cols == 1 && + std::is_same_v) { + F::addCreator("hypre", [](const O& op, const P& prm, const std::function&, std::size_t) { + return std::make_shared>(op.getmat(), prm); + }); + } +#endif } if constexpr (std::is_same_v>) { F::addCreator( diff --git a/opm/simulators/linalg/PressureBhpTransferPolicy.hpp b/opm/simulators/linalg/PressureBhpTransferPolicy.hpp index b80f3c749..f8fedcb25 100644 --- a/opm/simulators/linalg/PressureBhpTransferPolicy.hpp +++ b/opm/simulators/linalg/PressureBhpTransferPolicy.hpp @@ -139,7 +139,7 @@ namespace Opm } } else { coarseLevelMatrix_.reset( - new CoarseMatrix(fineLevelMatrix.N(), fineLevelMatrix.M(), CoarseMatrix::row_wise)); + new CoarseMatrix(fineLevelMatrix.N(), fineLevelMatrix.M(), fineLevelMatrix.nonzeroes(), CoarseMatrix::row_wise)); auto createIter = coarseLevelMatrix_->createbegin(); for (const auto& row : fineLevelMatrix) { for (auto col = row.begin(), cend = row.end(); col != cend; ++col) { diff --git a/opm/simulators/linalg/PressureSolverPolicy.hpp b/opm/simulators/linalg/PressureSolverPolicy.hpp index 684ed1c77..0669a725a 100644 --- a/opm/simulators/linalg/PressureSolverPolicy.hpp +++ b/opm/simulators/linalg/PressureSolverPolicy.hpp @@ -88,6 +88,11 @@ namespace Amg return linsolver_->category(); } + bool hasPerfectUpdate() const + { + return linsolver_->preconditioner().hasPerfectUpdate(); + } + void apply(X& x, X& b, double reduction, Dune::InverseOperatorResult& res) override { linsolver_->apply(x, b, reduction, res); diff --git a/opm/simulators/linalg/PressureTransferPolicy.hpp b/opm/simulators/linalg/PressureTransferPolicy.hpp index db52ba576..60d19904a 100644 --- a/opm/simulators/linalg/PressureTransferPolicy.hpp +++ b/opm/simulators/linalg/PressureTransferPolicy.hpp @@ -74,7 +74,7 @@ public: { using CoarseMatrix = typename CoarseOperator::matrix_type; const auto& fineLevelMatrix = fineOperator.getmat(); - coarseLevelMatrix_.reset(new CoarseMatrix(fineLevelMatrix.N(), fineLevelMatrix.M(), CoarseMatrix::row_wise)); + coarseLevelMatrix_.reset(new CoarseMatrix(fineLevelMatrix.N(), fineLevelMatrix.M(), fineLevelMatrix.nonzeroes(), CoarseMatrix::row_wise)); auto createIter = coarseLevelMatrix_->createbegin(); for (const auto& row : fineLevelMatrix) { diff --git a/opm/simulators/linalg/twolevelmethodcpr.hh b/opm/simulators/linalg/twolevelmethodcpr.hh index 97435c82c..fbf8638ed 100644 --- a/opm/simulators/linalg/twolevelmethodcpr.hh +++ b/opm/simulators/linalg/twolevelmethodcpr.hh @@ -495,6 +495,12 @@ public: { } + bool hasPerfectUpdate() const + { + // The two-level method has perfect update if both the finesmoother and coarse solver do. + return smoother_->hasPerfectUpdate() && coarseSolver_->hasPerfectUpdate(); + } + void apply(FineDomainType& v, const FineRangeType& d) { FineDomainType u(v); diff --git a/tests/HyprePreconditionerTestHelper.hpp b/tests/HyprePreconditionerTestHelper.hpp new file mode 100644 index 000000000..af817b4bf --- /dev/null +++ b/tests/HyprePreconditionerTestHelper.hpp @@ -0,0 +1,123 @@ +/* + Copyright 2024 SINTEF AS + Copyright 2024 Equinor ASA + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#ifndef TEST_HYPREPRECONDITIONER_HELPER_HPP +#define TEST_HYPREPRECONDITIONER_HELPER_HPP + +#include + +#include +#include +#include +#include +#include + +#include +#include + + +template +void setupLaplace2d(int N, Matrix& mat) +{ + const int nonZeroes = N*N * 5; // max 5 entries per row (diagonal + 4 neighbors) + mat.setBuildMode(Matrix::row_wise); + mat.setSize(N*N, N*N, nonZeroes); + + // Set up sparsity pattern + for (auto row = mat.createbegin(); row != mat.createend(); ++row) { + const int i = row.index(); + int x = i % N; + int y = i / N; + + row.insert(i); // diagonal + if (x > 0) // left neighbor + row.insert(i-1); + if (x < N-1) // right neighbor + row.insert(i+1); + if (y > 0) // upper neighbor + row.insert(i-N); + if (y < N-1) // lower neighbor + row.insert(i+N); + } + + // Fill the matrix with values + for (auto row = mat.begin(); row != mat.end(); ++row) { + const int i = row.index(); + int x = i % N; + int y = i / N; + + // Set diagonal + (*row)[i] = 4.0; + + // Set off-diagonal entries + if (x > 0) // left neighbor + (*row)[i-1] = -1.0; + if (x < N-1) // right neighbor + (*row)[i+1] = -1.0; + if (y > 0) // upper neighbor + (*row)[i-N] = -1.0; + if (y < N-1) // lower neighbor + (*row)[i+N] = -1.0; + } +} + +inline void testHyprePreconditioner(bool use_gpu) +{ + constexpr int N = 100; // 100x100 grid + using Matrix = Dune::BCRSMatrix>; + using Vector = Dune::BlockVector>; + + // Create matrix + Matrix matrix; + setupLaplace2d(N, matrix); + + // Create vectors + Vector x(N*N), b(N*N); + x = 100.0; // Initial guess + b = 0.0; // RHS + + // Create operator + using Operator = Dune::MatrixAdapter; + Operator op(matrix); + + // Set up HYPRE parameters + Opm::PropertyTree prm; + prm.put("use_gpu", use_gpu ? 1 : 0); + + // Create preconditioner + auto prec = std::make_shared>(matrix, prm); + + // Create solver + double reduction = 1e-8; + int maxit = 300; + int verbosity = 0; + Dune::LoopSolver solver(op, *prec, reduction, maxit, verbosity); + + // Solve + Dune::InverseOperatorResult res; + solver.apply(x, b, res); + + // Check convergence + BOOST_CHECK(res.converged); + BOOST_CHECK_LT(res.reduction, 1e-8); +} + + +#endif // TEST_HYPREPRECONDITIONER_HELPER_HPP diff --git a/tests/test_HyprePreconditionerCPU.cpp b/tests/test_HyprePreconditionerCPU.cpp new file mode 100644 index 000000000..0a029da9b --- /dev/null +++ b/tests/test_HyprePreconditionerCPU.cpp @@ -0,0 +1,68 @@ +/* + Copyright 2024 SINTEF AS + Copyright 2024 Equinor ASA + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#include + +#define BOOST_TEST_MODULE TestHyprePreconditionerCPU +#define BOOST_TEST_NO_MAIN + +#include +#include "MpiFixture.hpp" +#include "HyprePreconditionerTestHelper.hpp" + +#include + +#include +#include +#include +#include +#include + +#include +#include + +BOOST_GLOBAL_FIXTURE(MPIFixture); + +BOOST_AUTO_TEST_CASE(TestHyprePreconditionerCPU) +{ + testHyprePreconditioner(false); +} + +bool init_unit_test_func() +{ + return true; +} + +int main(int argc, char** argv) +{ + Dune::MPIHelper::instance(argc, argv); + +#if HYPRE_RELEASE_NUMBER >= 22900 + HYPRE_Initialize(); +#else + HYPRE_Init(); +#endif + + int result = boost::unit_test::unit_test_main(&init_unit_test_func, argc, argv); + + HYPRE_Finalize(); + + return result; +} diff --git a/tests/test_HyprePreconditionerGPU.cpp b/tests/test_HyprePreconditionerGPU.cpp new file mode 100644 index 000000000..e6630f86d --- /dev/null +++ b/tests/test_HyprePreconditionerGPU.cpp @@ -0,0 +1,60 @@ +/* + Copyright 2024 SINTEF AS + Copyright 2024 Equinor ASA + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#include + +#define BOOST_TEST_MODULE TestHyprePreconditionerGPU +#define BOOST_TEST_NO_MAIN +#include + +#include +#include + +#include "MpiFixture.hpp" +#include "HyprePreconditionerTestHelper.hpp" + +BOOST_GLOBAL_FIXTURE(MPIFixture); + +BOOST_AUTO_TEST_CASE(TestHyprePreconditionerGPU) +{ + testHyprePreconditioner(true); +} + +bool init_unit_test_func() +{ + return true; +} + +int main(int argc, char** argv) +{ + Dune::MPIHelper::instance(argc, argv); + +#if HYPRE_RELEASE_NUMBER >= 22900 + HYPRE_Initialize(); +#else + HYPRE_Init(); +#endif + + int result = boost::unit_test::unit_test_main(&init_unit_test_func, argc, argv); + + HYPRE_Finalize(); + + return result; +}