Merge pull request #5762 from jakobtorben/hypre_integration

Add Hypre BoomerAMG Support

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()

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()

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"

opm/simulators/flow/Main.cpp

@@ -35,6 +35,11 @@
 #include <opm/simulators/utils/DamarisOutputModule.hpp>
 #endif
 
+#if HAVE_HYPRE
+#include <HYPRE_config.h>
+#include <HYPRE_utilities.h>
+#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,

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 <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef OPM_HYPRE_PRECONDITIONER_HEADER_INCLUDED
+#define OPM_HYPRE_PRECONDITIONER_HEADER_INCLUDED
+
+#include <opm/common/ErrorMacros.hpp>
+#include <opm/common/TimingMacros.hpp>
+#include <opm/simulators/linalg/PreconditionerWithUpdate.hpp>
+#include <opm/simulators/linalg/PropertyTree.hpp>
+
+#include <dune/common/fmatrix.hh>
+#include <dune/istl/bcrsmatrix.hh>
+
+#include <HYPRE.h>
+#include <HYPRE_parcsr_ls.h>
+#include <HYPRE_krylov.h>
+#include <_hypre_utilities.h>
+
+#include <vector>
+#include <numeric>
+
+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 M, class X, class Y>
+class HyprePreconditioner : public Dune::PreconditionerWithUpdate<X,Y> {
+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<bool>("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<int>("print_level", 0));
+        HYPRE_BoomerAMGSetMaxIter(solver_, prm.get<int>("max_iter", 1));
+        HYPRE_BoomerAMGSetStrongThreshold(solver_, prm.get<double>("strong_threshold", 0.5));
+        HYPRE_BoomerAMGSetAggTruncFactor(solver_, prm.get<double>("agg_trunc_factor", 0.3));
+        HYPRE_BoomerAMGSetInterpType(solver_, prm.get<int>("interp_type", 6));
+        HYPRE_BoomerAMGSetMaxLevels(solver_, prm.get<int>("max_levels", 15));
+        HYPRE_BoomerAMGSetTol(solver_, prm.get<double>("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<int>("relax_type", 13));
+            HYPRE_BoomerAMGSetCoarsenType(solver_, prm.get<int>("coarsen_type", 10));
+            HYPRE_BoomerAMGSetAggNumLevels(solver_, prm.get<int>("agg_num_levels", 1));
+            HYPRE_BoomerAMGSetAggInterpType(solver_, prm.get<int>("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<HYPRE_Int> ncols_; //!< Number of columns per row.
+    std::vector<HYPRE_BigInt> rows_; //!< Row indices.
+    std::vector<HYPRE_BigInt> 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<int> 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

opm/simulators/linalg/OwningTwoLevelPreconditioner.hpp

@@ -174,7 +174,7 @@ public:
     }
 
     virtual bool hasPerfectUpdate() const override {
-        return false;
+        return twolevel_method_.hasPerfectUpdate();
     }
 
 private:

opm/simulators/linalg/PreconditionerFactory_impl.hpp

@@ -49,6 +49,10 @@
 // Include all cuistl/GPU preconditioners inside of this headerfile
 #include <opm/simulators/linalg/PreconditionerFactoryGPUIncludeWrapper.hpp>
 
+#if HAVE_HYPRE
+#include <opm/simulators/linalg/HyprePreconditioner.hpp>
+#endif
+
 
 namespace Opm {
 
@@ -543,6 +547,15 @@ struct StandardPreconditioners<Operator, Dune::Amg::SequentialInformation> {
                     return getRebuildOnUpdateWrapper<Dune::Amg::FastAMG<O, V>>(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<HYPRE_Real, typename V::field_type>) {
+                F::addCreator("hypre", [](const O& op, const P& prm, const std::function<V()>&, std::size_t) {
+                    return std::make_shared<Hypre::HyprePreconditioner<M, V, V>>(op.getmat(), prm);
+                });
+            }
+#endif
         }
         if constexpr (std::is_same_v<O, WellModelMatrixAdapter<M, V, V, false>>) {
             F::addCreator(

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) {

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);

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) {

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);

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 <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef TEST_HYPREPRECONDITIONER_HELPER_HPP
+#define TEST_HYPREPRECONDITIONER_HELPER_HPP
+
+#include <boost/test/unit_test.hpp>
+
+#include <dune/common/fmatrix.hh>
+#include <dune/istl/bcrsmatrix.hh>
+#include <dune/istl/bvector.hh>
+#include <dune/istl/operators.hh>
+#include <dune/istl/solvers.hh>
+
+#include <opm/simulators/linalg/HyprePreconditioner.hpp>
+#include <opm/simulators/linalg/PropertyTree.hpp>
+
+
+template<class Matrix>
+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<Dune::FieldMatrix<double, 1, 1>>;
+    using Vector = Dune::BlockVector<Dune::FieldVector<double, 1>>;
+
+    // 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<Matrix, Vector, Vector>;
+    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<Hypre::HyprePreconditioner<Matrix, Vector, Vector>>(matrix, prm);
+
+    // Create solver
+    double reduction = 1e-8;
+    int maxit = 300;
+    int verbosity = 0;
+    Dune::LoopSolver<Vector> 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

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 <http://www.gnu.org/licenses/>.
+*/
+
+#include <config.h>
+
+#define BOOST_TEST_MODULE TestHyprePreconditionerCPU
+#define BOOST_TEST_NO_MAIN
+
+#include <dune/common/parallel/mpihelper.hh>
+#include "MpiFixture.hpp"
+#include "HyprePreconditionerTestHelper.hpp"
+
+#include <boost/test/unit_test.hpp>
+
+#include <dune/common/fmatrix.hh>
+#include <dune/istl/bcrsmatrix.hh>
+#include <dune/istl/bvector.hh>
+#include <dune/istl/operators.hh>
+#include <dune/istl/solvers.hh>
+
+#include <opm/simulators/linalg/HyprePreconditioner.hpp>
+#include <opm/simulators/linalg/PropertyTree.hpp>
+
+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;
+}

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 <http://www.gnu.org/licenses/>.
+*/
+
+#include <config.h>
+
+#define BOOST_TEST_MODULE TestHyprePreconditionerGPU
+#define BOOST_TEST_NO_MAIN
+#include <boost/test/unit_test.hpp>
+
+#include <opm/simulators/linalg/HyprePreconditioner.hpp>
+#include <dune/common/parallel/mpihelper.hh>
+
+#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;
+}