changes to make true impes work with flexible solver

Cherry-picked commit 7da5ce4fec33. Which got cleaned up
(deleted trailing whitespace, tabs) and made compilable with recent changes.

opm/simulators/linalg/ISTLSolverEbosFlexible.hpp

@@ -65,6 +65,19 @@ class ISTLSolverEbosFlexible
 #endif
     using SolverType = Dune::FlexibleSolver<MatrixType, VectorType>;
 
+    // for quasiImpesWeights
+    typedef typename GET_PROP_TYPE(TypeTag, GlobalEqVector) Vector;
+    typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+    //typedef typename GET_PROP_TYPE(TypeTag, EclWellModel) WellModel;
+    //typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
+    typedef typename SparseMatrixAdapter::IstlMatrix Matrix;
+    typedef typename SparseMatrixAdapter::MatrixBlock MatrixBlockType;
+    typedef typename Vector::block_type BlockVector;
+    typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
+    typedef typename GET_PROP_TYPE(TypeTag, ThreadManager) ThreadManager;
+    typedef typename GridView::template Codim<0>::Entity Element;
+    typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
+
 
 public:
     static void registerParameters()
@@ -135,6 +148,30 @@ public:
             assert(recreate_solver == false);
             // Never recreate solver.
         }
+        if( prm_.get<std::string>("preconditioner.type") == "cpr" ||
+            prm_.get<std::string>("preconditioner.type") == "cprt"
+            )
+        {
+            bool transpose = false;
+            if(prm_.get<std::string>("preconditioner.type") == "cprt"){
+                transpose = true;
+            }
+
+
+            if(prm_.get<std::string>("preconditioner.weight_type") == "quasiimpes") {
+                VectorType weights = Opm::Amg::getQuasiImpesWeights<MatrixType, VectorType>(
+                    mat.istlMatrix(),
+                    prm_.get<int>("preconditioner.pressure_var_index"), transpose);
+                prm_.put("preconditioner.weights",weights);
+
+            }else if(prm_.get<std::string>("preconditioner.weight_type") == "trueimpes"  ){
+                VectorType weights =
+                    this->getTrueImpesWeights(b, prm_.get<int>("preconditioner.pressure_var_index"));
+            prm_.put("preconditioner.weights",weights);
+            }else{
+                throw std::runtime_error("no such weights implemented for cpr");
+            }
+        }
 
         if (recreate_solver || !solver_) {
             if (isParallel()) {
@@ -146,7 +183,7 @@ public:
             }
             rhs_ = b;
         } else {
-            solver_->preconditioner().update();
+            solver_->preconditioner().update(prm_);
             rhs_ = b;
         }
     }
@@ -204,6 +241,47 @@ protected:
         }
     }
 
+    VectorType getTrueImpesWeights(const VectorType& b,const int pressureVarIndex)
+    {
+        VectorType weights(b.size());
+        BlockVector rhs(0.0);
+        rhs[pressureVarIndex] = 1.0;
+        int index = 0;
+        ElementContext elemCtx(simulator_);
+        const auto& vanguard = simulator_.vanguard();
+        auto elemIt = vanguard.gridView().template begin</*codim=*/0>();
+        const auto& elemEndIt = vanguard.gridView().template end</*codim=*/0>();
+        for (; elemIt != elemEndIt; ++elemIt) {
+            const Element& elem = *elemIt;
+            elemCtx.updatePrimaryStencil(elem);
+            elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
+            const int numEq = MatrixType::block_type::rows;
+            Dune::FieldVector<Evaluation, numEq> storage;
+            unsigned threadId = ThreadManager::threadId();
+            simulator_.model().localLinearizer(threadId).localResidual().computeStorage(storage,elemCtx,/*spaceIdx=*/0, /*timeIdx=*/0);
+            Scalar extrusionFactor = elemCtx.intensiveQuantities(0, /*timeIdx=*/0).extrusionFactor();
+            Scalar scvVolume = elemCtx.stencil(/*timeIdx=*/0).subControlVolume(0).volume() * extrusionFactor;
+            Scalar storage_scale = scvVolume / elemCtx.simulator().timeStepSize();
+            MatrixBlockType block;
+            double pressure_scale = 50e5;
+            for (int ii = 0; ii < numEq; ++ii) {
+                for (int jj = 0; jj < numEq; ++jj) {
+                    block[ii][jj] = storage[ii].derivative(jj)/storage_scale;
+                    if (jj == pressureVarIndex) {
+                        block[ii][jj] *= pressure_scale;
+                    }
+                }
+            }
+            BlockVector bweights;
+            MatrixBlockType block_transpose = Opm::transposeDenseMatrix(block);
+            block_transpose.solve(bweights, rhs);
+            bweights /= 1000.0; // given normal densities this scales weights to about 1.
+            weights[index] = bweights;
+            ++index;
+        }
+        return weights;
+    }
+
 
     const Simulator& simulator_;
 

opm/simulators/linalg/OwningBlockPreconditioner.hpp

@@ -63,9 +63,9 @@ public:
     }
 
     // The update() function does nothing for a wrapped preconditioner.
-    virtual void update() override
+    virtual void update(const boost::property_tree::ptree& p) override
     {
-        orig_precond_.update();
+        orig_precond_.update(p);
     }
 
 private:

opm/simulators/linalg/OwningTwoLevelPreconditioner.hpp

@@ -53,6 +53,21 @@ namespace Dune
 template <class MatrixTypeT, class VectorTypeT>
 class FlexibleSolver;
 
+template <typename T, typename A, int i>
+std::ostream& operator<<(std::ostream& out,
+                         const BlockVector< FieldVector< T, i >, A >&  vector)
+{
+   Dune::writeMatrixMarket(vector, out);
+   return out;
+}
+
+    template <typename T, typename A, int i>
+    std::istream& operator>>(std::istream& input,
+                             BlockVector< FieldVector< T, i >, A >&  vector)
+{
+   Dune::readMatrixMarket(vector, input);
+   return input;
+}
 
 /// A version of the two-level preconditioner that is:
 /// - Self-contained, because it owns its policy components.
@@ -73,8 +88,7 @@ public:
         : linear_operator_(linearoperator)
         , finesmoother_(PrecFactory::create(linearoperator, prm.get_child("finesmoother")))
         , comm_(nullptr)
-        , weights_(Opm::Amg::getQuasiImpesWeights<MatrixType, VectorType>(
+        , weights_(prm.get<VectorType>("weights"))
-              linearoperator.getmat(), prm.get<int>("pressure_var_index"), transpose))
         , levelTransferPolicy_(dummy_comm_, weights_, prm.get<int>("pressure_var_index"))
         , coarseSolverPolicy_(prm.get_child("coarsesolver"))
         , twolevel_method_(linearoperator,
@@ -98,8 +112,7 @@ public:
         : linear_operator_(linearoperator)
         , finesmoother_(PrecFactory::create(linearoperator, prm.get_child("finesmoother"), comm))
         , comm_(&comm)
-        , weights_(Opm::Amg::getQuasiImpesWeights<MatrixType, VectorType>(
+        , weights_(prm.get<VectorType>("weights"))
-              linearoperator.getmat(), prm.get<int>("pressure_var_index"), transpose))
         , levelTransferPolicy_(*comm_, weights_, prm.get<int>("pressure_var_index"))
         , coarseSolverPolicy_(prm.get_child("coarsesolver"))
         , twolevel_method_(linearoperator,
@@ -110,6 +123,8 @@ public:
                            transpose ? 0 : 1)
         , prm_(prm)
     {
+        //Opm::Amg::getQuasiImpesWeights<MatrixType, VectorType>(
+        //      linearoperator.getmat(), prm.get<int>("pressure_var_index"), transpose))
         if (prm.get<int>("verbosity") > 10) {
             std::ofstream outfile(prm.get<std::string>("weights_filename"));
             if (!outfile) {
@@ -134,10 +149,11 @@ public:
         twolevel_method_.post(x);
     }
 
-    virtual void update() override
+    virtual void update(const boost::property_tree::ptree& prm) override
     {
-        Opm::Amg::getQuasiImpesWeights<MatrixType, VectorType>(
+        //Opm::Amg::getQuasiImpesWeights<MatrixType, VectorType>(
-            linear_operator_.getmat(), prm_.get<int>("pressure_var_index"), transpose, weights_);
+        //    linear_operator_.getmat(), prm_.get<int>("pressure_var_index"), transpose, weights_);
+        weights_ = prm.get<VectorType>("weights");
         updateImpl(comm_);
     }
 

opm/simulators/linalg/ParallelOverlappingILU0.hpp

@@ -921,7 +921,7 @@ public:
         DUNE_UNUSED_PARAMETER(x);
     }
 
-    virtual void update() override
+    virtual void update(const boost::property_tree::ptree& = boost::property_tree::ptree()) override
     {
         // (For older DUNE versions the communicator might be
         // invalid if redistribution in AMG happened on the coarset level.

opm/simulators/linalg/PreconditionerWithUpdate.hpp

@@ -22,7 +22,7 @@
 
 #include <dune/istl/preconditioner.hh>
 #include <memory>
-
+#include <boost/property_tree/ptree.hpp>
 namespace Dune
 {
 
@@ -31,7 +31,7 @@ template <class X, class Y>
 class PreconditionerWithUpdate : public Preconditioner<X, Y>
 {
 public:
-    virtual void update() = 0;
+    virtual void update(const boost::property_tree::ptree& pt) = 0;
 };
 
 template <class OriginalPreconditioner>
@@ -69,7 +69,7 @@ public:
     }
 
     // The update() function does nothing for a wrapped preconditioner.
-    virtual void update() override
+    virtual void update(const boost::property_tree::ptree& /* pt */) override
     {
     }
 

opm/simulators/linalg/PressureSolverPolicy.hpp

@@ -84,7 +84,8 @@ namespace Amg
 
             void updatePreconditioner()
             {
-                linsolver_->preconditioner().update();
+                pt::ptree prm;
+                linsolver_->preconditioner().update(prm);
             }
 
         private:

opm/simulators/linalg/amgcpr.hh

@@ -239,8 +239,9 @@ namespace Dune
       /**
        * @brief Update the coarse solver and the hierarchies.
        */
-      virtual void update();
+      virtual void update(const boost::property_tree::ptree& prm);
 
+      virtual void update();
       /**
        * @brief Check whether the coarse solver used is a direct solver.
        * @return True if the coarse level solver is a direct solver.
@@ -468,6 +469,11 @@ namespace Dune
       update();
     }
 
+    template<class M, class X, class S, class PI, class A>
+    void AMGCPR<M,X,S,PI,A>::update(const boost::property_tree::ptree& /* prm */)
+    {
+      update();
+    }
     template<class M, class X, class S, class PI, class A>
     void AMGCPR<M,X,S,PI,A>::update()
     {