diff --git a/opm/simulators/linalg/ISTLSolverEbos.hpp b/opm/simulators/linalg/ISTLSolverEbos.hpp index a42cead5c..64b11c9a2 100644 --- a/opm/simulators/linalg/ISTLSolverEbos.hpp +++ b/opm/simulators/linalg/ISTLSolverEbos.hpp @@ -240,7 +240,16 @@ protected: { parameters_.template init(); extractParallelGridInformationToISTL(simulator_.vanguard().grid(), parallelInformation_); - detail::findOverlapRowsAndColumns(simulator_.vanguard().grid(),overlapRowAndColumns_); + auto gridForConn = simulator_.vanguard().grid(); + const auto wellsForConn = simulator_.vanguard().schedule().getWellsatEnd(); + const bool useWellConn = EWOMS_GET_PARAM(TypeTag, bool, MatrixAddWellContributions); + + detail::setWellConnections(gridForConn, wellsForConn, useWellConn, wellConnectionsGraph_); + detail::findOverlapAndInterior(gridForConn, overlapRows_, interiorRows_); + if (gridForConn.comm().size() > 1) { + noGhostAdjacency(); + setGhostsInNoGhost(*noGhostMat_); + } } // nothing to clean here @@ -322,13 +331,8 @@ protected: { typedef WellModelMatrixAdapter< Matrix, Vector, Vector, WellModel, true > Operator; - auto ebosJacIgnoreOverlap = Matrix(*matrix_); - //remove ghost rows in local matrix - makeOverlapRowsInvalid(ebosJacIgnoreOverlap); - - //Not sure what actual_mat_for_prec is, so put ebosJacIgnoreOverlap as both variables - //to be certain that correct matrix is used for preconditioning. - Operator opA(ebosJacIgnoreOverlap, ebosJacIgnoreOverlap, wellModel, + copyJacToNoGhost(*matrix_, *noGhostMat_); + Operator opA(*noGhostMat_, *noGhostMat_, wellModel, parallelInformation_ ); assert( opA.comm() ); solve( opA, x, *rhs_, *(opA.comm()) ); @@ -384,9 +388,6 @@ protected: SPPointer sp(ScalarProductChooser::construct(parallelInformation_arg)); #endif - // Communicate if parallel. - parallelInformation_arg.copyOwnerToAll(istlb, istlb); - #if FLOW_SUPPORT_AMG // activate AMG if either flow_ebos is used or UMFPack is not available if( parameters_.linear_solver_use_amg_ || parameters_.use_cpr_) { @@ -644,29 +645,91 @@ protected: #endif } - /// Zero out off-diagonal blocks on rows corresponding to overlap cells - /// Diagonal blocks on ovelap rows are set to diag(1e100). - void makeOverlapRowsInvalid(Matrix& ebosJacIgnoreOverlap) const + /// Create sparsity pattern of matrix without off-diagonal ghost entries. + void noGhostAdjacency() { - //value to set on diagonal - Dune::FieldMatrix diag_block(0.0); - for (int eq = 0; eq < numEq; ++eq) - diag_block[eq][eq] = 1.0e100; + auto grid = simulator_.vanguard().grid(); + typedef typename Matrix::size_type size_type; + size_type numCells = grid.numCells(); + noGhostMat_.reset(new Matrix(numCells, numCells, Matrix::random)); - //loop over precalculated overlap rows and columns - for (auto row = overlapRowAndColumns_.begin(); row != overlapRowAndColumns_.end(); row++ ) + std::vector> pattern; + pattern.resize(grid.numCells()); + + auto lid = grid.localIdSet(); + const auto& gridView = grid.leafGridView(); + auto elemIt = gridView.template begin<0>(); + const auto& elemEndIt = gridView.template end<0>(); + + //Loop over cells + for (; elemIt != elemEndIt; ++elemIt) { - int lcell = row->first; - //diagonal block set to large value diagonal - ebosJacIgnoreOverlap[lcell][lcell] = diag_block; + const auto& elem = *elemIt; + size_type idx = lid.id(elem); + pattern[idx].insert(idx); - //loop over off diagonal blocks in overlap row - for (auto col = row->second.begin(); col != row->second.end(); ++col) + // Add well non-zero connections + for (auto wc = wellConnectionsGraph_[idx].begin(); wc!=wellConnectionsGraph_[idx].end(); ++wc) + pattern[idx].insert(*wc); + + // Add just a single element to ghost rows + if (elem.partitionType() != Dune::InteriorEntity) { - int ncell = *col; - //zero out block - ebosJacIgnoreOverlap[lcell][ncell] = 0.0; + noGhostMat_->setrowsize(idx, pattern[idx].size()); } + else { + 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()); + pattern[idx].insert(nid); + } + } + noGhostMat_->setrowsize(idx, pattern[idx].size()); + } + } + noGhostMat_->endrowsizes(); + for (size_type dofId = 0; dofId < numCells; ++dofId) + { + auto nabIdx = pattern[dofId].begin(); + auto endNab = pattern[dofId].end(); + for (; nabIdx != endNab; ++nabIdx) + { + noGhostMat_->addindex(dofId, *nabIdx); + } + } + noGhostMat_->endindices(); + } + + /// Set the ghost diagonal in noGhost to diag(1.0) + void setGhostsInNoGhost(Matrix& ng) + { + ng=0; + typedef typename Matrix::block_type MatrixBlockType; + MatrixBlockType diag_block(0.0); + for (int eq = 0; eq < Matrix::block_type::rows; ++eq) + diag_block[eq][eq] = 1.0; + + //loop over precalculated ghost rows and columns + for (auto row = overlapRows_.begin(); row != overlapRows_.end(); row++ ) + { + int lcell = *row; + //diagonal block set to 1 + ng[lcell][lcell] = diag_block; + } + } + + /// Copy interior rows to noghost matrix + void copyJacToNoGhost(const Matrix& jac, Matrix& ng) + { + //Loop over precalculated interior rows. + for (auto row = interiorRows_.begin(); row != interiorRows_.end(); row++ ) + { + //Copy row + ng[*row] = jac[*row]; } } @@ -864,10 +927,13 @@ protected: boost::any parallelInformation_; std::unique_ptr matrix_; + std::unique_ptr noGhostMat_; Vector *rhs_; std::unique_ptr matrix_for_preconditioner_; - std::vector>> overlapRowAndColumns_; + std::vector overlapRows_; + std::vector interiorRows_; + std::vector> wellConnectionsGraph_; FlowLinearSolverParameters parameters_; Vector weights_; bool scale_variables_; diff --git a/opm/simulators/linalg/ISTLSolverEbosCpr.hpp b/opm/simulators/linalg/ISTLSolverEbosCpr.hpp index ad3359727..dd99c7ce7 100644 --- a/opm/simulators/linalg/ISTLSolverEbosCpr.hpp +++ b/opm/simulators/linalg/ISTLSolverEbosCpr.hpp @@ -99,7 +99,7 @@ namespace Opm : SuperClass(simulator), oldMat() { extractParallelGridInformationToISTL(this->simulator_.vanguard().grid(), this->parallelInformation_); - detail::findOverlapRowsAndColumns(this->simulator_.vanguard().grid(), this->overlapRowAndColumns_); + detail::findOverlapAndInterior(this->simulator_.vanguard().grid(), this->overlapRows_, this->interiorRows_); } void prepare(const SparseMatrixAdapter& M, Vector& b) diff --git a/opm/simulators/linalg/ISTLSolverEbosFlexible.hpp b/opm/simulators/linalg/ISTLSolverEbosFlexible.hpp index 3403046ca..289780c9c 100644 --- a/opm/simulators/linalg/ISTLSolverEbosFlexible.hpp +++ b/opm/simulators/linalg/ISTLSolverEbosFlexible.hpp @@ -77,7 +77,7 @@ public: parameters_.template init(); prm_ = setupPropertyTree(parameters_); extractParallelGridInformationToISTL(simulator_.vanguard().grid(), parallelInformation_); - detail::findOverlapRowsAndColumns(simulator_.vanguard().grid(), overlapRowAndColumns_); + detail::findOverlapAndInterior(simulator_.vanguard().grid(), overlapRows_, interiorRows_); #if HAVE_MPI if (parallelInformation_.type() == typeid(ParallelISTLInformation)) { // Parallel case. @@ -176,27 +176,24 @@ public: protected: /// Zero out off-diagonal blocks on rows corresponding to overlap cells - /// Diagonal blocks on ovelap rows are set to diag(1e100). + /// Diagonal blocks on ovelap rows are set to diag(1.0). void makeOverlapRowsInvalid(MatrixType& matrix) const { - // Value to set on diagonal + //value to set on diagonal const int numEq = MatrixType::block_type::rows; typename MatrixType::block_type diag_block(0.0); for (int eq = 0; eq < numEq; ++eq) - diag_block[eq][eq] = 1.0e100; + diag_block[eq][eq] = 1.0; - // loop over precalculated overlap rows and columns - for (auto row = overlapRowAndColumns_.begin(); row != overlapRowAndColumns_.end(); row++) { - int lcell = row->first; - // diagonal block set to large value diagonal + //loop over precalculated overlap rows and columns + for (auto row = overlapRows_.begin(); row != overlapRows_.end(); row++ ) + { + int lcell = *row; + // Zero out row. + matrix[lcell] = 0.0; + + //diagonal block set to diag(1.0). matrix[lcell][lcell] = diag_block; - - // loop over off diagonal blocks in overlap row - for (auto col = row->second.begin(); col != row->second.end(); ++col) { - int ncell = *col; - // zero out block - matrix[lcell][ncell] = 0.0; - } } } @@ -212,7 +209,8 @@ protected: #if HAVE_MPI std::unique_ptr comm_; #endif - std::vector>> overlapRowAndColumns_; + std::vector overlapRows_; + std::vector interiorRows_; }; // end ISTLSolverEbosFlexible } // namespace Opm diff --git a/opm/simulators/linalg/ParallelOverlappingILU0.hpp b/opm/simulators/linalg/ParallelOverlappingILU0.hpp index 82975b5a6..de3048c81 100644 --- a/opm/simulators/linalg/ParallelOverlappingILU0.hpp +++ b/opm/simulators/linalg/ParallelOverlappingILU0.hpp @@ -750,7 +750,6 @@ public: { Range& md = reorderD(d); Domain& mv = reorderV(v); - copyOwnerToAll( md ); // iterator types typedef typename Range ::block_type dblock; @@ -778,8 +777,6 @@ public: mv[ i ] = rhs; // Lii = I } - copyOwnerToAll( mv ); - for( size_type i=0; i #include +#include namespace Opm { namespace detail { + /// \brief Find cell IDs for wells contained in local grid. + /// + /// Cell IDs of wells stored in a graph, so it can be used to create + /// an adjacency pattern. Only relevant when the UseWellContribusion option is set to true + /// \tparam The type of the DUNE grid. + /// \tparam Well vector type + /// \param grid The grid where we look for overlap cells. + /// \param wells List of wells contained in grid. + /// \param useWellConn Boolean that is true when UseWellContribusion is true + /// \param wellGraph Cell IDs of well cells stored in a graph. + template + void setWellConnections(const Grid& grid, const W& wells, bool useWellConn, std::vector>& wellGraph) + { + if ( grid.comm().size() > 1) + { + wellGraph.resize(grid.numCells()); + + if (useWellConn) { + const auto& cpgdim = grid.logicalCartesianSize(); + + std::vector cart(cpgdim[0]*cpgdim[1]*cpgdim[2], -1); + + for( int i=0; i < grid.numCells(); ++i ) + cart[grid.globalCell()[i]] = i; + + Dune::cpgrid::WellConnections well_indices; + well_indices.init(wells, cpgdim, cart); + + for (auto& well : well_indices) + { + for (auto perf = well.begin(); perf != well.end(); ++perf) + { + auto perf2 = perf; + for (++perf2; perf2 != well.end(); ++perf2) + { + wellGraph[*perf].insert(*perf2); + wellGraph[*perf2].insert(*perf); + } + } + } + } + } + } /// \brief Find the rows corresponding to overlap cells /// - /// Loop over grid and store cell ids of row-column pairs + /// Loop over grid and store cell ids of rows /// corresponding to overlap cells. /// \tparam The type of the DUNE grid. /// \param grid The grid where we look for overlap cells. - /// \param overlapRowAndColumns List where overlap rows and columns are stored. - template - void findOverlapRowsAndColumns(const Grid& grid, - std::vector>>& overlapRowAndColumns) + /// \param overlapRows List where overlap rows are stored. + /// \param interiorRows List where overlap rows are stored. + template + void findOverlapAndInterior(const Grid& grid, std::vector& overlapRows, + std::vector& interiorRows) { - // only relevant in parallel case. - if (grid.comm().size() > 1) { - // Numbering of cells + //only relevant in parallel case. + if ( grid.comm().size() > 1) + { + //Numbering of cells auto lid = grid.localIdSet(); const auto& gridView = grid.leafGridView(); auto elemIt = gridView.template begin<0>(); const auto& elemEndIt = gridView.template end<0>(); - // loop over cells in mesh - for (; elemIt != elemEndIt; ++elemIt) { + //loop over cells in mesh + for (; elemIt != elemEndIt; ++elemIt) + { const auto& elem = *elemIt; - - // If cell has partition type not equal to interior save row - if (elem.partitionType() != Dune::InteriorEntity) { - // local id of overlap cell - int lcell = lid.id(elem); - - std::vector columns; - // loop over faces of cell - auto isend = gridView.iend(elem); - for (auto is = gridView.ibegin(elem); is != isend; ++is) { - // check if face has neighbor - if (is->neighbor()) { - // get index of neighbor cell - int ncell = lid.id(is->outside()); - columns.push_back(ncell); - } - } - // add row to list - overlapRowAndColumns.push_back(std::pair>(lcell, columns)); + int lcell = lid.id(elem); + + if (elem.partitionType() != Dune::InteriorEntity) + { + //add row to list + overlapRows.push_back(lcell); + } else { + interiorRows.push_back(lcell); } } } } - } // namespace detail } // namespace Opm