Ignoring overlap cells in parallel ILU factorization (#1610)

* fixed the issue of including ghost cells in preconditioning, by zeroing out ghost rows and setting the diagonal to a large value.

* move altering of matrix to a function

* blockwise modification of matrix

* add findOverlapRowsAndColumns in BlacoilDetail. Add call to findOverlapRowsAndColumns in constructor of BlacoilModelEbos. Change makeOverlapRowsInvalid, by looping over precalculated inddies instead of grid

* Better formatting

opm/autodiff/BlackoilDetails.hpp

@@ -196,7 +196,56 @@ namespace detail {
             return grid.comm().sum(count);
         }
 
+        /// \brief Find the rows corresponding to overlap cells
+        ///
+        /// Loop over grid and store cell ids of row-column pairs 
+        /// 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 )
+        {
+            //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) 
+                {		
+                    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<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));	
+                    }    
+                }	      
+            }
+        }				       
     } // namespace detail
 } // namespace Opm
 

opm/autodiff/BlackoilModelEbos.hpp

@@ -166,6 +166,8 @@ namespace Opm {
         {
             // compute global sum of number of cells
             global_nc_ = detail::countGlobalCells(grid_);
+            //find rows of matrix corresponding to overlap
+            detail::findOverlapRowsAndColumns(grid_,overlapRowAndColumns_);
             if (!istlSolver_)
             {
                 OPM_THROW(std::logic_error,"solver down cast to ISTLSolver failed");
@@ -472,6 +474,32 @@ namespace Opm {
             return istlSolver().iterations();
         }
 
+        /// Zero out off-diagonal blocks on rows corresponding to overlap cells
+        /// Diagonal blocks on ovelap rows are set to diag(1e100).
+        void makeOverlapRowsInvalid(Mat& ebosJacIgnoreOverlap) const
+        {	  	  
+            //value to set on diagonal
+            MatrixBlockType diag_block(0.0);
+            for (int eq = 0; eq < numEq; ++eq)	    
+                diag_block[eq][eq] = 1.0e100;
+	  	  
+            //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
+                ebosJacIgnoreOverlap[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
+                    ebosJacIgnoreOverlap[lcell][ncell] = 0.0;
+                }
+            }    	  
+        }
+
         /// Solve the Jacobian system Jx = r where J is the Jacobian and
         /// r is the residual.
         void solveJacobianSystem(BVector& x) const
@@ -497,9 +525,16 @@ namespace Opm {
             const Mat& actual_mat_for_prec = matrix_for_preconditioner_ ? *matrix_for_preconditioner_.get() : ebosJac;
             // Solve system.
             if( isParallel() )
-            {
+            {	      
                 typedef WellModelMatrixAdapter< Mat, BVector, BVector, BlackoilWellModel<TypeTag>, true > Operator;
-                Operator opA(ebosJac, actual_mat_for_prec, wellModel(),
+
+                auto ebosJacIgnoreOverlap = Mat(ebosJac);
+                //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(),
                              istlSolver().parallelInformation() );
                 assert( opA.comm() );
                 istlSolver().solve( opA, x, ebosResid, *(opA.comm()) );
@@ -956,8 +991,8 @@ namespace Opm {
         double current_relaxation_;
         BVector dx_old_;
 
-        std::unique_ptr<Mat> matrix_for_preconditioner_;
+        std::unique_ptr<Mat> matrix_for_preconditioner_;        
-
+        std::vector<std::pair<int,std::vector<int>>> overlapRowAndColumns_;
     public:
         /// return the StandardWells object
         BlackoilWellModel<TypeTag>&