From 63058559bc8f3776818d2128bab839725ecb1f64 Mon Sep 17 00:00:00 2001 From: Markus Blatt Date: Mon, 2 Jul 2018 12:29:55 +0200 Subject: [PATCH] Added various other variants of MILU. These versions are inspired by the ones used in SuperLU and the enums to choose them have simuilar names, but without leading S (MILU_1-MILU_3). The following variants are supported (chosen by the enum MILU_VARIANT): ILU: plain ILU MILU_1: lump diagonal with dropped row entries. MILU_2: lump diagonal with the sum of the absolute values of the dropped row entries. MILU_3: if diagonal is positive add sum of dropped row entrires. Otherwise substract them. MILU_4: if diagonal is positive add sum of dropped row entrires. Otherwise do nothing --- opm/autodiff/CPRPreconditioner.hpp | 21 +- opm/autodiff/ISTLSolver.hpp | 14 +- .../NewtonIterationBlackoilInterleaved.hpp | 8 +- opm/autodiff/ParallelOverlappingILU0.hpp | 184 ++++++++++++++---- 4 files changed, 170 insertions(+), 57 deletions(-) diff --git a/opm/autodiff/CPRPreconditioner.hpp b/opm/autodiff/CPRPreconditioner.hpp index 78f3c269e..d52df9c96 100644 --- a/opm/autodiff/CPRPreconditioner.hpp +++ b/opm/autodiff/CPRPreconditioner.hpp @@ -29,6 +29,7 @@ #include #include +#include #include #include @@ -64,7 +65,7 @@ struct CPRParameter double cpr_relax_; double cpr_solver_tol_; int cpr_ilu_n_; - bool cpr_ilu_milu_; + MILU_VARIANT cpr_ilu_milu_; int cpr_max_ell_iter_; bool cpr_use_amg_; bool cpr_use_bicgstab_; @@ -81,12 +82,14 @@ 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_); cpr_solver_verbose_ = param.getDefault("cpr_solver_verbose", cpr_solver_verbose_); cpr_pressure_aggregation_ = param.getDefault("cpr_pressure_aggregation", cpr_pressure_aggregation_); + + std::string milu("ILU"); + cpr_ilu_milu_ = convertString2Milu(param.getDefault("ilu_milu", milu)); } void reset() @@ -94,7 +97,7 @@ struct CPRParameter cpr_relax_ = 1.0; cpr_solver_tol_ = 1e-2; cpr_ilu_n_ = 0; - cpr_ilu_milu_ = false; + cpr_ilu_milu_ = MILU_VARIANT::ILU; cpr_max_ell_iter_ = 25; cpr_use_amg_ = true; cpr_use_bicgstab_ = true; @@ -117,7 +120,7 @@ void setILUParameters(Opm::ParallelOverlappingILU0Args& args, template void setILUParameters(Opm::ParallelOverlappingILU0Args& args, - bool milu, int n=0) + MILU_VARIANT milu, int n=0) { args.setN(n); args.setMilu(milu); @@ -195,7 +198,7 @@ struct CPRSelector //! \param relax The relaxation factor to use. template std::shared_ptr > -createILU0Ptr(const M& A, const C& comm, double relax, bool milu) +createILU0Ptr(const M& A, const C& comm, double relax, MILU_VARIANT milu) { return std::make_shared >(A, comm, relax, milu); } @@ -205,7 +208,7 @@ createILU0Ptr(const M& A, const C& comm, double relax, bool milu) //! \param relax The relaxation factor to use. template std::shared_ptr > -createILUnPtr(const M& A, const C& comm, int ilu_n, double relax, bool milu) +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 ); } @@ -219,7 +222,7 @@ createILUnPtr(const M& A, const C& comm, int ilu_n, double relax, bool milu) template typename CPRSelector::EllipticPreconditionerPointer createEllipticPreconditionerPointer(const M& Ae, double relax, - bool milu, const P& comm) + MILU_VARIANT milu, const P& comm) { typedef typename CPRSelector ::EllipticPreconditioner ParallelPreconditioner; @@ -258,7 +261,7 @@ createAMGPreconditionerPointer(Op& opA, const double relax, const P& comm, template < class C, class Op, class P, class AMG > inline void -createAMGPreconditionerPointer(Op& opA, const double relax, const bool milu, const P& comm, std::unique_ptr< AMG >& amgPtr) +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; @@ -287,7 +290,7 @@ createAMGPreconditionerPointer(Op& opA, const double relax, const bool milu, con // \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 bool milu, const P& comm, std::unique_ptr< AMG >& amgPtr ) +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; diff --git a/opm/autodiff/ISTLSolver.hpp b/opm/autodiff/ISTLSolver.hpp index e253f68e2..745b9e60e 100644 --- a/opm/autodiff/ISTLSolver.hpp +++ b/opm/autodiff/ISTLSolver.hpp @@ -440,7 +440,7 @@ namespace Detail } const double relax = parameters_.ilu_relaxation_; - const bool ilu_milu = parameters_.ilu_milu_; + const MILU_VARIANT ilu_milu = parameters_.ilu_milu_; if ( parameters_.use_cpr_ ) { using Matrix = typename MatrixOperator::matrix_type; @@ -496,7 +496,7 @@ namespace Detail { const double relax = parameters_.ilu_relaxation_; const int ilu_fillin = parameters_.ilu_fillin_level_; - const bool ilu_milu = parameters_.ilu_milu_; + const MILU_VARIANT ilu_milu = parameters_.ilu_milu_; std::unique_ptr precond(new SeqPreconditioner(opA.getmat(), ilu_fillin, relax, ilu_milu)); return precond; } @@ -519,14 +519,14 @@ namespace Detail { typedef std::unique_ptr Pointer; const double relax = parameters_.ilu_relaxation_; - const bool ilu_milu = parameters_.ilu_milu_; + const MILU_VARIANT 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 bool milu) const + 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 ); } @@ -535,7 +535,7 @@ namespace Detail template void constructAMGPrecond(MatrixOperator& opA, const POrComm& comm, std::unique_ptr< AMG >& amg, std::unique_ptr< MatrixOperator >&, const double relax, - const bool milu) const + const MILU_VARIANT milu) const { ISTLUtility::template createAMGPreconditionerPointer( opA, relax, milu, comm, amg ); @@ -544,7 +544,7 @@ namespace Detail template void constructAMGPrecond(LinearOperator& /* linearOperator */, const POrComm& comm, std::unique_ptr< AMG >& amg, std::unique_ptr< MatrixOperator >& opA, const double relax, - const bool milu ) const + const MILU_VARIANT milu ) const { ISTLUtility::template createAMGPreconditionerPointer( *opA, relax, comm, amg, parameters_ ); @@ -553,7 +553,7 @@ namespace Detail template void - constructAMGPrecond(MatrixOperator& opA, const POrComm& comm, std::unique_ptr< AMG >& amg, std::unique_ptr< MatrixOperator >&, const double relax, const bool milu ) const + constructAMGPrecond(MatrixOperator& opA, const POrComm& comm, std::unique_ptr< AMG >& amg, std::unique_ptr< MatrixOperator >&, const double relax, const MILU_VARIANT milu ) const { ISTLUtility::template createAMGPreconditionerPointer( opA, relax, milu, comm, amg, parameters_ ); diff --git a/opm/autodiff/NewtonIterationBlackoilInterleaved.hpp b/opm/autodiff/NewtonIterationBlackoilInterleaved.hpp index d50a2123d..668f2b4f8 100644 --- a/opm/autodiff/NewtonIterationBlackoilInterleaved.hpp +++ b/opm/autodiff/NewtonIterationBlackoilInterleaved.hpp @@ -28,6 +28,7 @@ #include #include #include +#include #include #include @@ -44,7 +45,7 @@ namespace Opm int linear_solver_restart_; int linear_solver_verbosity_; int ilu_fillin_level_; - bool ilu_milu_; + Opm::MILU_VARIANT ilu_milu_; bool newton_use_gmres_; bool require_full_sparsity_pattern_; bool ignoreConvergenceFailure_; @@ -70,7 +71,8 @@ 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_); + std::string milu("ILU"); + ilu_milu_ = convertString2Milu(param.getDefault("ilu_milu", milu)); // Check whether to use cpr approach const std::string cprSolver = "cpr"; @@ -91,7 +93,7 @@ namespace Opm linear_solver_use_amg_ = false; ilu_fillin_level_ = 0; ilu_relaxation_ = 0.9; - ilu_milu_ = false; + ilu_milu_ = MILU_VARIANT::ILU; } }; diff --git a/opm/autodiff/ParallelOverlappingILU0.hpp b/opm/autodiff/ParallelOverlappingILU0.hpp index 4cfd45223..2b5953a85 100644 --- a/opm/autodiff/ParallelOverlappingILU0.hpp +++ b/opm/autodiff/ParallelOverlappingILU0.hpp @@ -28,6 +28,7 @@ #include #include +#include namespace Opm { @@ -37,19 +38,49 @@ namespace Opm template class ParallelOverlappingILU0; +enum class MILU_VARIANT{ + /// \brief Do not perform modified ILU + ILU = 0, + /// \brief \f$U_{ii} = U_{ii} +\f$ sum(dropped entries) + MILU_1 = 1, + /// \brief \f$U_{ii} = U_{ii} + sign(U_{ii}) * \f$ sum(dropped entries) + MILU_2 = 2, + /// \brief \f$U_{ii} = U_{ii} sign(U_{ii}) * \f$ sum(|dropped entries|) + MILU_3 = 3, + /// \brief \f$U_{ii} = U_{ii} + (U_{ii}>0?1:0) * \f$ sum(dropped entries) + MILU_4 = 4 +}; + +inline MILU_VARIANT convertString2Milu(std::string milu) +{ + if( 0 == milu.compare("MILU_1") ) + { + return MILU_VARIANT::MILU_1; + } + if ( 0 == milu.compare("MILU_2") ) + { + return MILU_VARIANT::MILU_2; + } + if ( 0 == milu.compare("MILU_3") ) + { + return MILU_VARIANT::MILU_3; + } + return MILU_VARIANT::ILU; +} + template class ParallelOverlappingILU0Args : public Dune::Amg::DefaultSmootherArgs { public: - ParallelOverlappingILU0Args(bool milu = false) + ParallelOverlappingILU0Args(MILU_VARIANT milu = MILU_VARIANT::ILU ) : milu_(milu) {} - void setMilu(bool milu) + void setMilu(MILU_VARIANT milu) { milu_ = milu; } - bool getMilu() const + MILU_VARIANT getMilu() const { return milu_; } @@ -62,7 +93,7 @@ class ParallelOverlappingILU0Args return n_; } private: - bool milu_; + MILU_VARIANT milu_; int n_; }; } // end namespace Opm @@ -118,12 +149,68 @@ namespace Opm namespace detail { - template - void milu0_decomposition(M& A, + struct IdentityFunctor + { + template + T operator()(const T& t) + { + return t; + } + }; + + struct OneFunctor + { + template + T operator()(const T&) + { + return 1.0; + } + }; + struct SignFunctor + { + template + double operator()(const T& t) + { + if ( t < 0.0 ) + { + return -1; + } + else + { + return 1.0; + } + } + }; + + struct IsPositiveFunctor + { + template + double operator()(const T& t) + { + if ( t < 0.0 ) + { + return 0; + } + else + { + return 1; + } + } + }; + struct AbsFunctor + { + template + T operator()(const T& t) + { + using std::abs; + return abs(t); + } + }; + + template + void milu0_decomposition(M& A, F1 absFunctor = F1(), F2 signFunctor = F2(), std::vector* diagonal = nullptr) { - using block = typename M::block_type; - if( diagonal ) { diagonal->reserve(A.N()); @@ -134,8 +221,7 @@ namespace Opm auto a_i_end = irow->end(); auto a_ik = irow->begin(); - block sum_dropped; - sum_dropped = 0.0; + std::array sum_dropped{}; // Eliminate entries in lower triangular matrix // and store factors for L @@ -170,7 +256,15 @@ namespace Opm } else { - sum_dropped += modifier; + auto entry = sum_dropped.begin(); + for( const auto& row: modifier ) + { + for( const auto& colEntry: row ) + { + *entry += absFunctor(-colEntry); + } + ++entry; + } ++a_kj; } } @@ -180,12 +274,12 @@ namespace Opm OPM_THROW(std::logic_error, "Matrix is missing diagonal for row " << irow.index()); int index = 0; - for(const auto& row: sum_dropped) + for(const auto& entry: sum_dropped) { - for(const auto& val: row) - { - (*a_ik)[index][index]-=val; - } + auto& bdiag = (*a_ik)[index][index]; + if(entry<0) + std::cout << entry << std::endl; + bdiag += signFunctor(bdiag) * entry; ++index; } @@ -197,6 +291,13 @@ namespace Opm } } + template + void milu0_decomposition(M& A, + std::vector* diagonal) + { + milu0_decomposition(A, detail::IdentityFunctor(), detail::OneFunctor(), + diagonal); + } //! compute ILU decomposition of A. A is overwritten by its decomposition template @@ -271,6 +372,7 @@ namespace Opm } } // end namespace detail + /// \brief A two-step version of an overlapping Schwarz preconditioner using one step ILU0 as /// /// This preconditioner differs from a ParallelRestrictedOverlappingSchwarz with @@ -374,14 +476,12 @@ 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). + \param milu The modified ILU variant to use. 0 means traditional ILU. \see MILU_VARIANT. */ template ParallelOverlappingILU0 (const Dune::BCRSMatrix& A, const int n, const field_type w, - bool milu) + MILU_VARIANT milu) : lower_(), upper_(), inv_(), @@ -398,14 +498,12 @@ 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). + \param milu The modified ILU variant to use. 0 means traditional ILU. \see MILU_VARIANT. */ template ParallelOverlappingILU0 (const Dune::BCRSMatrix& A, const ParallelInfo& comm, const int n, const field_type w, - bool milu) + MILU_VARIANT milu) : lower_(), upper_(), inv_(), @@ -422,13 +520,11 @@ 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). + \param milu The modified ILU variant to use. 0 means traditional ILU. \see MILU_VARIANT. */ template ParallelOverlappingILU0 (const Dune::BCRSMatrix& A, - const field_type w, bool milu) + const field_type w, MILU_VARIANT milu) : ParallelOverlappingILU0( A, 0, w, milu ) { } @@ -439,14 +535,12 @@ 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). + \param milu The modified ILU variant to use. 0 means traditional ILU. \see MILU_VARIANT. */ template ParallelOverlappingILU0 (const Dune::BCRSMatrix& A, const ParallelInfo& comm, const field_type w, - bool milu) + MILU_VARIANT milu) : lower_(), upper_(), inv_(), @@ -549,7 +643,7 @@ public: } protected: - void init( const Matrix& A, const int iluIteration, bool milu ) + void init( const Matrix& A, const int iluIteration, MILU_VARIANT milu ) { // (For older DUNE versions the communicator might be // invalid if redistribution in AMG happened on the coarset level. @@ -578,18 +672,32 @@ protected: if( iluIteration == 0 ) { // create ILU-0 decomposition ILU.reset( new Matrix( A ) ); - if ( milu ) - { - detail::milu0_decomposition ( *ILU ); - }else + switch ( milu ) { + case MILU_VARIANT::MILU_1: + detail::milu0_decomposition ( *ILU); + break; + case MILU_VARIANT::MILU_2: + detail::milu0_decomposition ( *ILU, detail::IdentityFunctor(), + detail::SignFunctor() ); + break; + case MILU_VARIANT::MILU_3: + detail::milu0_decomposition ( *ILU, detail::AbsFunctor(), + detail::SignFunctor() ); + break; + case MILU_VARIANT::MILU_4: + detail::milu0_decomposition ( *ILU, detail::IdentityFunctor(), + detail::IsPositiveFunctor() ); + break; + default: bilu0_decomposition( *ILU ); + break; } } else { // create ILU-n decomposition ILU.reset( new Matrix( A.N(), A.M(), Matrix::row_wise) ); - if ( milu ) + if ( milu > MILU_VARIANT::ILU ) { OpmLog::warning("MILU_N_NOT_SUPPORTED", "MILU variant only supported for zero fill-in"); }