added writing of coarse level matrix in cpr and corrected coarsening criterio in flexible solver

opm/simulators/linalg/BlackoilAmg.hpp

@@ -35,6 +35,7 @@
 #include <dune/istl/scalarproducts.hh>
 #include <dune/common/fvector.hh>
 #include <dune/common/fmatrix.hh>
+#include "matrixmarket_ewoms.hh"
 namespace Dune
 {
 namespace Amg
@@ -553,6 +554,13 @@ public:
     CoarseLevelSolver* createCoarseLevelSolver(LTP& transferPolicy)
     {
         coarseOperator_=transferPolicy.getCoarseLevelOperator();
+	if( param_->cpr_solver_verbose_ >10 ){
+	    std::string filename = "course_matrix_istl_blackoil.txt";
+	    std::ofstream filem(filename);
+	    //const auto& matrix = *transferPolicy.getCoarseMatrix();
+	    const auto& matrix = coarseOperator_->getmat(); 
+	    Dune::writeMatrixMarket(matrix, filem);
+	}
         const LevelTransferPolicy& transfer =
             reinterpret_cast<const LevelTransferPolicy&>(transferPolicy);
         AMGInverseOperator* inv = new AMGInverseOperator(param_,

opm/simulators/linalg/FlexibleSolver.hpp

@@ -130,6 +130,14 @@ private:
                                                               tol, // desired residual reduction factor
                                                               maxiter, // maximum number of iterations
                                                               verbosity));
+	} else if (solver_type == "cg") {
+            linsolver_.reset(new Dune::CGSolver<VectorType>(*linearoperator_,
+							    *scalarproduct_,
+							    *preconditioner_,
+							    tol, // desired residual reduction factor
+							    maxiter, // maximum number of iterations
+							    verbosity));    
+	    
         } else if (solver_type == "gmres") {
             int restart = prm.get<int>("restart");
             linsolver_.reset(new Dune::RestartedGMResSolver<VectorType>(*linearoperator_,

opm/simulators/linalg/PreconditionerFactory.hpp

@@ -29,6 +29,7 @@
 #include <opm/simulators/linalg/amgcpr.hh>
 
 #include <dune/istl/paamg/amg.hh>
+#include <dune/istl/paamg/kamg.hh>
 #include <dune/istl/paamg/fastamg.hh>
 #include <dune/istl/preconditioners.hh>
 
@@ -104,8 +105,10 @@ public:
     }
 
 private:
+    //using CriterionBase
+    //    = Dune::Amg::AggregationCriterion<Dune::Amg::SymmetricMatrixDependency<Matrix, Dune::Amg::FirstDiagonal>>;
     using CriterionBase
-        = Dune::Amg::AggregationCriterion<Dune::Amg::SymmetricMatrixDependency<Matrix, Dune::Amg::FirstDiagonal>>;
+        = Dune::Amg::AggregationCriterion<Dune::Amg::SymmetricDependency<Matrix, Dune::Amg::FirstDiagonal>>;
     using Criterion = Dune::Amg::CoarsenCriterion<CriterionBase>;
 
     // Helpers for creation of AMG preconditioner.
@@ -116,7 +119,9 @@ private:
         criterion.setAlpha(prm.get<double>("alpha"));
         criterion.setBeta(prm.get<double>("beta"));
         criterion.setMaxLevel(prm.get<int>("maxlevel"));
-        criterion.setSkipIsolated(false);
+        criterion.setSkipIsolated(prm.get<bool>("skip_isolated"));
+	criterion.setNoPreSmoothSteps(prm.get<int>("pre_smooth"));
+	criterion.setNoPostSmoothSteps(prm.get<int>("post_smooth"));
         criterion.setDebugLevel(prm.get<int>("verbosity"));
         return criterion;
     }
@@ -135,11 +140,21 @@ private:
     }
 
     template <class Smoother>
-    static PrecPtr makeAmgPreconditioner(const Operator& op, const boost::property_tree::ptree& prm)
+    static PrecPtr makeAmgPreconditioner(const Operator& op, const boost::property_tree::ptree& prm,bool use_kamg = false)
     {
         auto crit = amgCriterion(prm);
         auto sargs = amgSmootherArgs<Smoother>(prm);
-        return std::make_shared<Dune::Amg::AMGCPR<Operator, Vector, Smoother>>(op, crit, sargs);
+	if(use_kamg){
+	    return std::make_shared<
+		DummyUpdatePreconditioner<
+		    Dune::Amg::KAMG< Operator, Vector, Smoother>
+		    >
+		>(op, crit, sargs,
+		  prm.get<size_t>("max_krylov"),
+		  prm.get<double>("min_reduction")  );
+	}else{
+	    return std::make_shared<Dune::Amg::AMGCPR<Operator, Vector, Smoother>>(op, crit, sargs);
+	}
     }
 
     // Add a useful default set of preconditioners to the factory.
@@ -191,7 +206,7 @@ private:
         });
         doAddCreator("amg", [](const O& op, const P& prm, const C& comm) {
             const std::string smoother = prm.get<std::string>("smoother");
-            if (smoother == "ILU0") {
+            if (smoother == "ParOverILU0") {
                 using Smoother = Opm::ParallelOverlappingILU0<M, V, V, C>;
                 auto crit = amgCriterion(prm);
                 auto sargs = amgSmootherArgs<Smoother>(prm);
@@ -275,6 +290,32 @@ private:
                 throw std::runtime_error(msg);
             }
         });
+	doAddCreator("kamg", [](const O& op, const P& prm) {
+            const std::string smoother = prm.get<std::string>("smoother");
+            if (smoother == "ILU0") {
+                using Smoother = SeqILU0<M, V, V>;
+                return makeAmgPreconditioner<Smoother>(op, prm,true);
+            } else if (smoother == "Jac") {
+                using Smoother = SeqJac<M, V, V>;
+                return makeAmgPreconditioner<Smoother>(op, prm,true);
+            } else if (smoother == "SOR") {
+                using Smoother = SeqSOR<M, V, V>;
+                return makeAmgPreconditioner<Smoother>(op, prm,true);
+	    // } else if (smoother == "GS") {
+            //     using Smoother = SeqGS<M, V, V>;
+            //     return makeAmgPreconditioner<Smoother>(op, prm,true);
+            } else if (smoother == "SSOR") {
+                using Smoother = SeqSSOR<M, V, V>;
+                return makeAmgPreconditioner<Smoother>(op, prm,true);
+            } else if (smoother == "ILUn") {
+                using Smoother = SeqILUn<M, V, V>;
+                return makeAmgPreconditioner<Smoother>(op, prm,true);
+            } else {
+                std::string msg("No such smoother: ");
+                msg += smoother;
+                throw std::runtime_error(msg);
+            }
+        });
         doAddCreator("famg", [](const O& op, const P& prm) {
             auto crit = amgCriterion(prm);
             Dune::Amg::Parameters parms;

opm/simulators/linalg/PressureSolverPolicy.hpp

@@ -113,6 +113,16 @@ namespace Amg
             auto& tp = dynamic_cast<LevelTransferPolicy&>(transferPolicy); // TODO: make this unnecessary.
             PressureInverseOperator* inv
                 = new PressureInverseOperator(*coarseOperator_, prm_, tp.getCoarseLevelCommunication());
+	    if (prm_.get<int>("verbosity") > 10) {         
+		std::string filename = "course_matrix_istl_pressure_flex.txt";
+		std::ofstream filem(filename);
+		if (!filem) {
+		    throw std::runtime_error("Could not write matrix");
+		}
+		//const auto& matrix = *levelTransferPolicy_.getCoarseMatrix();
+		const auto& matrix = coarseOperator_->getmat(); 
+		Dune::writeMatrixMarket(matrix, filem);
+	    }
             return inv;
         }