diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 46887cc9885..4b900534985 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -154,7 +154,6 @@ list (APPEND PUBLIC_HEADER_FILES opm/simulators/linalg/OwningBlockPreconditioner.hpp opm/simulators/linalg/OwningTwoLevelPreconditioner.hpp opm/simulators/linalg/ParallelOverlappingILU0.hpp - opm/simulators/linalg/ParallelRestrictedAdditiveSchwarz.hpp opm/simulators/linalg/ParallelIstlInformation.hpp opm/simulators/linalg/PressureSolverPolicy.hpp opm/simulators/linalg/PressureTransferPolicy.hpp diff --git a/ebos/eclbasevanguard.hh b/ebos/eclbasevanguard.hh index 9bd3fa268ad..5c452ed50b1 100644 --- a/ebos/eclbasevanguard.hh +++ b/ebos/eclbasevanguard.hh @@ -77,6 +77,7 @@ NEW_PROP_TAG(EclOutputInterval); NEW_PROP_TAG(IgnoreKeywords); NEW_PROP_TAG(EnableExperiments); NEW_PROP_TAG(EdgeWeightsMethod); +NEW_PROP_TAG(OverlapLayers); SET_STRING_PROP(EclBaseVanguard, IgnoreKeywords, ""); SET_STRING_PROP(EclBaseVanguard, EclDeckFileName, ""); @@ -84,6 +85,7 @@ SET_INT_PROP(EclBaseVanguard, EclOutputInterval, -1); // use the deck-provided v SET_BOOL_PROP(EclBaseVanguard, EnableOpmRstFile, false); SET_BOOL_PROP(EclBaseVanguard, EclStrictParsing, false); SET_INT_PROP(EclBaseVanguard, EdgeWeightsMethod, 1); +SET_INT_PROP(EclBaseVanguard, OverlapLayers, 1); END_PROPERTIES @@ -129,6 +131,8 @@ public: "Use strict mode for parsing - all errors are collected before the applicaton exists."); EWOMS_REGISTER_PARAM(TypeTag, int, EdgeWeightsMethod, "Choose edge-weighing strategy: 0=uniform, 1=trans, 2=log(trans)."); + EWOMS_REGISTER_PARAM(TypeTag, int, OverlapLayers, + "Number of overlap layers added between subdomains in the parallel grid. If set to larger then 1, Restricted Additive Schwarz (RAS) will be used as parallel preconditioner instead of the default Block-Jacobi method. Expect lower number of linear solver iterations at the cost of increased parallel overhead."); } /*! @@ -261,6 +265,7 @@ public: std::string fileName = EWOMS_GET_PARAM(TypeTag, std::string, EclDeckFileName); edgeWeightsMethod_ = Dune::EdgeWeightMethod(EWOMS_GET_PARAM(TypeTag, int, EdgeWeightsMethod)); + overlapLayers_ = EWOMS_GET_PARAM(TypeTag, int, OverlapLayers); // Skip processing of filename if external deck already exists. if (!externalDeck_) @@ -449,6 +454,13 @@ public: */ Dune::EdgeWeightMethod edgeWeightsMethod() const { return edgeWeightsMethod_; } + + /* + * \brief Parameter deciding the number of overlap layers between subdomains in parallel runs. + */ + int overlapLayers() const + {return overlapLayers_;} + /*! * \brief Returns the name of the case. * @@ -608,6 +620,7 @@ private: Opm::SummaryConfig* eclSummaryConfig_; Dune::EdgeWeightMethod edgeWeightsMethod_; + int overlapLayers_; }; template diff --git a/ebos/eclcpgridvanguard.hh b/ebos/eclcpgridvanguard.hh index 6e26ba43db4..6b012f38c80 100644 --- a/ebos/eclcpgridvanguard.hh +++ b/ebos/eclcpgridvanguard.hh @@ -155,6 +155,7 @@ public: } Dune::EdgeWeightMethod edgeWeightsMethod = this->edgeWeightsMethod(); + int overlapLayers = this->overlapLayers(); // convert to transmissibility for faces // TODO: grid_->numFaces() is not generic. use grid_->size(1) instead? (might @@ -187,7 +188,7 @@ public: //distribute the grid and switch to the distributed view. { const auto wells = this->schedule().getWellsatEnd(); - defunctWellNames_ = std::get<1>(grid_->loadBalance(edgeWeightsMethod, &wells, faceTrans.data())); + defunctWellNames_ = std::get<1>(grid_->loadBalance(edgeWeightsMethod, &wells, faceTrans.data(), overlapLayers)); } grid_->switchToDistributedView(); diff --git a/opm/simulators/linalg/CPRPreconditioner.hpp b/opm/simulators/linalg/CPRPreconditioner.hpp index df1432c013c..fa330ca86c4 100644 --- a/opm/simulators/linalg/CPRPreconditioner.hpp +++ b/opm/simulators/linalg/CPRPreconditioner.hpp @@ -49,7 +49,6 @@ #include #include -#include #include namespace Opm { diff --git a/opm/simulators/linalg/ISTLSolverEbos.hpp b/opm/simulators/linalg/ISTLSolverEbos.hpp index b3ce39320ad..c7e5ee2726d 100644 --- a/opm/simulators/linalg/ISTLSolverEbos.hpp +++ b/opm/simulators/linalg/ISTLSolverEbos.hpp @@ -24,7 +24,6 @@ #include #include #include -#include #include #include #include @@ -257,6 +256,8 @@ class WellModelMatrixAdapter : public Dune::AssembledLinearOperator const auto wellsForConn = simulator_.vanguard().schedule().getWellsatEnd(); const bool useWellConn = EWOMS_GET_PARAM(TypeTag, bool, MatrixAddWellContributions); + overlapLayers_ = EWOMS_GET_PARAM(TypeTag, int, OverlapLayers); + detail::setWellConnections(gridForConn, wellsForConn, useWellConn, wellConnectionsGraph_); detail::findOverlapAndInterior(gridForConn, overlapRows_, interiorRows_); if (gridForConn.comm().size() > 1) { @@ -391,6 +392,11 @@ class WellModelMatrixAdapter : public Dune::AssembledLinearOperator // Construct scalar product. auto sp = Dune::createScalarProduct(parallelInformation_arg, category); + // If overlapLayers_ > 1 we are using Restricted Additive Schwarz and therefore need + // the residual at overlap DoFs. + if (overlapLayers_ > 1) + parallelInformation_arg.copyOwnerToAll(istlb, istlb); + #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_) { @@ -688,7 +694,7 @@ class WellModelMatrixAdapter : public Dune::AssembledLinearOperator pattern[idx].insert(*wc); // Add just a single element to ghost rows - if (elem.partitionType() != Dune::InteriorEntity) + if (elem.partitionType() == Dune::GhostEntity) { noGhostMat_->setrowsize(idx, pattern[idx].size()); } @@ -949,6 +955,7 @@ class WellModelMatrixAdapter : public Dune::AssembledLinearOperator std::vector overlapRows_; std::vector interiorRows_; std::vector> wellConnectionsGraph_; + int overlapLayers_; FlowLinearSolverParameters parameters_; Vector weights_; bool scale_variables_; diff --git a/opm/simulators/linalg/ParallelOverlappingILU0.hpp b/opm/simulators/linalg/ParallelOverlappingILU0.hpp index b990322e6b7..fb65c7ad777 100644 --- a/opm/simulators/linalg/ParallelOverlappingILU0.hpp +++ b/opm/simulators/linalg/ParallelOverlappingILU0.hpp @@ -509,7 +509,7 @@ namespace Opm if( j.index() == iIndex ) { inv[ row ] = (*j); - break; + break; } else if ( j.index() >= i.index() ) { diff --git a/opm/simulators/linalg/ParallelRestrictedAdditiveSchwarz.hpp b/opm/simulators/linalg/ParallelRestrictedAdditiveSchwarz.hpp deleted file mode 100644 index 12b6f133a2f..00000000000 --- a/opm/simulators/linalg/ParallelRestrictedAdditiveSchwarz.hpp +++ /dev/null @@ -1,228 +0,0 @@ -/* - Copyright 2015 Dr. Blatt - HPC-Simulation-Software & Services - Copyright 2015 Statoil AS - - This file is part of the Open Porous Media project (OPM). - - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ -#ifndef OPM_PARALLELRESTRICTEDADDITIVESCHWARZ_HEADER_INCLUDED -#define OPM_PARALLELRESTRICTEDADDITIVESCHWARZ_HEADER_INCLUDED - -#include -#include -#include -#include - -namespace Opm -{ - -template -class ParallelRestrictedOverlappingSchwarz; - -} // end namespace Opm - -namespace Dune -{ - -namespace Amg -{ - -/// \brief Tells AMG how to construct the Opm::ParallelOverlappingILU0 smoother -/// \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 -/// \tparam SeqPreconditioner The underlying sequential preconditioner to use. -template -struct ConstructionTraits > -{ - typedef DefaultParallelConstructionArgs Arguments; - typedef ConstructionTraits SeqConstructionTraits; - - /// \brief Construct a parallel restricted overlapping schwarz preconditioner. -#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 7) - typedef std::shared_ptr< Opm::ParallelRestrictedOverlappingSchwarz > ParallelRestrictedOverlappingSchwarzPointer; -#else - typedef Opm::ParallelRestrictedOverlappingSchwarz* ParallelRestrictedOverlappingSchwarzPointer; -#endif - - static inline ParallelRestrictedOverlappingSchwarzPointer - construct(Arguments& args) - { - return -#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 7) - std::make_shared( -#endif - new Opm::ParallelRestrictedOverlappingSchwarz - (*SeqConstructionTraits ::construct(args), - args.getComm()) -#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 7) - ); -#else - ; -#endif - } - - /// \brief Deconstruct and free a parallel restricted overlapping schwarz preconditioner. - static inline void deconstruct(Opm::ParallelRestrictedOverlappingSchwarz - * bp) - { - SeqConstructionTraits - ::deconstruct(static_cast(&bp->preconditioner)); - delete bp; - } - -}; - -/// \brief Tells AMG how to use Opm::ParallelOverlappingILU0 smoother -/// \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 -/// \tparam SeqPreconditioner The underlying sequential preconditioner to use. -template -struct SmootherTraits > -{ - typedef DefaultSmootherArgs Arguments; - -}; - -} // end namespace Amg - -} // end namespace Dune - -namespace Opm{ - -/// \brief Block parallel preconditioner. -/// -/// This is essentially a wrapper that takes a sequential -/// preconditioner. In each step the sequential preconditioner -/// is applied to the whole subdomain and then all owner data -/// points are updated on all other processes from the processor -/// that knows the complete matrix row for this data point (in dune-istl -/// speak that is the one that owns the data). -/// -/// Note that this is different from the usual approach in dune-istl where -/// the application of the sequential preconditioner only takes place on -/// the (owner) partition of the process disregarding any overlap/ghost region. -/// -/// For more information see https://www.cs.colorado.edu/~cai/papers/rash.pdf -/// -/// \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 -/// \tparam SeqPreconditioner The underlying sequential preconditioner to use. -template > -class ParallelRestrictedOverlappingSchwarz - : public Dune::Preconditioner { - friend class Dune::Amg - ::ConstructionTraits >; -public: - //! \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; - //! \brief The type of the communication object. - typedef ParallelInfo communication_type; - - // define the category - enum { - //! \brief The category the precondtioner is part of. - category=Dune::SolverCategory::overlapping - }; - - /*! \brief Constructor. - - constructor gets all parameters to operate the prec. - \param p The sequential preconditioner. - \param c The communication object for syncing overlap and copy - data points. (E.~g. OwnerOverlapCommunication ) - */ - ParallelRestrictedOverlappingSchwarz (SeqPreconditioner& p, const communication_type& c) - : preconditioner_(p), communication_(c) - { } - - /*! - \brief Prepare the preconditioner. - - \copydoc Preconditioner::pre(X&,Y&) - */ - virtual void pre (Domain& x, Range& b) - { - communication_.copyOwnerToAll(x,x); // make dirichlet values consistent - preconditioner_.pre(x,b); - } - - /*! - \brief Apply the preconditioner - - \copydoc Preconditioner::apply(X&,const Y&) - */ - virtual void apply (Domain& v, const Range& d) - { - apply(v, d); - } - - template - void apply (Domain& v, const Range& d) - { - // hack us a mutable d to prevent copying. - Range& md = const_cast(d); - communication_.copyOwnerToAll(md,md); - preconditioner_.template apply(v,d); - communication_.copyOwnerToAll(v,v); - // Make sure that d is the same as at the beginning of apply. - communication_.project(md); - } - - /*! - \brief Clean up. - - \copydoc Preconditioner::post(X&) - */ - virtual void post (Range& x) - { - preconditioner_.post(x); - } - -private: - //! \brief a sequential preconditioner - SeqPreconditioner& preconditioner_; - - //! \brief the communication object - const communication_type& communication_; -}; - - -} // end namespace OPM -#endif diff --git a/opm/simulators/linalg/findOverlapRowsAndColumns.hpp b/opm/simulators/linalg/findOverlapRowsAndColumns.hpp index 676e5278b4e..f7e0d6221a2 100644 --- a/opm/simulators/linalg/findOverlapRowsAndColumns.hpp +++ b/opm/simulators/linalg/findOverlapRowsAndColumns.hpp @@ -41,7 +41,7 @@ namespace detail template void setWellConnections(const Grid& grid, const W& wells, bool useWellConn, std::vector>& wellGraph) { - if ( grid.comm().size() > 1) + if ( grid.comm().size() > 1 ) { Dune::CartesianIndexMapper< Grid > cartMapper( grid ); const int numCells = cartMapper.compressedSize(); // grid.numCells() @@ -102,7 +102,7 @@ namespace detail const auto& elem = *elemIt; int lcell = lid.id(elem); - if (elem.partitionType() != Dune::InteriorEntity) + if (elem.partitionType() == Dune::GhostEntity) { //add row to list overlapRows.push_back(lcell);