New owners-first based linear algebra operations (SoMV, SP and ILU.apply).

ebos/eclthresholdpressure.hh

@@ -310,7 +310,7 @@ private:
                 const auto& intersection = *isIt;
 
                 // ignore boundary intersections for now (TODO?)
-                if (intersection.boundary())
+                if (!intersection.neighbor())
                     continue;
 
                 const auto& inside = intersection.inside();

opm/simulators/linalg/ISTLSolverEbos.hpp

@@ -185,6 +185,198 @@ 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
+
+#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6)
+    Dune::SolverCategory::Category category() const override
+    {
+        return overlapping ?
+            Dune::SolverCategory::overlapping : Dune::SolverCategory::sequential;
+    }
+#else
+    enum {
+        //! \brief The solver category.
+        category = overlapping ?
+        Dune::SolverCategory::overlapping :
+        Dune::SolverCategory::sequential
+    };
+#endif
+    
+    //! 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 boost::any& parallelInformation = boost::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 =
+                boost::any_cast<const ParallelISTLInformation&>( parallelInformation);
+            comm_.reset( new communication_type( info.communicator() ) );
+        }
+#endif
+    }
+    
+    virtual void apply( const X& x, Y& y ) const
+    {
+        unsigned row_count = 0;
+        for (auto row = A_.begin(); row_count < interiorSize_; ++row, ++row_count)
+        {
+            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
+    {
+        //auto first_row = A_.begin();
+        unsigned row_count = 0;
+        for (auto row = A_.begin(); row_count < interiorSize_; ++row, ++row_count)
+        {
+            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_;
+};
+
+/*!
+  \brief parallel ScalarProduct, that assumes ghost cells are ordered last 
+  Implements the Dune ScalarProduct interface
+ */
+template<class X, class C>
+class GhostLastScalarProduct : public Dune::ScalarProduct<X>
+{
+public:
+    typedef X domain_type;
+    typedef typename X::field_type field_type;
+    typedef typename Dune::FieldTraits<field_type>::real_type real_type;
+    typedef typename X::size_type size_type;
+
+    typedef C communication_type;
+    
+    enum {
+        category = Dune::SolverCategory::overlapping
+    };
+
+    GhostLastScalarProduct (const communication_type& com, size_type interiorSize)
+        : interiorSize_(interiorSize)
+    {        
+        comm_.reset(new communication_type(com.communicator()) );
+    }
+    virtual field_type dot (const X& x, const X& y)
+    {
+        field_type result = field_type(0.0);
+        for (size_type i = 0; i < interiorSize_; ++i)
+            result += x[i]*(y[i]);
+        return comm_->communicator().sum(result);
+    }
+    
+    virtual real_type norm (const X& x)
+    {
+        field_type result = field_type(0.0);
+        for (size_type i = 0; i < interiorSize_; ++i)
+            result += x[i].two_norm2();
+        return static_cast<double>(std::sqrt(comm_->communicator().sum(result)));
+    }
+    
+private:
+    std::unique_ptr< communication_type > comm_;
+    size_type interiorSize_;
+};
+
+template <class  X, class C, int c>
+struct GhostLastSPChooser
+{
+    typedef C communication_type;
+
+    enum { solverCategory = c};
+};
+
+template <class  X, class C>
+struct GhostLastSPChooser<X,C,Dune::SolverCategory::sequential>
+{
+    typedef Dune::SeqScalarProduct<X> ScalarProduct;
+    enum { category = Dune::SolverCategory::sequential};
+
+    static ScalarProduct* construct(const C& comm, size_t is)
+    {
+        return new ScalarProduct();
+    }
+};
+
+template <class  X, class C>
+struct GhostLastSPChooser<X,C,Dune::SolverCategory::overlapping>
+{
+    typedef  GhostLastScalarProduct<X,C> ScalarProduct;
+    typedef C communication_type;
+
+    enum {category = Dune::SolverCategory::overlapping};
+
+    static ScalarProduct* construct(const communication_type& comm, size_t is)
+    {
+        return new ScalarProduct(comm, is);
+    }
+};
+
     /// 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
@@ -253,15 +445,24 @@ protected:
             }
 #endif
             extractParallelGridInformationToISTL(simulator_.vanguard().grid(), parallelInformation_);
+
             const 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_);
+            ownersFirst_ = EWOMS_GET_PARAM(TypeTag, bool, OwnersFirst);
+            interiorCellNum_ = detail::numInteriorCells(simulator_.vanguard().grid(), ownersFirst_);
+
+            if (!ownersFirst_) {
+                detail::setWellConnections(gridForConn, wellsForConn, useWellConn, wellConnectionsGraph_);
+                detail::findOverlapAndInterior(gridForConn, overlapRows_, interiorRows_);
+
+                if (gridForConn.comm().size() > 1) {
+
+                    noGhostAdjacency();
+                    setGhostsInNoGhost(*noGhostMat_);
+                }
             }
         }
 
@@ -342,13 +543,24 @@ protected:
 
             if( isParallel() )
             {
-                typedef WellModelMatrixAdapter< Matrix, Vector, Vector, WellModel, true > Operator;
+                if ( ownersFirst_ ) {
+                    typedef WellModelGhostLastMatrixAdapter< Matrix, Vector, Vector, WellModel, true > Operator;
+                    Operator opA(*matrix_, *matrix_, wellModel, interiorCellNum_,
+                                 parallelInformation_ );
+                    
+                    assert( opA.comm() );
+                    solve( opA, x, *rhs_, *(opA.comm()) );
+                }
+                else {
+                    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()) );
+                    copyJacToNoGhost(*matrix_, *noGhostMat_);
+                    Operator opA(*noGhostMat_, *noGhostMat_, wellModel,
+                                 parallelInformation_ );
+                    assert( opA.comm() );
+                    solve( opA, x, *rhs_, *(opA.comm()) );
+                    
+                }
             }
             else
             {
@@ -388,12 +600,13 @@ protected:
                                              const POrComm& parallelInformation_arg,
                                              Dune::InverseOperatorResult& result) const
         {
+
             // Construct scalar product.
             auto sp = Dune::createScalarProduct<Vector,POrComm>(parallelInformation_arg, category);
 
 #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_)
-            {
+            {                
                 typedef ISTLUtility::CPRSelector< Matrix, Vector, Vector, POrComm>  CPRSelectorType;
                 typedef typename CPRSelectorType::Operator MatrixOperator;
 
@@ -457,12 +670,13 @@ protected:
                     solve(linearOperator, x, istlb, *sp, *precond, result);
                 } // end Dune call
 #else
-                // Construct preconditioner.
-                auto precond = constructPrecond(linearOperator, parallelInformation_arg);
+		// Construct preconditioner.
+		auto precond = constructPrecond(linearOperator, parallelInformation_arg);
 
-                // Solve.
-                solve(linearOperator, x, istlb, *sp, *precond, result);
+		// Solve.
+		solve(linearOperator, x, istlb, *sp, *precond, result);
 #endif
+
             }
         }
 
@@ -492,7 +706,9 @@ protected:
         typedef Dune::OwnerOverlapCopyCommunication<int, int> Comm;
         // 3x3 matrix block inversion was unstable from at least 2.3 until and
         // including 2.5.0
-        typedef ParallelOverlappingILU0<Matrix,Vector,Vector,Comm> ParPreconditioner;
+
+        typedef GhostLastParallelOverlappingILU0<Matrix,Vector,Vector,Comm> ParPreconditioner;
+
         template <class Operator>
         std::unique_ptr<ParPreconditioner>
         constructPrecond(Operator& opA, const Comm& comm) const
@@ -502,7 +718,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
 
@@ -946,9 +1162,14 @@ protected:
         Vector *rhs_;
         std::unique_ptr<Matrix> matrix_for_preconditioner_;
 
+
         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,65 @@ 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
+        size_t row_count = 0;
+        for (rowiterator i = A.begin(); row_count < interiorSize; ++i, ++row_count)
+        {
+            // 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 +568,7 @@ namespace Opm
             if( j.index() == iIndex )
             {
               inv[ row ] = (*j);
-	      break;
+              break;
             }
             else if ( j.index() >= i.index() )
             {
@@ -1026,5 +1085,365 @@ protected:
 
 };
 
+
+/// \brief A two-step version of an overlapping Schwarz preconditioner using one step ILU0 as
+///
+/// This preconditioner differs from a ParallelRestrictedOverlappingSchwarz with
+/// Dune:SeqILU0 in the follwing way:
+/// During apply we make sure that the current residual is consistent (i.e.
+/// each process knows the same value for each index. The we solve
+/// Ly= d for y and make y consistent again. Last we solve Ux = y and
+/// make sure that x is consistent.
+/// In contrast for ParallelRestrictedOverlappingSchwarz we solve (LU)x = d for x
+/// without forcing consistency between the two steps.
+/// \tparam Matrix The type of the Matrix.
+/// \tparam Domain The type of the Vector representing the domain.
+/// \tparam Range The type of the Vector representing the range.
+/// \tparam ParallelInfo The type of the parallel information object
+///         used, e.g. Dune::OwnerOverlapCommunication
+template<class Matrix, class Domain, class Range, class ParallelInfoT>
+class GhostLastParallelOverlappingILU0
+    : public Dune::Preconditioner<Domain,Range>
+{
+    typedef ParallelInfoT ParallelInfo;
+
+
+public:
+    //! \brief The matrix type the preconditioner is for.
+    typedef typename std::remove_const<Matrix>::type matrix_type;
+    //! \brief The domain type of the preconditioner.
+    typedef Domain domain_type;
+    //! \brief The range type of the preconditioner.
+    typedef Range range_type;
+    //! \brief The field type of the preconditioner.
+    typedef typename Domain::field_type field_type;
+
+    typedef typename matrix_type::block_type  block_type;
+    typedef typename matrix_type::size_type   size_type;
+
+protected:
+    struct CRS
+    {
+      CRS() : nRows_( 0 ) {}
+
+      size_type rows() const { return nRows_; }
+
+      size_type nonZeros() const
+      {
+        assert( rows_[ rows() ] != size_type(-1) );
+        return rows_[ rows() ];
+      }
+
+      void resize( const size_type nRows )
+      {
+          if( nRows_ != nRows )
+          {
+            nRows_ = nRows ;
+            rows_.resize( nRows_+1, size_type(-1) );
+          }
+      }
+
+      void reserveAdditional( const size_type nonZeros )
+      {
+          const size_type needed = values_.size() + nonZeros ;
+          if( values_.capacity() < needed )
+          {
+              const size_type estimate = needed * 1.1;
+              values_.reserve( estimate );
+              cols_.reserve( estimate );
+          }
+      }
+
+      void push_back( const block_type& value, const size_type index )
+      {
+          values_.push_back( value );
+          cols_.push_back( index );
+      }
+
+      std::vector< size_type  > rows_;
+      std::vector< block_type > values_;
+      std::vector< size_type  > cols_;
+      size_type nRows_;
+    };
+
+public:
+#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6)
+    Dune::SolverCategory::Category category() const override
+    {
+      return std::is_same<ParallelInfoT, Dune::Amg::SequentialInformation>::value ?
+              Dune::SolverCategory::sequential : Dune::SolverCategory::overlapping;
+    }
+
+#else
+    // define the category
+    enum {
+        //! \brief The category the preconditioner is part of.
+        category = std::is_same<ParallelInfoT, Dune::Amg::SequentialInformation>::value ?
+            Dune::SolverCategory::sequential : Dune::SolverCategory::overlapping
+    };
+#endif
+
+    /*! \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 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>
+    GhostLastParallelOverlappingILU0 (const Dune::BCRSMatrix<BlockType,Alloc>& A,
+                                      const int n, const field_type w,
+                                      MILU_VARIANT milu, unsigned interiorSize, bool redblack=false,
+                                      bool reorder_sphere=true)
+        : lower_(),
+          upper_(),
+          inv_(),
+          comm_(nullptr), 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), n, milu, redblack,
+              reorder_sphere  );
+    }
+
+    /*! \brief Constructor gets all parameters to operate the prec.
+      \param A The matrix to operate on.
+      \param comm   communication object, e.g. Dune::OwnerOverlapCopyCommunication
+      \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 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>
+    GhostLastParallelOverlappingILU0 (const Dune::BCRSMatrix<BlockType,Alloc>& A,
+                                      const ParallelInfo& comm, const int n, const field_type w,
+                                      MILU_VARIANT milu, unsigned 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), n, milu, redblack,
+              reorder_sphere );
+    }
+
+    /*! \brief Constructor.
+      Constructor gets all parameters to operate the prec.
+      \param A The matrix to operate on.
+      \param w The relaxation factor.
+      \param milu The modified ILU variant to use. 0 means traditional ILU. \see MILU_VARIANT.
+      \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>
+    GhostLastParallelOverlappingILU0 (const Dune::BCRSMatrix<BlockType,Alloc>& A,
+                                      const field_type w, MILU_VARIANT milu, unsigned interiorSize,
+                                      bool redblack=false,
+                                      bool reorder_sphere=true)
+        : GhostLastParallelOverlappingILU0( A, 0, w, milu, interiorSize, redblack, reorder_sphere )
+    {
+    }
+
+    /*! \brief Constructor.
+      Constructor gets all parameters to operate the prec.
+      \param A      The matrix to operate on.
+      \param comm   communication object, e.g. Dune::OwnerOverlapCopyCommunication
+      \param w      The relaxation factor.
+      \param milu   The modified ILU variant to use. 0 means traditional ILU. \see MILU_VARIANT.
+      \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>
+    GhostLastParallelOverlappingILU0 (const Dune::BCRSMatrix<BlockType,Alloc>& A,
+                                      const ParallelInfo& comm, const field_type w,
+                                      MILU_VARIANT milu, unsigned 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.
+      \copydoc Preconditioner::pre(X&,Y&)
+    */
+    virtual void pre (Domain& x, Range& b) override
+    {
+        DUNE_UNUSED_PARAMETER(x);
+        DUNE_UNUSED_PARAMETER(b);
+    }
+
+    /*!
+      \brief Apply the preconditoner.
+      \copydoc Preconditioner::apply(X&,const Y&)
+    */
+    virtual void apply (Domain& v, const Range& d) override
+    {
+        Range& md = const_cast<Range&>(d);;
+
+        // iterator types
+        typedef typename Range ::block_type  dblock;
+        typedef typename Domain::block_type  vblock;
+
+        const size_type iEnd = lower_.rows();
+        const size_type lastRow = iEnd - 1;
+        size_type interiorStart = 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 < lowerLoopEnd; ++ i )
+        {
+          dblock rhs( md[ i ] );
+          const size_type rowI     = lower_.rows_[ i ];
+          const size_type rowINext = lower_.rows_[ i+1 ];
+
+          for( size_type col = rowI; col < rowINext; ++ col )
+          {
+            lower_.values_[ col ].mmv( v[ lower_.cols_[ col ] ], rhs );
+          }
+
+          v[ i ] = rhs;  // Lii = I
+        }
+
+        for( size_type i = interiorStart; i < iEnd; ++ i )
+        {
+            vblock& vBlock = v[ lastRow - i ];
+            vblock rhs ( vBlock );
+            const size_type rowI     = upper_.rows_[ i ];
+            const size_type rowINext = upper_.rows_[ i+1 ];
+
+            for( size_type col = rowI; col < rowINext; ++ col )
+            {
+                upper_.values_[ col ].mmv( v[ upper_.cols_[ col ] ], rhs );
+            }
+
+            // apply inverse and store result
+            inv_[ i ].mv( rhs, vBlock);
+        }
+
+        copyOwnerToAll( v );
+
+        if( relaxation_ ) {
+            v *= w_;
+        }
+    }
+
+    template <class V>
+    void copyOwnerToAll( V& v ) const
+    {
+        if( comm_ ) {
+            comm_->copyOwnerToAll(v, v);
+        }
+    }
+
+    /*!
+      \brief Clean up.
+      \copydoc Preconditioner::post(X&)
+    */
+    virtual void post (Range& x) override
+    {
+        DUNE_UNUSED_PARAMETER(x);
+    }
+
+protected:
+    void init( const Matrix& A, const int iluIteration, MILU_VARIANT milu, bool redBlack, bool reorderSpheres )
+    {
+        // (For older DUNE versions the communicator might be
+        // invalid if redistribution in AMG happened on the coarset level.
+        // Therefore we check for nonzero size
+        if ( comm_ && comm_->communicator().size()<=0 )
+        {
+            if ( A.N() > 0 )
+            {
+                OPM_THROW(std::logic_error, "Expected a matrix with zero rows for an invalid communicator.");
+            }
+            else
+            {
+                // simply set the communicator to null
+                comm_ = nullptr;
+            }
+        }
+
+        int ilu_setup_successful = 1;
+        std::string message;
+        const int rank = ( comm_ ) ? comm_->communicator().rank() : 0;
+
+        std::unique_ptr< Matrix > ILU;
+
+        try
+        {
+            ILU.reset( new Matrix( A ) );
+            detail::ghost_last_bilu0_decomposition(*ILU, interiorSize_);
+        }
+        catch (const Dune::MatrixBlockError& error)
+        {
+            message = error.what();
+            std::cerr<<"Exception occured on process " << rank << " during " <<
+                "setup of ILU0 preconditioner with message: " <<
+                message<<std::endl;
+            ilu_setup_successful = 0;
+        }
+
+        // Check whether there was a problem on some process
+        const bool parallel_failure = comm_ && comm_->communicator().min(ilu_setup_successful) == 0;
+        const bool local_failure = ilu_setup_successful == 0;
+        if ( local_failure || parallel_failure )
+        {
+            throw Dune::MatrixBlockError();
+        }
+
+        // store ILU in simple CRS format
+        detail::convertToCRS( *ILU, lower_, upper_, inv_ );
+    }
+    
+protected:
+    //! \brief The ILU0 decomposition of the matrix.
+    CRS lower_;
+    CRS upper_;
+    std::vector< block_type > inv_;
+
+    const ParallelInfo* comm_;
+    //! \brief The relaxation factor to use.
+    const field_type w_;
+    const bool relaxation_;
+    size_type interiorSize_;
+};
+
+
 } // end namespace Opm
 #endif

opm/simulators/linalg/findOverlapRowsAndColumns.hpp

@@ -112,6 +112,28 @@ namespace detail
             }
         }
     }
+    template <class Grid>
+    size_t numInteriorCells(const Grid& grid, bool ownerFirst)
+    {
+        size_t numInterior = 0;
+        if (!ownerFirst)
+            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) {
+            const auto& elem = *elemIt;
+
+            // If cell has partition type not equal to interior save row
+            if (elem.partitionType() == Dune::InteriorEntity) {
+                numInterior++;
+            }
+        }
+
+        return numInterior;
+    }
 } // namespace detail
 } // namespace Opm