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

@@ -29,6 +29,7 @@
 
 #include <opm/common/utility/platform_dependent/disable_warnings.h>
 #include <opm/common/utility/parameters/ParameterGroup.hpp>
+#include <opm/autodiff/ParallelOverlappingILU0.hpp>
 
 #include <dune/istl/bvector.hh>
 #include <dune/istl/bcrsmatrix.hh>
@@ -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<T>& args,
 
 template<class T>
 void setILUParameters(Opm::ParallelOverlappingILU0Args<T>& args,
-                      bool milu, int n=0)
+                      MILU_VARIANT milu, int n=0)
 {
     args.setN(n);
     args.setMilu(milu);
@@ -195,7 +198,7 @@ struct CPRSelector<M,X,Y,Dune::Amg::SequentialInformation>
 //! \param relax The relaxation factor to use.
 template<class M, class X, class C>
 std::shared_ptr<ParallelOverlappingILU0<M,X,X,C> >
-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<ParallelOverlappingILU0<M,X,X,C> >(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<class M, class X, class C>
 std::shared_ptr<ParallelOverlappingILU0<M,X,X,C> >
-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<ParallelOverlappingILU0<M,X,X,C> >( 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<class M, class X=typename M::range_type, class P>
 typename CPRSelector<M,X,X,P>::EllipticPreconditionerPointer
 createEllipticPreconditionerPointer(const M& Ae, double relax,
-                                    bool milu, const P& comm)
+                                    MILU_VARIANT milu, const P& comm)
 {
     typedef typename CPRSelector<M,X,X,P >
         ::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;

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<SeqPreconditioner> precond(new SeqPreconditioner(opA.getmat(), ilu_fillin, relax, ilu_milu));
             return precond;
         }
@@ -519,14 +519,14 @@ namespace Detail
         {
             typedef std::unique_ptr<ParPreconditioner> 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 <class LinearOperator, class MatrixOperator, class POrComm, class AMG >
         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<pressureIndex>( *opA, relax, milu, comm, amg );
         }
@@ -535,7 +535,7 @@ namespace Detail
         template <class MatrixOperator, class POrComm, class AMG >
         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<pressureIndex>( opA, relax,
                                                                                  milu, comm, amg );
@@ -544,7 +544,7 @@ namespace Detail
         template <class C, class LinearOperator, class MatrixOperator, class POrComm, class AMG >
         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<C>( *opA, relax,
                                                                      comm, amg, parameters_ );
@@ -553,7 +553,7 @@ namespace Detail
 
         template <class C, class MatrixOperator, class POrComm, class AMG >
         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<C>( opA, relax, milu,
                                                                      comm, amg, parameters_ );

opm/autodiff/NewtonIterationBlackoilInterleaved.hpp

@@ -28,6 +28,7 @@
 #include <opm/autodiff/CPRPreconditioner.hpp>
 #include <opm/autodiff/NewtonIterationBlackoilInterface.hpp>
 #include <opm/common/utility/parameters/ParameterGroup.hpp>
+#include <opm/autodiff/ParallelOverlappingILU0.hpp>
 
 #include <array>
 #include <memory>
@@ -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;
         }
     };
 

opm/autodiff/ParallelOverlappingILU0.hpp

@@ -28,6 +28,7 @@
 #include <dune/istl/paamg/pinfo.hh>
 
 #include <type_traits>
+#include <numeric>
 
 namespace Opm
 {
@@ -37,19 +38,49 @@ namespace Opm
 template<class Matrix, class Domain, class Range, class ParallelInfo = Dune::Amg::SequentialInformation>
 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 F>
 class ParallelOverlappingILU0Args
     : public Dune::Amg::DefaultSmootherArgs<F>
 {
  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<class M>
+    struct IdentityFunctor
-    void milu0_decomposition(M& A,
+    {
+        template<class T>
+        T operator()(const T& t)
+        {
+            return t;
+        }
+    };
+
+    struct OneFunctor
+    {
+        template<class T>
+        T operator()(const T&)
+        {
+            return 1.0;
+        }
+    };
+    struct SignFunctor
+    {
+        template<class T>
+        double operator()(const T& t)
+        {
+            if ( t < 0.0 )
+            {
+                return -1;
+            }
+            else
+            {
+                return 1.0;
+            }
+        }
+    };
+
+    struct IsPositiveFunctor
+    {
+        template<class T>
+        double operator()(const T& t)
+        {
+            if ( t < 0.0 )
+            {
+                return 0;
+            }
+            else
+            {
+                return 1;
+            }
+        }
+    };
+    struct AbsFunctor
+    {
+        template<class T>
+        T operator()(const T& t)
+        {
+            using std::abs;
+            return abs(t);
+        }
+    };
+
+    template<class M, class F1=detail::IdentityFunctor, class F2=detail::OneFunctor >
+    void milu0_decomposition(M& A, F1 absFunctor = F1(), F2 signFunctor = F2(),
                              std::vector<typename M::block_type>* 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;
+            std::array<typename M::field_type, M::block_type::rows> sum_dropped{};
-            sum_dropped = 0.0;
 
             // 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)
+                auto& bdiag = (*a_ik)[index][index];
-                {
+                if(entry<0)
-                    (*a_ik)[index][index]-=val;
+                    std::cout << entry << std::endl;
-                }
+                bdiag += signFunctor(bdiag) * entry;
                 ++index;
             }
 
@@ -197,6 +291,13 @@ namespace Opm
         }
     }
 
+    template<class M>
+    void milu0_decomposition(M& A,
+                             std::vector<typename M::block_type>* diagonal)
+    {
+        milu0_decomposition(A, detail::IdentityFunctor(), detail::OneFunctor(),
+                            diagonal);
+    }
 
       //! compute ILU decomposition of A. A is overwritten by its decomposition
       template<class M, class CRS, class InvVector>
@@ -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
+      \param milu The modified ILU variant to use. 0 means traditional ILU. \see MILU_VARIANT.
-                  will get added to the diagonal of U to preserve the row sum
-                  for constant vectors (Ae = LUe).
     */
     template<class BlockType, class Alloc>
     ParallelOverlappingILU0 (const Dune::BCRSMatrix<BlockType,Alloc>& 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
+      \param milu The modified ILU variant to use. 0 means traditional ILU. \see MILU_VARIANT.
-                  will get added to the diagonal of U to preserve the row sum
-                  for constant vectors (Ae = LUe).
     */
     template<class BlockType, class Alloc>
     ParallelOverlappingILU0 (const Dune::BCRSMatrix<BlockType,Alloc>& 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
+      \param milu The modified ILU variant to use. 0 means traditional ILU. \see MILU_VARIANT.
-                  will get added to the diagonal of U to preserve the row sum
-                  for constant vectors (Ae = LUe).
     */
     template<class BlockType, class Alloc>
     ParallelOverlappingILU0 (const Dune::BCRSMatrix<BlockType,Alloc>& 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
+      \param milu   The modified ILU variant to use. 0 means traditional ILU. \see MILU_VARIANT.
-                  will get added to the diagonal of U to preserve the row sum
-                  for constant vectors (Ae = LUe).
     */
     template<class BlockType, class Alloc>
     ParallelOverlappingILU0 (const Dune::BCRSMatrix<BlockType,Alloc>& 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 )
+                switch ( milu )
-                {
-                    detail::milu0_decomposition ( *ILU );
-                }else
                 {
+                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");
                 }