From 0bae240a42199a937e9a0cee228882e673f7a193 Mon Sep 17 00:00:00 2001 From: Markus Blatt Date: Thu, 17 May 2018 12:01:02 +0200 Subject: [PATCH] Allow user to choose modified ILU0 decomposition. Using the parameter ilu_milu=true|false (default=false) the user can now choose to use the modified ILU0 decomposition. If selected values will not be dropped for nonzero entries but added to the diagonal of U. This approach will result in A*e = L*U*e for vector e with all entries beging 1. --- opm/autodiff/CPRPreconditioner.hpp | 59 +++++++++--- opm/autodiff/ISTLSolver.hpp | 32 ++++--- .../NewtonIterationBlackoilInterleaved.hpp | 3 + opm/autodiff/ParallelOverlappingILU0.hpp | 93 +++++++++++++++---- 4 files changed, 146 insertions(+), 41 deletions(-) diff --git a/opm/autodiff/CPRPreconditioner.hpp b/opm/autodiff/CPRPreconditioner.hpp index 2a425f7b5..78f3c269e 100644 --- a/opm/autodiff/CPRPreconditioner.hpp +++ b/opm/autodiff/CPRPreconditioner.hpp @@ -64,6 +64,7 @@ struct CPRParameter double cpr_relax_; double cpr_solver_tol_; int cpr_ilu_n_; + bool cpr_ilu_milu_; int cpr_max_ell_iter_; bool cpr_use_amg_; bool cpr_use_bicgstab_; @@ -80,6 +81,7 @@ struct CPRParameter 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_ilu_milu_ = param.getDefault("ilu_milu", cpr_ilu_milu_); 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_); @@ -92,6 +94,7 @@ struct CPRParameter cpr_relax_ = 1.0; cpr_solver_tol_ = 1e-2; cpr_ilu_n_ = 0; + cpr_ilu_milu_ = false; cpr_max_ell_iter_ = 25; cpr_use_amg_ = true; cpr_use_bicgstab_ = true; @@ -103,6 +106,31 @@ struct CPRParameter 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, + bool 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. /// @@ -167,9 +195,9 @@ struct CPRSelector //! \param relax The relaxation factor to use. template std::shared_ptr > -createILU0Ptr(const M& A, const C& comm, double relax) +createILU0Ptr(const M& A, const C& comm, double relax, bool milu) { - return std::make_shared >(A, comm, relax); + 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. @@ -177,25 +205,28 @@ createILU0Ptr(const M& A, const C& comm, double relax) //! \param relax The relaxation factor to use. template std::shared_ptr > -createILUnPtr(const M& A, const C& comm, int ilu_n, double relax) +createILUnPtr(const M& A, const C& comm, int ilu_n, double relax, bool milu) { - return std::make_shared >( A, comm, ilu_n, relax); + 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, - const P& comm) + bool milu, const P& comm) { typedef typename CPRSelector ::EllipticPreconditioner ParallelPreconditioner; typedef typename CPRSelector ::EllipticPreconditionerPointer EllipticPreconditionerPointer; - return EllipticPreconditionerPointer(new ParallelPreconditioner(Ae, comm, relax)); + return EllipticPreconditionerPointer(new ParallelPreconditioner(Ae, comm, relax, milu)); } template < class C, class Op, class P, class S, std::size_t index > @@ -220,13 +251,14 @@ createAMGPreconditionerPointer(Op& opA, const double relax, const P& comm, SmootherArgs smootherArgs; smootherArgs.iterations = 1; smootherArgs.relaxationFactor = relax; + setILUParameters(smootherArgs, params); 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 bool milu, const P& comm, std::unique_ptr< AMG >& amgPtr) { // TODO: revise choice of parameters int coarsenTarget=1200; @@ -243,6 +275,7 @@ createAMGPreconditionerPointer(Op& opA, const double relax, const P& comm, std:: SmootherArgs smootherArgs; smootherArgs.iterations = 1; smootherArgs.relaxationFactor = relax; + setILUParameters(smootherArgs, milu); amgPtr.reset( new AMG(opA, criterion, smootherArgs, comm ) ); } @@ -254,7 +287,7 @@ createAMGPreconditionerPointer(Op& opA, const double relax, const P& comm, std:: // \param amgPtr The unique_ptr to be filled (return) template < int pressureIndex=0, 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 bool milu, const P& comm, std::unique_ptr< AMG >& amgPtr ) { // type of matrix typedef typename Op::matrix_type M; @@ -268,7 +301,7 @@ createAMGPreconditionerPointer( Op& opA, const double relax, const P& comm, std: // The coarsening criterion used in the AMG typedef Dune::Amg::CoarsenCriterion Criterion; - createAMGPreconditionerPointer(opA, relax, comm, amgPtr); + createAMGPreconditionerPointer(opA, relax, milu, comm, amgPtr); } } // end namespace ISTLUtility @@ -370,10 +403,10 @@ createAMGPreconditionerPointer( Op& opA, const double relax, const P& comm, std: // create the preconditioner for the whole system. if( param_.cpr_ilu_n_ == 0 ) { - pre_ = ISTLUtility::createILU0Ptr( A_, comm, param_.cpr_relax_ ); + 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_); + pre_ = ISTLUtility::createILUnPtr( A_, comm, param_.cpr_ilu_n_, param_.cpr_relax_, param_.cpr_ilu_milu_); } } @@ -531,11 +564,11 @@ createAMGPreconditionerPointer( Op& opA, const double relax, const P& comm, std: { if( amg ) { - ISTLUtility::createAMGPreconditionerPointer( *opAe_ , param_.cpr_relax_, comm, amg_ ); + ISTLUtility::createAMGPreconditionerPointer( *opAe_ , param_.cpr_relax_, param_.cpr_ilu_milu_, comm, amg_ ); } else { - precond_ = ISTLUtility::createEllipticPreconditionerPointer( Ae_, param_.cpr_relax_, comm); + precond_ = ISTLUtility::createEllipticPreconditionerPointer( Ae_, param_.cpr_relax_, param_.cpr_ilu_milu_, comm); } } }; diff --git a/opm/autodiff/ISTLSolver.hpp b/opm/autodiff/ISTLSolver.hpp index dc31183fb..e253f68e2 100644 --- a/opm/autodiff/ISTLSolver.hpp +++ b/opm/autodiff/ISTLSolver.hpp @@ -440,6 +440,7 @@ namespace Detail } const double relax = parameters_.ilu_relaxation_; + const bool ilu_milu = parameters_.ilu_milu_; if ( parameters_.use_cpr_ ) { using Matrix = typename MatrixOperator::matrix_type; @@ -452,7 +453,7 @@ namespace Detail std::unique_ptr< AMG > amg; // Construct preconditioner. Criterion crit(15, 2000); - constructAMGPrecond( linearOperator, parallelInformation_arg, amg, opA, relax ); + constructAMGPrecond( linearOperator, parallelInformation_arg, amg, opA, relax, ilu_milu ); // Solve. solve(linearOperator, x, istlb, *sp, *amg, result); @@ -463,7 +464,7 @@ namespace Detail std::unique_ptr< AMG > amg; // Construct preconditioner. - constructAMGPrecond( linearOperator, parallelInformation_arg, amg, opA, relax ); + constructAMGPrecond( linearOperator, parallelInformation_arg, amg, opA, relax, ilu_milu ); // Solve. solve(linearOperator, x, istlb, *sp, *amg, result); @@ -495,7 +496,8 @@ namespace Detail { const double relax = parameters_.ilu_relaxation_; const int ilu_fillin = parameters_.ilu_fillin_level_; - std::unique_ptr precond(new SeqPreconditioner(opA.getmat(), ilu_fillin, relax)); + const bool ilu_milu = parameters_.ilu_milu_; + std::unique_ptr precond(new SeqPreconditioner(opA.getmat(), ilu_fillin, relax, ilu_milu)); return precond; } @@ -517,38 +519,44 @@ namespace Detail { typedef std::unique_ptr Pointer; const double relax = parameters_.ilu_relaxation_; - return Pointer(new ParPreconditioner(opA.getmat(), comm, relax)); + const bool ilu_milu = parameters_.ilu_milu_; + return Pointer(new ParPreconditioner(opA.getmat(), comm, relax, ilu_milu)); } #endif template void - constructAMGPrecond(LinearOperator& /* linearOperator */, const POrComm& comm, std::unique_ptr< AMG >& amg, std::unique_ptr< MatrixOperator >& opA, const double relax ) const + constructAMGPrecond(LinearOperator& /* linearOperator */, const POrComm& comm, std::unique_ptr< AMG >& amg, std::unique_ptr< MatrixOperator >& opA, const double relax, const bool milu) const { - ISTLUtility::template createAMGPreconditionerPointer( *opA, relax, comm, amg ); + ISTLUtility::template createAMGPreconditionerPointer( *opA, relax, milu, comm, amg ); } template void - constructAMGPrecond(MatrixOperator& opA, const POrComm& comm, std::unique_ptr< AMG >& amg, std::unique_ptr< MatrixOperator >&, const double relax ) const + constructAMGPrecond(MatrixOperator& opA, const POrComm& comm, std::unique_ptr< AMG >& amg, std::unique_ptr< MatrixOperator >&, const double relax, + const bool milu) const { - ISTLUtility::template createAMGPreconditionerPointer( opA, relax, comm, amg ); + 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 + constructAMGPrecond(LinearOperator& /* linearOperator */, const POrComm& comm, std::unique_ptr< AMG >& amg, std::unique_ptr< MatrixOperator >& opA, const double relax, + const bool milu ) const { - ISTLUtility::template createAMGPreconditionerPointer( *opA, relax, comm, amg, parameters_ ); + ISTLUtility::template createAMGPreconditionerPointer( *opA, relax, + comm, amg, parameters_ ); } template void - constructAMGPrecond(MatrixOperator& opA, const POrComm& comm, std::unique_ptr< AMG >& amg, std::unique_ptr< MatrixOperator >&, const double relax ) const + constructAMGPrecond(MatrixOperator& opA, const POrComm& comm, std::unique_ptr< AMG >& amg, std::unique_ptr< MatrixOperator >&, const double relax, const bool milu ) const { - ISTLUtility::template createAMGPreconditionerPointer( opA, relax, comm, amg, parameters_ ); + ISTLUtility::template createAMGPreconditionerPointer( opA, relax, milu, + 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 377aecd1d..d50a2123d 100644 --- a/opm/autodiff/NewtonIterationBlackoilInterleaved.hpp +++ b/opm/autodiff/NewtonIterationBlackoilInterleaved.hpp @@ -44,6 +44,7 @@ namespace Opm int linear_solver_restart_; int linear_solver_verbosity_; int ilu_fillin_level_; + bool ilu_milu_; bool newton_use_gmres_; bool require_full_sparsity_pattern_; bool ignoreConvergenceFailure_; @@ -69,6 +70,7 @@ namespace Opm linear_solver_use_amg_ = param.getDefault("linear_solver_use_amg", linear_solver_use_amg_ ); ilu_relaxation_ = param.getDefault("ilu_relaxation", ilu_relaxation_ ); ilu_fillin_level_ = param.getDefault("ilu_fillin_level", ilu_fillin_level_ ); + ilu_milu_ = param.getDefault("ilu_milu", ilu_milu_); // Check whether to use cpr approach const std::string cprSolver = "cpr"; @@ -89,6 +91,7 @@ namespace Opm linear_solver_use_amg_ = false; ilu_fillin_level_ = 0; ilu_relaxation_ = 0.9; + ilu_milu_ = false; } }; diff --git a/opm/autodiff/ParallelOverlappingILU0.hpp b/opm/autodiff/ParallelOverlappingILU0.hpp index 43339d9b5..cbfc6088c 100644 --- a/opm/autodiff/ParallelOverlappingILU0.hpp +++ b/opm/autodiff/ParallelOverlappingILU0.hpp @@ -37,6 +37,34 @@ namespace Opm template class ParallelOverlappingILU0; +template +class ParallelOverlappingILU0Args + : public Dune::Amg::DefaultSmootherArgs +{ + public: + ParallelOverlappingILU0Args(bool milu = false) + : milu_(milu) + {} + void setMilu(bool milu) + { + milu_ = milu; + } + bool getMilu() const + { + return milu_; + } + void setN(int n) + { + n_ = n; + } + int getN() const + { + return n_; + } + private: + bool milu_; + int n_; +}; } // end namespace Opm namespace Dune @@ -45,6 +73,13 @@ namespace Dune namespace Amg { + +template +struct SmootherTraits > +{ + using Arguments = Opm::ParallelOverlappingILU0Args; +}; + /// \brief Tells AMG how to construct the Opm::ParallelOverlappingILU0 smoother /// \tparam Matrix The type of the Matrix. /// \tparam Domain The type of the Vector representing the domain. @@ -54,17 +89,18 @@ namespace Amg template struct ConstructionTraits > { - typedef Dune::SeqILU0 T; + typedef Opm::ParallelOverlappingILU0 T; typedef DefaultParallelConstructionArgs Arguments; - typedef ConstructionTraits SeqConstructionTraits; static inline Opm::ParallelOverlappingILU0* construct(Arguments& args) { - return new Opm::ParallelOverlappingILU0(args.getMatrix(), - args.getComm(), - args.getArgs().relaxationFactor); + return new T(args.getMatrix(), + args.getComm(), + args.getArgs().getN(), + args.getArgs().relaxationFactor, + args.getArgs().getMilu()); } - static inline void deconstruct(Opm::ParallelOverlappingILU0* bp) + static inline void deconstruct(T* bp) { delete bp; } @@ -330,10 +366,14 @@ public: \param A The matrix to operate on. \param n ILU fill in level (for testing). This does not work in parallel. \param w The relaxation factor. + \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). */ template ParallelOverlappingILU0 (const Dune::BCRSMatrix& A, - const int n, const field_type w ) + const int n, const field_type w, + bool milu) : lower_(), upper_(), inv_(), @@ -342,7 +382,7 @@ public: { // BlockMatrix is a Subclass of FieldMatrix that just adds // methods. Therefore this cast should be safe. - init( reinterpret_cast(A), n ); + init( reinterpret_cast(A), n, milu ); } /*! \brief Constructor gets all parameters to operate the prec. @@ -350,10 +390,14 @@ public: \param comm communication object, e.g. Dune::OwnerOverlapCopyCommunication \param n ILU fill in level (for testing). This does not work in parallel. \param w The relaxation factor. + \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). */ template ParallelOverlappingILU0 (const Dune::BCRSMatrix& A, - const ParallelInfo& comm, const int n, const field_type w ) + const ParallelInfo& comm, const int n, const field_type w, + bool milu) : lower_(), upper_(), inv_(), @@ -362,7 +406,7 @@ public: { // BlockMatrix is a Subclass of FieldMatrix that just adds // methods. Therefore this cast should be safe. - init( reinterpret_cast(A), n ); + init( reinterpret_cast(A), n, milu ); } /*! \brief Constructor. @@ -370,11 +414,14 @@ public: Constructor gets all parameters to operate the prec. \param A The matrix to operate on. \param w The relaxation factor. + \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). */ template ParallelOverlappingILU0 (const Dune::BCRSMatrix& A, - const field_type w) - : ParallelOverlappingILU0( A, 0, w ) + const field_type w, bool milu) + : ParallelOverlappingILU0( A, 0, w, milu ) { } @@ -384,10 +431,14 @@ public: \param A The matrix to operate on. \param comm communication object, e.g. Dune::OwnerOverlapCopyCommunication \param w The relaxation factor. + \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). */ template ParallelOverlappingILU0 (const Dune::BCRSMatrix& A, - const ParallelInfo& comm, const field_type w) + const ParallelInfo& comm, const field_type w, + bool milu) : lower_(), upper_(), inv_(), @@ -396,7 +447,7 @@ public: { // BlockMatrix is a Subclass of FieldMatrix that just adds // methods. Therefore this cast should be safe. - init( reinterpret_cast(A), 0 ); + init( reinterpret_cast(A), 0, milu ); } /*! @@ -490,7 +541,7 @@ public: } protected: - void init( const Matrix& A, const int iluIteration ) + void init( const Matrix& A, const int iluIteration, bool milu ) { // (For older DUNE versions the communicator might be // invalid if redistribution in AMG happened on the coarset level. @@ -519,11 +570,21 @@ protected: if( iluIteration == 0 ) { // create ILU-0 decomposition ILU.reset( new Matrix( A ) ); - bilu0_decomposition( *ILU ); + if ( milu ) + { + detail::milu0_decomposition ( *ILU ); + }else + { + bilu0_decomposition( *ILU ); + } } else { // create ILU-n decomposition ILU.reset( new Matrix( A.N(), A.M(), Matrix::row_wise) ); + if ( milu ) + { + OpmLog::warning("MILU_N_NOT_SUPPORTED", "MILU variant only supported for zero fill-in"); + } bilu_decomposition( A, iluIteration, *ILU ); } }