From 866a6612559606abf0bf766339f101099f62b1f2 Mon Sep 17 00:00:00 2001 From: andrthu Date: Sat, 7 Dec 2019 15:38:19 +0100 Subject: [PATCH] New owners-first based linear algebra operations (SoMV, SP and ILU.apply). --- ebos/eclthresholdpressure.hh | 2 +- opm/simulators/linalg/ISTLSolverEbos.hpp | 257 ++++++++++- .../linalg/ParallelOverlappingILU0.hpp | 421 +++++++++++++++++- .../linalg/findOverlapRowsAndColumns.hpp | 22 + 4 files changed, 682 insertions(+), 20 deletions(-) diff --git a/ebos/eclthresholdpressure.hh b/ebos/eclthresholdpressure.hh index 1df8ec5c7..b5026aa58 100644 --- a/ebos/eclthresholdpressure.hh +++ b/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(); diff --git a/opm/simulators/linalg/ISTLSolverEbos.hpp b/opm/simulators/linalg/ISTLSolverEbos.hpp index b3ce39320..cb5691adf 100644 --- a/opm/simulators/linalg/ISTLSolverEbos.hpp +++ b/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 WellModelGhostLastMatrixAdapter : public Dune::AssembledLinearOperator +{ + typedef Dune::AssembledLinearOperator 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 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( 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 GhostLastScalarProduct : public Dune::ScalarProduct +{ +public: + typedef X domain_type; + typedef typename X::field_type field_type; + typedef typename Dune::FieldTraits::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(std::sqrt(comm_->communicator().sum(result))); + } + +private: + std::unique_ptr< communication_type > comm_; + size_type interiorSize_; +}; + +template +struct GhostLastSPChooser +{ + typedef C communication_type; + + enum { solverCategory = c}; +}; + +template +struct GhostLastSPChooser +{ + typedef Dune::SeqScalarProduct ScalarProduct; + enum { category = Dune::SolverCategory::sequential}; + + static ScalarProduct* construct(const C& comm, size_t is) + { + return new ScalarProduct(); + } +}; + +template +struct GhostLastSPChooser +{ + typedef GhostLastScalarProduct 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(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 Comm; // 3x3 matrix block inversion was unstable from at least 2.3 until and // including 2.5.0 - typedef ParallelOverlappingILU0 ParPreconditioner; + + typedef GhostLastParallelOverlappingILU0 ParPreconditioner; + template std::unique_ptr 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_for_preconditioner_; + std::vector overlapRows_; std::vector interiorRows_; std::vector> wellConnectionsGraph_; + + bool ownersFirst_; + size_t interiorCellNum_; + FlowLinearSolverParameters parameters_; Vector weights_; bool scale_variables_; diff --git a/opm/simulators/linalg/ParallelOverlappingILU0.hpp b/opm/simulators/linalg/ParallelOverlappingILU0.hpp index b990322e6..4e25df28a 100644 --- a/opm/simulators/linalg/ParallelOverlappingILU0.hpp +++ b/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 + 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() 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 GhostLastParallelOverlappingILU0 + : public Dune::Preconditioner +{ + typedef ParallelInfoT ParallelInfo; + + +public: + //! \brief The matrix type the preconditioner is for. + typedef typename std::remove_const::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::value ? + Dune::SolverCategory::sequential : Dune::SolverCategory::overlapping; + } + +#else + // define the category + enum { + //! \brief The category the preconditioner is part of. + category = std::is_same::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 + GhostLastParallelOverlappingILU0 (const Dune::BCRSMatrix& 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(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 + GhostLastParallelOverlappingILU0 (const Dune::BCRSMatrix& 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(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 + GhostLastParallelOverlappingILU0 (const Dune::BCRSMatrix& 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 + GhostLastParallelOverlappingILU0 (const Dune::BCRSMatrix& 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(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(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 + 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<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 diff --git a/opm/simulators/linalg/findOverlapRowsAndColumns.hpp b/opm/simulators/linalg/findOverlapRowsAndColumns.hpp index 676e5278b..4665fca08 100644 --- a/opm/simulators/linalg/findOverlapRowsAndColumns.hpp +++ b/opm/simulators/linalg/findOverlapRowsAndColumns.hpp @@ -112,6 +112,28 @@ namespace detail } } } + template + 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