diff --git a/opm/autodiff/BlackoilAmgCpr.hpp b/opm/autodiff/BlackoilAmgCpr.hpp index e35c53589..25adc1389 100644 --- a/opm/autodiff/BlackoilAmgCpr.hpp +++ b/opm/autodiff/BlackoilAmgCpr.hpp @@ -35,563 +35,511 @@ #include #include #include -namespace Dune -{ - namespace Amg - { - template - class UnSymmetricCriterion; - } -} - -namespace Dune -{ - - template - class MatrixBlock; - -} namespace Opm { - namespace Detail - { - template - std::unique_ptr scaleMatrixDRSPtr(const Operator& op, - const Communication& comm, - 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; - std::unique_ptr matrix(new Matrix(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;//, createOperator(op, *matrix, comm)); - } - } - - template - class OneComponentAggregationLevelTransferPolicyCpr; - - template + + template class OneComponentAggregationLevelTransferPolicyCpr - : public Dune::Amg::LevelTransferPolicyCpr::value> + : public Dune::Amg::LevelTransferPolicyCpr::value> { - typedef Dune::Amg::AggregatesMap AggregatesMap; + typedef Dune::Amg::AggregatesMap AggregatesMap; public: - using CoarseOperator = typename Detail::ScalarType::value; - typedef Dune::Amg::LevelTransferPolicy FatherType; - typedef Communication ParallelInformation; + using CoarseOperator = typename Detail::ScalarType::value; + typedef Dune::Amg::LevelTransferPolicy FatherType; + typedef Communication ParallelInformation; public: - OneComponentAggregationLevelTransferPolicyCpr(const Criterion& crit, const Communication& comm) - : criterion_(crit), communication_(&const_cast(comm)) - {} + OneComponentAggregationLevelTransferPolicyCpr(const Criterion& crit, const Communication& comm) + : criterion_(crit), communication_(&const_cast(comm)) + {} - void createCoarseLevelSystem(const Operator& fineOperator) - { - prolongDamp_ = 1; + void createCoarseLevelSystem(const Operator& fineOperator) + { + prolongDamp_ = 1; - 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(); + 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; - } + 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(); + 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*){}); + 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; - OperatorArgs oargs(*coarseLevelMatrix_, *coarseLevelCommunication_); - this->operator_.reset(Dune::Amg::ConstructionTraits::construct(oargs)); - } + this->lhs_.resize(this->coarseLevelMatrix_->M()); + this->rhs_.resize(this->coarseLevelMatrix_->N()); + using OperatorArgs = typename Dune::Amg::ConstructionTraits::Arguments; + OperatorArgs oargs(*coarseLevelMatrix_, *coarseLevelCommunication_); + this->operator_.reset(Dune::Amg::ConstructionTraits::construct(oargs)); + } - // compleately unsafe!!!!!! - void calculateCoarseEntries(const Operator& fineOperator)//const M& fineMatrix) - { - 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]; - } - } - } + // compleately unsafe!!!!!! + void calculateCoarseEntries(const Operator& fineOperator)//const M& fineMatrix) + { + 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]; + } + } + } - //template - // void calculateCoarseEntriesOld(const Operator& fineOperator)//const M& fineMatrix) - // { - // 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][COMPONENT_INDEX]; - // } - // } - // } - // } - // } + //template + // void calculateCoarseEntriesOld(const Operator& fineOperator)//const M& fineMatrix) + // { + // 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][COMPONENT_INDEX]; + // } + // } + // } + // } + // } - void moveToCoarseLevel(const typename FatherType::FineRangeType& fine) - { - // Set coarse vector to zero - this->rhs_=0; + void moveToCoarseLevel(const typename FatherType::FineRangeType& fine) + { + // Set coarse vector to zero + this->rhs_=0; - auto end = fine.end(), begin=fine.begin(); + auto end = fine.end(), begin=fine.begin(); - for(auto block=begin; block != end; ++block) - { - this->rhs_[block-begin] = (*block)[COMPONENT_INDEX]; - } + for(auto block=begin; block != end; ++block) + { + this->rhs_[block-begin] = (*block)[COMPONENT_INDEX]; + } - this->lhs_=0; - } + this->lhs_=0; + } - void moveToFineLevel(typename FatherType::FineDomainType& fine) - { + void moveToFineLevel(typename FatherType::FineDomainType& fine) + { - auto end=fine.end(), begin=fine.begin(); + auto end=fine.end(), begin=fine.begin(); - for(auto block=begin; block != end; ++block) - { - (*block)[COMPONENT_INDEX] = this->lhs_[block-begin]; - } + for(auto block=begin; block != end; ++block) + { + (*block)[COMPONENT_INDEX] = this->lhs_[block-begin]; + } - } + } - OneComponentAggregationLevelTransferPolicyCpr* clone() const - { - return new OneComponentAggregationLevelTransferPolicyCpr(*this); - } + OneComponentAggregationLevelTransferPolicyCpr* clone() const + { + return new OneComponentAggregationLevelTransferPolicyCpr(*this); + } - const Communication& getCoarseLevelCommunication() const - { - return *coarseLevelCommunication_; - } + 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_; + typename Operator::matrix_type::field_type prolongDamp_; + //std::shared_ptr aggregatesMap_; + Criterion criterion_; + Communication* communication_; + std::shared_ptr coarseLevelCommunication_; + std::shared_ptr coarseLevelMatrix_; }; - namespace Detail - { + namespace Detail + { + /** + * @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 OneStepAMGCoarseSolverPolicyNoSolve + { + 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. + */ + OneStepAMGCoarseSolverPolicyNoSolve(const CPRParameter* param, const SmootherArgs& args, const Criterion& c) + : param_(param), smootherArgs_(args), criterion_(c) + {} + /** @brief Copy constructor. */ + OneStepAMGCoarseSolverPolicyNoSolve(const OneStepAMGCoarseSolverPolicyNoSolve& 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, + typename AMGType::Operator& op, + const Criterion& crit, + const typename AMGType::SmootherArgs& args, + const Communication& comm) + : param_(param), amg_(), op_(op), crit_(crit), args_(args), comm_(comm) + { + amg_.reset(new AMGType(op, crit,args, comm)); + } + + void updateAmgPreconditioner(typename AMGType::Operator& op){ + //op_ = op; + //amg_->recalculateHierarchy(); + amg_->updateSolver(crit_, op, comm_); + //amg_.reset(new AMGType(op, crit_,args_, comm_)); + //amg_->recalculateGalerkin(); + } +#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6) + Dune::SolverCategory::Category category() const override + { + return std::is_same::value ? + Dune::SolverCategory::sequential : Dune::SolverCategory::overlapping; + } +#endif + + + void apply(X& x, X& b, double reduction, Dune::InverseOperatorResult& res) override + { + DUNE_UNUSED_PARAMETER(reduction); + DUNE_UNUSED_PARAMETER(res); +#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6) + auto sp = Dune::createScalarProduct(comm_, op_.category()); +#else + using Chooser = Dune::ScalarProductChooser; + auto sp = Chooser::construct(comm_); +#endif + Dune::Preconditioner* prec = amg_.get(); + // 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; + if ( param_->cpr_ell_solvetype_ == 0 ) + { + // Category of preconditioner will be checked at compile time. Therefore we need + // to cast to the derived class +#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6) + Dune::BiCGSTABSolver solver(const_cast(op_), *sp, *prec, + tolerance, maxit, verbosity); +#else + Dune::BiCGSTABSolver solver(const_cast(op_), *sp, + reinterpret_cast(*prec), + tolerance, maxit, verbosity); +#endif + solver.apply(x,b,res); + + } + else if (param_->cpr_ell_solvetype_ == 1) + { + // Category of preconditioner will be checked at compile time. Therefore we need + // to cast to the derived class +#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6) + Dune::CGSolver solver(const_cast(op_), *sp, *prec, + tolerance, maxit, verbosity); +#else + Dune::CGSolver solver(const_cast(op_), *sp, + reinterpret_cast(*prec), + tolerance, maxit, verbosity); +#endif + solver.apply(x,b,res); + } + else + { + // X v(x); + // prec->pre(x,b); + // op_->applyscaleadd(-1,x,b); + // v = 0; + // prec->apply(v,b); + // x += v; + // op_->applyscaleadd(-1,x,b); + // prec->post(x); +#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6) + Dune::LoopSolver solver(const_cast(op_), *sp, *prec, + tolerance, maxit, verbosity); +#else + Dune::LoopSolver solver(const_cast(op_), *sp, + reinterpret_cast(*prec), + tolerance, maxit, verbosity); +#endif + solver.apply(x,b,res); + } +#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6) + +#else + delete sp; +#endif + } + + 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::unique_ptr op_; + typename AMGType::Operator& op_; + Criterion crit_; + typename AMGType::SmootherArgs args_; + 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 + void setCoarseOperator(LTP& transferPolicy){ + coarseOperator_= transferPolicy.getCoarseLevelOperator(); + } + 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; //std::shared_ptr >(inv); + + } + //void recalculateGalerkin(){ + // coarseOperator_.recalculateHierarchy(); + //} + 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_; + }; + + + } // end namespace Detail + + /** - * @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. + * \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 OneStepAMGCoarseSolverPolicyNoSolve + template + class BlackoilAmgCpr + : public Dune::Preconditioner { 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. - */ - OneStepAMGCoarseSolverPolicyNoSolve(const CPRParameter* param, const SmootherArgs& args, const Criterion& c) - : param_(param), smootherArgs_(args), criterion_(c) - {} - /** @brief Copy constructor. */ - OneStepAMGCoarseSolverPolicyNoSolve(const OneStepAMGCoarseSolverPolicyNoSolve& 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, - typename AMGType::Operator& op, - const Criterion& crit, - const typename AMGType::SmootherArgs& args, - const Communication& comm) - : param_(param), amg_(), op_(op), crit_(crit), args_(args), comm_(comm) - { - amg_.reset(new AMGType(op, crit,args, comm)); - } + /** \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; - void updateAmgPreconditioner(typename AMGType::Operator& op){ - //op_ = op; - //amg_->recalculateHierarchy(); - amg_->updateSolver(crit_, op, comm_); - //amg_.reset(new AMGType(op, crit_,args_, comm_)); - //amg_->recalculateGalerkin(); - } + 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 = + OneComponentAggregationLevelTransferPolicyCpr; + using CoarseSolverPolicy = + Detail::OneStepAMGCoarseSolverPolicyNoSolve; + using TwoLevelMethod = + Dune::Amg::TwoLevelMethodCpr; + public: #if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6) Dune::SolverCategory::Category category() const override { return std::is_same::value ? Dune::SolverCategory::sequential : Dune::SolverCategory::overlapping; } -#endif - - - void apply(X& x, X& b, double reduction, Dune::InverseOperatorResult& res) override - { - DUNE_UNUSED_PARAMETER(reduction); - DUNE_UNUSED_PARAMETER(res); -#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6) - auto sp = Dune::createScalarProduct(comm_, op_.category()); -#else - using Chooser = Dune::ScalarProductChooser; - auto sp = Chooser::construct(comm_); -#endif - Dune::Preconditioner* prec = amg_.get(); - // 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; - if ( param_->cpr_ell_solvetype_ == 0 ) - { - // Category of preconditioner will be checked at compile time. Therefore we need - // to cast to the derived class -#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6) - Dune::BiCGSTABSolver solver(const_cast(op_), *sp, *prec, - tolerance, maxit, verbosity); -#else - Dune::BiCGSTABSolver solver(const_cast(op_), *sp, - reinterpret_cast(*prec), - tolerance, maxit, verbosity); -#endif - solver.apply(x,b,res); - - } - else if (param_->cpr_ell_solvetype_ == 1) - { - // Category of preconditioner will be checked at compile time. Therefore we need - // to cast to the derived class -#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6) - Dune::CGSolver solver(const_cast(op_), *sp, *prec, - tolerance, maxit, verbosity); -#else - Dune::CGSolver solver(const_cast(op_), *sp, - reinterpret_cast(*prec), - tolerance, maxit, verbosity); -#endif - solver.apply(x,b,res); - } - else - { - // X v(x); - // prec->pre(x,b); - // op_->applyscaleadd(-1,x,b); - // v = 0; - // prec->apply(v,b); - // x += v; - // op_->applyscaleadd(-1,x,b); - // prec->post(x); -#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6) - Dune::LoopSolver solver(const_cast(op_), *sp, *prec, - tolerance, maxit, verbosity); -#else - Dune::LoopSolver solver(const_cast(op_), *sp, - reinterpret_cast(*prec), - tolerance, maxit, verbosity); -#endif - solver.apply(x,b,res); - } -#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6) - #else - delete sp; -#endif - } - - 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::unique_ptr op_; - typename AMGType::Operator& op_; - Criterion crit_; - typename AMGType::SmootherArgs args_; - 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 - void setCoarseOperator(LTP& transferPolicy){ - coarseOperator_= transferPolicy.getCoarseLevelOperator(); - } - 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; //std::shared_ptr >(inv); - - } - //void recalculateGalerkin(){ - // coarseOperator_.recalculateHierarchy(); - //} - 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_; - }; - - - } // end namespace Detail - - - /** - * \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 first - * coarse level system, either by simply extracting the coupling between the components at COMPONENT_INDEX - * in the matrix blocks or by extracting them and 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 the pressure). - */ - template - class BlackoilAmgCpr - : 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 CoarseSmoother = typename Detail::ScalarType::value; - using FineCriterion = - typename Detail::OneComponentCriterionType::value; - using CoarseCriterion = typename Detail::ScalarType::value; - using LevelTransferPolicy = - OneComponentAggregationLevelTransferPolicyCpr; - using CoarseSolverPolicy = - Detail::OneStepAMGCoarseSolverPolicyNoSolve; - using TwoLevelMethod = - Dune::Amg::TwoLevelMethodCpr; - public: -#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6) - Dune::SolverCategory::Category category() const override - { - return std::is_same::value ? - Dune::SolverCategory::sequential : Dune::SolverCategory::overlapping; - } -#else - // define the category - enum { - //! \brief The category the precondtioner is part of. - category = Operator::category - }; + // define the category + enum { + //! \brief The category the precondtioner is part of. + category = Operator::category + }; #endif - /** - * \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. - */ - BlackoilAmgCpr(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::scaleMatrixDRSPtr(fineOperator, comm, - COMPONENT_INDEX, weights_, param)), - scaledMatrixOperator_(Detail::createOperatorPtr(fineOperator, *scaledMatrix_, comm)), - smoother_(Detail::constructSmoother(*scaledMatrixOperator_, - smargs, comm)), - levelTransferPolicy_(criterion, comm), - coarseSolverPolicy_(¶m, smargs, criterion), - twoLevelMethod_(*scaledMatrixOperator_, - smoother_, - levelTransferPolicy_, - coarseSolverPolicy_, 0, 1) - { - } - void updatePreconditioner(const typename TwoLevelMethod::FineDomainType& weights, - const Operator& fineOperator, - const SmootherArgs& smargs, - const Communication& comm){ - weights_ = weights; - *scaledMatrix_ = *Detail::scaleMatrixDRSPtr(fineOperator, comm, - COMPONENT_INDEX, weights_, param_); - //*scaledMatrixOperator_ = *Detail::createOperatorPtr(fineOperator,*scaledMatrix_,comm); - smoother_ .reset(Detail::constructSmoother(*scaledMatrixOperator_, - smargs, comm)); - twoLevelMethod_.updatePreconditioner(*scaledMatrixOperator_, - smoother_, - coarseSolverPolicy_); - } + /** + * \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. + */ + BlackoilAmgCpr(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), + coarseSolverPolicy_(¶m, smargs, criterion), + twoLevelMethod_(*scaledMatrixOperator_, + smoother_, + levelTransferPolicy_, + coarseSolverPolicy_, 0, 1) + { + } + + void updatePreconditioner(const typename TwoLevelMethod::FineDomainType& weights, + const Operator& fineOperator, + const SmootherArgs& smargs, + const Communication& comm) + { + weights_ = weights; + *scaledMatrix_ = *Detail::scaleMatrixDRS(fineOperator, COMPONENT_INDEX, weights_, param_); + smoother_.reset(Detail::constructSmoother(*scaledMatrixOperator_, smargs, comm)); + twoLevelMethod_.updatePreconditioner(*scaledMatrixOperator_, + smoother_, + coarseSolverPolicy_); + } - void pre(typename TwoLevelMethod::FineDomainType& x, - typename TwoLevelMethod::FineRangeType& b) override - { - twoLevelMethod_.pre(x,b); - } + 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 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_; - typename TwoLevelMethod::FineDomainType weights_;//make copy - std::unique_ptr scaledMatrix_; - std::unique_ptr scaledMatrixOperator_; - //Operator scaledMatrixOperator_; - //std::tuple, Operator> - std::shared_ptr smoother_; - LevelTransferPolicy levelTransferPolicy_; - CoarseSolverPolicy coarseSolverPolicy_; - TwoLevelMethod twoLevelMethod_; - //BlockVector weights_; - }; + 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_; + typename TwoLevelMethod::FineDomainType weights_; //make copy + std::unique_ptr scaledMatrix_; + std::unique_ptr scaledMatrixOperator_; + std::shared_ptr smoother_; + LevelTransferPolicy levelTransferPolicy_; + CoarseSolverPolicy coarseSolverPolicy_; + TwoLevelMethod twoLevelMethod_; + }; } // end namespace Opm #endif