Merge pull request #2220 from andrthu/owners-first

[For testing] Enable "owners first" ordering. Altered SpMV, SP and ILU.

ebos/eclbasevanguard.hh

@@ -78,6 +78,7 @@ NEW_PROP_TAG(EclOutputInterval);
 NEW_PROP_TAG(IgnoreKeywords);
 NEW_PROP_TAG(EnableExperiments);
 NEW_PROP_TAG(EdgeWeightsMethod);
+NEW_PROP_TAG(OwnerCellsFirst);
 
 SET_STRING_PROP(EclBaseVanguard, IgnoreKeywords, "");
 SET_STRING_PROP(EclBaseVanguard, EclDeckFileName, "");
@@ -86,6 +87,7 @@ SET_BOOL_PROP(EclBaseVanguard, EnableOpmRstFile, false);
 SET_BOOL_PROP(EclBaseVanguard, EclStrictParsing, false);
 SET_BOOL_PROP(EclBaseVanguard, SchedRestart, true);
 SET_INT_PROP(EclBaseVanguard, EdgeWeightsMethod, 1);
+SET_BOOL_PROP(EclBaseVanguard, OwnerCellsFirst, false);
 
 END_PROPERTIES
 
@@ -133,6 +135,8 @@ public:
                              "When restarting: should we try to initialize wells and groups from historical SCHEDULE section.");
         EWOMS_REGISTER_PARAM(TypeTag, int, EdgeWeightsMethod,
                              "Choose edge-weighing strategy: 0=uniform, 1=trans, 2=log(trans).");
+        EWOMS_REGISTER_PARAM(TypeTag, bool, OwnerCellsFirst,
+                             "Order cells owned by rank before ghost/overlap cells.");
     }
 
     /*!
@@ -265,6 +269,7 @@ public:
 
         std::string fileName = EWOMS_GET_PARAM(TypeTag, std::string, EclDeckFileName);
         edgeWeightsMethod_   = Dune::EdgeWeightMethod(EWOMS_GET_PARAM(TypeTag, int, EdgeWeightsMethod));
+        ownersFirst_ = EWOMS_GET_PARAM(TypeTag, bool, OwnerCellsFirst);
 
         // Skip processing of filename if external deck already exists.
         if (!externalDeck_)
@@ -451,6 +456,13 @@ public:
      */
     Dune::EdgeWeightMethod edgeWeightsMethod() const
     { return edgeWeightsMethod_; }
+
+    /*!
+     * \brief Parameter that decide if cells owned by rank are ordered before ghost cells.
+     */
+    bool ownersFirst() const
+    { return ownersFirst_; }
+    
     /*!
      * \brief Returns the name of the case.
      *
@@ -619,6 +631,7 @@ private:
     Opm::SummaryConfig* eclSummaryConfig_;
 
     Dune::EdgeWeightMethod edgeWeightsMethod_;
+    bool ownersFirst_;
 
 protected:
     /*! \brief The cell centroids after loadbalance was called.

ebos/eclcpgridvanguard.hh

@@ -159,6 +159,7 @@ public:
             }
 
             Dune::EdgeWeightMethod edgeWeightsMethod = this->edgeWeightsMethod();
+            bool ownersFirst = this->ownersFirst();
 
             // convert to transmissibility for faces
             // TODO: grid_->numFaces() is not generic. use grid_->size(1) instead? (might
@@ -204,7 +205,7 @@ public:
 
                     PropsCentroidsDataHandle<Dune::CpGrid> handle(*grid_, eclState, eclGrid, this->centroids_,
                                                                   cartesianIndexMapper());
-                    defunctWellNames_ = std::get<1>(grid_->loadBalance(handle, edgeWeightsMethod, &wells, faceTrans.data()));
+                    defunctWellNames_ = std::get<1>(grid_->loadBalance(handle, edgeWeightsMethod, &wells, faceTrans.data(), ownersFirst));
                 }
                 catch(const std::bad_cast& e)
                 {

ebos/eclthresholdpressure.hh

@@ -309,8 +309,9 @@ private:
                 // store intersection, this might be costly
                 const auto& intersection = *isIt;
 
-                // ignore boundary intersections for now (TODO?)
                 if (intersection.boundary())
+                    continue; // ignore boundary intersections for now (TODO?)
+                else if (!intersection.neighbor()) //processor boundary but not domain boundary
                     continue;
 
                 const auto& inside = intersection.inside();

opm/simulators/linalg/ISTLSolverEbos.hpp

@@ -185,6 +185,108 @@ protected:
   std::shared_ptr< communication_type > comm_;
 };
 
+
+/*!
+   \brief Adapter to turn a matrix into a linear operator.
+   Adapts a matrix to the assembled linear operator interface.
+   We assume parallel ordering, where ghost rows are located after interior rows
+ */
+template<class M, class X, class Y, class WellModel, bool overlapping >
+class WellModelGhostLastMatrixAdapter : public Dune::AssembledLinearOperator<M,X,Y>
+{
+    typedef Dune::AssembledLinearOperator<M,X,Y> BaseType;
+
+public:
+    typedef M matrix_type;
+    typedef X domain_type;
+    typedef Y range_type;
+    typedef typename X::field_type field_type;
+
+#if HAVE_MPI
+    typedef Dune::OwnerOverlapCopyCommunication<int,int> communication_type;
+#else
+    typedef Dune::CollectiveCommunication< int > communication_type;
+#endif
+
+
+    Dune::SolverCategory::Category category() const override
+    {
+        return overlapping ?
+            Dune::SolverCategory::overlapping : Dune::SolverCategory::sequential;
+    }
+    
+    //! constructor: just store a reference to a matrix
+    WellModelGhostLastMatrixAdapter (const M& A,
+                                     const M& A_for_precond,
+                                     const WellModel& wellMod,
+                                     const size_t interiorSize,
+                                     const std::any& parallelInformation OPM_UNUSED_NOMPI = std::any() )
+        : A_( A ), A_for_precond_(A_for_precond), wellMod_( wellMod ), interiorSize_(interiorSize), comm_()
+    {
+#if HAVE_MPI
+        if( parallelInformation.type() == typeid(ParallelISTLInformation) )
+        {
+            const ParallelISTLInformation& info =
+                std::any_cast<const ParallelISTLInformation&>( parallelInformation);
+            comm_.reset( new communication_type( info.communicator() ) );
+        }
+#endif
+    }
+    
+    virtual void apply( const X& x, Y& y ) const
+    {
+        for (auto row = A_.begin(); row.index() < interiorSize_; ++row)
+        {
+            y[row.index()]=0;
+            auto endc = (*row).end();
+            for (auto col = (*row).begin(); col != endc; ++col)
+                (*col).umv(x[col.index()], y[row.index()]);
+        }
+        
+        // add well model modification to y
+        wellMod_.apply(x, y );
+
+        ghostLastProject( y );
+    }
+
+    // y += \alpha * A * x
+    virtual void applyscaleadd (field_type alpha, const X& x, Y& y) const
+    {
+        for (auto row = A_.begin(); row.index() < interiorSize_; ++row)
+        {
+            auto endc = (*row).end();
+            for (auto col = (*row).begin(); col != endc; ++col)
+                (*col).usmv(alpha, x[col.index()], y[row.index()]);
+        }
+        // add scaled well model modification to y
+        wellMod_.applyScaleAdd( alpha, x, y );
+        
+        ghostLastProject( y );
+    }
+
+    virtual const matrix_type& getmat() const { return A_for_precond_; }
+
+    communication_type* comm()
+    {
+        return comm_.operator->();
+    }
+    
+protected:
+    void ghostLastProject(Y& y) const
+    {
+        size_t end = y.size();
+        for (size_t i = interiorSize_; i < end; ++i)
+            y[i] = 0;
+    }
+    
+    const matrix_type& A_ ;
+    const matrix_type& A_for_precond_ ;
+    const WellModel& wellMod_;
+    size_t interiorSize_;
+    
+    std::unique_ptr< communication_type > comm_;
+};
+
     /// This class solves the fully implicit black-oil system by
     /// solving the reduced system (after eliminating well variables)
     /// as a block-structured matrix (one block for all cell variables) for a fixed
@@ -257,11 +359,19 @@ protected:
             const auto wellsForConn = simulator_.vanguard().schedule().getWellsatEnd();
             const bool useWellConn = EWOMS_GET_PARAM(TypeTag, bool, MatrixAddWellContributions);
 
-            detail::setWellConnections(gridForConn, wellsForConn, useWellConn, wellConnectionsGraph_);
+            ownersFirst_ = EWOMS_GET_PARAM(TypeTag, bool, OwnerCellsFirst);
-            detail::findOverlapAndInterior(gridForConn, overlapRows_, interiorRows_);
+            interiorCellNum_ = detail::numMatrixRowsToUseInSolver(simulator_.vanguard().grid(), ownersFirst_);
-            if (gridForConn.comm().size() > 1) {
+
-                noGhostAdjacency();
+            if (!ownersFirst_ || parameters_.linear_solver_use_amg_  || parameters_.use_cpr_ ) {
-                setGhostsInNoGhost(*noGhostMat_);
+                detail::setWellConnections(gridForConn, wellsForConn, useWellConn, wellConnectionsGraph_);
+                detail::findOverlapAndInterior(gridForConn, overlapRows_, interiorRows_);
+                if (gridForConn.comm().size() > 1) {
+
+                    noGhostAdjacency();
+                    setGhostsInNoGhost(*noGhostMat_);
+                }
+                if (ownersFirst_)
+                    OpmLog::warning("OwnerCellsFirst option is true, but ignored.");
             }
         }
 
@@ -342,13 +452,26 @@ protected:
 
             if( isParallel() )
             {
-                typedef WellModelMatrixAdapter< Matrix, Vector, Vector, WellModel, true > Operator;
+                if ( ownersFirst_ && (!parameters_.linear_solver_use_amg_ || !parameters_.use_cpr_) ) {
+                    typedef WellModelGhostLastMatrixAdapter< Matrix, Vector, Vector, WellModel, true > Operator;
+                    Operator opA(*matrix_, *matrix_, wellModel, interiorCellNum_,
+                                 parallelInformation_ );
                     
-                copyJacToNoGhost(*matrix_, *noGhostMat_);
+                    assert( opA.comm() );
-                Operator opA(*noGhostMat_, *noGhostMat_, wellModel,
+                    solve( opA, x, *rhs_, *(opA.comm()) );
-                             parallelInformation_ );
+                }
-                assert( opA.comm() );
+                else {
-                solve( opA, x, *rhs_, *(opA.comm()) );
+
+                    typedef WellModelMatrixAdapter< Matrix, Vector, Vector, WellModel, true > Operator;
+
+                    copyJacToNoGhost(*matrix_, *noGhostMat_);
+                    Operator opA(*noGhostMat_, *noGhostMat_, wellModel,
+                                 parallelInformation_ );
+    
+                    assert( opA.comm() );
+                    solve( opA, x, *rhs_, *(opA.comm()) );
+                    
+                }
             }
             else
             {
@@ -502,7 +625,7 @@ protected:
             const MILU_VARIANT ilu_milu  = parameters_.ilu_milu_;
             const bool ilu_redblack = parameters_.ilu_redblack_;
             const bool ilu_reorder_spheres = parameters_.ilu_reorder_sphere_;
-            return Pointer(new ParPreconditioner(opA.getmat(), comm, relax, ilu_milu, ilu_redblack, ilu_reorder_spheres));
+            return Pointer(new ParPreconditioner(opA.getmat(), comm, relax, ilu_milu, interiorCellNum_, ilu_redblack, ilu_reorder_spheres));
         }
 #endif
 
@@ -949,6 +1072,10 @@ protected:
         std::vector<int> overlapRows_;
         std::vector<int> interiorRows_;
         std::vector<std::set<int>> wellConnectionsGraph_;
+
+        bool ownersFirst_;
+        size_t interiorCellNum_;
+
         FlowLinearSolverParameters parameters_;
         Vector weights_;
         bool scale_variables_;

opm/simulators/linalg/ParallelOverlappingILU0.hpp

@@ -435,6 +435,64 @@ namespace Opm
         }
     }
 
+    //! Compute Blocked ILU0 decomposition, when we know junk ghost rows are located at the end of A
+    template<class M>
+    void ghost_last_bilu0_decomposition (M& A, size_t interiorSize)
+    {
+        // iterator types
+        typedef typename M::RowIterator rowiterator;
+        typedef typename M::ColIterator coliterator;
+        typedef typename M::block_type block;
+
+        // implement left looking variant with stored inverse
+        for (rowiterator i = A.begin(); i.index() < interiorSize; ++i)
+        {
+            // coliterator is diagonal after the following loop
+            coliterator endij=(*i).end();           // end of row i
+            coliterator ij;
+
+            // eliminate entries left of diagonal; store L factor
+            for (ij=(*i).begin(); ij.index()<i.index(); ++ij)
+            {
+                // find A_jj which eliminates A_ij
+                coliterator jj = A[ij.index()].find(ij.index());
+                
+                // compute L_ij = A_jj^-1 * A_ij
+                (*ij).rightmultiply(*jj);
+
+                // modify row
+                coliterator endjk=A[ij.index()].end();    // end of row j
+                coliterator jk=jj; ++jk;
+                coliterator ik=ij; ++ik;
+                while (ik!=endij && jk!=endjk)
+                    if (ik.index()==jk.index())
+                    {
+                        block B(*jk);
+                        B.leftmultiply(*ij);
+                        *ik -= B;
+                        ++ik; ++jk;
+                    }
+                    else
+                    {
+                        if (ik.index()<jk.index())
+                            ++ik;
+                        else
+                            ++jk;
+                    }
+            }
+
+            // invert pivot and store it in A
+            if (ij.index()!=i.index())
+                DUNE_THROW(Dune::ISTLError,"diagonal entry missing");
+            try {
+                (*ij).invert();   // compute inverse of diagonal block
+            }
+            catch (Dune::FMatrixError & e) {
+                DUNE_THROW(Dune::ISTLError,"ILU failed to invert matrix block");
+            }
+        }
+    }
+
       //! compute ILU decomposition of A. A is overwritten by its decomposition
       template<class M, class CRS, class InvVector>
       void convertToCRS(const M& A, CRS& lower, CRS& upper, InvVector& inv )
@@ -509,7 +567,7 @@ namespace Opm
             if( j.index() == iIndex )
             {
               inv[ row ] = (*j);
-	      break;
+              break;
             }
             else if ( j.index() >= i.index() )
             {
@@ -635,6 +693,7 @@ public:
           comm_(nullptr), w_(w),
           relaxation_( std::abs( w - 1.0 ) > 1e-15 )
     {
+        interiorSize_ = A.N();
         // BlockMatrix is a Subclass of FieldMatrix that just adds
         // methods. Therefore this cast should be safe.
         init( reinterpret_cast<const Matrix&>(A), n, milu, redblack,
@@ -664,6 +723,7 @@ public:
           comm_(&comm), w_(w),
           relaxation_( std::abs( w - 1.0 ) > 1e-15 )
     {
+        interiorSize_ = A.N();
         // BlockMatrix is a Subclass of FieldMatrix that just adds
         // methods. Therefore this cast should be safe.
         init( reinterpret_cast<const Matrix&>(A), n, milu, redblack,
@@ -714,12 +774,46 @@ public:
           comm_(&comm), w_(w),
           relaxation_( std::abs( w - 1.0 ) > 1e-15 )
     {
+        interiorSize_ = A.N();
         // BlockMatrix is a Subclass of FieldMatrix that just adds
         // methods. Therefore this cast should be safe.
         init( reinterpret_cast<const Matrix&>(A), 0, milu, redblack,
               reorder_sphere );
     }
 
+    /*! \brief Constructor.
+
+      Constructor gets all parameters to operate the prec.
+      \param A The matrix to operate on.
+      \param n ILU fill in level (for testing). This does not work in parallel.
+      \param w The relaxation factor.
+      \param milu The modified ILU variant to use. 0 means traditional ILU. \see MILU_VARIANT.
+      \param interiorSize The number of interior/owner rows in the matrix.
+      \param redblack Whether to use a red-black ordering.
+      \param reorder_sphere If true, we start the reordering at a root node.
+                            The vertices on each layer aound it (same distance) are
+                            ordered consecutivly. If false, we preserver the order of
+                            the vertices with the same color.
+    */
+    template<class BlockType, class Alloc>
+    ParallelOverlappingILU0 (const Dune::BCRSMatrix<BlockType,Alloc>& A,
+                             const ParallelInfo& comm,
+                             const field_type w, MILU_VARIANT milu,
+                             size_type interiorSize, bool redblack=false,
+                             bool reorder_sphere=true)
+        : lower_(),
+          upper_(),
+          inv_(),
+          comm_(&comm), w_(w),
+          relaxation_( std::abs( w - 1.0 ) > 1e-15 ),
+          interiorSize_(interiorSize)
+    {
+        // BlockMatrix is a Subclass of FieldMatrix that just adds
+        // methods. Therefore this cast should be safe.
+        init( reinterpret_cast<const Matrix&>(A), 0, milu, redblack,
+              reorder_sphere  );
+    }
+
     /*!
       \brief Prepare the preconditioner.
 
@@ -747,13 +841,15 @@ public:
 
         const size_type iEnd = lower_.rows();
         const size_type lastRow = iEnd - 1;
+        size_type upperLoppStart = iEnd - interiorSize_;
+        size_type lowerLoopEnd = interiorSize_;
         if( iEnd != upper_.rows() )
         {
             OPM_THROW(std::logic_error,"ILU: number of lower and upper rows must be the same");
         }
 
         // lower triangular solve
-        for( size_type i=0; i<iEnd; ++ i )
+        for( size_type i=0; i<lowerLoopEnd; ++ i )
         {
           dblock rhs( md[ i ] );
           const size_type rowI     = lower_.rows_[ i ];
@@ -767,7 +863,7 @@ public:
           mv[ i ] = rhs;  // Lii = I
         }
 
-        for( size_type i=0; i<iEnd; ++ i )
+        for( size_type i=upperLoppStart; i<iEnd; ++ i )
         {
             vblock& vBlock = mv[ lastRow - i ];
             vblock rhs ( vBlock );
@@ -911,7 +1007,10 @@ protected:
                                                   detail::IsPositiveFunctor() );
                     break;
                 default:
-                    bilu0_decomposition( *ILU );
+                    if (interiorSize_ == A.N())
+                        bilu0_decomposition( *ILU );
+                    else
+                        detail::ghost_last_bilu0_decomposition(*ILU, interiorSize_);
                     break;
                 }
             }
@@ -1023,7 +1122,7 @@ protected:
     //! \brief The relaxation factor to use.
     const field_type w_;
     const bool relaxation_;
-
+    size_type interiorSize_;
 };
 
 } // end namespace Opm

opm/simulators/linalg/findOverlapRowsAndColumns.hpp

@@ -112,6 +112,33 @@ namespace detail
             }
         }
     }
+
+    /// \brief If ownerFirst=true, returns the number of interior cells in grid, else just numCells().
+    ///
+    /// If cells in grid is ordered so that interior/owner cells come before overlap/copy cells, the method
+    /// returns the number of interior cells numInterior. In the linear solver only the first numInterior rows of
+    /// the matrix are needed.
+    template <class Grid>
+    size_t numMatrixRowsToUseInSolver(const Grid& grid, bool ownerFirst)
+    {
+        size_t numInterior = 0;
+        if (!ownerFirst || grid.comm().size()==1)
+            return grid.numCells();
+        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) {
+
+            // Count only the interior cells.
+            if (elemIt->partitionType() == Dune::InteriorEntity) {
+                numInterior++;
+            }
+        }
+
+        return numInterior;
+    }
 } // namespace detail
 } // namespace Opm