From d623695d2a17a6e98c91c8f959d385aaaebe7091 Mon Sep 17 00:00:00 2001 From: hnil Date: Wed, 9 Aug 2023 10:22:33 +0200 Subject: [PATCH] - moded all bda spesific tings to separete class --- opm/simulators/linalg/ISTLSolverEbos.cpp | 197 ----------------------- opm/simulators/linalg/ISTLSolverEbos.hpp | 130 +-------------- 2 files changed, 3 insertions(+), 324 deletions(-) diff --git a/opm/simulators/linalg/ISTLSolverEbos.cpp b/opm/simulators/linalg/ISTLSolverEbos.cpp index 89f2465bf..6aeb1932d 100644 --- a/opm/simulators/linalg/ISTLSolverEbos.cpp +++ b/opm/simulators/linalg/ISTLSolverEbos.cpp @@ -202,169 +202,6 @@ void FlexibleSolverInfo::create(const Matrix& matrix, } } -#if COMPILE_BDA_BRIDGE -template -BdaSolverInfo:: -BdaSolverInfo(const std::string& accelerator_mode, - const int linear_solver_verbosity, - const int maxit, - const double tolerance, - const int platformID, - const int deviceID, - const bool opencl_ilu_parallel, - const std::string& linsolver) - : bridge_(std::make_unique(accelerator_mode, - linear_solver_verbosity, maxit, - tolerance, platformID, deviceID, - opencl_ilu_parallel, linsolver)) - , accelerator_mode_(accelerator_mode) -{} - -template -BdaSolverInfo::~BdaSolverInfo() = default; - -template -template -void BdaSolverInfo:: -prepare(const Grid& grid, - const Dune::CartesianIndexMapper& cartMapper, - const std::vector& wellsForConn, - const std::vector& cellPartition, - const size_t nonzeroes, - const bool useWellConn) -{ - if (numJacobiBlocks_ > 1) { - detail::setWellConnections(grid, cartMapper, wellsForConn, - useWellConn, - wellConnectionsGraph_, - numJacobiBlocks_); - this->blockJacobiAdjacency(grid, cellPartition, nonzeroes); - } -} - -template -bool BdaSolverInfo:: -apply(Vector& rhs, - const bool useWellConn, - WellContribFunc getContribs, - const int rank, - Matrix& matrix, - Vector& x, - Dune::InverseOperatorResult& result) -{ - bool use_gpu = bridge_->getUseGpu(); - if (use_gpu) { - 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 amgcl or rocalution -#if HAVE_CUDA || HAVE_OPENCL - if (!useWellConn) { - getContribs(*wellContribs); - } -#endif - - if (numJacobiBlocks_ > 1) { - this->copyMatToBlockJac(matrix, *blockJacobiForGPUILU0_); - // Const_cast needed since the CUDA stuff overwrites values for better matrix condition.. - bridge_->solve_system(&matrix, blockJacobiForGPUILU0_.get(), - numJacobiBlocks_, rhs, *wellContribs, result); - } - else - bridge_->solve_system(&matrix, &matrix, - numJacobiBlocks_, rhs, *wellContribs, result); - if (result.converged) { - // get result vector x from non-Dune backend, iff solve was successful - bridge_->get_result(x); - return true; - } else { - // warn about CPU fallback - // BdaBridge might have disabled its BdaSolver for this simulation due to some error - // in that case the BdaBridge is disabled and flexibleSolver is always used - // or maybe the BdaSolver did not converge in time, then it will be used next linear solve - if (rank == 0) { - OpmLog::warning(bridge_->getAccleratorName() + " did not converge, now trying Dune to solve current linear system..."); - } - } - } - - return false; -} - -template -bool BdaSolverInfo:: -gpuActive() -{ - return bridge_->getUseGpu(); -} - -template -template -void BdaSolverInfo:: -blockJacobiAdjacency(const Grid& grid, - const std::vector& cell_part, - size_t nonzeroes) -{ - using size_type = typename Matrix::size_type; - using Iter = typename Matrix::CreateIterator; - size_type numCells = grid.size(0); - blockJacobiForGPUILU0_ = std::make_unique(numCells, numCells, - nonzeroes, Matrix::row_wise); - - const auto& lid = grid.localIdSet(); - const auto& gridView = grid.leafGridView(); - auto elemIt = gridView.template begin<0>(); // should never overrun, since blockJacobiForGPUILU0_ is initialized with numCells rows - - //Loop over cells - for (Iter row = blockJacobiForGPUILU0_->createbegin(); row != blockJacobiForGPUILU0_->createend(); ++elemIt, ++row) - { - const auto& elem = *elemIt; - size_type idx = lid.id(elem); - row.insert(idx); - - // Add well non-zero connections - for (const auto wc : wellConnectionsGraph_[idx]) { - row.insert(wc); - } - - int locPart = cell_part[idx]; - - //Add neighbor if it is on the same part - auto isend = gridView.iend(elem); - for (auto is = gridView.ibegin(elem); is!=isend; ++is) - { - //check if face has neighbor - if (is->neighbor()) - { - size_type nid = lid.id(is->outside()); - int nabPart = cell_part[nid]; - if (locPart == nabPart) { - row.insert(nid); - } - } - } - } -} - -template -void BdaSolverInfo:: -copyMatToBlockJac(const Matrix& mat, Matrix& blockJac) -{ - auto rbegin = blockJac.begin(); - auto rend = blockJac.end(); - auto outerRow = mat.begin(); - for (auto row = rbegin; row != rend; ++row, ++outerRow) { - auto outerCol = (*outerRow).begin(); - for (auto col = (*row).begin(); col != (*row).end(); ++col) { - // outerRow is guaranteed to have all column entries that row has! - while(outerCol.index() < col.index()) ++outerCol; - assert(outerCol.index() == col.index()); - *col = *outerCol; // copy nonzero block - } - } -} -#endif // COMPILE_BDA_BRIDGE - template using BM = Dune::BCRSMatrix>; template @@ -380,40 +217,6 @@ using CommunicationType = Dune::CollectiveCommunication; template void makeOverlapRowsInvalid>(BM&, const std::vector&); \ template struct FlexibleSolverInfo,BV,CommunicationType>; -#if COMPILE_BDA_BRIDGE - -#define INSTANCE_GRID(Dim, Grid) \ - template void BdaSolverInfo,BV>:: \ - prepare(const Grid&, \ - const Dune::CartesianIndexMapper&, \ - const std::vector&, \ - const std::vector&, \ - const size_t, const bool); -using PolyHedralGrid3D = Dune::PolyhedralGrid<3, 3>; -#if HAVE_DUNE_ALUGRID -#if HAVE_MPI - using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridMPIComm>; -#else - using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridNoComm>; -#endif //HAVE_MPI -#define INSTANCE(Dim) \ - template struct BdaSolverInfo,BV>; \ - INSTANCE_GRID(Dim,Dune::CpGrid) \ - INSTANCE_GRID(Dim,ALUGrid3CN) \ - INSTANCE_GRID(Dim,PolyHedralGrid3D) \ - INSTANCE_FLEX(Dim) -#else -#define INSTANCE(Dim) \ - template struct BdaSolverInfo,BV>; \ - INSTANCE_GRID(Dim,Dune::CpGrid) \ - INSTANCE_GRID(Dim,PolyHedralGrid3D) \ - INSTANCE_FLEX(Dim) -#endif -#else -#define INSTANCE(Dim) \ - INSTANCE_FLEX(Dim) -#endif // COMPILE_BDA_BRIDGE - INSTANCE(1) INSTANCE(2) INSTANCE(3) diff --git a/opm/simulators/linalg/ISTLSolverEbos.hpp b/opm/simulators/linalg/ISTLSolverEbos.hpp index 4b2a6db76..7c99052e2 100644 --- a/opm/simulators/linalg/ISTLSolverEbos.hpp +++ b/opm/simulators/linalg/ISTLSolverEbos.hpp @@ -87,10 +87,6 @@ public: namespace Opm { -#if COMPILE_BDA_BRIDGE -template class BdaBridge; -class WellContributions; -#endif namespace detail { @@ -116,60 +112,6 @@ struct FlexibleSolverInfo size_t interiorCellNum_ = 0; }; -#if COMPILE_BDA_BRIDGE -template -struct BdaSolverInfo -{ - using WellContribFunc = std::function; - using Bridge = BdaBridge; - - BdaSolverInfo(const std::string& accelerator_mode, - const int linear_solver_verbosity, - const int maxit, - const double tolerance, - const int platformID, - const int deviceID, - const bool opencl_ilu_parallel, - const std::string& linsolver); - - ~BdaSolverInfo(); - - template - void prepare(const Grid& grid, - const Dune::CartesianIndexMapper& cartMapper, - const std::vector& wellsForConn, - const std::vector& cellPartition, - const size_t nonzeroes, - const bool useWellConn); - - bool apply(Vector& rhs, - const bool useWellConn, - WellContribFunc getContribs, - const int rank, - Matrix& matrix, - Vector& x, - Dune::InverseOperatorResult& result); - - bool gpuActive(); - - int numJacobiBlocks_ = 0; - -private: - /// Create sparsity pattern for block-Jacobi matrix based on partitioning of grid. - /// Do not initialize the values, that is done in copyMatToBlockJac() - template - void blockJacobiAdjacency(const Grid& grid, - const std::vector& cell_part, - size_t nonzeroes); - - void copyMatToBlockJac(const Matrix& mat, Matrix& blockJac); - - std::unique_ptr bridge_; - std::string accelerator_mode_; - std::unique_ptr blockJacobiForGPUILU0_; - std::vector> wellConnectionsGraph_; -}; -#endif #ifdef HAVE_MPI /// Copy values in parallel. @@ -270,37 +212,6 @@ std::unique_ptr blockJacobiAdjacency(const Grid& grid, #if HAVE_MPI comm_.reset( new CommunicationType( simulator_.vanguard().grid().comm() ) ); #endif - -#if COMPILE_BDA_BRIDGE - { - std::string accelerator_mode = EWOMS_GET_PARAM(TypeTag, std::string, AcceleratorMode); - if ((simulator_.vanguard().grid().comm().size() > 1) && (accelerator_mode != "none")) { - if (on_io_rank) { - OpmLog::warning("Cannot use GPU with MPI, GPU are disabled"); - } - accelerator_mode = "none"; - } - const int platformID = EWOMS_GET_PARAM(TypeTag, int, OpenclPlatformId); - const int deviceID = EWOMS_GET_PARAM(TypeTag, int, BdaDeviceId); - const int maxit = EWOMS_GET_PARAM(TypeTag, int, LinearSolverMaxIter); - const double tolerance = EWOMS_GET_PARAM(TypeTag, double, LinearSolverReduction); - const bool opencl_ilu_parallel = EWOMS_GET_PARAM(TypeTag, bool, OpenclIluParallel); - const int linear_solver_verbosity = parameters_.linear_solver_verbosity_; - std::string linsolver = EWOMS_GET_PARAM(TypeTag, std::string, LinearSolver); - bdaBridge = std::make_unique>(accelerator_mode, - linear_solver_verbosity, - maxit, - tolerance, - platformID, - deviceID, - opencl_ilu_parallel, - linsolver); - } -#else - if (EWOMS_GET_PARAM(TypeTag, std::string, AcceleratorMode) != "none") { - OPM_THROW(std::logic_error,"Cannot use accelerated solver since CUDA, OpenCL and amgcl were not found by cmake"); - } -#endif extractParallelGridInformationToISTL(simulator_.vanguard().grid(), parallelInformation_); // For some reason simulator_.model().elementMapper() is not initialized at this stage @@ -341,7 +252,7 @@ std::unique_ptr blockJacobiAdjacency(const Grid& grid, void prepare(const Matrix& M, Vector& b) { - OPM_TIMEBLOCK(istlSolverEbosPrepare); + OPM_TIMEBLOCK(istlSolverEbosPrepare); const bool firstcall = (matrix_ == nullptr); #if HAVE_MPI if (firstcall) { @@ -359,14 +270,6 @@ std::unique_ptr blockJacobiAdjacency(const Grid& grid, useWellConn_ = EWOMS_GET_PARAM(TypeTag, bool, MatrixAddWellContributions); // setup sparsity pattern for jacobi matrix for preconditioner (only used for openclSolver) -#if HAVE_OPENCL - bdaBridge->numJacobiBlocks_ = EWOMS_GET_PARAM(TypeTag, int, NumJacobiBlocks); - bdaBridge->prepare(simulator_.vanguard().grid(), - simulator_.vanguard().cartesianIndexMapper(), - simulator_.vanguard().schedule().getWellsatEnd(), - simulator_.vanguard().cellPartition(), - getMatrix().nonzeroes(), useWellConn_); -#endif } else { // Pointers should not change if ( &M != matrix_ ) { @@ -379,13 +282,7 @@ std::unique_ptr blockJacobiAdjacency(const Grid& grid, if (isParallel() && prm_.get("preconditioner.type") != "ParOverILU0") { detail::makeOverlapRowsInvalid(getMatrix(), overlapRows_); } -#if COMPILE_BDA_BRIDGE - if(!bdaBridge->gpuActive()){ - prepareFlexibleSolver(); - } -#else prepareFlexibleSolver(); -#endif } @@ -406,7 +303,7 @@ std::unique_ptr blockJacobiAdjacency(const Grid& grid, bool solve(Vector& x) { - OPM_TIMEBLOCK(istlSolverEbosSolve); + OPM_TIMEBLOCK(istlSolverEbosSolve); calls_ += 1; // Write linear system if asked for. const int verbosity = prm_.get("verbosity", 0); @@ -420,25 +317,8 @@ std::unique_ptr blockJacobiAdjacency(const Grid& grid, // Solve system. Dune::InverseOperatorResult result; - -#if COMPILE_BDA_BRIDGE - std::function getContribs = - [this](WellContributions& w) - { - this->simulator_.problem().wellModel().getWellContributions(w); - }; - if (!bdaBridge->apply(*rhs_, useWellConn_, getContribs, - simulator_.gridView().comm().rank(), - const_cast(this->getMatrix()), - x, result)) -#endif { OPM_TIMEBLOCK(flexibleSolverApply); -#if COMPILE_BDA_BRIDGE - if(bdaBridge->gpuActive()){ - prepareFlexibleSolver(); - } -#endif assert(flexibleSolver_.solver_); flexibleSolver_.solver_->apply(x, *rhs_, result); } @@ -521,7 +401,7 @@ std::unique_ptr blockJacobiAdjacency(const Grid& grid, } else { - OPM_TIMEBLOCK(flexibleSolverUpdate); + OPM_TIMEBLOCK(flexibleSolverUpdate); flexibleSolver_.pre_->update(); } } @@ -609,10 +489,6 @@ std::unique_ptr blockJacobiAdjacency(const Grid& grid, Matrix* matrix_; Vector *rhs_; -#if COMPILE_BDA_BRIDGE - std::unique_ptr> bdaBridge; -#endif - detail::FlexibleSolverInfo flexibleSolver_; std::vector overlapRows_; std::vector interiorRows_;