Merge pull request #2205 from andrthu/adjacency-change

Remove non-contributing computations in parallel solver by changing matrix sparsity pattern.

opm/simulators/linalg/ISTLSolverEbos.hpp

@@ -240,7 +240,16 @@ protected:
         {
             parameters_.template init<TypeTag>();
             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<Scalar, numEq, numEq> 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<std::set<size_type>> 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> matrix_;
+        std::unique_ptr<Matrix> noGhostMat_;
         Vector *rhs_;
         std::unique_ptr<Matrix> matrix_for_preconditioner_;
 
-        std::vector<std::pair<int,std::vector<int>>> overlapRowAndColumns_;
+        std::vector<int> overlapRows_;
+        std::vector<int> interiorRows_;
+        std::vector<std::set<int>> wellConnectionsGraph_;
         FlowLinearSolverParameters parameters_;
         Vector weights_;
         bool scale_variables_;

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)

opm/simulators/linalg/ISTLSolverEbosFlexible.hpp

@@ -77,7 +77,7 @@ public:
         parameters_.template init<TypeTag>();
         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<Communication> comm_;
 #endif
-    std::vector<std::pair<int, std::vector<int>>> overlapRowAndColumns_;
+    std::vector<int> overlapRows_;
+    std::vector<int> interiorRows_;
 }; // end ISTLSolverEbosFlexible
 
 } // namespace Opm

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<iEnd; ++ i )
         {
             vblock& vBlock = mv[ lastRow - i ];

opm/simulators/linalg/findOverlapRowsAndColumns.hpp

@@ -22,59 +22,94 @@
 
 #include <vector>
 #include <utility>
+#include <opm/grid/common/WellConnections.hpp>
 
 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<class Grid, class W>
+    void setWellConnections(const Grid& grid, const W& wells, bool useWellConn, std::vector<std::set<int>>& wellGraph)
+    {
+        if ( grid.comm().size() > 1) 
+        {
+            wellGraph.resize(grid.numCells());
+
+            if (useWellConn) {
+                const auto& cpgdim = grid.logicalCartesianSize();
+
+                std::vector<int> 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 <class Grid>
-    void findOverlapRowsAndColumns(const Grid& grid,
-                                   std::vector<std::pair<int, std::vector<int>>>& overlapRowAndColumns)
+    /// \param overlapRows List where overlap rows are stored.
+    /// \param interiorRows List where overlap rows are stored.
+    template<class Grid>
+    void findOverlapAndInterior(const Grid& grid, std::vector<int>& overlapRows, 
+                                std::vector<int>& 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;
+                int lcell = lid.id(elem);
                 
-                // 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<int> 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<int, std::vector<int>>(lcell, columns));
+                if (elem.partitionType() != Dune::InteriorEntity)
+                {  
+                    //add row to list
+                    overlapRows.push_back(lcell);
+                } else {
+                    interiorRows.push_back(lcell);
                 }
             }
         }
     }
-
 } // namespace detail
 } // namespace Opm