Merge pull request #252 from dr-robertk/PR/cprpreconditioner-more-flexible-choice-of-ilu

This commits allows for flexible choice of either ILU(0) or ILU(n).

opm/autodiff/CPRPreconditioner.hpp

@@ -78,11 +78,14 @@ namespace Opm
         //! \brief Elliptic Operator
         typedef Dune::MatrixAdapter<M,X,X> Operator;
 
+        //! \brief preconditioner for the whole system (here either ILU(0) or ILU(n)
+        typedef Dune::Preconditioner<X,X> WholeSystemPreconditioner;
+
         //! \brief ilu-0 preconditioner for the elliptic system
-        typedef Dune::SeqILU0<M,X,X> Preconditioner;
+        typedef Dune::SeqILU0<M,X,X> EllipticPreconditioner;
 
         //! \brief amg preconditioner for the elliptic system
-        typedef Preconditioner Smoother;
+        typedef EllipticPreconditioner Smoother;
         typedef Dune::Amg::AMG<Operator, X, Smoother> AMG;
 
         /*! \brief Constructor.
@@ -95,8 +98,9 @@ namespace Opm
           \param useBiCG if true, BiCG solver is used (default), otherwise CG solver
         */
         CPRPreconditioner (const M& A, const M& Ae, const field_type relax,
-                           const bool useAMG  = false,
+                           const unsigned int ilu_n,
-                           const bool useBiCG = true )
+                           const bool useAMG,
+                           const bool useBiCG )
             : A_(A),
               Ae_(Ae),
               de_( Ae_.N() ),
@@ -105,15 +109,20 @@ namespace Opm
               opAe_( Ae_ ),
               precond_(), // ilu0 preconditioner for elliptic system
               amg_(),     // amg  preconditioner for elliptic system
-              ILU_(A),    // copy A (will be overwritten by ILU decomp)
+              pre_(), // copy A will be made be the preconditioner
-              vilu_( ILU_.N() ),
+              vilu_( A_.N() ),
               relax_(relax),
               use_bicg_solver_( useBiCG )
         {
             // create appropriate preconditioner for elliptic system
             createPreconditioner( useAMG );
 
-            Dune::bilu0_decomposition(ILU_);
+            if( ilu_n == 0 ) {
+                pre_.reset( new Dune::SeqILU0<M,X,X>( A_, relax_) );
+            }
+            else {
+                pre_.reset( new Dune::SeqILUn<M,X,X>( A_, ilu_n, relax_) );
+            }
         }
 
         /*!
@@ -151,14 +160,17 @@ namespace Opm
             dmodified_ = d;
             A_.mmv(v, dmodified_);
 
-            // Apply ILU0.
+            // Apply Preconditioner for whole system (relax will be applied already)
-            Dune::bilu_backsolve(ILU_, vilu_, dmodified_);
+            pre_->apply( vilu_, dmodified_);
-            v += vilu_;
 
             // don't apply relaxation if relax_ == 1
-            if( std::abs( relax_ - 1.0 ) < 1e-12 ) return;
+            if( std::abs( relax_ - 1.0 ) < 1e-12 ) {
-
+                v += vilu_;
+            }
+            else {
                 v *= relax_;
+                v += vilu_;
+            }
         }
 
         /*!
@@ -227,12 +239,12 @@ namespace Opm
         Operator opAe_;
 
         //! \brief ILU0 preconditioner for the elliptic system
-        std::unique_ptr< Preconditioner > precond_;
+        std::unique_ptr< EllipticPreconditioner > precond_;
         //! \brief AMG preconditioner with ILU0 smoother
         std::unique_ptr< AMG > amg_;
 
-        //! \brief The ILU0 decomposition of the matrix.
+        //! \brief The preconditioner for the whole system
-        matrix_type ILU_;
+        std::unique_ptr< WholeSystemPreconditioner > pre_;
 
         //! \brief temporary variables for ILU solve
         Y vilu_;
@@ -266,7 +278,7 @@ namespace Opm
               amg_ = std::unique_ptr< AMG > (new AMG(opAe_, criterion, smootherArgs));
             }
             else
-              precond_ = std::unique_ptr< Preconditioner > (new Preconditioner( Ae_, relax_ ));
+              precond_ = std::unique_ptr< EllipticPreconditioner > (new EllipticPreconditioner( Ae_, relax_ ));
        }
     };
 

opm/autodiff/NewtonIterationBlackoilCPR.cpp

@@ -117,8 +117,10 @@ namespace Opm
     NewtonIterationBlackoilCPR::NewtonIterationBlackoilCPR(const parameter::ParameterGroup& param)
         : iterations_( 0 )
     {
-        use_amg_ = param.getDefault("cpr_use_amg", false);
+        cpr_relax_        = param.getDefault("cpr_relax", 1.0);
-        use_bicgstab_ = param.getDefault("cpr_use_bicgstab", true);
+        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);
     }
 
 
@@ -193,8 +195,7 @@ namespace Opm
         // Construct preconditioner.
         // typedef Dune::SeqILU0<Mat,Vector,Vector> Preconditioner;
         typedef Opm::CPRPreconditioner<Mat,Vector,Vector> Preconditioner;
-        const double relax = 1.0;
+        Preconditioner precond(istlA, istlAe, cpr_relax_, cpr_ilu_n_, cpr_use_amg_, cpr_use_bicgstab_);
-        Preconditioner precond(istlA, istlAe, relax, use_amg_, use_bicgstab_);
 
         // Construct linear solver.
         const double tolerance = 1e-3;

opm/autodiff/NewtonIterationBlackoilCPR.hpp

@@ -43,6 +43,8 @@ namespace Opm
         ///                    the preconditioning and choice of
         ///                    linear solvers.
         ///                    Parameters:
+        ///                        cpr_relax        (default 1.0) relaxation for the preconditioner
+        ///                        cpr_ilu_n        (default 0) use ILU(n) for preconditioning of the linear system
         ///                        cpr_use_amg      (default false) if true, use AMG preconditioner for elliptic part
         ///                        cpr_use_bicgstab (default true)  if true, use BiCGStab (else use CG) for elliptic part
         NewtonIterationBlackoilCPR(const parameter::ParameterGroup& param);
@@ -58,8 +60,10 @@ namespace Opm
         virtual int iterations () const { return iterations_; }
     private:
         mutable int iterations_;
-        bool use_amg_;
+        double cpr_relax_;
-        bool use_bicgstab_;
+        unsigned int cpr_ilu_n_;
+        bool cpr_use_amg_;
+        bool cpr_use_bicgstab_;
     };
 
 } // namespace Opm