diff --git a/opm/autodiff/BlackoilAmg.hpp b/opm/autodiff/BlackoilAmg.hpp index 72ac125c2..d1030b949 100644 --- a/opm/autodiff/BlackoilAmg.hpp +++ b/opm/autodiff/BlackoilAmg.hpp @@ -316,13 +316,13 @@ public: * @param args The arguments used for constructing the smoother. * @param c The crition used for the aggregation within AMG. */ - OneStepAMGCoarseSolverPolicy(const SmootherArgs& args, const Criterion& c) - : smootherArgs_(args), criterion_(c), use_amg_(true) + OneStepAMGCoarseSolverPolicy(const CPRParameter* param, const SmootherArgs& args, const Criterion& c) + : param_(param), smootherArgs_(args), criterion_(c) {} /** @brief Copy constructor. */ OneStepAMGCoarseSolverPolicy(const OneStepAMGCoarseSolverPolicy& other) - : coarseOperator_(other.coarseOperator_), smootherArgs_(other.smootherArgs_), - criterion_(other.criterion_), use_amg_(other.use_amg_) + : param_(other.param_), coarseOperator_(other.coarseOperator_), smootherArgs_(other.smootherArgs_), + criterion_(other.criterion_) {} private: /** @@ -333,14 +333,14 @@ private: */ struct AMGInverseOperator : public Dune::InverseOperator { - AMGInverseOperator(const typename AMGType::Operator& op, + AMGInverseOperator(const CPRParameter* param, + const typename AMGType::Operator& op, const Criterion& crit, const typename AMGType::SmootherArgs& args, - const Communication& comm, - bool use_amg) - : amg_(), smoother_(), op_(op), comm_(comm), first_(true) + const Communication& comm) + : param_(param), amg_(), smoother_(), op_(op), comm_(comm), first_(true) { - if ( use_amg ) + if ( param_->cpr_use_amg_ ) { amg_.reset(new AMGType(op, crit,args, comm)); } @@ -386,8 +386,24 @@ private: { prec = smoother_.get(); } - Dune::BiCGSTABSolver solver(const_cast(op_), *sp, *prec , 1e-2, 25, 0); - solver.apply(x,b,res); + // 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_use_bicgstab_ ) + { + Dune::BiCGSTABSolver solver(const_cast(op_), *sp, *prec, + tolerance, maxit, verbosity); + solver.apply(x,b,res); + } + else + { + Dune::CGSolver solver(const_cast(op_), *sp, *prec, + tolerance, maxit, verbosity); + solver.apply(x,b,res); + } + #if ! DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6) delete sp; #endif @@ -410,6 +426,7 @@ private: { } private: + const CPRParameter* param_; X x_; std::unique_ptr amg_; std::unique_ptr smoother_; @@ -434,11 +451,11 @@ public: coarseOperator_=transferPolicy.getCoarseLevelOperator(); const LevelTransferPolicy& transfer = reinterpret_cast(transferPolicy); - AMGInverseOperator* inv = new AMGInverseOperator(*coarseOperator_, + AMGInverseOperator* inv = new AMGInverseOperator(param_, + *coarseOperator_, criterion_, smootherArgs_, - transfer.getCoarseLevelCommunication(), - use_amg_); + transfer.getCoarseLevelCommunication()); return inv; //std::shared_ptr >(inv); @@ -447,11 +464,12 @@ public: 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_; - bool use_amg_; }; template @@ -856,14 +874,16 @@ public: }; #endif - BlackoilAmg(const Operator& fineOperator, const Criterion& criterion, + BlackoilAmg(const CPRParameter& param, + const Operator& fineOperator, const Criterion& criterion, const SmootherArgs& smargs, const Communication& comm) - : scaledMatrixOperator_(Detail::scaleMatrixQuasiImpes(fineOperator, comm, + : param_(param), + scaledMatrixOperator_(Detail::scaleMatrixQuasiImpes(fineOperator, comm, COMPONENT_INDEX)), smoother_(Detail::constructSmoother(std::get<1>(scaledMatrixOperator_), smargs, comm)), levelTransferPolicy_(criterion, comm), - coarseSolverPolicy_(smargs, criterion), + coarseSolverPolicy_(¶m, smargs, criterion), twoLevelMethod_(std::get<1>(scaledMatrixOperator_), smoother_, levelTransferPolicy_, coarseSolverPolicy_, 0, 1) @@ -888,6 +908,7 @@ public: twoLevelMethod_.apply(v, scaledD); } private: + const CPRParameter& param_; std::tuple, Operator> scaledMatrixOperator_; std::shared_ptr smoother_; LevelTransferPolicy levelTransferPolicy_; diff --git a/opm/autodiff/CPRPreconditioner.hpp b/opm/autodiff/CPRPreconditioner.hpp index de1f92ded..7df676d4f 100644 --- a/opm/autodiff/CPRPreconditioner.hpp +++ b/opm/autodiff/CPRPreconditioner.hpp @@ -52,6 +52,12 @@ namespace Opm { +struct CPRParameter; + +template +class BlackoilAmg; + namespace ISTLUtility { /// @@ -149,9 +155,36 @@ createEllipticPreconditionerPointer(const M& Ae, double relax, return EllipticPreconditionerPointer(new ParallelPreconditioner(Ae, comm, relax)); } +template < class C, class Op, class P, class S, std::size_t index > +inline void +createAMGPreconditionerPointer(Op& opA, const double relax, const P& comm, + std::unique_ptr< BlackoilAmg >& amgPtr, + const CPRParameter& params = CPRParameter()) +{ + using AMG = BlackoilAmg; + // 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; + + amgPtr.reset( new AMG( params, opA, criterion, smootherArgs, comm ) ); +} + template < class C, class Op, class P, class AMG > inline void -createAMGPreconditionerPointer(Op& opA, const double relax, const P& comm, std::unique_ptr< AMG >& amgPtr ) +createAMGPreconditionerPointer(Op& opA, const double relax, const P& comm, std::unique_ptr< AMG >& amgPtr, + const CPRParameter& params = CPRParameter()) { // TODO: revise choice of parameters int coarsenTarget=1200; diff --git a/opm/autodiff/FlowMainEbos.hpp b/opm/autodiff/FlowMainEbos.hpp index 1f3f4f3bd..8578aa6ee 100755 --- a/opm/autodiff/FlowMainEbos.hpp +++ b/opm/autodiff/FlowMainEbos.hpp @@ -610,7 +610,14 @@ namespace Opm void setupLinearSolver() { typedef typename BlackoilModelEbos :: ISTLSolverType ISTLSolverType; - + const std::string cprSolver = "cpr"; + if (!param_.has("solver_approach") ) + { + if ( eclState().getSimulationConfig().useCPR() ) + { + param_.insertParameter("solver_approach", cprSolver); + } + } extractParallelGridInformationToISTL(grid(), parallel_information_); fis_solver_.reset( new ISTLSolverType( param_, parallel_information_ ) ); } diff --git a/opm/autodiff/ISTLSolver.hpp b/opm/autodiff/ISTLSolver.hpp index cc8089e05..ca2a12e3d 100644 --- a/opm/autodiff/ISTLSolver.hpp +++ b/opm/autodiff/ISTLSolver.hpp @@ -404,7 +404,7 @@ namespace Opm parallelInformation_arg.copyOwnerToAll(istlb, istlb); #if FLOW_SUPPORT_AMG // activate AMG if either flow_ebos is used or UMFPack is not available - if( parameters_.linear_solver_use_amg_ ) + if( parameters_.linear_solver_use_amg_ || parameters_.use_cpr_) { typedef ISTLUtility::CPRSelector< Matrix, Vector, Vector, POrComm> CPRSelectorType; typedef typename CPRSelectorType::Operator MatrixOperator; @@ -418,7 +418,7 @@ namespace Opm } const double relax = parameters_.ilu_relaxation_; - if ( ! parameters_.amg_blackoil_system_ ) + if ( ! parameters_.use_cpr_ ) { typedef typename CPRSelectorType::AMG AMG; std::unique_ptr< AMG > amg; @@ -519,7 +519,7 @@ namespace Opm void constructAMGPrecond(LinearOperator& /* linearOperator */, const POrComm& comm, std::unique_ptr< AMG >& amg, std::unique_ptr< MatrixOperator >& opA, const double relax ) const { - ISTLUtility::template createAMGPreconditionerPointer( *opA, relax, comm, amg ); + ISTLUtility::template createAMGPreconditionerPointer( *opA, relax, comm, amg, parameters_ ); } @@ -527,7 +527,7 @@ namespace Opm void constructAMGPrecond(MatrixOperator& opA, const POrComm& comm, std::unique_ptr< AMG >& amg, std::unique_ptr< MatrixOperator >&, const double relax ) const { - ISTLUtility::template createAMGPreconditionerPointer( opA, relax, comm, amg ); + ISTLUtility::template createAMGPreconditionerPointer( opA, relax, comm, amg, parameters_ ); } /// \brief Solve the system using the given preconditioner and scalar product. template diff --git a/opm/autodiff/NewtonIterationBlackoilInterleaved.hpp b/opm/autodiff/NewtonIterationBlackoilInterleaved.hpp index ab6a71304..377aecd1d 100644 --- a/opm/autodiff/NewtonIterationBlackoilInterleaved.hpp +++ b/opm/autodiff/NewtonIterationBlackoilInterleaved.hpp @@ -25,6 +25,7 @@ #ifndef OPM_NEWTONITERATIONBLACKOILINTERLEAVED_HEADER_INCLUDED #define OPM_NEWTONITERATIONBLACKOILINTERLEAVED_HEADER_INCLUDED +#include #include #include @@ -35,6 +36,7 @@ namespace Opm { /// This class carries all parameters for the NewtonIterationBlackoilInterleaved class struct NewtonIterationBlackoilInterleavedParameters + : public CPRParameter { double linear_solver_reduction_; double ilu_relaxation_; @@ -46,11 +48,12 @@ namespace Opm bool require_full_sparsity_pattern_; bool ignoreConvergenceFailure_; bool linear_solver_use_amg_; - bool amg_blackoil_system_; + bool use_cpr_; NewtonIterationBlackoilInterleavedParameters() { reset(); } // read values from parameter class NewtonIterationBlackoilInterleavedParameters( const ParameterGroup& param ) + : CPRParameter(param) { // set default parameters reset(); @@ -64,15 +67,18 @@ namespace Opm require_full_sparsity_pattern_ = param.getDefault("require_full_sparsity_pattern", require_full_sparsity_pattern_); ignoreConvergenceFailure_ = param.getDefault("linear_solver_ignoreconvergencefailure", ignoreConvergenceFailure_); linear_solver_use_amg_ = param.getDefault("linear_solver_use_amg", linear_solver_use_amg_ ); - amg_blackoil_system_ = param.getDefault("amg_blackoil_system", amg_blackoil_system_); ilu_relaxation_ = param.getDefault("ilu_relaxation", ilu_relaxation_ ); ilu_fillin_level_ = param.getDefault("ilu_fillin_level", ilu_fillin_level_ ); + + // Check whether to use cpr approach + const std::string cprSolver = "cpr"; + use_cpr_ = ( param.getDefault("solver_approach", std::string()) == cprSolver ); } // set default values void reset() { - amg_blackoil_system_ = false; + use_cpr_ = false; newton_use_gmres_ = false; linear_solver_reduction_ = 1e-2; linear_solver_maxiter_ = 150;