Resort to support the usual parameters (as for flow_legacy)

These are solver_approach=cpr cpr_use_amg=(true|false), etc.

opm/autodiff/BlackoilAmg.hpp

@@ -316,13 +316,13 @@ public:
      * @param args The arguments used for constructing the smoother.
      * @param c The crition used for the aggregation within AMG.
      */
-    OneStepAMGCoarseSolverPolicy(const SmootherArgs& args, const Criterion& c)
+    OneStepAMGCoarseSolverPolicy(const CPRParameter* param, const SmootherArgs& args, const Criterion& c)
-        : smootherArgs_(args), criterion_(c), use_amg_(true)
+        : param_(param), smootherArgs_(args), criterion_(c)
     {}
     /** @brief Copy constructor. */
     OneStepAMGCoarseSolverPolicy(const OneStepAMGCoarseSolverPolicy& other)
-        : coarseOperator_(other.coarseOperator_), smootherArgs_(other.smootherArgs_),
+        : param_(other.param_), coarseOperator_(other.coarseOperator_), smootherArgs_(other.smootherArgs_),
-          criterion_(other.criterion_), use_amg_(other.use_amg_)
+          criterion_(other.criterion_)
     {}
 private:
     /**
@@ -333,14 +333,14 @@ private:
      */
     struct AMGInverseOperator : public Dune::InverseOperator<X,X>
     {
-        AMGInverseOperator(const typename AMGType::Operator& op,
+        AMGInverseOperator(const CPRParameter* param,
+                           const typename AMGType::Operator& op,
                            const Criterion& crit,
                            const typename AMGType::SmootherArgs& args,
-                           const Communication& comm,
+                           const Communication& comm)
-                           bool use_amg)
+            : param_(param), amg_(), smoother_(), op_(op), comm_(comm), first_(true)
-            : amg_(), smoother_(), op_(op), comm_(comm), first_(true)
         {
-            if ( use_amg )
+            if ( param_->cpr_use_amg_ )
             {
                 amg_.reset(new AMGType(op, crit,args, comm));
             }
@@ -386,8 +386,24 @@ private:
             {
                 prec = smoother_.get();
             }
-            Dune::BiCGSTABSolver<X> solver(const_cast<typename AMGType::Operator&>(op_), *sp, *prec , 1e-2, 25, 0);
+            // Linear solver parameters
-            solver.apply(x,b,res);
+            const double tolerance = param_->cpr_solver_tol_;
+            const int maxit        = param_->cpr_max_ell_iter_;
+            const int verbosity    = ( param_->cpr_solver_verbose_ &&
+                                       comm_.communicator().rank()==0 ) ? 1 : 0;
+            if ( param_->cpr_use_bicgstab_ )
+            {
+                Dune::BiCGSTABSolver<X> solver(const_cast<typename AMGType::Operator&>(op_), *sp, *prec,
+                                               tolerance, maxit, verbosity);
+                solver.apply(x,b,res);
+            }
+            else
+            {
+                Dune::CGSolver<X> solver(const_cast<typename AMGType::Operator&>(op_), *sp, *prec,
+                                         tolerance, maxit, verbosity);
+                solver.apply(x,b,res);
+            }
+
 #if ! DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6)
             delete sp;
 #endif
@@ -410,6 +426,7 @@ private:
         {
         }
     private:
+        const CPRParameter* param_;
         X x_;
         std::unique_ptr<AMGType> amg_;
         std::unique_ptr<Smoother> smoother_;
@@ -434,11 +451,11 @@ public:
         coarseOperator_=transferPolicy.getCoarseLevelOperator();
         const LevelTransferPolicy& transfer =
             reinterpret_cast<const LevelTransferPolicy&>(transferPolicy);
-        AMGInverseOperator* inv = new AMGInverseOperator(*coarseOperator_,
+        AMGInverseOperator* inv = new AMGInverseOperator(param_,
+                                                         *coarseOperator_,
                                                          criterion_,
                                                          smootherArgs_,
-                                                         transfer.getCoarseLevelCommunication(),
+                                                         transfer.getCoarseLevelCommunication());
-                                                         use_amg_);
 
         return inv; //std::shared_ptr<InverseOperator<X,X> >(inv);
 
@@ -447,11 +464,12 @@ public:
 private:
     /** @brief The coarse level operator. */
     std::shared_ptr<Operator> coarseOperator_;
+    /** @brief The parameters for the CPR preconditioner. */
+    const CPRParameter* param_;
     /** @brief The arguments used to construct the smoother. */
     SmootherArgs smootherArgs_;
     /** @brief The coarsening criterion. */
     Criterion criterion_;
-    bool use_amg_;
 };
 
 template<class Smoother, class Operator, class Communication>
@@ -856,14 +874,16 @@ public:
     };
 #endif
 
-    BlackoilAmg(const Operator& fineOperator, const Criterion& criterion,
+    BlackoilAmg(const CPRParameter& param,
+                const Operator& fineOperator, const Criterion& criterion,
                 const SmootherArgs& smargs, const Communication& comm)
-        : scaledMatrixOperator_(Detail::scaleMatrixQuasiImpes(fineOperator, comm,
+        : param_(param),
+          scaledMatrixOperator_(Detail::scaleMatrixQuasiImpes(fineOperator, comm,
                                                               COMPONENT_INDEX)),
           smoother_(Detail::constructSmoother<Smoother>(std::get<1>(scaledMatrixOperator_),
                                                         smargs, comm)),
           levelTransferPolicy_(criterion, comm),
-          coarseSolverPolicy_(smargs, criterion),
+          coarseSolverPolicy_(&param, smargs, criterion),
           twoLevelMethod_(std::get<1>(scaledMatrixOperator_), smoother_,
                           levelTransferPolicy_,
                           coarseSolverPolicy_, 0, 1)
@@ -888,6 +908,7 @@ public:
         twoLevelMethod_.apply(v, scaledD);
     }
 private:
+    const CPRParameter& param_;
     std::tuple<std::unique_ptr<Matrix>, Operator> scaledMatrixOperator_;
     std::shared_ptr<Smoother> smoother_;
     LevelTransferPolicy levelTransferPolicy_;

opm/autodiff/CPRPreconditioner.hpp

@@ -52,6 +52,12 @@
 namespace Opm
 {
 
+struct CPRParameter;
+
+template<typename O, typename S, typename C,
+         typename P, std::size_t COMPONENT_INDEX>
+class BlackoilAmg;
+
 namespace ISTLUtility
 {
 ///
@@ -149,9 +155,36 @@ createEllipticPreconditionerPointer(const M& Ae, double relax,
     return EllipticPreconditionerPointer(new ParallelPreconditioner(Ae, comm, relax));
 }
 
+template < class C, class Op, class P, class S, std::size_t index >
+inline void
+createAMGPreconditionerPointer(Op& opA, const double relax, const P& comm,
+                               std::unique_ptr< BlackoilAmg<Op,S,C,P,index> >& amgPtr,
+                               const CPRParameter& params = CPRParameter())
+{
+    using AMG = BlackoilAmg<Op,S,C,P,index>;
+    // TODO: revise choice of parameters
+    int coarsenTarget=1200;
+    using Criterion = C;
+    Criterion criterion(15, coarsenTarget);
+    criterion.setDebugLevel( 0 ); // no debug information, 1 for printing hierarchy information
+    criterion.setDefaultValuesIsotropic(2);
+    criterion.setNoPostSmoothSteps( 1 );
+    criterion.setNoPreSmoothSteps( 1 );
+
+    // for DUNE 2.2 we also need to pass the smoother args
+    typedef typename AMG::Smoother Smoother;
+    typedef typename Dune::Amg::SmootherTraits<Smoother>::Arguments  SmootherArgs;
+    SmootherArgs  smootherArgs;
+    smootherArgs.iterations = 1;
+    smootherArgs.relaxationFactor = relax;
+
+    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 P& comm, std::unique_ptr< AMG >& amgPtr,
+                               const CPRParameter& params = CPRParameter())
 {
     // TODO: revise choice of parameters
     int coarsenTarget=1200;

opm/autodiff/FlowMainEbos.hpp

@@ -610,7 +610,14 @@ namespace Opm
         void setupLinearSolver()
         {
             typedef typename BlackoilModelEbos<TypeTag> :: ISTLSolverType ISTLSolverType;
-
+            const std::string cprSolver = "cpr";
+            if (!param_.has("solver_approach") )
+            {
+                if ( eclState().getSimulationConfig().useCPR() )
+                {
+                    param_.insertParameter("solver_approach", cprSolver);
+                }
+            }
             extractParallelGridInformationToISTL(grid(), parallel_information_);
             fis_solver_.reset( new ISTLSolverType( param_, parallel_information_ ) );
         }

opm/autodiff/ISTLSolver.hpp

@@ -404,7 +404,7 @@ namespace Opm
             parallelInformation_arg.copyOwnerToAll(istlb, istlb);
 
 #if FLOW_SUPPORT_AMG // activate AMG if either flow_ebos is used or UMFPack is not available
-            if( parameters_.linear_solver_use_amg_ )
+            if( parameters_.linear_solver_use_amg_ || parameters_.use_cpr_)
             {
                 typedef ISTLUtility::CPRSelector< Matrix, Vector, Vector, POrComm>  CPRSelectorType;
                 typedef typename CPRSelectorType::Operator MatrixOperator;
@@ -418,7 +418,7 @@ namespace Opm
                 }
 
                 const double relax = parameters_.ilu_relaxation_;
-                if ( ! parameters_.amg_blackoil_system_ )
+                if ( ! parameters_.use_cpr_ )
                 {
                     typedef typename CPRSelectorType::AMG AMG;
                     std::unique_ptr< AMG > amg;
@@ -519,7 +519,7 @@ namespace Opm
         void
         constructAMGPrecond(LinearOperator& /* linearOperator */, const POrComm& comm, std::unique_ptr< AMG >& amg, std::unique_ptr< MatrixOperator >& opA, const double relax ) const
         {
-            ISTLUtility::template createAMGPreconditionerPointer<C>( *opA, relax, comm, amg );
+            ISTLUtility::template createAMGPreconditionerPointer<C>( *opA, relax, comm, amg, parameters_ );
         }
 
 
@@ -527,7 +527,7 @@ namespace Opm
         void
         constructAMGPrecond(MatrixOperator& opA, const POrComm& comm, std::unique_ptr< AMG >& amg, std::unique_ptr< MatrixOperator >&, const double relax ) const
         {
-            ISTLUtility::template createAMGPreconditionerPointer<C>( opA, relax, comm, amg );
+            ISTLUtility::template createAMGPreconditionerPointer<C>( opA, relax, comm, amg, parameters_ );
         }
         /// \brief Solve the system using the given preconditioner and scalar product.
         template <class Operator, class ScalarProd, class Precond>

opm/autodiff/NewtonIterationBlackoilInterleaved.hpp

@@ -25,6 +25,7 @@
 #ifndef OPM_NEWTONITERATIONBLACKOILINTERLEAVED_HEADER_INCLUDED
 #define OPM_NEWTONITERATIONBLACKOILINTERLEAVED_HEADER_INCLUDED
 
+#include <opm/autodiff/CPRPreconditioner.hpp>
 #include <opm/autodiff/NewtonIterationBlackoilInterface.hpp>
 #include <opm/common/utility/parameters/ParameterGroup.hpp>
 
@@ -35,6 +36,7 @@ namespace Opm
 {
     /// This class carries all parameters for the NewtonIterationBlackoilInterleaved class
     struct NewtonIterationBlackoilInterleavedParameters
+        : public CPRParameter
     {
         double linear_solver_reduction_;
         double ilu_relaxation_;
@@ -46,11 +48,12 @@ namespace Opm
         bool   require_full_sparsity_pattern_;
         bool   ignoreConvergenceFailure_;
         bool   linear_solver_use_amg_;
-        bool   amg_blackoil_system_;
+        bool   use_cpr_;
 
         NewtonIterationBlackoilInterleavedParameters() { reset(); }
         // read values from parameter class
         NewtonIterationBlackoilInterleavedParameters( const ParameterGroup& param )
+            : CPRParameter(param)
         {
             // set default parameters
             reset();
@@ -64,15 +67,18 @@ namespace Opm
             require_full_sparsity_pattern_ = param.getDefault("require_full_sparsity_pattern", require_full_sparsity_pattern_);
             ignoreConvergenceFailure_ = param.getDefault("linear_solver_ignoreconvergencefailure", ignoreConvergenceFailure_);
             linear_solver_use_amg_    = param.getDefault("linear_solver_use_amg", linear_solver_use_amg_ );
-            amg_blackoil_system_      = param.getDefault("amg_blackoil_system", amg_blackoil_system_);
             ilu_relaxation_           = param.getDefault("ilu_relaxation", ilu_relaxation_ );
             ilu_fillin_level_         = param.getDefault("ilu_fillin_level",  ilu_fillin_level_ );
+
+            // Check whether to use cpr approach
+            const std::string cprSolver = "cpr";
+            use_cpr_ = ( param.getDefault("solver_approach", std::string()) == cprSolver );
         }
 
         // set default values
         void reset()
         {
-            amg_blackoil_system_     = false;
+            use_cpr_     = false;
             newton_use_gmres_        = false;
             linear_solver_reduction_ = 1e-2;
             linear_solver_maxiter_   = 150;