Merge pull request #4178 from Tongdongq/rocalution

Add rocalutionSolver

CMakeLists.txt

@@ -297,6 +297,12 @@ if(OpenCL_FOUND)
   endif()
 endif()
 
+find_package(rocalution)
+if(ROCALUTION_FOUND)
+  set(HAVE_ROCALUTION 1)
+endif()
+
+
 # read the list of components from this file (in the project directory);
 # it should set various lists with the names of the files to include
 include (CMakeLists_files.cmake)
@@ -538,6 +544,10 @@ if(OpenCL_FOUND)
   target_link_libraries( opmsimulators PUBLIC ${OpenCL_LIBRARIES} )
 endif()
 
+if(ROCALUTION_FOUND)
+  target_include_directories(opmsimulators PUBLIC ${rocalution_INCLUDE_DIR}/rocalution)
+endif()
+
 if(VexCL_FOUND)
   target_link_libraries( opmsimulators PUBLIC OPM::VexCL::OpenCL )
 endif()

CMakeLists_files.cmake

@@ -143,7 +143,10 @@ if(OPENCL_FOUND)
   list (APPEND MAIN_SOURCE_FILES opm/simulators/linalg/bda/opencl/openclSolverBackend.cpp)
   list (APPEND MAIN_SOURCE_FILES opm/simulators/linalg/bda/opencl/openclWellContributions.cpp)
 endif()
-if(CUDA_FOUND OR OPENCL_FOUND OR HAVE_FPGA OR HAVE_AMGCL)
+if(ROCALUTION_FOUND)
+  list (APPEND MAIN_SOURCE_FILES opm/simulators/linalg/bda/rocalutionSolverBackend.cpp)
+endif()
+if(CUDA_FOUND OR OPENCL_FOUND OR HAVE_FPGA OR amgcl_FOUND OR ROCALUTION_FOUND)
   list (APPEND MAIN_SOURCE_FILES opm/simulators/linalg/bda/BdaBridge.cpp)
 endif()
 if(HAVE_FPGA)
@@ -151,7 +154,7 @@ if(HAVE_FPGA)
   list (APPEND MAIN_SOURCE_FILES opm/simulators/linalg/bda/FPGASolverBackend.cpp)
   list (APPEND MAIN_SOURCE_FILES opm/simulators/linalg/bda/FPGAUtils.cpp)
 endif()
-if(HAVE_AMGCL)
+if(amgcl_FOUND)
   list (APPEND MAIN_SOURCE_FILES opm/simulators/linalg/bda/amgclSolverBackend.cpp)
   if(CUDA_FOUND)
     list (APPEND MAIN_SOURCE_FILES opm/simulators/linalg/bda/cuda/amgclSolverBackend.cu)
@@ -206,6 +209,9 @@ if(OPENCL_FOUND)
   list(APPEND TEST_SOURCE_FILES tests/test_solvetransposed3x3.cpp)
   list(APPEND TEST_SOURCE_FILES tests/test_csrToCscOffsetMap.cpp)
 endif()
+if(ROCALUTION_FOUND)
+  list(APPEND TEST_SOURCE_FILES tests/test_rocalutionSolver.cpp)
+endif()
 
 list (APPEND TEST_DATA_FILES
   tests/SUMMARY_DECK_NON_CONSTANT_POROSITY.DATA
@@ -310,6 +316,7 @@ list (APPEND PUBLIC_HEADER_FILES
   opm/simulators/linalg/bda/opencl/openclWellContributions.hpp
   opm/simulators/linalg/bda/Matrix.hpp
   opm/simulators/linalg/bda/MultisegmentWellContribution.hpp
+  opm/simulators/linalg/bda/rocalutionSolverBackend.hpp
   opm/simulators/linalg/bda/WellContributions.hpp
   opm/simulators/linalg/amgcpr.hh
   opm/simulators/linalg/twolevelmethodcpr.hh

opm-simulators-prereqs.cmake

@@ -11,6 +11,7 @@ set (opm-simulators_CONFIG_VAR
   HAVE_FPGA
   HAVE_AMGCL
   HAVE_VEXCL
+  HAVE_ROCALUTION
   HAVE_SUITESPARSE_UMFPACK_H
   HAVE_DUNE_ISTL
   DUNE_ISTL_WITH_CHECKING
@@ -40,6 +41,8 @@ set (opm-simulators_DEPS
   "SuiteSparse REQUIRED COMPONENTS umfpack"
   # SuperLU direct solver
   "SuperLU"
+  # ROCALUTION from ROCM framework
+  "rocalution"
   # OPM dependency
   "opm-common REQUIRED"
   "opm-material REQUIRED"

opm/simulators/linalg/ISTLSolverEbos.cpp

@@ -31,7 +31,7 @@
 #include <opm/simulators/linalg/ParallelIstlInformation.hpp>
 #include <opm/simulators/utils/ParallelCommunication.hpp>
 
-#if HAVE_CUDA || HAVE_OPENCL || HAVE_FPGA || HAVE_AMGCL
+#if HAVE_CUDA || HAVE_OPENCL || HAVE_FPGA || HAVE_AMGCL || HAVE_ROCALUTION
 #include <opm/simulators/linalg/bda/BdaBridge.hpp>
 #include <opm/simulators/linalg/bda/WellContributions.hpp>
 #include <iostream>
@@ -171,7 +171,7 @@ void FlexibleSolverInfo<Matrix,Vector,Comm>::create(const Matrix& matrix,
     }
 }
 
-#if HAVE_CUDA || HAVE_OPENCL || HAVE_FPGA || HAVE_AMGCL
+#if HAVE_CUDA || HAVE_OPENCL || HAVE_FPGA || HAVE_AMGCL || HAVE_ROCALUTION
 template<class Matrix, class Vector>
 BdaSolverInfo<Matrix,Vector>::
 BdaSolverInfo(const std::string& accelerator_mode,
@@ -229,7 +229,7 @@ apply(Vector& rhs,
         auto wellContribs = WellContributions::create(accelerator_mode_, useWellConn);
         bridge_->initWellContributions(*wellContribs, x.N() * x[0].N());
 
-        // the WellContributions can only be applied separately with CUDA or OpenCL, not with an FPGA or amgcl
+        // the WellContributions can only be applied separately with CUDA or OpenCL, not with an FPGA, amgcl or rocalution
 #if HAVE_CUDA || HAVE_OPENCL
         if (!useWellConn) {
             getContribs(*wellContribs);
@@ -345,7 +345,7 @@ using CommunicationType = Dune::CollectiveCommunication<int>;
     template void makeOverlapRowsInvalid<BM<Dim>>(BM<Dim>&, const std::vector<int>&); \
     template struct FlexibleSolverInfo<BM<Dim>,BV<Dim>,CommunicationType>;
 
-#if HAVE_CUDA || HAVE_OPENCL || HAVE_FPGA || HAVE_AMGCL
+#if HAVE_CUDA || HAVE_OPENCL || HAVE_FPGA || HAVE_AMGCL || HAVE_ROCALUTION
 #define INSTANCE(Dim) \
     template struct BdaSolverInfo<BM<Dim>,BV<Dim>>; \
     template void BdaSolverInfo<BM<Dim>,BV<Dim>>:: \

opm/simulators/linalg/ISTLSolverEbos.hpp

@@ -84,7 +84,7 @@ public:
 namespace Opm
 {
 
-#if HAVE_CUDA || HAVE_OPENCL || HAVE_FPGA || HAVE_AMGCL
+#if HAVE_CUDA || HAVE_OPENCL || HAVE_FPGA || HAVE_AMGCL || HAVE_ROCALUTION
 template<class Matrix, class Vector, int block_size> class BdaBridge;
 class WellContributions;
 #endif
@@ -113,7 +113,7 @@ struct FlexibleSolverInfo
     size_t interiorCellNum_ = 0;
 };
 
-#if HAVE_CUDA || HAVE_OPENCL || HAVE_FPGA || HAVE_AMGCL
+#if HAVE_CUDA || HAVE_OPENCL || HAVE_FPGA || HAVE_AMGCL || HAVE_ROCALUTION
 template<class Matrix, class Vector>
 struct BdaSolverInfo
 {
@@ -246,7 +246,7 @@ std::unique_ptr<Matrix> blockJacobiAdjacency(const Grid& grid,
                                      EWOMS_PARAM_IS_SET(TypeTag, int, LinearSolverMaxIter),
                                      EWOMS_PARAM_IS_SET(TypeTag, int, CprMaxEllIter));
 
-#if HAVE_CUDA || HAVE_OPENCL || HAVE_FPGA || HAVE_AMGCL
+#if HAVE_CUDA || HAVE_OPENCL || HAVE_FPGA || HAVE_AMGCL || HAVE_ROCALUTION
             {
                 std::string accelerator_mode = EWOMS_GET_PARAM(TypeTag, std::string, AcceleratorMode);
                 if ((simulator_.vanguard().grid().comm().size() > 1) && (accelerator_mode != "none")) {
@@ -394,7 +394,7 @@ std::unique_ptr<Matrix> blockJacobiAdjacency(const Grid& grid,
 
             // Use GPU if: available, chosen by user, and successful.
             // Use FPGA if: support compiled, chosen by user, and successful.
-#if HAVE_CUDA || HAVE_OPENCL || HAVE_FPGA || HAVE_AMGCL
+#if HAVE_CUDA || HAVE_OPENCL || HAVE_FPGA || HAVE_AMGCL || HAVE_ROCALUTION
             std::function<void(WellContributions&)> getContribs =
                 [this](WellContributions& w)
                 {
@@ -566,7 +566,7 @@ std::unique_ptr<Matrix> blockJacobiAdjacency(const Grid& grid,
         Matrix* matrix_;
         Vector *rhs_;
 
-#if HAVE_CUDA || HAVE_OPENCL || HAVE_FPGA || HAVE_AMGCL
+#if HAVE_CUDA || HAVE_OPENCL || HAVE_FPGA || HAVE_AMGCL || HAVE_ROCALUTION
         std::unique_ptr<detail::BdaSolverInfo<Matrix, Vector>> bdaBridge;
 #endif
 

opm/simulators/linalg/bda/BdaBridge.cpp

@@ -45,6 +45,10 @@
 #include <opm/simulators/linalg/bda/amgclSolverBackend.hpp>
 #endif
 
+#if HAVE_ROCALUTION
+#include <opm/simulators/linalg/bda/rocalutionSolverBackend.hpp>
+#endif
+
 typedef Dune::InverseOperatorResult InverseOperatorResult;
 
 namespace Opm
@@ -92,12 +96,19 @@ BdaBridge<BridgeMatrix, BridgeVector, block_size>::BdaBridge(std::string acceler
         backend.reset(new Opm::Accelerator::amgclSolverBackend<block_size>(linear_solver_verbosity, maxit, tolerance, platformID, deviceID));
 #else
         OPM_THROW(std::logic_error, "Error amgclSolver was chosen, but amgcl was not found by CMake");
+#endif
+    } else if (accelerator_mode.compare("rocalution") == 0) {
+#if HAVE_ROCALUTION
+        use_gpu = true; // should be replaced by a 'use_bridge' boolean
+        backend.reset(new Opm::Accelerator::rocalutionSolverBackend<block_size>(linear_solver_verbosity, maxit, tolerance));
+#else
+        OPM_THROW(std::logic_error, "Error rocalutionSolver was chosen, but rocalution was not found by CMake");
 #endif
     } else if (accelerator_mode.compare("none") == 0) {
         use_gpu = false;
         use_fpga = false;
     } else {
-        OPM_THROW(std::logic_error, "Error unknown value for parameter 'AcceleratorMode', should be passed like '--accelerator-mode=[none|cusparse|opencl|fpga|amgcl]");
+        OPM_THROW(std::logic_error, "Error unknown value for parameter 'AcceleratorMode', should be passed like '--accelerator-mode=[none|cusparse|opencl|fpga|amgcl|rocalution]");
     }
 }
 

opm/simulators/linalg/bda/WellContributions.cpp

@@ -63,6 +63,12 @@ WellContributions::create(const std::string& accelerator_mode, bool useWellConn)
         }
         return std::make_unique<WellContributions>();
     }
+    else if(accelerator_mode.compare("rocalution") == 0){
+        if (!useWellConn) {
+            OPM_THROW(std::logic_error, "Error rocalution requires --matrix-add-well-contributions=true");
+        }
+        return std::make_unique<WellContributions>();
+    }
     else{
         OPM_THROW(std::logic_error, "Invalid accelerator mode");
     }

opm/simulators/linalg/bda/rocalutionSolverBackend.cpp

@@ -0,0 +1,228 @@
+/*
+  Copyright 2022 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>
+#include <cmath>
+#include <sstream>
+#include <fmt/format.h>
+
+#include <opm/common/OpmLog/OpmLog.hpp>
+#include <opm/common/ErrorMacros.hpp>
+#include <dune/common/timer.hh>
+
+#include <opm/simulators/linalg/bda/rocalutionSolverBackend.hpp>
+
+#include <rocalution.hpp>
+#include <base/matrix_formats_ind.hpp> // check if blocks are interpreted as row-major or column-major
+
+namespace Opm
+{
+namespace Accelerator
+{
+
+using Opm::OpmLog;
+using Dune::Timer;
+
+template <unsigned int block_size>
+rocalutionSolverBackend<block_size>::rocalutionSolverBackend(int verbosity_, int maxit_, double tolerance_) : BdaSolver<block_size>(verbosity_, maxit_, tolerance_) {
+    rocalution::init_rocalution();
+    rocalution::info_rocalution();
+    roc_solver = std::make_unique<rocalution::BiCGStab<rocalution::LocalMatrix<double>, rocalution::LocalVector<double>, double> >();
+    roc_prec = std::make_unique<rocalution::ILU<rocalution::LocalMatrix<double>, rocalution::LocalVector<double>, double> >();
+    roc_solver->Verbose(0);
+    roc_solver->Init(/*abs_tol=*/1e-15, tolerance, /*divergence_tol=*/1e3, maxit);
+}
+
+
+template <unsigned int block_size>
+rocalutionSolverBackend<block_size>::~rocalutionSolverBackend() {
+    rocalution::stop_rocalution();
+}
+
+
+template <unsigned int block_size>
+void rocalutionSolverBackend<block_size>::initialize(BlockedMatrix *matrix) {
+    this->Nb = matrix->Nb;
+    this->N = Nb * block_size;
+    this->nnzb = matrix->nnzbs;
+    this->nnz = nnzb * block_size * block_size;
+
+    std::ostringstream out;
+    out << fmt::format("Initializing rocalution, matrix size: {} blockrows, nnzb: {}\n", Nb, nnzb);
+    out << fmt::format("Maxit: {}, tolerance: {:e}\n", maxit, tolerance);
+    OpmLog::info(out.str());
+
+    h_x.resize(Nb * block_size);
+
+    initialized = true;
+} // end initialize()
+
+
+template <unsigned int block_size>
+void rocalutionSolverBackend<block_size>::convert_matrix(BlockedMatrix *matrix) {
+    Timer t;
+
+    for(int i = 0; i < Nb+1; ++i){
+        tmp_rowpointers[i] = matrix->rowPointers[i];
+    }
+    for(int i = 0; i < nnzb; ++i){
+        tmp_colindices[i] = matrix->colIndices[i];
+    }
+
+    // convert values inside block from row major to col major
+    // this is the same as transposing a block
+    // when compiling rocm from scratch, it is possible to choose the direction, making this transposing unnecessary
+    // BCSR_IND_BASE == 0: rocalution expects column-major
+    // BCSR_IND_BASE == 1: rocalution expects row-major
+    if (BCSR_IND_BASE == 0) {
+        for(int i = 0; i < nnzb; ++i){
+            tmp_nnzvalues[i * block_size * block_size + 0] = matrix->nnzValues[i * block_size * block_size + 0];
+            tmp_nnzvalues[i * block_size * block_size + 1] = matrix->nnzValues[i * block_size * block_size + 3];
+            tmp_nnzvalues[i * block_size * block_size + 2] = matrix->nnzValues[i * block_size * block_size + 6];
+            tmp_nnzvalues[i * block_size * block_size + 3] = matrix->nnzValues[i * block_size * block_size + 1];
+            tmp_nnzvalues[i * block_size * block_size + 4] = matrix->nnzValues[i * block_size * block_size + 4];
+            tmp_nnzvalues[i * block_size * block_size + 5] = matrix->nnzValues[i * block_size * block_size + 7];
+            tmp_nnzvalues[i * block_size * block_size + 6] = matrix->nnzValues[i * block_size * block_size + 2];
+            tmp_nnzvalues[i * block_size * block_size + 7] = matrix->nnzValues[i * block_size * block_size + 5];
+            tmp_nnzvalues[i * block_size * block_size + 8] = matrix->nnzValues[i * block_size * block_size + 8];
+        }
+    }
+    if (verbosity >= 3) {
+        std::ostringstream out;
+        out << "rocalutionSolver::convert_matrix(): " << t.stop() << " s";
+        OpmLog::info(out.str());
+    }
+}
+
+
+// copy result to host memory
+// caller must be sure that x is a valid array
+template <unsigned int block_size>
+void rocalutionSolverBackend<block_size>::get_result(double *x) {
+    Timer t;
+
+    std::copy(h_x.begin(), h_x.end(), x);
+
+    if (verbosity >= 3) {
+        std::ostringstream out;
+        out << "rocalutionSolver::get_result(): " << t.stop() << " s";
+        OpmLog::info(out.str());
+    }
+} // end get_result()
+
+
+template <unsigned int block_size>
+SolverStatus rocalutionSolverBackend<block_size>::solve_system(std::shared_ptr<BlockedMatrix> matrix,
+                                                           double *b,
+                                                           [[maybe_unused]] std::shared_ptr<BlockedMatrix> jacMatrix,
+                                                           [[maybe_unused]] WellContributions& wellContribs,
+                                                           BdaResult &res)
+{
+    if (initialized == false) {
+        initialize(matrix.get());
+    }
+
+    tmp_rowpointers = new int[Nb+1];
+    tmp_colindices = new int[nnzb];
+    tmp_nnzvalues = new double[nnzb*block_size*block_size];
+
+    convert_matrix(matrix.get());
+
+    rocalution::LocalVector<double> roc_x;
+    rocalution::LocalVector<double> roc_rhs;
+    rocalution::LocalMatrix<double> roc_mat;
+
+    // this also transfers ownership to the allocated memory to rocalution
+    // and sets the tmp_* pointers to nullptr
+    roc_mat.SetDataPtrBCSR(
+        &tmp_rowpointers,
+        &tmp_colindices,
+        &tmp_nnzvalues,
+        "matrix A", nnzb, Nb, Nb, block_size);
+
+    roc_mat.MoveToAccelerator();
+    roc_x.MoveToAccelerator();
+    roc_rhs.MoveToAccelerator();
+
+    roc_x.Allocate("x", roc_mat.GetN());
+    roc_rhs.Allocate("rhs", roc_mat.GetN());
+
+    // initialize vectors
+    roc_rhs.CopyFromData(b);
+    roc_x.Zeros();
+
+    roc_solver->Clear();
+    roc_solver->SetOperator(roc_mat);
+    roc_solver->SetPreconditioner(*roc_prec);
+
+    // the implementation of ILU::ReBuildNumeric() does not exist at the time of writing
+    // so it just calls ILU::Build() everytime
+    roc_solver->ReBuildNumeric();
+
+    double norm_0 = roc_rhs.Norm(); // since the initial guess is a vector with 0s, initial error is norm(b)
+
+    // actually solve
+    Dune::Timer t_solve;
+    roc_solver->Solve(roc_rhs, &roc_x);
+
+    // roc_solver->GetSolverStatus() returns:
+    // 0, if no criteria has been reached yet
+    // 1, if absolute tolerance has been reached
+    // 2, if relative tolerance has been reached
+    // 3, if divergence tolerance has been reached
+    // 4, if maximum number of iteration has been reached
+
+    res.elapsed = t_solve.stop();
+    res.iterations = roc_solver->GetIterationCount();
+    res.reduction = roc_solver->GetCurrentResidual() / norm_0;
+    res.conv_rate  = static_cast<double>(pow(res.reduction, 1.0 / res.iterations));
+    res.converged = (roc_solver->GetSolverStatus() == 2);
+
+
+    // copy solution vector to host vector
+    // if roc_x could be reused, this should be removed here
+    // and roc_x should be directly copied into x in get_result()
+    roc_x.MoveToHost();
+    roc_x.CopyToData(h_x.data());
+
+    if (verbosity >= 1) {
+        std::ostringstream out;
+        out << "=== converged: " << res.converged << ", conv_rate: " << res.conv_rate << ", time: " << res.elapsed << \
+            ", time per iteration: " << res.elapsed / res.iterations << ", iterations: " << res.iterations;
+        OpmLog::info(out.str());
+    }
+
+    return SolverStatus::BDA_SOLVER_SUCCESS;
+}
+
+
+#define INSTANTIATE_BDA_FUNCTIONS(n) \
+template rocalutionSolverBackend<n>::rocalutionSolverBackend(int, int, double);
+
+INSTANTIATE_BDA_FUNCTIONS(1);
+INSTANTIATE_BDA_FUNCTIONS(2);
+INSTANTIATE_BDA_FUNCTIONS(3);
+INSTANTIATE_BDA_FUNCTIONS(4);
+INSTANTIATE_BDA_FUNCTIONS(5);
+INSTANTIATE_BDA_FUNCTIONS(6);
+
+#undef INSTANTIATE_BDA_FUNCTIONS
+
+} // namespace Accelerator
+} // namespace Opm

opm/simulators/linalg/bda/rocalutionSolverBackend.hpp

@@ -0,0 +1,109 @@
+/*
+  Copyright 2022 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_ROCALUTIONSOLVER_BACKEND_HEADER_INCLUDED
+#define OPM_ROCALUTIONSOLVER_BACKEND_HEADER_INCLUDED
+
+#include <opm/simulators/linalg/bda/BdaResult.hpp>
+#include <opm/simulators/linalg/bda/BdaSolver.hpp>
+#include <opm/simulators/linalg/bda/WellContributions.hpp>
+
+namespace rocalution {
+template<class Matrix, class Vector, class Scalar> class BiCGStab;
+template<class Matrix, class Vector, class Scalar> class ILU;
+template<class Scalar> class LocalMatrix;
+template<class Scalar> class LocalVector;
+}
+
+namespace Opm
+{
+namespace Accelerator
+{
+
+/// This class implements a rocalution based linear solver solver on GPU
+/// It uses ilu0-bicgstab
+template <unsigned int block_size>
+class rocalutionSolverBackend : public BdaSolver<block_size>
+{
+    typedef BdaSolver<block_size> Base;
+
+    using Base::N;
+    using Base::Nb;
+    using Base::nnz;
+    using Base::nnzb;
+    using Base::verbosity;
+    using Base::platformID;
+    using Base::deviceID;
+    using Base::maxit;
+    using Base::tolerance;
+    using Base::initialized;
+
+private:
+    std::vector<double> h_x; // store solution vector on host
+    int *tmp_rowpointers;    // store matrix on host, this pointer is given to and freed by rocalution
+    int *tmp_colindices;     // store matrix on host, this pointer is given to and freed by rocalution
+    double *tmp_nnzvalues;   // store matrix on host, this pointer is given to and freed by rocalution
+
+    // must be declared in this order, to prevent a segfault during the test
+    std::unique_ptr<rocalution::ILU<rocalution::LocalMatrix<double>, rocalution::LocalVector<double>, double> > roc_prec;
+    std::unique_ptr<rocalution::BiCGStab<rocalution::LocalMatrix<double>, rocalution::LocalVector<double>, double> > roc_solver;
+
+    /// Initialize sizes and allocate memory
+    /// \param[in] matrix     matrix A
+    void initialize(BlockedMatrix *matrix);
+
+    /// Convert matrix to rocalution format
+    /// copy matrix to raw pointers, which are given to and freed by rocalution
+    /// \param[in] matrix     matrix A
+    void convert_matrix(BlockedMatrix *matrix);
+
+public:
+
+    /// Construct a rocalutionSolver
+    /// also initialize rocalution library and rocalution variables
+    /// \param[in] linear_solver_verbosity    verbosity of rocalutionSolver
+    /// \param[in] maxit                      maximum number of iterations for rocalutionSolver
+    /// \param[in] tolerance                  required relative tolerance for rocalutionSolver
+    rocalutionSolverBackend(int linear_solver_verbosity, int maxit, double tolerance);
+
+    /// Destroy a rocalutionSolver, and free memory
+    ~rocalutionSolverBackend();
+
+    /// Solve linear system, A*x = b, matrix A must be in blocked-CSR format
+    /// \param[in] matrix         matrix A
+    /// \param[in] b              input vector, contains N values
+    /// \param[in] jacMatrix      matrix for preconditioner
+    /// \param[in] wellContribs   WellContributions, to apply them separately, instead of adding them to matrix A
+    /// \param[inout] res         summary of solver result
+    /// \return                   status code
+    SolverStatus solve_system(std::shared_ptr<BlockedMatrix> matrix, double *b,
+        std::shared_ptr<BlockedMatrix> jacMatrix, WellContributions& wellContribs, BdaResult &res) override;
+
+    /// Get result after linear solve, and peform postprocessing if necessary
+    /// \param[inout] x          resulting x vector, caller must guarantee that x points to a valid array
+    void get_result(double *x) override;
+    
+}; // end class rocalutionSolverBackend
+
+} // namespace Accelerator
+} // namespace Opm
+
+#endif
+
+

tests/test_cusparseSolver.cpp

@@ -22,12 +22,6 @@
 
 #define BOOST_TEST_MODULE OPM_test_cusparseSolver
 #include <boost/test/unit_test.hpp>
-#include <boost/version.hpp>
-
-#include <dune/common/version.hh>
-
-#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6) && \
-    BOOST_VERSION / 100 % 1000 > 48
 
 #include <opm/simulators/linalg/bda/BdaBridge.hpp>
 #include <opm/simulators/linalg/bda/WellContributions.hpp>
@@ -158,7 +152,7 @@ void test3(const pt::ptree& prm)
 }
 
 
-BOOST_AUTO_TEST_CASE(TestDefaultPreconditionerFactory)
+BOOST_AUTO_TEST_CASE(TestCusparseSolver)
 {
     pt::ptree prm;
 
@@ -175,14 +169,3 @@ BOOST_AUTO_TEST_CASE(TestDefaultPreconditionerFactory)
         BOOST_WARN_MESSAGE(true, "Problem with initializing a device. skipping test");
     }
 }
-
-
-#else
-
-// Do nothing if we do not have at least Dune 2.6.
-BOOST_AUTO_TEST_CASE(DummyTest)
-{
-    BOOST_REQUIRE(true);
-}
-
-#endif

tests/test_flexiblesolver.cpp

@@ -21,12 +21,6 @@
 
 #define BOOST_TEST_MODULE OPM_test_FlexibleSolver
 #include <boost/test/unit_test.hpp>
-#include <boost/version.hpp>
-
-#include <dune/common/version.hh>
-
-#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6) && \
-    BOOST_VERSION / 100 % 1000 > 48
 
 #include <opm/simulators/linalg/FlexibleSolver.hpp>
 #include <opm/simulators/linalg/getQuasiImpesWeights.hpp>
@@ -124,13 +118,3 @@ BOOST_AUTO_TEST_CASE(TestFlexibleSolver)
         }
     }
 }
-
-#else
-
-// Do nothing if we do not have at least Dune 2.6.
-BOOST_AUTO_TEST_CASE(DummyTest)
-{
-    BOOST_REQUIRE(true);
-}
-
-#endif

tests/test_openclSolver.cpp

@@ -22,12 +22,6 @@
 
 #define BOOST_TEST_MODULE OPM_test_openclSolver
 #include <boost/test/unit_test.hpp>
-#include <boost/version.hpp>
-
-#include <dune/common/version.hh>
-
-#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6) && \
-    BOOST_VERSION / 100 % 1000 > 48
 
 #include <opm/simulators/linalg/bda/BdaBridge.hpp>
 #include <opm/simulators/linalg/bda/WellContributions.hpp>
@@ -189,7 +183,7 @@ void test3(const pt::ptree& prm)
 }
 
 
-BOOST_AUTO_TEST_CASE(TestDefaultPreconditionerFactory)
+BOOST_AUTO_TEST_CASE(TestOpenclSolver)
 {
     pt::ptree prm;
 
@@ -206,14 +200,3 @@ BOOST_AUTO_TEST_CASE(TestDefaultPreconditionerFactory)
         BOOST_WARN_MESSAGE(true, "Problem with initializing Platform. skipping test");
     }
 }
-
-
-#else
-
-// Do nothing if we do not have at least Dune 2.6.
-BOOST_AUTO_TEST_CASE(DummyTest)
-{
-    BOOST_REQUIRE(true);
-}
-
-#endif

tests/test_preconditionerfactory.cpp

@@ -21,12 +21,6 @@
 
 #define BOOST_TEST_MODULE OPM_test_PreconditionerFactory
 #include <boost/test/unit_test.hpp>
-#include <boost/version.hpp>
-
-#include <dune/common/version.hh>
-
-#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6) && \
-    BOOST_VERSION / 100 % 1000 > 48
 
 #include <opm/simulators/linalg/matrixblock.hh>
 #include <opm/simulators/linalg/ilufirstelement.hh>
@@ -370,15 +364,3 @@ BOOST_AUTO_TEST_CASE(TestWithRepeatingOperator)
     // Test with 3x3 block solvers.
     test3rep(prm);
 }
-
-
-
-#else
-
-// Do nothing if we do not have at least Dune 2.6.
-BOOST_AUTO_TEST_CASE(DummyTest)
-{
-    BOOST_REQUIRE(true);
-}
-
-#endif

tests/test_rocalutionSolver.cpp

@@ -0,0 +1,163 @@
+/*
+  Copyright 2019 SINTEF Digital, Mathematics and Cybernetics.
+  Copyright 2022 Equinor
+
+  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 OPM_test_rocalutionSolver
+#include <boost/test/unit_test.hpp>
+
+#include <opm/simulators/linalg/bda/BdaBridge.hpp>
+#include <opm/simulators/linalg/bda/WellContributions.hpp>
+#include <rocalution.hpp>
+
+#include <dune/common/fvector.hh>
+#include <dune/istl/bvector.hh>
+#include <dune/istl/bcrsmatrix.hh>
+#include <dune/istl/matrixmarket.hh>
+#include <dune/istl/solvers.hh>
+#include <dune/istl/preconditioners.hh>
+
+#include <boost/property_tree/json_parser.hpp>
+#include <boost/property_tree/ptree.hpp>
+
+template <int bz>
+using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, bz, bz>>;
+template <int bz>
+using Vector = Dune::BlockVector<Dune::FieldVector<double, bz>>;
+
+template <int bz>
+void readLinearSystem(const std::string& matrix_filename, const std::string& rhs_filename, Matrix<bz>& matrix, Vector<bz>& rhs)
+{
+    {
+        std::ifstream mfile(matrix_filename);
+        if (!mfile) {
+            throw std::runtime_error("Could not read matrix file");
+        }
+        readMatrixMarket(matrix, mfile);
+    }
+    {
+        std::ifstream rhsfile(rhs_filename);
+        if (!rhsfile) {
+            throw std::runtime_error("Could not read rhs file");
+        }
+        readMatrixMarket(rhs, rhsfile);
+    }
+}
+
+template <int bz>
+Dune::BlockVector<Dune::FieldVector<double, bz>>
+getDuneSolution(Matrix<bz>& matrix, Vector<bz>& rhs)
+{
+    Dune::InverseOperatorResult result;
+
+    Vector<bz> x(rhs.size());
+
+    typedef Dune::MatrixAdapter<Matrix<bz>,Vector<bz>,Vector<bz> > Operator;
+    Operator fop(matrix);
+    double relaxation = 0.9;
+    Dune::SeqILU<Matrix<bz>,Vector<bz>,Vector<bz> > prec(matrix, relaxation);
+    double reduction = 1e-2;
+    int maxit = 10;
+    int verbosity = 0;
+    Dune::BiCGSTABSolver<Vector<bz> > solver(fop, prec, reduction, maxit, verbosity);
+    solver.apply(x, rhs, result);
+    return x;
+}
+
+template <int bz>
+Dune::BlockVector<Dune::FieldVector<double, bz>>
+testRocalutionSolver(const boost::property_tree::ptree& prm, Matrix<bz>& matrix, Vector<bz>& rhs)
+{
+    const int linear_solver_verbosity = prm.get<int>("verbosity");
+    const int maxit = prm.get<int>("maxiter");
+    const double tolerance = prm.get<double>("tol");
+    const bool opencl_ilu_parallel(true);
+    const int platformID = 0;
+    const int deviceID = 0;
+    const std::string accelerator_mode("rocalution");
+    const std::string fpga_bitstream("empty");    // unused
+    const std::string linsolver("ilu0");
+    Dune::InverseOperatorResult result;
+
+    Vector<bz> x(rhs.size());
+    auto wellContribs = Opm::WellContributions::create(accelerator_mode, true);
+    std::unique_ptr<Opm::BdaBridge<Matrix<bz>, Vector<bz>, bz> > bridge;
+    try {
+        bridge = std::make_unique<Opm::BdaBridge<Matrix<bz>, Vector<bz>, bz> >(accelerator_mode,
+                                                                               fpga_bitstream,
+                                                                               linear_solver_verbosity,
+                                                                               maxit,
+                                                                               tolerance,
+                                                                               platformID,
+                                                                               deviceID,
+                                                                               opencl_ilu_parallel,
+                                                                               linsolver);
+    } catch (const std::logic_error& error) {
+        BOOST_WARN_MESSAGE(true, error.what());
+    }
+    auto mat2 = matrix; // deep copy to make sure nnz values are in contiguous memory
+                        // matrix created by readMatrixMarket() did not have contiguous memory
+    bridge->solve_system(&mat2, &mat2, /*numJacobiBlocks=*/0, rhs, *wellContribs, result);
+    bridge->get_result(x);
+
+    return x;
+}
+
+namespace pt = boost::property_tree;
+
+void test3(const pt::ptree& prm)
+{
+    const int bz = 3;
+    Matrix<bz> matrix;
+    Vector<bz> rhs;
+    readLinearSystem("matr33.txt", "rhs3.txt", matrix, rhs);
+    Vector<bz> rhs2 = rhs; // deep copy, getDuneSolution() changes values in rhs vector
+    auto duneSolution = getDuneSolution<bz>(matrix, rhs);
+    auto sol = testRocalutionSolver<bz>(prm, matrix, rhs2);
+
+    BOOST_REQUIRE_EQUAL(sol.size(), duneSolution.size());
+    for (size_t i = 0; i < sol.size(); ++i) {
+        for (int row = 0; row < bz; ++row) {
+            BOOST_CHECK_CLOSE(sol[i][row], duneSolution[i][row], 1e-3);
+        }
+    }
+}
+
+
+BOOST_AUTO_TEST_CASE(TestRocalutionSolver)
+{
+    pt::ptree prm;
+
+    // Read parameters.
+    {
+        std::ifstream file("options_flexiblesolver.json");
+        pt::read_json(file, prm);
+    }
+
+    rocalution::init_rocalution();
+    auto rocalution_backend_descriptor = rocalution::_get_backend_descriptor();
+
+    if (rocalution_backend_descriptor->accelerator) {
+        // test rocalution with 3x3 blocks
+        test3(prm);
+    } else {
+        BOOST_WARN_MESSAGE(true, "Problem with initializing a device. skipping test");
+    }
+}