diff --git a/opm/autodiff/CPRPreconditioner.hpp b/opm/autodiff/CPRPreconditioner.hpp index 3a58b0013..344752810 100644 --- a/opm/autodiff/CPRPreconditioner.hpp +++ b/opm/autodiff/CPRPreconditioner.hpp @@ -27,6 +27,7 @@ #include #include +#include #include #include @@ -39,6 +40,7 @@ #include #include #include +#include #include @@ -93,6 +95,10 @@ struct CPRSelector typedef Dune::SeqILU0 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 Dune::Amg::FastAMG AMG; + /// \brief creates an Operator from the matrix /// \param M The matrix to use. /// \param p The parallel information to use. @@ -119,6 +125,9 @@ struct CPRSelector > ParallelPreconditionerDeleter > > 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. @@ -256,6 +265,41 @@ createEllipticPreconditionerPointer(const M& Ae, double relax, #endif } // end namespace + struct CPRParameter + { + double cpr_relax_; + double cpr_solver_tol_; + int cpr_ilu_n_; + int cpr_max_ell_iter_; + bool cpr_use_amg_; + bool cpr_use_bicgstab_; + + CPRParameter() { reset(); } + + CPRParameter( const parameter::ParameterGroup& param) + { + // reset values to default + reset(); + + cpr_relax_ = param.getDefault("cpr_relax", cpr_relax_); + cpr_solver_tol_ = param.getDefault("cpr_solver_tol", cpr_solver_tol_); + cpr_ilu_n_ = param.getDefault("cpr_ilu_n", cpr_ilu_n_); + cpr_max_ell_iter_ = param.getDefault("cpr_max_elliptic_iter",cpr_max_ell_iter_); + cpr_use_amg_ = param.getDefault("cpr_use_amg", cpr_use_amg_); + cpr_use_bicgstab_ = param.getDefault("cpr_use_bicgstab", cpr_use_bicgstab_); + } + + void reset() + { + cpr_relax_ = 1.0; + cpr_solver_tol_ = 1e-4; + cpr_ilu_n_ = 0; + cpr_max_ell_iter_ = 5000; + cpr_use_amg_ = false; + cpr_use_bicgstab_ = true; + } + }; + /*! \brief CPR preconditioner. @@ -313,8 +357,7 @@ createEllipticPreconditionerPointer(const M& Ae, double relax, EllipticPreconditionerPointer; //! \brief amg preconditioner for the elliptic system - typedef EllipticPreconditioner Smoother; - typedef Dune::Amg::AMG AMG; + typedef typename CPRSelector::AMG AMG; /*! \brief Constructor. @@ -327,12 +370,10 @@ createEllipticPreconditionerPointer(const M& Ae, double relax, \param paralleInformation The information about the parallelization, if this is a parallel run */ - CPRPreconditioner (const M& A, const M& Ae, const field_type relax, - const unsigned int ilu_n, - const bool useAMG, - const bool useBiCG, + CPRPreconditioner (const CPRParameter& param, const M& A, const M& Ae, const ParallelInformation& comm=ParallelInformation()) - : A_(A), + : param_( param ), + A_(A), Ae_(Ae), de_( Ae_.N() ), ve_( Ae_.M() ), @@ -342,18 +383,16 @@ createEllipticPreconditionerPointer(const M& Ae, double relax, amg_(), // amg preconditioner for elliptic system pre_(), // copy A will be made be the preconditioner vilu_( A_.N() ), - relax_(relax), - use_bicg_solver_( useBiCG ), comm_(comm) { // create appropriate preconditioner for elliptic system - createPreconditioner( useAMG, comm ); + createPreconditioner( param_.cpr_use_amg_, comm ); - if( ilu_n == 0 ) { - pre_ = createILU0Ptr( A_, relax_, comm ); + if( param_.cpr_ilu_n_ == 0 ) { + pre_ = createILU0Ptr( A_, param_.cpr_relax_, comm ); } else { - pre_ = createILUnPtr( A_, ilu_n, relax_, comm ); + pre_ = createILUnPtr( A_, param_.cpr_ilu_n_, param_.cpr_relax_, comm ); } } @@ -396,11 +435,11 @@ createEllipticPreconditionerPointer(const M& Ae, double relax, pre_->apply( vilu_, dmodified_); // don't apply relaxation if relax_ == 1 - if( std::abs( relax_ - 1.0 ) < 1e-12 ) { + if( std::abs( param_.cpr_relax_ - 1.0 ) < 1e-12 ) { v += vilu_; } else { - v *= relax_; + v *= param_.cpr_relax_; v += vilu_; } } @@ -418,8 +457,8 @@ createEllipticPreconditionerPointer(const M& Ae, double relax, void solveElliptic(Y& x, Y& de) { // Linear solver parameters - const double tolerance = 1e-4; - const int maxit = 5000; + const double tolerance = param_.cpr_solver_tol_; + const int maxit = param_.cpr_max_ell_iter_; const int verbosity = 0; // operator result containing iterations etc. @@ -435,7 +474,7 @@ createEllipticPreconditionerPointer(const M& Ae, double relax, if( amg_ ) { // Solve system with AMG - if( use_bicg_solver_ ) { + if( param_.cpr_use_bicgstab_ ) { Dune::BiCGSTABSolver linsolve(*opAe_, *sp, (*amg_), tolerance, maxit, verbosity); linsolve.apply(x, de, result); } @@ -448,7 +487,7 @@ createEllipticPreconditionerPointer(const M& Ae, double relax, { assert( precond_ ); // Solve system with ILU-0 - if( use_bicg_solver_ ) { + if( param_.cpr_use_bicgstab_ ) { Dune::BiCGSTABSolver linsolve(*opAe_, *sp, (*precond_), tolerance, maxit, verbosity); linsolve.apply(x, de, result); } @@ -464,6 +503,9 @@ createEllipticPreconditionerPointer(const M& Ae, double relax, } } + //! \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. @@ -491,12 +533,6 @@ createEllipticPreconditionerPointer(const M& Ae, double relax, //! \brief temporary variables for ILU solve Y vilu_; - //! \brief The relaxation factor to use. - field_type relax_; - - //! \brief true if ISTL BiCGSTABSolver is used, otherwise ISTL CGSolver is used - const bool use_bicg_solver_; - //! \brief The information about the parallelization const P& comm_; protected: @@ -504,30 +540,35 @@ createEllipticPreconditionerPointer(const M& Ae, double relax, { if( amg ) { - typedef Dune::Amg::CoarsenCriterion< Dune::Amg::SymmetricCriterion > Criterion; - typedef typename Dune::Amg::SmootherTraits::Arguments SmootherArgs; + //! \brief The coupling metric used in the AMG + typedef Dune::Amg::FirstDiagonal CouplingMetric; - SmootherArgs smootherArgs; + //! \brief The coupling criterion used in the AMG + typedef Dune::Amg::SymmetricCriterion CritBase; - smootherArgs.iterations = 1; - smootherArgs.relaxationFactor = relax_; + //! \brief The coarsening criterion used in the AMG + typedef Dune::Amg::CoarsenCriterion Criterion; int coarsenTarget=1200; Criterion criterion(15,coarsenTarget); criterion.setDebugLevel( 0 ); // no debug information, 1 for printing hierarchy information criterion.setDefaultValuesIsotropic(2); - criterion.setAlpha(.67); - criterion.setBeta(1.0e-6); - criterion.setMaxLevel(10); - amg_ = std::unique_ptr< AMG > (new AMG(*opAe_, criterion, smootherArgs)); + //criterion.setAlpha(.67); + //criterion.setGamma(1); + //criterion.setBeta(1.0e-6); + //criterion.setMaxLevel(10); + criterion.setNoPostSmoothSteps( 1 ); + criterion.setNoPreSmoothSteps( 1 ); + amg_ = std::unique_ptr< AMG > (new AMG(*opAe_, criterion));//, smootherArgs)); } else { - precond_ = createEllipticPreconditionerPointer( Ae_, relax_, comm); + precond_ = createEllipticPreconditionerPointer( Ae_, param_.cpr_relax_, comm); } } }; + } // namespace Opm #endif // OPM_CPRPRECONDITIONER_HEADER_INCLUDED diff --git a/opm/autodiff/NewtonIterationBlackoilCPR.cpp b/opm/autodiff/NewtonIterationBlackoilCPR.cpp index 229b48de0..56dc6df44 100644 --- a/opm/autodiff/NewtonIterationBlackoilCPR.cpp +++ b/opm/autodiff/NewtonIterationBlackoilCPR.cpp @@ -99,13 +99,10 @@ namespace Opm /// Construct a system solver. NewtonIterationBlackoilCPR::NewtonIterationBlackoilCPR(const parameter::ParameterGroup& param, const boost::any& parallelInformation) - : iterations_( 0 ), parallelInformation_(parallelInformation) + : cpr_param_( param ), + iterations_( 0 ), + parallelInformation_(parallelInformation) { - cpr_relax_ = param.getDefault("cpr_relax", 1.0); - cpr_ilu_n_ = param.getDefault("cpr_ilu_n", 0); - cpr_use_amg_ = param.getDefault("cpr_use_amg", false); - cpr_use_bicgstab_ = param.getDefault("cpr_use_bicgstab", true); - linear_solver_reduction_ = param.getDefault("linear_solver_reduction", 1e-3 ); linear_solver_maxiter_ = param.getDefault("linear_solver_maxiter", 150 ); linear_solver_restart_ = param.getDefault("linear_solver_restart", 40 ); diff --git a/opm/autodiff/NewtonIterationBlackoilCPR.hpp b/opm/autodiff/NewtonIterationBlackoilCPR.hpp index 3723ca8ba..80fb53954 100644 --- a/opm/autodiff/NewtonIterationBlackoilCPR.hpp +++ b/opm/autodiff/NewtonIterationBlackoilCPR.hpp @@ -46,6 +46,7 @@ namespace Opm typedef Dune::FieldMatrix MatrixBlockType; typedef Dune::BCRSMatrix Mat; typedef Dune::BlockVector Vector; + public: /// Construct a system solver. @@ -91,9 +92,9 @@ namespace Opm sp(ScalarProductChooser::construct(parallelInformation)); // Construct preconditioner. // typedef Dune::SeqILU0 Preconditioner; - typedef Opm::CPRPreconditioner Preconditioner; + typedef Opm::CPRPreconditioner Preconditioner; parallelInformation.copyOwnerToAll(istlb, istlb); - Preconditioner precond(opA.getmat(), istlAe, cpr_relax_, cpr_ilu_n_, cpr_use_amg_, cpr_use_bicgstab_, parallelInformation); + Preconditioner precond(cpr_param_, opA.getmat(), istlAe, parallelInformation); // TODO: Revise when linear solvers interface opm-core is done // Construct linear solver. @@ -112,11 +113,9 @@ namespace Opm } } + CPRParameter cpr_param_; + mutable int iterations_; - double cpr_relax_; - unsigned int cpr_ilu_n_; - bool cpr_use_amg_; - bool cpr_use_bicgstab_; bool newton_use_gmres_; boost::any parallelInformation_;