diff --git a/ebos/eclbasevanguard.hh b/ebos/eclbasevanguard.hh index fef9eae80..1bba603e3 100644 --- a/ebos/eclbasevanguard.hh +++ b/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. diff --git a/ebos/eclcpgridvanguard.hh b/ebos/eclcpgridvanguard.hh index 80abf8176..627fcda80 100644 --- a/ebos/eclcpgridvanguard.hh +++ b/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 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) { diff --git a/ebos/eclthresholdpressure.hh b/ebos/eclthresholdpressure.hh index 1df8ec5c7..2048921da 100644 --- a/ebos/eclthresholdpressure.hh +++ b/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(); diff --git a/opm/simulators/linalg/ISTLSolverEbos.hpp b/opm/simulators/linalg/ISTLSolverEbos.hpp index b3ce39320..b05fe6f90 100644 --- a/opm/simulators/linalg/ISTLSolverEbos.hpp +++ b/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 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 + + + 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( 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_); - detail::findOverlapAndInterior(gridForConn, overlapRows_, interiorRows_); - if (gridForConn.comm().size() > 1) { - noGhostAdjacency(); - setGhostsInNoGhost(*noGhostMat_); + ownersFirst_ = EWOMS_GET_PARAM(TypeTag, bool, OwnerCellsFirst); + interiorCellNum_ = detail::numMatrixRowsToUseInSolver(simulator_.vanguard().grid(), ownersFirst_); + + if (!ownersFirst_ || parameters_.linear_solver_use_amg_ || parameters_.use_cpr_ ) { + 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_ ); + + assert( opA.comm() ); + solve( opA, x, *rhs_, *(opA.comm()) ); + } + else { - copyJacToNoGhost(*matrix_, *noGhostMat_); - Operator opA(*noGhostMat_, *noGhostMat_, wellModel, - parallelInformation_ ); - assert( opA.comm() ); - 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 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..5739fdef0 100644 --- a/opm/simulators/linalg/ParallelOverlappingILU0.hpp +++ b/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 + 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() 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(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(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(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 + ParallelOverlappingILU0 (const Dune::BCRSMatrix& 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(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 + 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