diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 100b4fc21..8a0ac9aed 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -170,10 +170,8 @@ list (APPEND PUBLIC_HEADER_FILES opm/simulators/linalg/bda/openclSolverBackend.hpp opm/simulators/linalg/bda/MultisegmentWellContribution.hpp opm/simulators/linalg/bda/WellContributions.hpp - opm/simulators/linalg/BlackoilAmg.hpp opm/simulators/linalg/amgcpr.hh opm/simulators/linalg/twolevelmethodcpr.hh - opm/simulators/linalg/CPRPreconditioner.hpp opm/simulators/linalg/ExtractParallelGridInformationToISTL.hpp opm/simulators/linalg/FlexibleSolver.hpp opm/simulators/linalg/FlexibleSolver_impl.hpp diff --git a/flow/flow_blackoil_dunecpr.cpp b/flow/flow_blackoil_dunecpr.cpp index d17219ff2..3a4d0c3df 100644 --- a/flow/flow_blackoil_dunecpr.cpp +++ b/flow/flow_blackoil_dunecpr.cpp @@ -50,14 +50,6 @@ namespace Opm { static constexpr int value = 100; }; template - struct UseAmg { // probably not used - static constexpr bool value = true; - }; - template - struct UseCpr { - static constexpr bool value = true; - }; - template struct CprMaxEllIter { static constexpr int value = 1; }; @@ -70,17 +62,9 @@ namespace Opm { static constexpr int value = 3; }; template - struct CprSolverVerbose { - static constexpr int value = 0; - }; - template struct LinearSolverConfiguration { static constexpr auto value = "ilu0"; }; - template - struct SystemStrategy { - static constexpr auto value = "quasiimpes"; - }; template struct FluidSystem diff --git a/flow/flow_ebos_blackoil.cpp b/flow/flow_ebos_blackoil.cpp index b59e42556..ff6df4dd1 100644 --- a/flow/flow_ebos_blackoil.cpp +++ b/flow/flow_ebos_blackoil.cpp @@ -16,9 +16,6 @@ */ #include "config.h" -// Define making clear that the simulator supports AMG -#define FLOW_SUPPORT_AMG 1 - #include #include diff --git a/flow/flow_ebos_gasoil.cpp b/flow/flow_ebos_gasoil.cpp index 26547bc63..33dc8ae38 100644 --- a/flow/flow_ebos_gasoil.cpp +++ b/flow/flow_ebos_gasoil.cpp @@ -16,9 +16,6 @@ */ #include "config.h" -// Define making clear that the simulator supports AMG -#define FLOW_SUPPORT_AMG 1 - #include #include diff --git a/flow/flow_ebos_oilwater.cpp b/flow/flow_ebos_oilwater.cpp index a75c79a44..78a1f6b6d 100644 --- a/flow/flow_ebos_oilwater.cpp +++ b/flow/flow_ebos_oilwater.cpp @@ -16,9 +16,6 @@ */ #include "config.h" -// Define making clear that the simulator supports AMG -#define FLOW_SUPPORT_AMG 1 - #include #include diff --git a/flow/flow_ebos_oilwater_brine.cpp b/flow/flow_ebos_oilwater_brine.cpp index 6872b3f4f..5acdc9b28 100644 --- a/flow/flow_ebos_oilwater_brine.cpp +++ b/flow/flow_ebos_oilwater_brine.cpp @@ -16,9 +16,6 @@ */ #include "config.h" -// Define making clear that the simulator supports AMG -#define FLOW_SUPPORT_AMG 1 - #include #include diff --git a/flow/flow_ebos_oilwater_polymer.cpp b/flow/flow_ebos_oilwater_polymer.cpp index e4ed3e684..4c7a6bdf9 100644 --- a/flow/flow_ebos_oilwater_polymer.cpp +++ b/flow/flow_ebos_oilwater_polymer.cpp @@ -16,9 +16,6 @@ */ #include "config.h" -// Define making clear that the simulator supports AMG -#define FLOW_SUPPORT_AMG 1 - #include #include diff --git a/flow/flow_ebos_oilwater_polymer_injectivity.cpp b/flow/flow_ebos_oilwater_polymer_injectivity.cpp index fae09d421..fb6ef02fd 100644 --- a/flow/flow_ebos_oilwater_polymer_injectivity.cpp +++ b/flow/flow_ebos_oilwater_polymer_injectivity.cpp @@ -16,9 +16,6 @@ */ #include "config.h" -// Define making clear that the simulator supports AMG -#define FLOW_SUPPORT_AMG 1 - #include #include diff --git a/opm/simulators/linalg/BlackoilAmg.hpp b/opm/simulators/linalg/BlackoilAmg.hpp deleted file mode 100644 index 86235ccd7..000000000 --- a/opm/simulators/linalg/BlackoilAmg.hpp +++ /dev/null @@ -1,1043 +0,0 @@ -/* - Copyright 2017 Dr. Blatt - HPC-Simulation-Software & Services - - 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_AMG_HEADER_INCLUDED -#define OPM_AMG_HEADER_INCLUDED - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -namespace Dune -{ -namespace Amg -{ -template -class UnSymmetricCriterion; -} -} - -namespace Opm -{ - namespace Amg - { - template - class Element - { - public: - enum { /* @brief We preserve the sign.*/ - is_sign_preserving = true - }; - - template - typename M::field_type operator()(const M& m) const - { - return m[Row][Column]; - } - }; - } // namespace Amg -} // namespace Opm - -namespace Dune -{ - -template -class MatrixBlock; - -} - -namespace Opm -{ - -namespace Detail -{ - -/** - * \brief Creates a MatrixAdapter as an operator, storing it in a unique_ptr. - * - * The first argument is used to specify the return type using function overloading. - * \param matrix The matrix to wrap. - */ -template -std::unique_ptr< Dune::MatrixAdapter > -createOperator(const Dune::MatrixAdapter&, const M& matrix, const T&) -{ - return std::make_unique< Dune::MatrixAdapter >(matrix); -} - -/** - * \brief Creates an OverlappingSchwarzOperator as an operator, storing it in a unique_ptr. - * - * The first argument is used to specify the return type using function overloading. - * \param matrix The matrix to wrap. - * \param comm The object encapsulating the parallelization information. - */ -template -std::unique_ptr< Dune::OverlappingSchwarzOperator > -createOperator(const Dune::OverlappingSchwarzOperator&, const M& matrix, const T& comm) -{ - return std::make_unique< Dune::OverlappingSchwarzOperator >(matrix, comm); -} - -//! \brief Applies diagonal scaling to the discretization Matrix (Scheichl, 2003) -//! -//! See section 3.2.3 of Scheichl, Masson: Decoupling and Block Preconditioning for -//! Sedimentary Basin Simulations, 2003. -//! \param op The operator that stems from the discretization. -//! \param comm The communication objecte describing the data distribution. -//! \param pressureEqnIndex The index of the pressure in the matrix block -//! \retun A pair of the scaled matrix and the associated operator- -template -std::unique_ptr -scaleMatrixDRS(const Operator& op, std::size_t pressureEqnIndex, const Vector& weights, const Opm::CPRParameter& param) -{ - using Matrix = typename Operator::matrix_type; - using Block = typename Matrix::block_type; - using BlockVector = typename Vector::block_type; - auto matrix = std::make_unique(op.getmat()); - if (param.cpr_use_drs_) { - const auto endi = matrix->end(); - for (auto i = matrix->begin(); i != endi; ++i) { - const BlockVector& bw = weights[i.index()]; - const auto endj = (*i).end(); - for (auto j = (*i).begin(); j != endj; ++j) { - Block& block = *j; - BlockVector& bvec = block[pressureEqnIndex]; - // should introduce limits which also change the weights - block.mtv(bw, bvec); - } - } - } - return matrix; -} - -//! \brief Applies diagonal scaling to the discretization Matrix (Scheichl, 2003) -//! -//! See section 3.2.3 of Scheichl, Masson: Decoupling and Block Preconditioning for -//! Sedimentary Basin Simulations, 2003. -//! \param vector The vector to scale -//! \param pressureEqnIndex The index of the pressure in the matrix block -template -void scaleVectorDRS(Vector& vector, std::size_t pressureEqnIndex, const Opm::CPRParameter& param, const Vector& weights) -{ - using Block = typename Vector::block_type; - if (param.cpr_use_drs_) { - for (std::size_t j = 0; j < vector.size(); ++j) { - Block& block = vector[j]; - const Block& bw = weights[j]; - block[pressureEqnIndex] = bw.dot(block); - } - } -} - -//! \brief TMP to create the scalar pendant to a real block matrix, vector, smoother, etc. -//! -//! \code -//! using Scalar = ScalarType::value; -//! \endcode -//! will get the corresponding scalar type. -template -struct ScalarType -{ -}; - -template -struct ScalarType > -{ - typedef Dune::FieldVector value; -}; - -template -struct ScalarType > -{ - typedef Dune::FieldMatrix value; -}; - -template -struct ScalarType > -{ - typedef Dune::MatrixBlock value; -}; - -template -struct ScalarType > -{ - typedef Opm::MatrixBlock value; -}; - -template -struct ScalarType > -{ - using ScalarBlock = typename ScalarType::value; - using ScalarAllocator = typename Allocator::template rebind::other; - typedef Dune::BCRSMatrix value; -}; - -template -struct ScalarType > -{ - using ScalarBlock = typename ScalarType::value; - using ScalarAllocator = typename Allocator::template rebind::other; - typedef Dune::BlockVector value; -}; - -template -struct ScalarType > -{ - typedef Dune::SeqScalarProduct::value> value; -}; - -#define ComposeScalarTypeForSeqPrecond(PREC) \ - template \ - struct ScalarType > \ - { \ - typedef PREC::value, \ - typename ScalarType::value, \ - typename ScalarType::value, \ - l> value; \ - } - -ComposeScalarTypeForSeqPrecond(Dune::SeqJac); -ComposeScalarTypeForSeqPrecond(Dune::SeqSOR); -ComposeScalarTypeForSeqPrecond(Dune::SeqSSOR); -#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 7) -ComposeScalarTypeForSeqPrecond(Dune::SeqILU); -#else -// Dune::SeqGS and Dune::SeqSOR are the same in DUNE 2.7 -ComposeScalarTypeForSeqPrecond(Dune::SeqGS); -ComposeScalarTypeForSeqPrecond(Dune::SeqILU0); -ComposeScalarTypeForSeqPrecond(Dune::SeqILUn); -#endif - -#undef ComposeScalarTypeForSeqPrecond - -template -struct ScalarType > -{ - typedef Dune::Richardson::value, - typename ScalarType::value> value; -}; - -template -struct ScalarType > -{ - typedef Dune::OverlappingSchwarzOperator::value, - typename ScalarType::value, - typename ScalarType::value, - C> value; -}; - -template -struct ScalarType > -{ - typedef Dune::MatrixAdapter::value, - typename ScalarType::value, - typename ScalarType::value> value; -}; - -template -struct ScalarType > -{ - typedef Dune::OverlappingSchwarzScalarProduct::value, - C> value; -}; - -template -struct ScalarType > -{ - typedef Dune::NonoverlappingSchwarzScalarProduct::value, - C> value; -}; - -template -struct ScalarType > -{ - typedef Dune::BlockPreconditioner::value, - typename ScalarType::value, - C, - typename ScalarType::value> value; -}; - -template -struct ScalarType > -{ - typedef ParallelOverlappingILU0::value, - typename ScalarType::value, - typename ScalarType::value, - C> value; -}; - -template -struct ScalarType,N> > > -{ - using value = Dune::Amg::CoarsenCriterion::value>, Dune::Amg::FirstDiagonal> >; -}; - -template -struct ScalarType,N> > > -{ - using value = Dune::Amg::CoarsenCriterion::value>, Dune::Amg::FirstDiagonal> >; -}; - -template -struct OneComponentCriterionType -{}; - -template -struct OneComponentCriterionType,N> >, COMPONENT_INDEX, VARIABLE_INDEX> -{ - using value = Dune::Amg::CoarsenCriterion, Opm::Amg::Element > >; -}; - -template -struct OneComponentCriterionType,N> >, COMPONENT_INDEX, VARIABLE_INDEX> -{ - using value = Dune::Amg::CoarsenCriterion, Opm::Amg::Element > >; -}; - -template -class OneComponentAggregationLevelTransferPolicy; - - -/** - * @brief A policy class for solving the coarse level system using one step of AMG. - * @tparam O The type of the linear operator used. - * @tparam S The type of the smoother used in AMG. - * @tparam C The type of the crition used for the aggregation within AMG. - * @tparam C1 The type of the information about the communication. Either - * Dune::OwnerOverlapCopyCommunication or Dune::SequentialInformation. - */ -template -class OneStepAMGCoarseSolverPolicy -{ -public: - typedef P LevelTransferPolicy; - /** @brief The type of the linear operator used. */ - typedef O Operator; - /** @brief The type of the range and domain of the operator. */ - typedef typename O::range_type X; - /** @brief The type of the crition used for the aggregation within AMG.*/ - typedef C Criterion; - /** @brief The type of the communication used for AMG.*/ - typedef typename P::ParallelInformation Communication; - /** @brief The type of the smoother used in AMG. */ - typedef S Smoother; - /** @brief The type of the arguments used for constructing the smoother. */ - typedef typename Dune::Amg::SmootherTraits::Arguments SmootherArgs; - /** @brief The type of the AMG construct on the coarse level.*/ - typedef Dune::Amg::AMGCPR AMGType; - /** - * @brief Constructs the coarse solver policy. - * @param args The arguments used for constructing the smoother. - * @param c The crition used for the aggregation within AMG. - */ - OneStepAMGCoarseSolverPolicy(const CPRParameter* param, const SmootherArgs& args, const Criterion& c) - : param_(param), smootherArgs_(args), criterion_(c) - {} - /** @brief Copy constructor. */ - OneStepAMGCoarseSolverPolicy(const OneStepAMGCoarseSolverPolicy& other) - : param_(other.param_), coarseOperator_(other.coarseOperator_), smootherArgs_(other.smootherArgs_), - criterion_(other.criterion_) - {} -private: - /** - * @brief A wrapper that makes an inverse operator out of AMG. - * - * The operator will use one step of AMG to approximately solve - * the coarse level system. - */ - struct AMGInverseOperator : public Dune::InverseOperator - { - AMGInverseOperator(const CPRParameter* param, - const typename AMGType::Operator& op, - const Criterion& crit, - const typename AMGType::SmootherArgs& args, - const Communication& comm) - : param_(param), amg_(), smoother_(), op_(op), crit_(crit), comm_(comm) - { - if ( param_->cpr_use_amg_ ) - { - amg_.reset(new AMGType(op, crit,args, comm)); - } - else - { - typename Dune::Amg::ConstructionTraits::Arguments cargs; - cargs.setMatrix(op.getmat()); - cargs.setComm(comm); - cargs.setArgs(args); -#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 7) - smoother_ = Dune::Amg::ConstructionTraits::construct(cargs); -#else - smoother_.reset(Dune::Amg::ConstructionTraits::construct(cargs)); -#endif - } - } - - void updatePreconditioner() - { - amg_->updateSolver(crit_, op_, comm_); - } - - Dune::SolverCategory::Category category() const override - { - return std::is_same::value ? - Dune::SolverCategory::sequential : Dune::SolverCategory::overlapping; - } - - void apply(X& x, X& b, double reduction, Dune::InverseOperatorResult& res) override - { - DUNE_UNUSED_PARAMETER(reduction); - DUNE_UNUSED_PARAMETER(res); - - auto sp = Dune::createScalarProduct(comm_, op_.category()); - Dune::Preconditioner* prec = amg_.get(); - if ( ! amg_ ) - { - prec = smoother_.get(); - } - // Linear solver parameters - const double tolerance = param_->cpr_solver_tol_; - const int maxit = param_->cpr_max_ell_iter_; - int verbosity = 0; - if (comm_.communicator().rank() == 0) { - verbosity = param_->cpr_solver_verbose_; - } - - if ( param_->cpr_ell_solvetype_ == 0) - { - Dune::BiCGSTABSolver solver(const_cast(op_), *sp, *prec, - tolerance, maxit, verbosity); - solver.apply(x,b,res); - } - else if (param_->cpr_ell_solvetype_ == 1) - { - Dune::CGSolver solver(const_cast(op_), *sp, *prec, - tolerance, maxit, verbosity); - solver.apply(x,b,res); - } - else - { - Dune::LoopSolver solver(const_cast(op_), *sp, *prec, - tolerance, maxit, verbosity); - solver.apply(x,b,res); - } - - // Warn if unknown options. - if (param_->cpr_ell_solvetype_ > 2 && comm_.communicator().rank() == 0) { - OpmLog::warning("cpr_ell_solver_type_unknown", "Unknown CPR elliptic solver type specification, using LoopSolver."); - } - } - - void apply(X& x, X& b, Dune::InverseOperatorResult& res) override - { - return apply(x,b,1e-8,res); - } - - ~AMGInverseOperator() - {} - AMGInverseOperator(const AMGInverseOperator& other) - : x_(other.x_), amg_(other.amg_) - { - } - private: - const CPRParameter* param_; - X x_; - std::unique_ptr amg_; - std::shared_ptr smoother_; - const typename AMGType::Operator& op_; - Criterion crit_; - const Communication& comm_; - }; - -public: - /** @brief The type of solver constructed for the coarse level. */ - typedef AMGInverseOperator CoarseLevelSolver; - - /** - * @brief Constructs a coarse level solver. - * - * @param transferPolicy The policy describing the transfer between levels. - * @return A pointer to the constructed coarse level solver. - */ - template - CoarseLevelSolver* createCoarseLevelSolver(LTP& transferPolicy) - { - coarseOperator_=transferPolicy.getCoarseLevelOperator(); - const LevelTransferPolicy& transfer = - reinterpret_cast(transferPolicy); - AMGInverseOperator* inv = new AMGInverseOperator(param_, - *coarseOperator_, - criterion_, - smootherArgs_, - transfer.getCoarseLevelCommunication()); - - return inv; - } - - template - void setCoarseOperator(LTP& transferPolicy) - { - coarseOperator_= transferPolicy.getCoarseLevelOperator(); - } - -private: - /** @brief The coarse level operator. */ - std::shared_ptr coarseOperator_; - /** @brief The parameters for the CPR preconditioner. */ - const CPRParameter* param_; - /** @brief The arguments used to construct the smoother. */ - SmootherArgs smootherArgs_; - /** @brief The coarsening criterion. */ - Criterion criterion_; -}; - -template -std::shared_ptr< Smoother > -constructSmoother(const Operator& op, - const typename Dune::Amg::SmootherTraits::Arguments& smargs, - const Communication& comm) -{ - typename Dune::Amg::ConstructionTraits::Arguments args; - args.setMatrix(op.getmat()); - args.setComm(comm); - args.setArgs(smargs); - // for DUNE < 2.7 ConstructionTraits::construct returns a raw - // pointer, therefore std::make_shared cannot be used here - return std::shared_ptr< Smoother > (Dune::Amg::ConstructionTraits::construct(args)); -} - -template -const Dune::Amg::OverlapVertex* -buildOverlapVertices(const G& graph, const C& pinfo, - Dune::Amg::AggregatesMap& aggregates, - const S& overlap, - std::size_t& overlapCount) -{ - // count the overlap vertices. - overlapCount = 0; - - const auto& lookup=pinfo.globalLookup(); - - for ( const auto& vertex: graph ) { - const auto* pair = lookup.pair(vertex); - - if(pair!=0 && overlap.contains(pair->local().attribute())) - ++overlapCount; - } - // Allocate space - using Vertex = typename G::VertexDescriptor; - using OverlapVertex = Dune::Amg::OverlapVertex; - - auto* overlapVertices = new OverlapVertex[overlapCount==0 ? 1 : overlapCount]; - if(overlapCount==0) - return overlapVertices; - - // Initialize them - overlapCount=0; - for ( const auto& vertex: graph ) { - const auto* pair = lookup.pair(vertex); - - if(pair!=0 && overlap.contains(pair->local().attribute())) { - overlapVertices[overlapCount].aggregate = &aggregates[pair->local()]; - overlapVertices[overlapCount].vertex = pair->local(); - ++overlapCount; - } - } - std::sort(overlapVertices, overlapVertices+overlapCount, - [](const OverlapVertex& v1, const OverlapVertex& v2) - { - return *v1.aggregate < *v2.aggregate; - }); - // due to the sorting the isolated aggregates (to be skipped) are at the end. - - return overlapVertices; -} - -template -void buildCoarseSparseMatrix(M& coarseMatrix, G& fineGraph, - const V& visitedMap, - const C& pinfo, - Dune::Amg::AggregatesMap& aggregates, - const S& overlap) -{ - using OverlapVertex = Dune::Amg ::OverlapVertex; - std::size_t count; - - const OverlapVertex* overlapVertices = buildOverlapVertices(fineGraph, - pinfo, - aggregates, - overlap, - count); - - // Reset the visited flags of all vertices. - // As the isolated nodes will be skipped we simply mark them as visited -#ifndef NDEBUG - const auto UNAGGREGATED = Dune::Amg::AggregatesMap::UNAGGREGATED; -#endif - const auto ISOLATED = Dune::Amg::AggregatesMap::ISOLATED; - - for ( const auto& vertex: fineGraph ) { - assert(aggregates[vertex] != UNAGGREGATED); - put(visitedMap, vertex, aggregates[vertex]==ISOLATED); - } - - Dune::Amg::SparsityBuilder sparsityBuilder(coarseMatrix); - - Dune::Amg::ConnectivityConstructor::examine(fineGraph, visitedMap, pinfo, - aggregates, overlap, - overlapVertices, - overlapVertices+count, - sparsityBuilder); - delete[] overlapVertices; -} - -template -void buildCoarseSparseMatrix(M& coarseMatrix, G& fineGraph, const V& visitedMap, - const Dune::Amg::SequentialInformation& pinfo, - Dune::Amg::AggregatesMap& aggregates, - const S&) -{ - // Reset the visited flags of all vertices. - // As the isolated nodes will be skipped we simply mark them as visited -#ifndef NDEBUG - const auto UNAGGREGATED = Dune::Amg::AggregatesMap::UNAGGREGATED; -#endif - const auto ISOLATED = Dune::Amg::AggregatesMap::ISOLATED; - - for(const auto& vertex: fineGraph ) { - assert(aggregates[vertex] != UNAGGREGATED); - put(visitedMap, vertex, aggregates[vertex]==ISOLATED); - } - - Dune::Amg::SparsityBuilder sparsityBuilder(coarseMatrix); - - Dune::Amg::ConnectivityConstructor - ::examine(fineGraph, visitedMap, pinfo, aggregates, sparsityBuilder); -} - -} // end namespace Detail - -/** - * @brief A LevelTransferPolicy that uses aggregation to construct the coarse level system. - * @tparam Operator The type of the fine level operator. - * @tparam Criterion The criterion that describes the aggregation procedure. - * @tparam Communication The class that describes the communication pattern. - */ -template -class OneComponentAggregationLevelTransferPolicy - : public Dune::Amg::LevelTransferPolicyCpr::value> -{ - typedef Dune::Amg::AggregatesMap AggregatesMap; -public: - using CoarseOperator = typename Detail::ScalarType::value; - typedef Dune::Amg::LevelTransferPolicy FatherType; - typedef Communication ParallelInformation; - -public: - OneComponentAggregationLevelTransferPolicy(const Criterion& crit, const Communication& comm, - bool cpr_pressure_aggregation) - : criterion_(crit), communication_(&const_cast(comm)), - cpr_pressure_aggregation_(cpr_pressure_aggregation) - {} - - void createCoarseLevelSystem(const Operator& fineOperator) - { - prolongDamp_ = 1; - - if ( cpr_pressure_aggregation_ ) - { - typedef Dune::Amg::PropertiesGraphCreator GraphCreator; - typedef typename GraphCreator::PropertiesGraph PropertiesGraph; - typedef typename GraphCreator::GraphTuple GraphTuple; - - typedef typename PropertiesGraph::VertexDescriptor Vertex; - - std::vector excluded(fineOperator.getmat().N(), false); - - using OverlapFlags = Dune::NegateSet; - GraphTuple graphs = GraphCreator::create(fineOperator, excluded, - *communication_, OverlapFlags()); - - aggregatesMap_.reset(new AggregatesMap(std::get<1>(graphs)->maxVertex()+1)); - - int noAggregates, isoAggregates, oneAggregates, skippedAggregates; - using std::get; - std::tie(noAggregates, isoAggregates, oneAggregates, skippedAggregates) = - aggregatesMap_->buildAggregates(fineOperator.getmat(), *(get<1>(graphs)), - criterion_, true); - - using CommunicationArgs = typename Dune::Amg::ConstructionTraits::Arguments; -#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 7) - CommunicationArgs commArgs(communication_->communicator(), communication_->category()); -#else - CommunicationArgs commArgs(communication_->communicator(), communication_->getSolverCategory()); -#endif -#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 7) - coarseLevelCommunication_ = Dune::Amg::ConstructionTraits::construct(commArgs); -#else - coarseLevelCommunication_.reset(Dune::Amg::ConstructionTraits::construct(commArgs)); -#endif - using Iterator = typename std::vector::iterator; - using std::get; - auto visitedMap = get(Dune::Amg::VertexVisitedTag(), *(get<1>(graphs))); - communication_->buildGlobalLookup(fineOperator.getmat().N()); - std::size_t aggregates = - Dune::Amg::IndicesCoarsener - ::coarsen(*communication_, *get<1>(graphs), visitedMap, - *aggregatesMap_, *coarseLevelCommunication_, - noAggregates); - GraphCreator::free(graphs); - coarseLevelCommunication_->buildGlobalLookup(aggregates); - Dune::Amg::AggregatesPublisher - ::publish(*aggregatesMap_, - *communication_, - coarseLevelCommunication_->globalLookup()); - std::vector& visited=excluded; - - std::fill(visited.begin(), visited.end(), false); - - Dune::IteratorPropertyMap - visitedMap2(visited.begin(), Dune::IdentityMap()); - using CoarseMatrix = typename CoarseOperator::matrix_type; - coarseLevelMatrix_.reset(new CoarseMatrix(aggregates, aggregates, - CoarseMatrix::row_wise)); - Detail::buildCoarseSparseMatrix(*coarseLevelMatrix_, *get<0>(graphs), visitedMap2, - *communication_, - *aggregatesMap_, - OverlapFlags()); - delete get<0>(graphs); - communication_->freeGlobalLookup(); - if( static_cast(this->coarseLevelMatrix_->N()) - < criterion_.coarsenTarget()) - { - coarseLevelCommunication_->freeGlobalLookup(); - } - calculateCoarseEntriesWithAggregatesMap(fineOperator); - } - else - { - using CoarseMatrix = typename CoarseOperator::matrix_type; - const auto& fineLevelMatrix = fineOperator.getmat(); - coarseLevelMatrix_.reset(new CoarseMatrix(fineLevelMatrix.N(), fineLevelMatrix.M(), CoarseMatrix::row_wise)); - auto createIter = coarseLevelMatrix_->createbegin(); - - for ( const auto& row: fineLevelMatrix ) - { - for ( auto col = row.begin(), cend = row.end(); col != cend; ++col) - { - createIter.insert(col.index()); - } - ++createIter; - } - - auto coarseRow = coarseLevelMatrix_->begin(); - for ( const auto& row: fineLevelMatrix ) - { - auto coarseCol = coarseRow->begin(); - - for ( auto col = row.begin(), cend = row.end(); col != cend; ++col, ++coarseCol ) - { - assert( col.index() == coarseCol.index() ); - *coarseCol = (*col)[COMPONENT_INDEX][VARIABLE_INDEX]; - } - ++coarseRow; - } - coarseLevelCommunication_.reset(communication_, [](Communication*){}); - } - - this->lhs_.resize(this->coarseLevelMatrix_->M()); - this->rhs_.resize(this->coarseLevelMatrix_->N()); - using OperatorArgs = typename Dune::Amg::ConstructionTraits::Arguments; -#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 7) - OperatorArgs oargs(coarseLevelMatrix_, *coarseLevelCommunication_); - this->operator_ = Dune::Amg::ConstructionTraits::construct(oargs); -#else - OperatorArgs oargs(*coarseLevelMatrix_, *coarseLevelCommunication_); - this->operator_.reset(Dune::Amg::ConstructionTraits::construct(oargs)); -#endif - } - - void calculateCoarseEntriesWithAggregatesMap(const Operator& fineOperator) - { - const auto& fineMatrix = fineOperator.getmat(); - *coarseLevelMatrix_ = 0; - for(auto row = fineMatrix.begin(), rowEnd = fineMatrix.end(); - row != rowEnd; ++row) - { - const auto& i = (*aggregatesMap_)[row.index()]; - if(i != AggregatesMap::ISOLATED) - { - for(auto entry = row->begin(), entryEnd = row->end(); - entry != entryEnd; ++entry) - { - const auto& j = (*aggregatesMap_)[entry.index()]; - if ( j != AggregatesMap::ISOLATED ) - { - (*coarseLevelMatrix_)[i][j] += (*entry)[COMPONENT_INDEX][VARIABLE_INDEX]; - } - } - } - } - } - - virtual void calculateCoarseEntries(const Operator& fineOperator) - { - const auto& fineMatrix = fineOperator.getmat(); - *coarseLevelMatrix_ = 0; - for(auto row = fineMatrix.begin(), rowEnd = fineMatrix.end(); - row != rowEnd; ++row) - { - const auto& i = row.index(); - for(auto entry = row->begin(), entryEnd = row->end(); - entry != entryEnd; ++entry) - { - const auto& j = entry.index(); - (*coarseLevelMatrix_)[i][j] += (*entry)[COMPONENT_INDEX][VARIABLE_INDEX]; - } - } - } - - void moveToCoarseLevel(const typename FatherType::FineRangeType& fine) - { - // Set coarse vector to zero - this->rhs_=0; - - if ( cpr_pressure_aggregation_ ) - { - auto end = fine.end(), begin=fine.begin(); - - for(auto block=begin; block != end; ++block) - { - const auto& vertex = (*aggregatesMap_)[block-begin]; - if(vertex != AggregatesMap::ISOLATED) - { - this->rhs_[vertex] += (*block)[COMPONENT_INDEX]; - } - } - } - else - { - auto end = fine.end(), begin=fine.begin(); - - for(auto block=begin; block != end; ++block) - { - this->rhs_[block-begin] = (*block)[COMPONENT_INDEX]; - } - } - - this->lhs_=0; - } - - void moveToFineLevel(typename FatherType::FineDomainType& fine) - { - if( cpr_pressure_aggregation_ ) - { - this->lhs_ *= prolongDamp_; - auto end=fine.end(), begin=fine.begin(); - - for(auto block=begin; block != end; ++block) - { - const auto& vertex = (*aggregatesMap_)[block-begin]; - if(vertex != AggregatesMap::ISOLATED) - (*block)[COMPONENT_INDEX] += this->lhs_[vertex]; - } - communication_->copyOwnerToAll(fine,fine); - } - else - { - auto end=fine.end(), begin=fine.begin(); - - for(auto block=begin; block != end; ++block) - { - (*block)[COMPONENT_INDEX] = this->lhs_[block-begin]; - } - } - } - - OneComponentAggregationLevelTransferPolicy* clone() const - { - return new OneComponentAggregationLevelTransferPolicy(*this); - } - - const Communication& getCoarseLevelCommunication() const - { - return *coarseLevelCommunication_; - } -private: - typename Operator::matrix_type::field_type prolongDamp_; - std::shared_ptr aggregatesMap_; - Criterion criterion_; - Communication* communication_; - std::shared_ptr coarseLevelCommunication_; - std::shared_ptr coarseLevelMatrix_; - bool cpr_pressure_aggregation_; -}; - -/** - * \brief An algebraic twolevel or multigrid approach for solving blackoil (supports CPR with and without AMG) - * - * This preconditioner first decouples the component used for coarsening using a simple scaling - * approach (e.g. Scheichl, Masson 2013,\see scaleMatrixDRS). Then it constructs the - * coarse level system. The coupling is defined by the weights corresponding to the element located at - * (COMPONENT_INDEX, VARIABLE_INDEX) in the block matrix. Then the coarse level system is constructed - * either by extracting these elements, or by applying aggregation to them directly. This coarse level - * can be solved either by AMG or by ILU. The preconditioner is configured using CPRParameter. - * \tparam O The type of the operator (encapsulating a BCRSMatrix). - * \tparam S The type of the smoother. - * \tparam C The type of coarsening criterion to use. - * \tparam P The type of the class describing the parallelization. - * \tparam COMPONENT_INDEX The index of the component to use for coarsening (usually water). - * \tparam VARIABLE_INDEX The index of the variable to use for coarsening (usually pressure). - */ -template -class BlackoilAmg - : public Dune::Preconditioner -{ -public: - /** \brief The type of the operator (encapsulating a BCRSMatrix). */ - using Operator = O; - /** \brief The type of coarsening criterion to use. */ - using Criterion = C; - /** \brief The type of the class describing the parallelization. */ - using Communication = P; - /** \brief The type of the smoother. */ - using Smoother = S; - /** \brief The type of the smoother arguments for construction. */ - using SmootherArgs = typename Dune::Amg::SmootherTraits::Arguments; - -protected: - using Matrix = typename Operator::matrix_type; - using CoarseOperator = typename Detail::ScalarType::value; - using CoarseSmoother = typename Detail::ScalarType::value; - using FineCriterion = - typename Detail::OneComponentCriterionType::value; - using CoarseCriterion = typename Detail::ScalarType::value; - using LevelTransferPolicy = - OneComponentAggregationLevelTransferPolicy; - using CoarseSolverPolicy = - Detail::OneStepAMGCoarseSolverPolicy; - using TwoLevelMethod = - Dune::Amg::TwoLevelMethodCpr; -public: - Dune::SolverCategory::Category category() const override - { - return std::is_same::value ? - Dune::SolverCategory::sequential : Dune::SolverCategory::overlapping; - } - - /** - * \brief Constructor. - * \param param The parameters used for configuring the solver. - * \param fineOperator The operator of the fine level. - * \param criterion The criterion describing the coarsening approach. - * \param smargs The arguments for constructing the smoother. - * \param comm The information about the parallelization. - */ - BlackoilAmg(const CPRParameter& param, - const typename TwoLevelMethod::FineDomainType& weights, - const Operator& fineOperator, const Criterion& criterion, - const SmootherArgs& smargs, const Communication& comm) - : param_(param), - weights_(weights), - scaledMatrix_(Detail::scaleMatrixDRS(fineOperator, COMPONENT_INDEX, weights, param)), - scaledMatrixOperator_(Detail::createOperator(fineOperator, *scaledMatrix_, comm)), - smoother_( Detail::constructSmoother(*scaledMatrixOperator_, smargs, comm)), - levelTransferPolicy_(criterion, comm, param.cpr_pressure_aggregation_), - coarseSolverPolicy_(¶m, smargs, criterion), - twoLevelMethod_(*scaledMatrixOperator_, smoother_, - levelTransferPolicy_, coarseSolverPolicy_, 0, 1) - { - } - - void pre(typename TwoLevelMethod::FineDomainType& x, - typename TwoLevelMethod::FineRangeType& b) override - { - twoLevelMethod_.pre(x,b); - } - - void post(typename TwoLevelMethod::FineDomainType& x) override - { - twoLevelMethod_.post(x); - } - - void apply(typename TwoLevelMethod::FineDomainType& v, - const typename TwoLevelMethod::FineRangeType& d) override - { - auto scaledD = d; - Detail::scaleVectorDRS(scaledD, COMPONENT_INDEX, param_, weights_); - twoLevelMethod_.apply(v, scaledD); - } -private: - const CPRParameter& param_; - const typename TwoLevelMethod::FineDomainType& weights_; - std::unique_ptr scaledMatrix_; - std::unique_ptr scaledMatrixOperator_; - std::shared_ptr smoother_; - LevelTransferPolicy levelTransferPolicy_; - CoarseSolverPolicy coarseSolverPolicy_; - TwoLevelMethod twoLevelMethod_; -}; - -namespace ISTLUtility -{ -/// -/// \brief A traits class for selecting the types of the preconditioner. -/// -/// \tparam M The type of the matrix. -/// \tparam X The type of the domain of the linear problem. -/// \tparam Y The type of the range of the linear problem. -/// \tparam P The type of the parallel information. -/// \tparam C The type of the coarsening criterion to use. -/// \tparam index The pressure index. -//// -template -struct BlackoilAmgSelector -{ - using Criterion = C; - using Selector = CPRSelector; - using ParallelInformation = typename Selector::ParallelInformation; - using Operator = typename Selector::Operator; - using Smoother = typename Selector::EllipticPreconditioner; - using AMG = BlackoilAmg; -}; -} // end namespace ISTLUtility -} // end namespace Opm -#endif diff --git a/opm/simulators/linalg/CPRPreconditioner.hpp b/opm/simulators/linalg/CPRPreconditioner.hpp deleted file mode 100644 index df1432c01..000000000 --- a/opm/simulators/linalg/CPRPreconditioner.hpp +++ /dev/null @@ -1,528 +0,0 @@ -/* - Copyright 2014 SINTEF ICT, Applied Mathematics. - Copyright 2014 IRIS AS. - Copyright 2015 Dr. Blatt - HPC-Simulation-Software & Services - Copyright 2015 NTNU - 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_CPRPRECONDITIONER_HEADER_INCLUDED -#define OPM_CPRPRECONDITIONER_HEADER_INCLUDED - -#include -#include - -#include -#include -#include -#include - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include - -#include -#include - -#include -#include -namespace Opm -{ - -template -class BlackoilAmg; - -namespace Amg -{ - template - class Element; -} - -namespace ISTLUtility -{ - -template -void setILUParameters(Opm::ParallelOverlappingILU0Args& args, - const CPRParameter& params) -{ - args.setN(params.cpr_ilu_n_); - args.setMilu(params.cpr_ilu_milu_); -} - -template -void setILUParameters(Opm::ParallelOverlappingILU0Args& args, - MILU_VARIANT milu, int n=0) -{ - args.setN(n); - args.setMilu(milu); -} - -template -void setILUParameters(S&, const P&) -{} - -template -void setILUParameters(S&, bool, int) -{} - -/// -/// \brief A traits class for selecting the types of the preconditioner. -/// -/// \tparam M The type of the matrix. -/// \tparam X The type of the domain of the linear problem. -/// \tparam Y The type of the range of the linear problem. -/// \tparam P The type of the parallel information. -//// -template -struct CPRSelector -{ - /// \brief The information about the parallelization and communication - using ParallelInformation = P; - /// \brief The operator type; - typedef Dune::OverlappingSchwarzOperator Operator; - typedef ParallelOverlappingILU0 - EllipticPreconditioner; - /// \brief The type of the unique pointer to the preconditioner of the elliptic part. - typedef std::unique_ptr - EllipticPreconditionerPointer; - - typedef EllipticPreconditioner Smoother; - typedef Dune::Amg::AMG AMG; - - /// \brief creates an Operator from the matrix - /// \param M The matrix to use. - /// \param p The parallel information to use. - static Operator* makeOperator(const M& m, const ParallelInformation& p) - { - return new Operator(m, p); - } -}; - -template -struct CPRSelector -{ - /// \brief The information about the parallelization and communication - typedef Dune::Amg::SequentialInformation ParallelInformation; - /// \brief The operator type; - typedef Dune::MatrixAdapter Operator; - /// \brief The type of the preconditioner used for the elliptic part. - typedef ParallelOverlappingILU0 - EllipticPreconditioner; - /// \brief The type of the unique pointer to the preconditioner of the elliptic part. - typedef std::unique_ptr EllipticPreconditionerPointer; - - /// \brief type of AMG used to precondition the elliptic system. - typedef EllipticPreconditioner Smoother; - typedef Dune::Amg::AMG AMG; - - /// \brief creates an Operator from the matrix - /// \param M The matrix to use. - /// \param p The parallel information to use. - static Operator* makeOperator(const M& m, const ParallelInformation&) - { - return new Operator(m); - } - -}; -//! \brief Creates and initializes a unique pointer to an sequential ILU0 preconditioner. -//! \param A The matrix of the linear system to solve. -//! \param relax The relaxation factor to use. -template -std::shared_ptr > -createILU0Ptr(const M& A, const C& comm, double relax, MILU_VARIANT milu) -{ - return std::make_shared >(A, comm, relax, milu); -} -//! \brief Creates and initializes a shared pointer to an ILUn preconditioner. -//! \param A The matrix of the linear system to solve. -//! \param ilu_n The n parameter for the extension of the nonzero pattern. -//! \param relax The relaxation factor to use. -template -std::shared_ptr > -createILUnPtr(const M& A, const C& comm, int ilu_n, double relax, MILU_VARIANT milu) -{ - return std::make_shared >( A, comm, ilu_n, relax, milu ); -} -/// \brief Creates the elliptic preconditioner (ILU0) -/// \param Ae The matrix of the elliptic system. -/// \param relax The relaxation parameter for ILU0. -/// \param milu If true, the modified ilu approach is used. Dropped elements -/// will get added to the diagonal of U to preserve the row sum -/// for constant vectors (Ae = LUe). -/// \param comm The object describing the parallelization information and communication. -template -typename CPRSelector::EllipticPreconditionerPointer -createEllipticPreconditionerPointer(const M& Ae, double relax, - MILU_VARIANT milu, const P& comm) -{ - typedef typename CPRSelector - ::EllipticPreconditioner ParallelPreconditioner; - - typedef typename CPRSelector - ::EllipticPreconditionerPointer EllipticPreconditionerPointer; - return EllipticPreconditionerPointer(new ParallelPreconditioner(Ae, comm, relax, milu)); -} - -template < class C, class Op, class P, class S, std::size_t PressureEqnIndex, std::size_t PressureVarIndex, class Vector> -inline void -createAMGPreconditionerPointer(Op& opA, const double relax, const P& comm, - std::unique_ptr< BlackoilAmg >& amgPtr, - const CPRParameter& params, - const Vector& weights) -{ - using AMG = BlackoilAmg; - int verbosity = 0; - if (comm.communicator().rank() == 0) { - verbosity = params.cpr_solver_verbose_; - } - // TODO: revise choice of parameters - int coarsenTarget=1200; - using Criterion = C; - Criterion criterion(15, coarsenTarget); - criterion.setDebugLevel( verbosity ); // no debug information, 1 for printing hierarchy information - criterion.setDefaultValuesIsotropic(2); - criterion.setNoPostSmoothSteps( 1 ); - criterion.setNoPreSmoothSteps( 1 ); - - // Since DUNE 2.2 we also need to pass the smoother args instead of steps directly - typedef typename AMG::Smoother Smoother; - typedef typename Dune::Amg::SmootherTraits::Arguments SmootherArgs; - SmootherArgs smootherArgs; - smootherArgs.iterations = 1; - smootherArgs.relaxationFactor = relax; - setILUParameters(smootherArgs, params); - - amgPtr.reset( new AMG( params, weights, opA, criterion, smootherArgs, comm ) ); -} - -template < class C, class Op, class P, class AMG > -inline void -createAMGPreconditionerPointer(Op& opA, const double relax, const MILU_VARIANT milu, const P& comm, std::unique_ptr< AMG >& amgPtr) -{ - // TODO: revise choice of parameters - int coarsenTarget=1200; - using Criterion = C; - Criterion criterion(15, coarsenTarget); - criterion.setDebugLevel( 0 ); // no debug information, 1 for printing hierarchy information - criterion.setDefaultValuesIsotropic(2); - criterion.setNoPostSmoothSteps( 1 ); - criterion.setNoPreSmoothSteps( 1 ); - - // for DUNE 2.2 we also need to pass the smoother args - typedef typename AMG::Smoother Smoother; - typedef typename Dune::Amg::SmootherTraits::Arguments SmootherArgs; - SmootherArgs smootherArgs; - smootherArgs.iterations = 1; - smootherArgs.relaxationFactor = relax; - setILUParameters(smootherArgs, milu); - - amgPtr.reset( new AMG(opA, criterion, smootherArgs, comm ) ); -} - -/// \brief Creates the elliptic preconditioner (ILU0) -/// \param opA The operator representing the matrix of the system. -/// \param relax The relaxation parameter for ILU0. -/// \param comm The object describing the parallelization information and communication. -// \param amgPtr The unique_ptr to be filled (return) -template < int PressureEqnIndex, int PressureVarIndex, class Op, class P, class AMG > -inline void -createAMGPreconditionerPointer( Op& opA, const double relax, const MILU_VARIANT milu, const P& comm, std::unique_ptr< AMG >& amgPtr ) -{ - // type of matrix - typedef typename Op::matrix_type M; - - // The coupling metric used in the AMG - typedef Opm::Amg::Element CouplingMetric; - - // The coupling criterion used in the AMG - typedef Dune::Amg::SymmetricCriterion CritBase; - - // The coarsening criterion used in the AMG - typedef Dune::Amg::CoarsenCriterion Criterion; - - createAMGPreconditionerPointer(opA, relax, milu, comm, amgPtr); -} - -} // end namespace ISTLUtility - - /*! - \brief CPR preconditioner. - - This is a two-stage preconditioner, combining an elliptic-type - partial solution with ILU0 for the whole system. - - \tparam M The matrix type to operate on - \tparam X Type of the update - \tparam Y Type of the defect - \tparam P Type of the parallel information. If not provided - this will be Dune::Amg::SequentialInformation. - The preconditioner is parallel if this is - Dune::OwnerOverlapCopyCommunication - */ - template - class CPRPreconditioner : public Dune::Preconditioner - { - // prohibit copying for now - CPRPreconditioner( const CPRPreconditioner& ); - - public: - //! \brief The type describing the parallel information - typedef P ParallelInformation; - //! \brief The matrix type the preconditioner is for. - typedef typename std::remove_const::type matrix_type; - //! \brief The domain type of the preconditioner. - typedef X domain_type; - //! \brief The range type of the preconditioner. - typedef Y range_type; - //! \brief The field type of the preconditioner. - typedef typename X::field_type field_type; - - // define the category - Dune::SolverCategory::Category category() const override - { - return std::is_same::value ? - Dune::SolverCategory::sequential : Dune::SolverCategory::overlapping; - } - - typedef ISTLUtility::CPRSelector CPRSelectorType ; - - //! \brief Elliptic Operator - typedef typename CPRSelectorType::Operator Operator; - - //! \brief preconditioner for the whole system (here either ILU(0) or ILU(n) - typedef Dune::Preconditioner WholeSystemPreconditioner; - - //! \brief type of the unique pointer to the ilu-0 preconditioner - //! used the for the elliptic system - typedef typename CPRSelectorType::EllipticPreconditionerPointer - EllipticPreconditionerPointer; - - //! \brief amg preconditioner for the elliptic system - typedef typename CPRSelectorType::AMG AMG; - - /*! \brief Constructor. - - Constructor gets all parameters to operate the prec. - \param A The matrix to operate on. - \param Ae The top-left elliptic part of A. - \param relax The ILU0 relaxation factor. - \param useAMG if true, AMG is used as a preconditioner for the elliptic sub-system, otherwise ilu-0 (default) - \param useBiCG if true, BiCG solver is used (default), otherwise CG solver - \param paralleInformation The information about the parallelization, if this is a - parallel run - */ - CPRPreconditioner (const CPRParameter& param, const M& A, const M& Ae, - const ParallelInformation& comm=ParallelInformation(), - const ParallelInformation& commAe=ParallelInformation()) - : param_( param ), - A_(A), - Ae_(Ae), - de_( Ae_.N() ), - ve_( Ae_.M() ), - dmodified_( A_.N() ), - opAe_(CPRSelectorType::makeOperator(Ae_, commAe)), - precond_(), // ilu0 preconditioner for elliptic system - amg_(), // amg preconditioner for elliptic system - pre_(), // copy A will be made be the preconditioner - vilu_( A_.N() ), - comm_(comm), - commAe_(commAe) - { - // create appropriate preconditioner for elliptic system - createEllipticPreconditioner( param_.cpr_use_amg_, commAe_ ); - - // create the preconditioner for the whole system. - if( param_.cpr_ilu_n_ == 0 ) { - pre_ = ISTLUtility::createILU0Ptr( A_, comm, param_.cpr_relax_, param_.cpr_ilu_milu_ ); - } - else { - pre_ = ISTLUtility::createILUnPtr( A_, comm, param_.cpr_ilu_n_, param_.cpr_relax_, param_.cpr_ilu_milu_); - } - } - - /*! - \brief Prepare the preconditioner. - - \copydoc Preconditioner::pre(X&,Y&) - */ - virtual void pre (X& /*x*/, Y& /*b*/) - { - } - - /*! - \brief Apply the preconditoner. - - \copydoc Preconditioner::apply(X&,const Y&) - */ - virtual void apply (X& v, const Y& d) - { - // Extract part of d corresponding to elliptic part. - // Note: Assumes that the elliptic part comes first. - std::copy_n(d.begin(), de_.size(), de_.begin()); - - // Solve elliptic part, extend solution to full. - // reset result - ve_ = 0; - solveElliptic( ve_, de_ ); - - //reset return value - v = 0.0; - // Again assuming that the elliptic part comes first. - std::copy(ve_.begin(), ve_.end(), v.begin()); - - // Subtract elliptic residual from initial residual. - // dmodified = d - A * vfull - dmodified_ = d; - A_.mmv(v, dmodified_); - // A is not parallel, do communication manually. - comm_.copyOwnerToAll(dmodified_, dmodified_); - - // Apply Preconditioner for whole system (relax will be applied already) - pre_->apply( vilu_, dmodified_); - - // don't apply relaxation if relax_ == 1 - if( std::abs( param_.cpr_relax_ - 1.0 ) < 1e-12 ) { - v += vilu_; - } - else { - v *= param_.cpr_relax_; - v += vilu_; - } - } - - /*! - \brief Clean up. - - \copydoc Preconditioner::post(X&) - */ - virtual void post (X& /*x*/) - { - } - - protected: - void solveElliptic(Y& x, Y& de) - { - // Linear solver parameters - const double tolerance = param_.cpr_solver_tol_; - const int maxit = param_.cpr_max_ell_iter_; - const int verbosity = ( param_.cpr_solver_verbose_ && - comm_.communicator().rank()==0 ) ? 1 : 0; - - // operator result containing iterations etc. - Dune::InverseOperatorResult result; - - // the scalar product chooser - auto sp = Dune::createScalarProduct(commAe_, category()); - - if( amg_ ) - { - // Solve system with AMG - if( param_.cpr_use_bicgstab_ ) { - Dune::BiCGSTABSolver linsolve(*opAe_, *sp, (*amg_), tolerance, maxit, verbosity); - linsolve.apply(x, de, result); - } - else { - Dune::CGSolver linsolve(*opAe_, *sp, (*amg_), tolerance, maxit, verbosity); - linsolve.apply(x, de, result); - } - } - else - { - assert( precond_ ); - // Solve system with ILU-0 - if( param_.cpr_use_bicgstab_ ) { - Dune::BiCGSTABSolver linsolve(*opAe_, *sp, (*precond_), tolerance, maxit, verbosity); - linsolve.apply(x, de, result); - } - else { - Dune::CGSolver linsolve(*opAe_, *sp, (*precond_), tolerance, maxit, verbosity); - linsolve.apply(x, de, result); - } - - } - - if (!result.converged) { - OPM_THROW(LinearSolverProblem, "CPRPreconditioner failed to solve elliptic subsystem."); - } - } - - //! \brief Parameter collection for CPR - const CPRParameter& param_; - - //! \brief The matrix for the full linear problem. - const matrix_type& A_; - //! \brief The elliptic part of the matrix. - const matrix_type& Ae_; - - //! \brief temporary variables for elliptic solve - Y de_, ve_, dmodified_; - - //! \brief elliptic operator - std::unique_ptr opAe_; - - //! \brief ILU0 preconditioner for the elliptic system - EllipticPreconditionerPointer precond_; - //! \brief AMG preconditioner with ILU0 smoother - std::unique_ptr< AMG > amg_; - - //! \brief The preconditioner for the whole system - //! - //! We have to use a shared_ptr instead of a unique_ptr - //! as we need to use a custom allocator based on dynamic - //! information. But for unique_ptr the type of this deleter - //! has to be available at coompile time. - std::shared_ptr< WholeSystemPreconditioner > pre_; - - //! \brief temporary variables for ILU solve - Y vilu_; - - //! \brief The information about the parallelization of the whole system. - const P& comm_; - //! \brief The information about the parallelization of the elliptic part - //! of the system - const P& commAe_; - protected: - void createEllipticPreconditioner( const bool amg, const P& comm ) - { - if( amg ) - { - ISTLUtility::createAMGPreconditionerPointer( *opAe_ , param_.cpr_relax_, param_.cpr_ilu_milu_, comm, amg_ ); - } - else - { - precond_ = ISTLUtility::createEllipticPreconditionerPointer( Ae_, param_.cpr_relax_, param_.cpr_ilu_milu_, comm); - } - } - }; - - -} // namespace Opm - -#endif // OPM_CPRPRECONDITIONER_HEADER_INCLUDED diff --git a/opm/simulators/linalg/FlowLinearSolverParameters.hpp b/opm/simulators/linalg/FlowLinearSolverParameters.hpp index 21f73102b..80b5f1d61 100644 --- a/opm/simulators/linalg/FlowLinearSolverParameters.hpp +++ b/opm/simulators/linalg/FlowLinearSolverParameters.hpp @@ -94,34 +94,14 @@ struct LinearSolverIgnoreConvergenceFailure{ using type = UndefinedProperty; }; template -struct UseAmg { - using type = UndefinedProperty; -}; -template -struct UseCpr { - using type = UndefinedProperty; -}; -template struct PreconditionerAddWellContributions { using type = UndefinedProperty; }; template -struct SystemStrategy { - using type = UndefinedProperty; -}; -template struct ScaleLinearSystem { using type = UndefinedProperty; }; template -struct CprSolverVerbose { - using type = UndefinedProperty; -}; -template -struct CprUseDrs { - using type = UndefinedProperty; -}; -template struct CprMaxEllIter { using type = UndefinedProperty; }; @@ -205,14 +185,6 @@ struct LinearSolverIgnoreConvergenceFailure static constexpr bool value = false; }; template -struct UseAmg { - static constexpr bool value = false; -}; -template -struct UseCpr { - static constexpr bool value = false; -}; -template struct LinearSolverBackend { using type = Opm::ISTLSolverEbos; }; @@ -221,22 +193,10 @@ struct PreconditionerAddWellContributions { static constexpr bool value = false; }; template -struct SystemStrategy { - static constexpr auto value = "none"; -}; -template struct ScaleLinearSystem { static constexpr bool value = false; }; template -struct CprSolverVerbose { - static constexpr int value = 0; -}; -template -struct CprUseDrs { - static constexpr bool value = false; -}; -template struct CprMaxEllIter { static constexpr int value = 20; }; @@ -277,48 +237,48 @@ namespace Opm /** * \brief Parameters used to configure the CPRPreconditioner. */ - struct CPRParameter - { - double cpr_relax_; - double cpr_solver_tol_; - int cpr_ilu_n_; - MILU_VARIANT cpr_ilu_milu_; - bool cpr_ilu_redblack_; - bool cpr_ilu_reorder_sphere_; - bool cpr_use_drs_; - int cpr_max_ell_iter_; - int cpr_ell_solvetype_; - bool cpr_use_amg_; - bool cpr_use_bicgstab_; - int cpr_solver_verbose_; - bool cpr_pressure_aggregation_; - int cpr_reuse_setup_; - CPRParameter() { reset(); } + // struct CPRParameter + // { + // double cpr_relax_; + // double cpr_solver_tol_; + // int cpr_ilu_n_; + // MILU_VARIANT cpr_ilu_milu_; + // bool cpr_ilu_redblack_; + // bool cpr_ilu_reorder_sphere_; + // bool cpr_use_drs_; + // int cpr_max_ell_iter_; + // int cpr_ell_solvetype_; + // bool cpr_use_amg_; + // bool cpr_use_bicgstab_; + // int cpr_solver_verbose_; + // bool cpr_pressure_aggregation_; + // int cpr_reuse_setup_; + // CPRParameter() { reset(); } - void reset() - { - cpr_solver_tol_ = 1e-2; - cpr_ilu_n_ = 0; - cpr_ilu_milu_ = MILU_VARIANT::ILU; - cpr_ilu_redblack_ = false; - cpr_ilu_reorder_sphere_ = true; - cpr_max_ell_iter_ = 25; - cpr_ell_solvetype_ = 0; - cpr_use_drs_ = false; - cpr_use_amg_ = true; - cpr_use_bicgstab_ = true; - cpr_solver_verbose_ = 0; - cpr_pressure_aggregation_ = false; - cpr_reuse_setup_ = 0; - } - }; + // void reset() + // { + // cpr_solver_tol_ = 1e-2; + // cpr_ilu_n_ = 0; + // cpr_ilu_milu_ = MILU_VARIANT::ILU; + // cpr_ilu_redblack_ = false; + // cpr_ilu_reorder_sphere_ = true; + // cpr_max_ell_iter_ = 25; + // cpr_ell_solvetype_ = 0; + // cpr_use_drs_ = false; + // cpr_use_amg_ = true; + // cpr_use_bicgstab_ = true; + // cpr_solver_verbose_ = 0; + // cpr_pressure_aggregation_ = false; + // cpr_reuse_setup_ = 0; + // } + // }; /// This class carries all parameters for the NewtonIterationBlackoilInterleaved class. struct FlowLinearSolverParameters - : public CPRParameter + //: public CPRParameter { double linear_solver_reduction_; double ilu_relaxation_; @@ -332,15 +292,18 @@ namespace Opm bool newton_use_gmres_; bool require_full_sparsity_pattern_; bool ignoreConvergenceFailure_; - bool linear_solver_use_amg_; - bool use_cpr_; - std::string system_strategy_; + // bool linear_solver_use_amg_; + // bool use_cpr_; + // std::string system_strategy_; bool scale_linear_system_; std::string linear_solver_configuration_; std::string linear_solver_configuration_json_file_; std::string gpu_mode_; int bda_device_id_; int opencl_platform_id_; + int cpr_max_ell_iter_ = 20; + int cpr_reuse_setup_ = 0; + bool use_gpu_; template void init() @@ -358,14 +321,14 @@ namespace Opm newton_use_gmres_ = EWOMS_GET_PARAM(TypeTag, bool, UseGmres); require_full_sparsity_pattern_ = EWOMS_GET_PARAM(TypeTag, bool, LinearSolverRequireFullSparsityPattern); ignoreConvergenceFailure_ = EWOMS_GET_PARAM(TypeTag, bool, LinearSolverIgnoreConvergenceFailure); - linear_solver_use_amg_ = EWOMS_GET_PARAM(TypeTag, bool, UseAmg); - use_cpr_ = EWOMS_GET_PARAM(TypeTag, bool, UseCpr); - system_strategy_ = EWOMS_GET_PARAM(TypeTag, std::string, SystemStrategy); + // linear_solver_use_amg_ = EWOMS_GET_PARAM(TypeTag, bool, UseAmg); + // use_cpr_ = EWOMS_GET_PARAM(TypeTag, bool, UseCpr); + // system_strategy_ = EWOMS_GET_PARAM(TypeTag, std::string, SystemStrategy); scale_linear_system_ = EWOMS_GET_PARAM(TypeTag, bool, ScaleLinearSystem); - cpr_solver_verbose_ = EWOMS_GET_PARAM(TypeTag, int, CprSolverVerbose); - cpr_use_drs_ = EWOMS_GET_PARAM(TypeTag, bool, CprUseDrs); + // cpr_solver_verbose_ = EWOMS_GET_PARAM(TypeTag, int, CprSolverVerbose); + // cpr_use_drs_ = EWOMS_GET_PARAM(TypeTag, bool, CprUseDrs); cpr_max_ell_iter_ = EWOMS_GET_PARAM(TypeTag, int, CprMaxEllIter); - cpr_ell_solvetype_ = EWOMS_GET_PARAM(TypeTag, int, CprEllSolvetype); + // cpr_ell_solvetype_ = EWOMS_GET_PARAM(TypeTag, int, CprEllSolvetype); cpr_reuse_setup_ = EWOMS_GET_PARAM(TypeTag, int, CprReuseSetup); linear_solver_configuration_ = EWOMS_GET_PARAM(TypeTag, std::string, LinearSolverConfiguration); linear_solver_configuration_json_file_ = EWOMS_GET_PARAM(TypeTag, std::string, LinearSolverConfigurationJsonFile); @@ -389,14 +352,14 @@ namespace Opm EWOMS_REGISTER_PARAM(TypeTag, bool, UseGmres, "Use GMRES as the linear solver"); EWOMS_REGISTER_PARAM(TypeTag, bool, LinearSolverRequireFullSparsityPattern, "Produce the full sparsity pattern for the linear solver"); EWOMS_REGISTER_PARAM(TypeTag, bool, LinearSolverIgnoreConvergenceFailure, "Continue with the simulation like nothing happened after the linear solver did not converge"); - EWOMS_REGISTER_PARAM(TypeTag, bool, UseAmg, "Use AMG as the linear solver's preconditioner"); - EWOMS_REGISTER_PARAM(TypeTag, bool, UseCpr, "Use CPR as the linear solver's preconditioner"); - EWOMS_REGISTER_PARAM(TypeTag, std::string, SystemStrategy, "Strategy for reformulating and scaling linear system (none: no scaling -- should not be used with CPR, original: use weights that are equivalent to no scaling -- should not be used with CPR, simple: form pressure equation as simple sum of conservation equations, quasiimpes: form pressure equation based on diagonal block, trueimpes: form pressure equation based on linearization of accumulation term)"); + // EWOMS_REGISTER_PARAM(TypeTag, bool, UseAmg, "Use AMG as the linear solver's preconditioner"); + // EWOMS_REGISTER_PARAM(TypeTag, bool, UseCpr, "Use CPR as the linear solver's preconditioner"); + // EWOMS_REGISTER_PARAM(TypeTag, std::string, SystemStrategy, "Strategy for reformulating and scaling linear system (none: no scaling -- should not be used with CPR, original: use weights that are equivalent to no scaling -- should not be used with CPR, simple: form pressure equation as simple sum of conservation equations, quasiimpes: form pressure equation based on diagonal block, trueimpes: form pressure equation based on linearization of accumulation term)"); EWOMS_REGISTER_PARAM(TypeTag, bool, ScaleLinearSystem, "Scale linear system according to equation scale and primary variable types"); - EWOMS_REGISTER_PARAM(TypeTag, int, CprSolverVerbose, "Verbosity of cpr solver (0: silent, 1: print summary of inner linear solver, 2: print extensive information about inner linear solve, including setup information)"); - EWOMS_REGISTER_PARAM(TypeTag, bool, CprUseDrs, "Use dynamic row sum using weights"); + // EWOMS_REGISTER_PARAM(TypeTag, int, CprSolverVerbose, "Verbosity of cpr solver (0: silent, 1: print summary of inner linear solver, 2: print extensive information about inner linear solve, including setup information)"); + // EWOMS_REGISTER_PARAM(TypeTag, bool, CprUseDrs, "Use dynamic row sum using weights"); EWOMS_REGISTER_PARAM(TypeTag, int, CprMaxEllIter, "MaxIterations of the elliptic pressure part of the cpr solver"); - EWOMS_REGISTER_PARAM(TypeTag, int, CprEllSolvetype, "Solver type of elliptic pressure solve (0: bicgstab, 1: cg, 2: only amg preconditioner)"); + // EWOMS_REGISTER_PARAM(TypeTag, int, CprEllSolvetype, "Solver type of elliptic pressure solve (0: bicgstab, 1: cg, 2: only amg preconditioner)"); EWOMS_REGISTER_PARAM(TypeTag, int, CprReuseSetup, "Reuse Amg Setup"); EWOMS_REGISTER_PARAM(TypeTag, std::string, LinearSolverConfiguration, "Configuration of solver valid is: ilu0 (default), cpr_quasiimpes, cpr_trueimpes or file (specified in LinearSolverConfigurationJsonFile) "); EWOMS_REGISTER_PARAM(TypeTag, std::string, LinearSolverConfigurationJsonFile, "Filename of JSON configuration for flexible linear solver system."); @@ -410,7 +373,7 @@ namespace Opm // set default values void reset() { - use_cpr_ = false; + // use_cpr_ = false; newton_use_gmres_ = false; linear_solver_reduction_ = 1e-2; linear_solver_maxiter_ = 150; @@ -418,7 +381,7 @@ namespace Opm linear_solver_verbosity_ = 0; require_full_sparsity_pattern_ = false; ignoreConvergenceFailure_ = false; - linear_solver_use_amg_ = false; + // linear_solver_use_amg_ = false; ilu_fillin_level_ = 0; ilu_relaxation_ = 0.9; ilu_milu_ = MILU_VARIANT::ILU; diff --git a/opm/simulators/linalg/ISTLSolverEbos.hpp b/opm/simulators/linalg/ISTLSolverEbos.hpp index 8efd1defb..d33af33b6 100644 --- a/opm/simulators/linalg/ISTLSolverEbos.hpp +++ b/opm/simulators/linalg/ISTLSolverEbos.hpp @@ -21,33 +21,18 @@ #ifndef OPM_ISTLSOLVER_EBOS_HEADER_INCLUDED #define OPM_ISTLSOLVER_EBOS_HEADER_INCLUDED -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - #include #include +#include +#include +#include +#include +#include +#include +#include +#include +#include -#include -#include -#include -#include -#include -#include - -#include #if HAVE_CUDA || HAVE_OPENCL #include @@ -82,16 +67,6 @@ public: namespace Opm { -template -DenseMatrix transposeDenseMatrix(const DenseMatrix& M) -{ - DenseMatrix tmp; - for (int i = 0; i < M.rows; ++i) - for (int j = 0; j < M.cols; ++j) - tmp[j][i] = M[i][j]; - - return tmp; -} /// This class solves the fully implicit black-oil system by /// solving the reduced system (after eliminating well variables) @@ -108,23 +83,12 @@ DenseMatrix transposeDenseMatrix(const DenseMatrix& M) using Indices = GetPropType; using WellModel = GetPropType; using Simulator = GetPropType; - typedef typename SparseMatrixAdapter::IstlMatrix Matrix; - - typedef typename SparseMatrixAdapter::MatrixBlock MatrixBlockType; - typedef typename Vector::block_type BlockVector; - using Evaluation = GetPropType; + using Matrix = typename SparseMatrixAdapter::IstlMatrix; using ThreadManager = GetPropType; - typedef typename GridView::template Codim<0>::Entity Element; using ElementContext = GetPropType; using FlexibleSolverType = Dune::FlexibleSolver; using AbstractOperatorType = Dune::AssembledLinearOperator; using WellModelOperator = WellModelAsLinearOperator; - // Due to miscibility oil <-> gas the water eqn is the one we can replace with a pressure equation. - static const bool waterEnabled = Indices::waterEnabled; - static const int pindex = (waterEnabled) ? BlackOilDefaultIndexTraits::waterCompIdx : BlackOilDefaultIndexTraits::oilCompIdx; - enum { pressureEqnIndex = pindex }; - enum { pressureVarIndex = Indices::pressureSwitchIdx }; - static const int numEq = Indices::numEq; #if HAVE_CUDA || HAVE_OPENCL static const unsigned int block_size = Matrix::block_type::rows; @@ -132,13 +96,13 @@ DenseMatrix transposeDenseMatrix(const DenseMatrix& M) #endif #if HAVE_MPI - typedef Dune::OwnerOverlapCopyCommunication communication_type; + using CommunicationType = Dune::OwnerOverlapCopyCommunication; #else - typedef Dune::CollectiveCommunication< int > communication_type; + using CommunicationType = Dune::CollectiveCommunication< int >; #endif public: - typedef Dune::AssembledLinearOperator< Matrix, Vector, Vector > AssembledLinearOperatorType; + using AssembledLinearOperatorType = Dune::AssembledLinearOperator< Matrix, Vector, Vector >; static void registerParameters() { @@ -148,29 +112,27 @@ DenseMatrix transposeDenseMatrix(const DenseMatrix& M) /// Construct a system solver. /// \param[in] parallelInformation In the case of a parallel run /// with dune-istl the information about the parallelization. - ISTLSolverEbos(const Simulator& simulator) + explicit ISTLSolverEbos(const Simulator& simulator) : simulator_(simulator), iterations_( 0 ), converged_(false), matrix_() { + const bool on_io_rank = (simulator.gridView().comm().rank() == 0); #if HAVE_MPI - comm_.reset( new communication_type( simulator_.vanguard().grid().comm() ) ); + comm_.reset( new CommunicationType( simulator_.vanguard().grid().comm() ) ); #endif parameters_.template init(); - useFlexible_ = parameters_.use_cpr_ || EWOMS_PARAM_IS_SET(TypeTag, std::string, LinearSolverConfiguration); - - if (useFlexible_) - { - prm_ = setupPropertyTree(parameters_); - } + prm_ = setupPropertyTree(parameters_); const auto& gridForConn = simulator_.vanguard().grid(); #if HAVE_CUDA || HAVE_OPENCL std::string gpu_mode = EWOMS_GET_PARAM(TypeTag, std::string, GpuMode); int platformID = EWOMS_GET_PARAM(TypeTag, int, OpenclPlatformId); int deviceID = EWOMS_GET_PARAM(TypeTag, int, BdaDeviceId); if (gridForConn.comm().size() > 1 && gpu_mode.compare("none") != 0) { - OpmLog::warning("Warning cannot use GPU with MPI, GPU is disabled"); + if (on_io_rank) { + OpmLog::warning("Warning cannot use GPU with MPI, GPU is disabled"); + } gpu_mode = "none"; } const int maxit = EWOMS_GET_PARAM(TypeTag, int, LinearSolverMaxIter); @@ -184,34 +146,33 @@ DenseMatrix transposeDenseMatrix(const DenseMatrix& M) } #endif extractParallelGridInformationToISTL(simulator_.vanguard().grid(), parallelInformation_); - useWellConn_ = EWOMS_GET_PARAM(TypeTag, bool, MatrixAddWellContributions); - ownersFirst_ = EWOMS_GET_PARAM(TypeTag, bool, OwnerCellsFirst); - interiorCellNum_ = detail::numMatrixRowsToUseInSolver(simulator_.vanguard().grid(), ownersFirst_); - if ( isParallel() && (!ownersFirst_ || parameters_.linear_solver_use_amg_) ) { - detail::setWellConnections(gridForConn, simulator_.vanguard().schedule().getWellsatEnd(), useWellConn_, wellConnectionsGraph_); - // For some reason simulator_.model().elementMapper() is not initialized at this stage - // Hence const auto& elemMapper = simulator_.model().elementMapper(); does not work. - // Set it up manually - using ElementMapper = - Dune::MultipleCodimMultipleGeomTypeMapper; - ElementMapper elemMapper(simulator_.vanguard().gridView(), Dune::mcmgElementLayout()); - detail::findOverlapAndInterior(gridForConn, elemMapper, overlapRows_, interiorRows_); - noGhostAdjacency(); - setGhostsInNoGhost(*noGhostMat_); - if (ownersFirst_) - OpmLog::warning("OwnerCellsFirst option is true, but ignored."); + // For some reason simulator_.model().elementMapper() is not initialized at this stage + // Hence const auto& elemMapper = simulator_.model().elementMapper(); does not work. + // Set it up manually + using ElementMapper = + Dune::MultipleCodimMultipleGeomTypeMapper; + ElementMapper elemMapper(simulator_.vanguard().grid().leafGridView(), Dune::mcmgElementLayout()); + detail::findOverlapAndInterior(simulator_.vanguard().grid(), elemMapper, overlapRows_, interiorRows_); + + useWellConn_ = EWOMS_GET_PARAM(TypeTag, bool, MatrixAddWellContributions); + const bool ownersFirst = EWOMS_GET_PARAM(TypeTag, bool, OwnerCellsFirst); + if (!ownersFirst) { + const std::string msg = "The linear solver no longer supports --owner-cells-first=false."; + if (on_io_rank) { + OpmLog::error(msg); + } + OPM_THROW_NOLOG(std::runtime_error, msg); } - if (useFlexible_) - { - // Print parameters to PRT/DBG logs. - if (simulator.gridView().comm().rank() == 0) { - std::ostringstream os; - os << "Property tree for linear solver:\n"; - boost::property_tree::write_json(os, prm_, true); - OpmLog::note(os.str()); - } + interiorCellNum_ = detail::numMatrixRowsToUseInSolver(simulator_.vanguard().grid(), true); + + // Print parameters to PRT/DBG logs. + if (on_io_rank) { + std::ostringstream os; + os << "Property tree for linear solver:\n"; + boost::property_tree::write_json(os, prm_, true); + OpmLog::note(os.str()); } } @@ -233,85 +194,27 @@ DenseMatrix transposeDenseMatrix(const DenseMatrix& M) #endif // update matrix entries for solvers. - if (noGhostMat_) - { - copyJacToNoGhost(M.istlMatrix(), *noGhostMat_); - } - else - { - if (firstcall) - { - // ebos will not change the matrix object. Hence simply store a pointer - // to the original one with a deleter that does nothing. - // Outch! We need to be able to scale the linear system! Hence const_cast - matrix_ = const_cast(&M.istlMatrix()); - } - else - { - // Pointers should not change - if ( &(M.istlMatrix()) != matrix_ ) - { + if (firstcall) { + // ebos will not change the matrix object. Hence simply store a pointer + // to the original one with a deleter that does nothing. + // Outch! We need to be able to scale the linear system! Hence const_cast + matrix_ = const_cast(&M.istlMatrix()); + } else { + // Pointers should not change + if ( &(M.istlMatrix()) != matrix_ ) { OPM_THROW(std::logic_error, "Matrix objects are expected to be reused when reassembling!" <<" old pointer was " << matrix_ << ", new one is " << (&M.istlMatrix()) ); - } } } rhs_ = &b; - if (useFlexible_) - { - prepareFlexibleSolver(); - } - else - { - this->scaleSystem(); + if (isParallel()) { + makeOverlapRowsInvalid(getMatrix()); } + prepareFlexibleSolver(); firstcall = false; } - void scaleSystem() - { - if (useWellConn_) { - bool form_cpr = true; - if (parameters_.system_strategy_ == "quasiimpes") { - weights_ = getQuasiImpesWeights(); - } else if (parameters_.system_strategy_ == "trueimpes") { - weights_ = getStorageWeights(); - } else if (parameters_.system_strategy_ == "simple") { - BlockVector bvec(1.0); - weights_ = getSimpleWeights(bvec); - } else if (parameters_.system_strategy_ == "original") { - BlockVector bvec(0.0); - bvec[pressureEqnIndex] = 1; - weights_ = getSimpleWeights(bvec); - } else { - if (parameters_.system_strategy_ != "none") { - OpmLog::warning("unknown_system_strategy", "Unknown linear solver system strategy: '" + parameters_.system_strategy_ + "', applying 'none' strategy."); - } - form_cpr = false; - } - if (parameters_.scale_linear_system_) { - // also scale weights - this->scaleEquationsAndVariables(weights_); - } - if (form_cpr && !(parameters_.cpr_use_drs_)) { - scaleMatrixAndRhs(weights_); - } - if (weights_.size() == 0) { - // if weights are not set cpr_use_drs_=false; - parameters_.cpr_use_drs_ = false; - } - } else { - if (parameters_.use_cpr_ && parameters_.cpr_use_drs_) { - OpmLog::warning("DRS_DISABLE", "Disabling DRS as matrix does not contain well contributions"); - } - parameters_.cpr_use_drs_ = false; - if (parameters_.scale_linear_system_) { - // also scale weights - this->scaleEquationsAndVariables(weights_); - } - } - } void setResidual(Vector& /* b */) { // rhs_ = &b; // Must be handled in prepare() instead. @@ -326,55 +229,9 @@ DenseMatrix transposeDenseMatrix(const DenseMatrix& M) } bool solve(Vector& x) { - const int verbosity = useFlexible_ ? prm_.get("verbosity", 0) : parameters_.linear_solver_verbosity_; + // Write linear system if asked for. + const int verbosity = prm_.get("verbosity", 0); const bool write_matrix = verbosity > 10; - // Solve system. - - if (useFlexible_) - { - Dune::InverseOperatorResult res; - assert(flexibleSolver_); - flexibleSolver_->apply(x, *rhs_, res); - iterations_ = res.iterations; - if (write_matrix) { - Opm::Helper::writeSystem(simulator_, //simulator is only used to get names - getMatrix(), - *rhs_, - comm_.get()); - } - - return converged_ = res.converged; - } - - const WellModel& wellModel = simulator_.problem().wellModel(); - const WellModelAsLinearOperator wellOp(wellModel); - - if( isParallel() ) - { - if ( ownersFirst_ && !parameters_.linear_solver_use_amg_ && !useFlexible_) { - typedef WellModelGhostLastMatrixAdapter< Matrix, Vector, Vector, true > Operator; - assert(matrix_); - Operator opA(getMatrix(), wellOp, interiorCellNum_); - solve( opA, x, *rhs_, *comm_ ); - } - else { - typedef WellModelMatrixAdapter< Matrix, Vector, Vector, true > Operator; - assert (noGhostMat_); - Operator opA(getMatrix(), wellOp, comm_ ); - solve( opA, x, *rhs_, *comm_ ); - } - } - else - { - typedef WellModelMatrixAdapter< Matrix, Vector, Vector, false > Operator; - Operator opA(getMatrix(), wellOp); - solve( opA, x, *rhs_ ); - } - - if (parameters_.scale_linear_system_) { - scaleSolution(x); - } - if (write_matrix) { Opm::Helper::writeSystem(simulator_, //simulator is only used to get names getMatrix(), @@ -382,6 +239,49 @@ DenseMatrix transposeDenseMatrix(const DenseMatrix& M) comm_.get()); } + // Solve system. + Dune::InverseOperatorResult result; + bool gpu_was_used = false; + + // Use GPU if: available, chosen by user, and successful. +#if HAVE_CUDA || HAVE_OPENCL + bool use_gpu = bdaBridge->getUseGpu(); + if (use_gpu) { + const std::string gpu_mode = EWOMS_GET_PARAM(TypeTag, std::string, GpuMode); + WellContributions wellContribs(gpu_mode); + if (!useWellConn_) { + simulator_.problem().wellModel().getWellContributions(wellContribs); + } + // Const_cast needed since the CUDA stuff overwrites values for better matrix condition.. + bdaBridge->solve_system(const_cast(&getMatrix()), *rhs_, wellContribs, result); + if (result.converged) { + // get result vector x from non-Dune backend, iff solve was successful + bdaBridge->get_result(x); + gpu_was_used = true; + } else { + // CPU fallback + use_gpu = bdaBridge->getUseGpu(); // update value, BdaBridge might have disabled cusparseSolver + if (use_gpu && simulator_.gridView().comm().rank() == 0) { + if (gpu_mode.compare("cusparse") == 0) { + OpmLog::warning("cusparseSolver did not converge, now trying Dune to solve current linear system..."); + } + if (gpu_mode.compare("opencl") == 0) { + OpmLog::warning("openclSolver did not converge, now trying Dune to solve current linear system..."); + } + } + } + } +#endif + + // Otherwise, use flexible istl solver. + if (!gpu_was_used) { + assert(flexibleSolver_); + flexibleSolver_->apply(x, *rhs_, result); + } + + // Check convergence, iterations etc. + checkConvergence(result); + return converged_; } @@ -399,112 +299,6 @@ DenseMatrix transposeDenseMatrix(const DenseMatrix& M) const std::any& parallelInformation() const { return parallelInformation_; } protected: - /// \brief construct the CPR preconditioner and the solver. - /// \tparam P The type of the parallel information. - /// \param parallelInformation the information about the parallelization. - template - void constructPreconditionerAndSolve(LinearOperator& linearOperator, - Vector& x, Vector& istlb, - 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; - - std::unique_ptr< MatrixOperator > opA; - - if( ! std::is_same< LinearOperator, MatrixOperator > :: value ) - { - // create new operator in case linear operator and matrix operator differ - opA.reset( CPRSelectorType::makeOperator( linearOperator.getmat(), parallelInformation_arg ) ); - } - - const double relax = parameters_.ilu_relaxation_; - const MILU_VARIANT ilu_milu = parameters_.ilu_milu_; - if ( parameters_.use_cpr_ ) - { - using MatrixType = typename MatrixOperator::matrix_type; - using CouplingMetric = Opm::Amg::Element; - using CritBase = Dune::Amg::SymmetricCriterion; - using Criterion = Dune::Amg::CoarsenCriterion; - using AMG = typename ISTLUtility - ::BlackoilAmgSelector< MatrixType, Vector, Vector,POrComm, Criterion, pressureEqnIndex, pressureVarIndex >::AMG; - - std::unique_ptr< AMG > amg; - // Construct preconditioner. - Criterion crit(15, 2000); - constructAMGPrecond( linearOperator, parallelInformation_arg, amg, opA, relax, ilu_milu ); - - // Solve. - solve(linearOperator, x, istlb, *sp, *amg, result); - } - else - { - typedef typename CPRSelectorType::AMG AMG; - std::unique_ptr< AMG > amg; - - // Construct preconditioner. - constructAMGPrecond( linearOperator, parallelInformation_arg, amg, opA, relax, ilu_milu ); - - // Solve. - solve(linearOperator, x, istlb, *sp, *amg, result); - } - } - else -#endif - { - // tries to solve linear system -#if HAVE_CUDA || HAVE_OPENCL - bool use_gpu = bdaBridge->getUseGpu(); - if (use_gpu) { - const std::string gpu_mode = EWOMS_GET_PARAM(TypeTag, std::string, GpuMode); - WellContributions wellContribs(gpu_mode); - if (!useWellConn_) { - simulator_.problem().wellModel().getWellContributions(wellContribs); - } - // Const_cast needed since the CUDA stuff overwrites values for better matrix condition.. - bdaBridge->solve_system(const_cast(&getMatrix()), istlb, wellContribs, result); - if (result.converged) { - // get result vector x from non-Dune backend, iff solve was successful - bdaBridge->get_result(x); - } else { - // CPU fallback - use_gpu = bdaBridge->getUseGpu(); // update value, BdaBridge might have disabled cusparseSolver - if (use_gpu) { - if(gpu_mode.compare("cusparse") == 0){ - OpmLog::warning("cusparseSolver did not converge, now trying Dune to solve current linear system..."); - } - - if(gpu_mode.compare("opencl") == 0){ - OpmLog::warning("openclSolver did not converge, now trying Dune to solve current linear system..."); - } - } - - // call Dune - auto precond = constructPrecond(linearOperator, parallelInformation_arg); - solve(linearOperator, x, istlb, *sp, *precond, result); - } - } else { // gpu is not selected or disabled - auto precond = constructPrecond(linearOperator, parallelInformation_arg); - solve(linearOperator, x, istlb, *sp, *precond, result); - } -#else - // Construct preconditioner. - auto precond = constructPrecond(linearOperator, parallelInformation_arg); - - // Solve. - solve(linearOperator, x, istlb, *sp, *precond, result); -#endif - } - } - // 3x3 matrix block inversion was unstable at least 2.3 until and including // 2.5.0. There may still be some issue with the 4x4 matrix block inversion @@ -514,157 +308,13 @@ DenseMatrix transposeDenseMatrix(const DenseMatrix& M) Matrix::block_type::cols> >, Vector, Vector> SeqPreconditioner; - - template - std::unique_ptr constructPrecond(Operator& opA, const Dune::Amg::SequentialInformation&) const - { - const double relax = parameters_.ilu_relaxation_; - const int ilu_fillin = parameters_.ilu_fillin_level_; - const MILU_VARIANT ilu_milu = parameters_.ilu_milu_; - const bool ilu_redblack = parameters_.ilu_redblack_; - const bool ilu_reorder_spheres = parameters_.ilu_reorder_sphere_; - auto precond = std::make_unique(opA.getmat(), ilu_fillin, relax, ilu_milu, ilu_redblack, ilu_reorder_spheres); - return precond; - } - #if HAVE_MPI typedef Dune::OwnerOverlapCopyCommunication Comm; // 3x3 matrix block inversion was unstable from at least 2.3 until and // including 2.5.0 typedef ParallelOverlappingILU0 ParPreconditioner; - template - std::unique_ptr - constructPrecond(Operator& opA, const Comm& comm) const - { - typedef std::unique_ptr Pointer; - const double relax = parameters_.ilu_relaxation_; - 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, interiorCellNum_, ilu_redblack, ilu_reorder_spheres)); - } #endif - template - void - constructAMGPrecond(LinearOperator& /* linearOperator */, const POrComm& comm, std::unique_ptr< AMG >& amg, std::unique_ptr< MatrixOperator >& opA, const double relax, const MILU_VARIANT milu) const - { - ISTLUtility::template createAMGPreconditionerPointer( *opA, relax, milu, comm, amg ); - } - - - template - void - constructAMGPrecond(LinearOperator& /* linearOperator */, const POrComm& comm, std::unique_ptr< AMG >& amg, std::unique_ptr< MatrixOperator >& opA, const double relax, - const MILU_VARIANT /* milu */ ) const - { - ISTLUtility::template createAMGPreconditionerPointer( *opA, relax, - comm, amg, parameters_, weights_ ); - } - - - /// \brief Solve the system using the given preconditioner and scalar product. - template - void solve(Operator& opA, Vector& x, Vector& istlb, ScalarProd& sp, Precond& precond, Dune::InverseOperatorResult& result) const - { - // TODO: Revise when linear solvers interface opm-core is done - // Construct linear solver. - // GMRes solver - int verbosity = 0; - if (simulator_.gridView().comm().rank() == 0) - verbosity = parameters_.linear_solver_verbosity_; - - if ( parameters_.newton_use_gmres_ ) { - Dune::RestartedGMResSolver linsolve(opA, sp, precond, - parameters_.linear_solver_reduction_, - parameters_.linear_solver_restart_, - parameters_.linear_solver_maxiter_, - verbosity); - // Solve system. - linsolve.apply(x, istlb, result); - } - else { // BiCGstab solver - Dune::BiCGSTABSolver linsolve(opA, sp, precond, - parameters_.linear_solver_reduction_, - parameters_.linear_solver_maxiter_, - verbosity); - // Solve system. - linsolve.apply(x, istlb, result); - } - } - - - /// Solve the linear system Ax = b, with A being the - /// combined derivative matrix of the residual and b - /// being the residual itself. - /// \param[in] A matrix A - /// \param[inout] x solution to be computed x - /// \param[in] b right hand side b - void solve(Matrix& A, Vector& x, Vector& b ) const - { - // Parallel version is deactivated until we figure out how to do it properly. -#if HAVE_MPI - if (parallelInformation_.type() == typeid(ParallelISTLInformation)) - { - // Construct operator, scalar product and vectors needed. - typedef Dune::OverlappingSchwarzOperator Operator; - Operator opA(A, *comm_); - solve( opA, x, b, *comm_ ); - } - else -#endif - { - // Construct operator, scalar product and vectors needed. - Dune::MatrixAdapter< Matrix, Vector, Vector> opA( A ); - solve( opA, x, b ); - } - } - - /// Solve the linear system Ax = b, with A being the - /// combined derivative matrix of the residual and b - /// being the residual itself. - /// \param[in] A matrix A - /// \param[inout] x solution to be computed x - /// \param[in] b right hand side b - template - void solve(Operator& opA OPM_UNUSED_NOMPI, - Vector& x OPM_UNUSED_NOMPI, - Vector& b OPM_UNUSED_NOMPI, - Comm& comm OPM_UNUSED_NOMPI) const - { - Dune::InverseOperatorResult result; - // Parallel version is deactivated until we figure out how to do it properly. -#if HAVE_MPI - if (parallelInformation_.type() == typeid(ParallelISTLInformation)) - { - // Construct operator, scalar product and vectors needed. - constructPreconditionerAndSolve(opA, x, b, comm, result); - } - else -#endif - { - OPM_THROW(std::logic_error,"this method if for parallel solve only"); - } - - checkConvergence( result ); - } - - /// Solve the linear system Ax = b, with A being the - /// combined derivative matrix of the residual and b - /// being the residual itself. - /// \param[in] A matrix A - /// \param[inout] x solution to be computed x - /// \param[in] b right hand side b - template - void solve(Operator& opA, Vector& x, Vector& b ) const - { - Dune::InverseOperatorResult result; - // Construct operator, scalar product and vectors needed. - Dune::Amg::SequentialInformation info; - constructPreconditionerAndSolve(opA, x, b, info, result); - checkConvergence( result ); - } - void checkConvergence( const Dune::InverseOperatorResult& result ) const { // store number of iterations @@ -688,9 +338,53 @@ DenseMatrix transposeDenseMatrix(const DenseMatrix& M) } void prepareFlexibleSolver() + { + + std::function weightsCalculator = getWeightsCalculator(); + + if (shouldCreateSolver()) { + if (isParallel()) { +#if HAVE_MPI + if (useWellConn_) { + using ParOperatorType = Dune::OverlappingSchwarzOperator; + linearOperatorForFlexibleSolver_ = std::make_unique(getMatrix(), *comm_); + flexibleSolver_ = std::make_unique(*linearOperatorForFlexibleSolver_, *comm_, prm_, weightsCalculator); + } else { + using ParOperatorType = WellModelGhostLastMatrixAdapter; + wellOperator_ = std::make_unique(simulator_.problem().wellModel()); + linearOperatorForFlexibleSolver_ = std::make_unique(getMatrix(), *wellOperator_, interiorCellNum_); + flexibleSolver_ = std::make_unique(*linearOperatorForFlexibleSolver_, *comm_, prm_, weightsCalculator); + } +#endif + } else { + if (useWellConn_) { + using SeqOperatorType = Dune::MatrixAdapter; + linearOperatorForFlexibleSolver_ = std::make_unique(getMatrix()); + flexibleSolver_ = std::make_unique(*linearOperatorForFlexibleSolver_, prm_, weightsCalculator); + } else { + using SeqOperatorType = WellModelMatrixAdapter; + wellOperator_ = std::make_unique(simulator_.problem().wellModel()); + linearOperatorForFlexibleSolver_ = std::make_unique(getMatrix(), *wellOperator_); + flexibleSolver_ = std::make_unique(*linearOperatorForFlexibleSolver_, prm_, weightsCalculator); + } + } + } + else + { + flexibleSolver_->preconditioner().update(); + } + } + + + /// Return true if we should (re)create the whole solver, + /// instead of just calling update() on the preconditioner. + bool shouldCreateSolver() const { // Decide if we should recreate the solver or just do // a minimal preconditioner update. + if (!flexibleSolver_) { + return true; + } const int newton_iteration = this->simulator_.model().newtonMethod().numIterations(); bool recreate_solver = false; if (this->parameters_.cpr_reuse_setup_ == 0) { @@ -711,177 +405,43 @@ DenseMatrix transposeDenseMatrix(const DenseMatrix& M) assert(recreate_solver == false); // Never recreate solver. } + return recreate_solver; + } + + /// Return an appropriate weight function if a cpr preconditioner is asked for. + std::function getWeightsCalculator() const + { std::function weightsCalculator; auto preconditionerType = prm_.get("preconditioner.type", "cpr"); - if( preconditionerType == "cpr" || - preconditionerType == "cprt" - ) - { - bool transpose = false; - if(preconditionerType == "cprt"){ - transpose = true; - } - - auto weightsType = prm_.get("preconditioner.weight_type", "quasiimpes"); - auto pressureIndex = this->prm_.get("preconditioner.pressure_var_index", 1); - if(weightsType == "quasiimpes") { + if (preconditionerType == "cpr" || preconditionerType == "cprt") { + const bool transpose = preconditionerType == "cprt"; + const auto weightsType = prm_.get("preconditioner.weight_type", "quasiimpes"); + const auto pressureIndex = this->prm_.get("preconditioner.pressure_var_index", 1); + if (weightsType == "quasiimpes") { // weighs will be created as default in the solver - weightsCalculator = - [this, transpose, pressureIndex](){ - return Opm::Amg::getQuasiImpesWeights(getMatrix(), - pressureIndex, - transpose); - }; - - }else if(weightsType == "trueimpes" ){ - weightsCalculator = - [this](){ - return this->getStorageWeights(); - }; - }else{ - OPM_THROW(std::invalid_argument, "Weights type " << weightsType << "not implemented for cpr." + weightsCalculator = [this, transpose, pressureIndex]() { + return Opm::Amg::getQuasiImpesWeights(this->getMatrix(), pressureIndex, transpose); + }; + } else if (weightsType == "trueimpes") { + weightsCalculator = [this, pressureIndex]() { + return this->getTrueImpesWeights(pressureIndex); + }; + } else { + OPM_THROW(std::invalid_argument, + "Weights type " << weightsType << "not implemented for cpr." << " Please use quasiimpes or trueimpes."); } } - - if (recreate_solver || !flexibleSolver_) { - if (isParallel()) { -#if HAVE_MPI - if (useWellConn_) { - assert(noGhostMat_); - using ParOperatorType = Dune::OverlappingSchwarzOperator; - linearOperatorForFlexibleSolver_ = std::make_unique(getMatrix(), *comm_); - flexibleSolver_ = std::make_unique(*linearOperatorForFlexibleSolver_, *comm_, prm_, weightsCalculator); - } else { - if (!ownersFirst_) { - OPM_THROW(std::runtime_error, "In parallel, the flexible solver requires " - "--owner-cells-first=true when --matrix-add-well-contributions=false is used."); - } - using ParOperatorType = WellModelGhostLastMatrixAdapter; - wellOperator_ = std::make_unique(simulator_.problem().wellModel()); - linearOperatorForFlexibleSolver_ = std::make_unique(getMatrix(), *wellOperator_, interiorCellNum_); - flexibleSolver_ = std::make_unique(*linearOperatorForFlexibleSolver_, *comm_, prm_, weightsCalculator); - } -#endif - } else { - if (useWellConn_) { - using SeqLinearOperator = Dune::MatrixAdapter; - linearOperatorForFlexibleSolver_ = std::make_unique(getMatrix()); - flexibleSolver_ = std::make_unique(*linearOperatorForFlexibleSolver_, prm_, weightsCalculator); - } else { - using SeqLinearOperator = WellModelMatrixAdapter; - wellOperator_ = std::make_unique(simulator_.problem().wellModel()); - linearOperatorForFlexibleSolver_ = std::make_unique(getMatrix(), *wellOperator_); - flexibleSolver_ = std::make_unique(*linearOperatorForFlexibleSolver_, prm_, weightsCalculator); - } - } - } - else - { - flexibleSolver_->preconditioner().update(); - } + return weightsCalculator; } - /// Create sparsity pattern of matrix without off-diagonal ghost entries. - void noGhostAdjacency() - { - const auto& grid = simulator_.vanguard().grid(); - const auto& gridView = simulator_.vanguard().gridView(); - // For some reason simulator_.model().elementMapper() is not initialized at this stage. - // Hence const auto& elemMapper = simulator_.model().elementMapper(); does not work. - // Set it up manually - using ElementMapper = - Dune::MultipleCodimMultipleGeomTypeMapper; - ElementMapper elemMapper(gridView, Dune::mcmgElementLayout()); - typedef typename Matrix::size_type size_type; - size_type numCells = grid.size( 0 ); - noGhostMat_.reset(new Matrix(numCells, numCells, Matrix::random)); - - std::vector> pattern; - pattern.resize(numCells); - - auto elemIt = gridView.template begin<0>(); - const auto& elemEndIt = gridView.template end<0>(); - - //Loop over cells - for (; elemIt != elemEndIt; ++elemIt) - { - const auto& elem = *elemIt; - size_type idx = elemMapper.index(elem); - pattern[idx].insert(idx); - - // Add well non-zero connections - for (auto wc = wellConnectionsGraph_[idx].begin(); wc!=wellConnectionsGraph_[idx].end(); ++wc) - pattern[idx].insert(*wc); - - // Add just a single element to ghost rows - if (elem.partitionType() != Dune::InteriorEntity) - { - noGhostMat_->setrowsize(idx, pattern[idx].size()); - } - else { - auto isend = gridView.iend(elem); - for (auto is = gridView.ibegin(elem); is!=isend; ++is) - { - //check if face has neighbor - if (is->neighbor()) - { - size_type nid = elemMapper.index(is->outside()); - pattern[idx].insert(nid); - } - } - noGhostMat_->setrowsize(idx, pattern[idx].size()); - } - } - noGhostMat_->endrowsizes(); - for (size_type dofId = 0; dofId < numCells; ++dofId) - { - auto nabIdx = pattern[dofId].begin(); - auto endNab = pattern[dofId].end(); - for (; nabIdx != endNab; ++nabIdx) - { - noGhostMat_->addindex(dofId, *nabIdx); - } - } - noGhostMat_->endindices(); - } - - /// Set the ghost diagonal in noGhost to diag(1.0) - void setGhostsInNoGhost(Matrix& ng) - { - ng=0; - typedef typename Matrix::block_type MatrixBlockTypeT; - MatrixBlockTypeT diag_block(0.0); - for (int eq = 0; eq < Matrix::block_type::rows; ++eq) - diag_block[eq][eq] = 1.0; - - //loop over precalculated ghost rows and columns - for (auto row = overlapRows_.begin(); row != overlapRows_.end(); row++ ) - { - int lcell = *row; - //diagonal block set to 1 - ng[lcell][lcell] = diag_block; - } - } - - /// Copy interior rows to noghost matrix - void copyJacToNoGhost(const Matrix& jac, Matrix& ng) - { - //Loop over precalculated interior rows. - for (auto row = interiorRows_.begin(); row != interiorRows_.end(); row++ ) - { - //Copy row - ng[*row] = jac[*row]; - } - } // Weights to make approximate pressure equations. // Calculated from the storage terms (only) of the // conservation equations, ignoring all other terms. - Vector getStorageWeights() const + Vector getTrueImpesWeights(int pressureVarIndex) const { Vector weights(rhs_->size()); ElementContext elemCtx(simulator_); @@ -891,137 +451,38 @@ DenseMatrix transposeDenseMatrix(const DenseMatrix& M) return weights; } - // Interaction between the CPR weights (the function argument 'weights') - // and the variable and equation weights from - // simulator_.model().primaryVarWeight() and - // simulator_.model().eqWeight() is nontrivial and does not work - // at the moment. Possibly refactoring of ewoms weight treatment - // is needed. In the meantime this function shows what needs to be - // done to integrate the weights properly. - void scaleEquationsAndVariables(Vector& weights) + + /// Zero out off-diagonal blocks on rows corresponding to overlap cells + /// Diagonal blocks on ovelap rows are set to diag(1.0). + void makeOverlapRowsInvalid(Matrix& matrix) const { - // loop over primary variables - const auto endi = getMatrix().end(); - for (auto i = getMatrix().begin(); i != endi; ++i) { - const auto endj = (*i).end(); - BlockVector& brhs = (*rhs_)[i.index()]; - for (auto j = (*i).begin(); j != endj; ++j) { - MatrixBlockType& block = *j; - for (std::size_t ii = 0; ii < block.rows; ii++ ) { - for (std::size_t jj = 0; jj < block.cols; jj++) { - double var_scale = simulator_.model().primaryVarWeight(i.index(),jj); - block[ii][jj] /= var_scale; - block[ii][jj] *= simulator_.model().eqWeight(i.index(), ii); - } - } + //value to set on diagonal + const int numEq = Matrix::block_type::rows; + typename Matrix::block_type diag_block(0.0); + for (int eq = 0; eq < numEq; ++eq) + diag_block[eq][eq] = 1.0; + + //loop over precalculated overlap rows and columns + for (auto row = overlapRows_.begin(); row != overlapRows_.end(); row++ ) + { + int lcell = *row; + // Zero out row. + matrix[lcell] = 0.0; + + //diagonal block set to diag(1.0). + matrix[lcell][lcell] = diag_block; } - for (std::size_t ii = 0; ii < brhs.size(); ii++) { - brhs[ii] *= simulator_.model().eqWeight(i.index(), ii); - } - if (weights.size() == getMatrix().N()) { - BlockVector& bw = weights[i.index()]; - for (std::size_t ii = 0; ii < brhs.size(); ii++) { - bw[ii] /= simulator_.model().eqWeight(i.index(), ii); - } - double abs_max = - *std::max_element(bw.begin(), bw.end(), [](double a, double b){ return std::abs(a) < std::abs(b); } ); - bw /= abs_max; - } - } } - void scaleSolution(Vector& x) - { - for (std::size_t i = 0; i < x.size(); ++i) { - auto& bx = x[i]; - for (std::size_t jj = 0; jj < bx.size(); jj++) { - double var_scale = simulator_.model().primaryVarWeight(i,jj); - bx[jj] /= var_scale; - } - } - } - - Vector getQuasiImpesWeights() - { - return Amg::getQuasiImpesWeights(getMatrix(), pressureVarIndex, /* transpose=*/ true); - } - - Vector getSimpleWeights(const BlockVector& rhs) - { - Vector weights(rhs_->size(), 0); - for (auto& bw : weights) { - bw = rhs; - } - return weights; - } - - void scaleMatrixAndRhs(const Vector& weights) - { - using Block = typename Matrix::block_type; - const auto endi = getMatrix().end(); - for (auto i = getMatrix().begin(); i !=endi; ++i) { - const BlockVector& bweights = weights[i.index()]; - BlockVector& brhs = (*rhs_)[i.index()]; - const auto endj = (*i).end(); - for (auto j = (*i).begin(); j != endj; ++j) { - // assume it is something on all rows - Block& block = (*j); - BlockVector neweq(0.0); - for (std::size_t ii = 0; ii < block.rows; ii++) { - for (std::size_t jj = 0; jj < block.cols; jj++) { - neweq[jj] += bweights[ii]*block[ii][jj]; - } - } - block[pressureEqnIndex] = neweq; - } - Scalar newrhs(0.0); - for (std::size_t ii = 0; ii < brhs.size(); ii++) { - newrhs += bweights[ii]*brhs[ii]; - } - brhs[pressureEqnIndex] = newrhs; - } - } - - static void multBlocksInMatrix(Matrix& ebosJac, const MatrixBlockType& trans, const bool left = true) - { - const int n = ebosJac.N(); - for (int row_index = 0; row_index < n; ++row_index) { - auto& row = ebosJac[row_index]; - auto* dataptr = row.getptr(); - for (int elem = 0; elem < row.N(); ++elem) { - auto& block = dataptr[elem]; - if (left) { - block = block.leftmultiply(trans); - } else { - block = block.rightmultiply(trans); - } - } - } - } - - static void multBlocksVector(Vector& ebosResid_cp, const MatrixBlockType& leftTrans) - { - for (auto& bvec : ebosResid_cp) { - auto bvec_new = bvec; - leftTrans.mv(bvec, bvec_new); - bvec = bvec_new; - } - } - - static void scaleCPRSystem(Matrix& M_cp, Vector& b_cp, const MatrixBlockType& leftTrans) - { - multBlocksInMatrix(M_cp, leftTrans, true); - multBlocksVector(b_cp, leftTrans); - } Matrix& getMatrix() { - return noGhostMat_ ? *noGhostMat_ : *matrix_; + return *matrix_; } const Matrix& getMatrix() const { - return noGhostMat_ ? *noGhostMat_ : *matrix_; + return *matrix_; } const Simulator& simulator_; @@ -1031,7 +492,6 @@ DenseMatrix transposeDenseMatrix(const DenseMatrix& M) // non-const to be able to scale the linear system Matrix* matrix_; - std::unique_ptr noGhostMat_; Vector *rhs_; std::unique_ptr flexibleSolver_; @@ -1041,17 +501,14 @@ DenseMatrix transposeDenseMatrix(const DenseMatrix& M) std::vector interiorRows_; std::vector> wellConnectionsGraph_; - bool ownersFirst_; bool useWellConn_; - bool useFlexible_; size_t interiorCellNum_; FlowLinearSolverParameters parameters_; boost::property_tree::ptree prm_; - Vector weights_; bool scale_variables_; - std::shared_ptr< communication_type > comm_; + std::shared_ptr< CommunicationType > comm_; }; // end ISTLSolver } // namespace Opm diff --git a/opm/simulators/linalg/setupPropertyTree_impl.hpp b/opm/simulators/linalg/setupPropertyTree_impl.hpp index 25f85b4a8..3e0f6ec57 100644 --- a/opm/simulators/linalg/setupPropertyTree_impl.hpp +++ b/opm/simulators/linalg/setupPropertyTree_impl.hpp @@ -50,18 +50,13 @@ setupPropertyTree(const FlowLinearSolverParameters& p) #else OPM_THROW(std::invalid_argument, "--linear-solver-configuration=file not supported with " - << "boost version. Needs versoin > 1.48."); + << "boost version. Needs version > 1.48."); #endif } else { std::string conf = p.linear_solver_configuration_; // Support old UseCpr if not configuration was set - if (!EWOMS_PARAM_IS_SET(TypeTag, std::string, LinearSolverConfiguration) && p.use_cpr_) - { - conf = "cpr_trueimpes"; - } - if((conf == "cpr_trueimpes") || (conf == "cpr_quasiimpes")){ prm.put("tol", p.linear_solver_reduction_); if (EWOMS_PARAM_IS_SET(TypeTag, int, LinearSolverMaxIter)) diff --git a/tests/test_blackoil_amg.cpp b/tests/test_blackoil_amg.cpp index 490f205ce..8ef0f09fa 100644 --- a/tests/test_blackoil_amg.cpp +++ b/tests/test_blackoil_amg.cpp @@ -21,13 +21,15 @@ #define BOOST_TEST_MODULE BlackoilAmgTest #define BOOST_TEST_NO_MAIN #include -#include +#include +#include #include #include #include #include #include +#include #if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 7) #include #else @@ -35,6 +37,7 @@ #endif #include #include +#include class MPIError { public: @@ -260,6 +263,7 @@ M setupAnisotropic2d(int N, Dune::ParallelIndexSet& indices, const Dune:: return mat; } + //BOOST_AUTO_TEST_CASE(runBlackoilAmgLaplace) void runBlackoilAmgLaplace() { @@ -289,30 +293,18 @@ void runBlackoilAmgLaplace() setBoundary(x, b, N, comm.indexSet()); Operator fop(mat, comm); - typedef Dune::Amg::CoarsenCriterion > - Criterion; - typedef Dune::SeqSSOR Smoother; - //typedef Dune::SeqJac Smoother; - //typedef Dune::SeqILU0 Smoother; - //typedef Dune::SeqILUn Smoother; - typedef Dune::BlockPreconditioner ParSmoother; - typedef typename Dune::Amg::SmootherTraits::Arguments SmootherArgs; - Dune::OverlappingSchwarzScalarProduct sp(comm); - Dune::InverseOperatorResult r; - SmootherArgs smootherArgs; - Criterion criterion; - smootherArgs.iterations = 1; - Opm::CPRParameter param; + boost::property_tree::ptree prm; + prm.put("type", "amg"); + std::function weights = [&mat]() { + return Opm::Amg::getQuasiImpesWeights(mat, 0, false); + }; + auto amg = Opm::PreconditionerFactory::create(fop, prm, weights, comm); + + Dune::CGSolver amgCG(fop, sp, *amg, 10e-8, 300, (ccomm.rank()==0) ? 2 : 0); - Opm::BlackoilAmg amg(param, - {}, - fop, criterion, - smootherArgs, - comm); - Dune::CGSolver amgCG(fop, sp, amg, 10e-8, 300, (ccomm.rank()==0) ? 2 : 0); amgCG.apply(x,b,r); }