Allow well operators with FlexibleSolver.

opm/simulators/linalg/FlexibleSolver.hpp

@@ -43,54 +43,56 @@ public:
     using MatrixType = MatrixTypeT;
     using VectorType = VectorTypeT;
 
+    /// Base class type of the operator passed to the solver.
+    using AbstractOperatorType = Dune::AssembledLinearOperator<MatrixType, VectorType, VectorType>;
+    /// Base class type of the contained preconditioner.
+    using AbstractPrecondType = Dune::PreconditionerWithUpdate<VectorType, VectorType>;
+
     /// Create a sequential solver.
-    FlexibleSolver(const MatrixType& matrix,
+    FlexibleSolver(AbstractOperatorType& op,
                    const boost::property_tree::ptree& prm,
-                   const std::function<VectorTypeT()>& weightsCalculator = std::function<VectorTypeT()>());
+                   const std::function<VectorType()>& weightsCalculator = std::function<VectorType()>());
 
     /// Create a parallel solver (if Comm is e.g. OwnerOverlapCommunication).
     template <class Comm>
-    FlexibleSolver(const MatrixType& matrix,
+    FlexibleSolver(AbstractOperatorType& op,
                    const Comm& comm,
                    const boost::property_tree::ptree& prm,
-                   const std::function<VectorTypeT()>& weightsCalculator = std::function<VectorTypeT()>());
+                   const std::function<VectorType()>& weightsCalculator = std::function<VectorType()>());
 
     virtual void apply(VectorType& x, VectorType& rhs, Dune::InverseOperatorResult& res) override;
 
     virtual void apply(VectorType& x, VectorType& rhs, double reduction, Dune::InverseOperatorResult& res) override;
 
-    /// Type of the contained preconditioner.
-    using AbstractPrecondType = Dune::PreconditionerWithUpdate<VectorType, VectorType>;
-
     /// Access the contained preconditioner.
     AbstractPrecondType& preconditioner();
 
     virtual Dune::SolverCategory::Category category() const override;
 
 private:
-    using AbstractOperatorType = Dune::AssembledLinearOperator<MatrixType, VectorType, VectorType>;
     using AbstractScalarProductType = Dune::ScalarProduct<VectorType>;
     using AbstractSolverType = Dune::InverseOperator<VectorType, VectorType>;
 
     // Machinery for making sequential or parallel operators/preconditioners/scalar products.
     template <class Comm>
-    void initOpPrecSp(const MatrixType& matrix, const boost::property_tree::ptree& prm,
+    void initOpPrecSp(AbstractOperatorType& op, const boost::property_tree::ptree& prm,
-                      const std::function<VectorTypeT()> weightsCalculator, const Comm& comm);
+                      const std::function<VectorType()> weightsCalculator, const Comm& comm);
 
-    void initOpPrecSp(const MatrixType& matrix, const boost::property_tree::ptree& prm,
+    void initOpPrecSp(AbstractOperatorType& op, const boost::property_tree::ptree& prm,
-                      const std::function<VectorTypeT()> weightsCalculator, const Dune::Amg::SequentialInformation&);
+                      const std::function<VectorType()> weightsCalculator, const Dune::Amg::SequentialInformation&);
 
-    void initSolver(const boost::property_tree::ptree& prm, bool isMaster);
+    void initSolver(const boost::property_tree::ptree& prm, const bool is_iorank);
 
     // Main initialization routine.
     // Call with Comm == Dune::Amg::SequentialInformation to get a serial solver.
     template <class Comm>
-    void init(const MatrixType& matrix,
+    void init(AbstractOperatorType& op,
               const Comm& comm,
               const boost::property_tree::ptree& prm,
-              const std::function<VectorTypeT()> weightsCalculator);
+              const std::function<VectorType()> weightsCalculator);
 
-    std::shared_ptr<AbstractOperatorType> linearoperator_;
+    AbstractOperatorType* linearoperator_for_solver_;
+    std::shared_ptr<AbstractOperatorType> linearoperator_for_precond_;
     std::shared_ptr<AbstractPrecondType> preconditioner_;
     std::shared_ptr<AbstractScalarProductType> scalarproduct_;
     std::shared_ptr<AbstractSolverType> linsolver_;

opm/simulators/linalg/FlexibleSolver_impl.hpp

@@ -39,23 +39,23 @@ namespace Dune
     /// Create a sequential solver.
     template <class MatrixType, class VectorType>
     FlexibleSolver<MatrixType, VectorType>::
-    FlexibleSolver(const MatrixType& matrix,
+    FlexibleSolver(AbstractOperatorType& op,
                    const boost::property_tree::ptree& prm,
                    const std::function<VectorType()>& weightsCalculator)
     {
-        init(matrix, Dune::Amg::SequentialInformation(), prm, weightsCalculator);
+        init(op, Dune::Amg::SequentialInformation(), prm, weightsCalculator);
     }
 
     /// Create a parallel solver (if Comm is e.g. OwnerOverlapCommunication).
     template <class MatrixType, class VectorType>
     template <class Comm>
     FlexibleSolver<MatrixType, VectorType>::
-    FlexibleSolver(const MatrixType& matrix,
+    FlexibleSolver(AbstractOperatorType& op,
                    const Comm& comm,
                    const boost::property_tree::ptree& prm,
                    const std::function<VectorType()>& weightsCalculator)
     {
-        init(matrix, comm, prm, weightsCalculator);
+        init(op, comm, prm, weightsCalculator);
     }
 
     template <class MatrixType, class VectorType>
@@ -88,7 +88,7 @@ namespace Dune
     FlexibleSolver<MatrixType, VectorType>::
     category() const
     {
-        return linearoperator_->category();
+        return linearoperator_for_solver_->category();
     }
 
     // Machinery for making sequential or parallel operators/preconditioners/scalar products.
@@ -96,56 +96,64 @@ namespace Dune
     template <class Comm>
     void
     FlexibleSolver<MatrixType, VectorType>::
-    initOpPrecSp(const MatrixType& matrix, const boost::property_tree::ptree& prm,
+    initOpPrecSp(AbstractOperatorType& op,
-                 const std::function<VectorType()> weightsCalculator, const Comm& comm)
+                 const boost::property_tree::ptree& prm,
+                 const std::function<VectorType()> weightsCalculator,
+                 const Comm& comm)
     {
         // Parallel case.
-        using ParOperatorType = Dune::OverlappingSchwarzOperator<MatrixType, VectorType, VectorType, Comm>;
         using pt = const boost::property_tree::ptree;
-        auto linop = std::make_shared<ParOperatorType>(matrix, comm);
+        using ParOperatorType = Dune::OverlappingSchwarzOperator<MatrixType, VectorType, VectorType, Comm>;
-        linearoperator_ = linop;
+        linearoperator_for_solver_ = &op;
+        auto op_prec = std::make_shared<ParOperatorType>(op.getmat(), comm);
         auto child = prm.get_child_optional("preconditioner");
-        preconditioner_
+        preconditioner_ = Opm::PreconditionerFactory<ParOperatorType, Comm>::create(*op_prec,
-            = Opm::PreconditionerFactory<ParOperatorType, Comm>::create(*linop, child? *child : pt(),
+                                                                                    child ? *child : pt(),
-                                                                        weightsCalculator, comm);
+                                                                                    weightsCalculator,
-        scalarproduct_ = Dune::createScalarProduct<VectorType, Comm>(comm, linearoperator_->category());
+                                                                                    comm);
+        scalarproduct_ = Dune::createScalarProduct<VectorType, Comm>(comm, op.category());
+        linearoperator_for_precond_ = op_prec;
     }
 
     template <class MatrixType, class VectorType>
     void
     FlexibleSolver<MatrixType, VectorType>::
-    initOpPrecSp(const MatrixType& matrix, const boost::property_tree::ptree& prm,
+    initOpPrecSp(AbstractOperatorType& op,
-                 const std::function<VectorType()> weightsCalculator, const Dune::Amg::SequentialInformation&)
+                 const boost::property_tree::ptree& prm,
+                 const std::function<VectorType()> weightsCalculator,
+                 const Dune::Amg::SequentialInformation&)
     {
         // Sequential case.
-        using SeqOperatorType = Dune::MatrixAdapter<MatrixType, VectorType, VectorType>;
         using pt = const boost::property_tree::ptree;
-        auto linop = std::make_shared<SeqOperatorType>(matrix);
+        using SeqOperatorType = Dune::MatrixAdapter<MatrixType, VectorType, VectorType>;
-        linearoperator_ = linop;
+        linearoperator_for_solver_ = &op;
+        auto op_prec = std::make_shared<SeqOperatorType>(op.getmat());
         auto child = prm.get_child_optional("preconditioner");
-        preconditioner_ = Opm::PreconditionerFactory<SeqOperatorType>::create(*linop, child? *child : pt(),
+        preconditioner_ = Opm::PreconditionerFactory<SeqOperatorType>::create(*op_prec,
+                                                                              child ? *child : pt(),
                                                                               weightsCalculator);
         scalarproduct_ = std::make_shared<Dune::SeqScalarProduct<VectorType>>();
+        linearoperator_for_precond_ = op_prec;
     }
 
     template <class MatrixType, class VectorType>
     void
     FlexibleSolver<MatrixType, VectorType>::
-    initSolver(const boost::property_tree::ptree& prm, bool isMaster)
+    initSolver(const boost::property_tree::ptree& prm, const bool is_iorank)
     {
         const double tol = prm.get<double>("tol", 1e-2);
         const int maxiter = prm.get<int>("maxiter", 200);
-        const int verbosity = isMaster? prm.get<int>("verbosity", 0) : 0;
+        const int verbosity = is_iorank ? prm.get<int>("verbosity", 0) : 0;
         const std::string solver_type = prm.get<std::string>("solver", "bicgstab");
         if (solver_type == "bicgstab") {
-            linsolver_.reset(new Dune::BiCGSTABSolver<VectorType>(*linearoperator_,
+            linsolver_.reset(new Dune::BiCGSTABSolver<VectorType>(*linearoperator_for_solver_,
                                                                   *scalarproduct_,
                                                                   *preconditioner_,
                                                                   tol, // desired residual reduction factor
                                                                   maxiter, // maximum number of iterations
                                                                   verbosity));
         } else if (solver_type == "loopsolver") {
-            linsolver_.reset(new Dune::LoopSolver<VectorType>(*linearoperator_,
+            linsolver_.reset(new Dune::LoopSolver<VectorType>(*linearoperator_for_solver_,
                                                               *scalarproduct_,
                                                               *preconditioner_,
                                                               tol, // desired residual reduction factor
@@ -153,7 +161,7 @@ namespace Dune
                                                               verbosity));
         } else if (solver_type == "gmres") {
             int restart = prm.get<int>("restart", 15);
-            linsolver_.reset(new Dune::RestartedGMResSolver<VectorType>(*linearoperator_,
+            linsolver_.reset(new Dune::RestartedGMResSolver<VectorType>(*linearoperator_for_solver_,
                                                                         *scalarproduct_,
                                                                         *preconditioner_,
                                                                         tol,
@@ -163,7 +171,7 @@ namespace Dune
 #if HAVE_SUITESPARSE_UMFPACK
         } else if (solver_type == "umfpack") {
             bool dummy = false;
-            linsolver_.reset(new Dune::UMFPack<MatrixType>(linearoperator_->getmat(), verbosity, dummy));
+            linsolver_.reset(new Dune::UMFPack<MatrixType>(linearoperator_for_solver_->getmat(), verbosity, dummy));
 #endif
         } else {
             OPM_THROW(std::invalid_argument, "Properties: Solver " << solver_type << " not known.");
@@ -177,13 +185,13 @@ namespace Dune
     template <class Comm>
     void
     FlexibleSolver<MatrixType, VectorType>::
-    init(const MatrixType& matrix,
+    init(AbstractOperatorType& op,
          const Comm& comm,
          const boost::property_tree::ptree& prm,
          const std::function<VectorType()> weightsCalculator)
     {
-        initOpPrecSp(matrix, prm, weightsCalculator, comm);
+        initOpPrecSp(op, prm, weightsCalculator, comm);
-        initSolver(prm, comm.communicator().rank()==0);
+        initSolver(prm, comm.communicator().rank() == 0);
     }
 
 } // namespace Dune
@@ -198,28 +206,33 @@ using BM = Dune::BCRSMatrix<Dune::FieldMatrix<double, N, N>>;
 template <int N>
 using OBM = Dune::BCRSMatrix<Opm::MatrixBlock<double, N, N>>;
 
-// INSTANTIATE_CONSTRUCTOR instantiates the constructor that is a template,
-// this is only needed in the MPI case, since otherwise the Comm type is
-// not a template argument but always SequentialInformation.
 #if HAVE_MPI
+
 using Comm = Dune::OwnerOverlapCopyCommunication<int, int>;
-#define INSTANTIATE_FLEXIBLESOLVER_CONSTRUCTOR(n) \
+
-template Dune::FlexibleSolver<OBM<n>, BV<n>>::FlexibleSolver(const MatrixType& matrix,                        \
+// Note: we must instantiate the constructor that is a template.
+// This is only needed in the parallel case, since otherwise the Comm type is
+// not a template argument but always SequentialInformation.
+
+#define INSTANTIATE_FLEXIBLESOLVER(N)                  \
+template class Dune::FlexibleSolver<BM<N>, BV<N>>;     \
+template class Dune::FlexibleSolver<OBM<N>, BV<N>>;    \
+template Dune::FlexibleSolver<BM<N>, BV<N>>::FlexibleSolver(AbstractOperatorType& op,                         \
+                                                            const Comm& comm,                                 \
+                                                            const boost::property_tree::ptree& prm,           \
+                                                            const std::function<BV<N>()>& weightsCalculator); \
+template Dune::FlexibleSolver<OBM<N>, BV<N>>::FlexibleSolver(AbstractOperatorType& op,                        \
                                                              const Comm& comm,                                \
                                                              const boost::property_tree::ptree& prm,          \
-                                                             const std::function<BV<n>()>& weightsCalculator);
+                                                             const std::function<BV<N>()>& weightsCalculator);
-#else
-#define INSTANTIATE_FLEXIBLESOLVER_CONSTRUCTOR(n)
-#endif
 
-// INSTANTIATE instantiates the class including any templated constructors if necessary.
+#else // HAVE_MPI
-#define INSTANTIATE_FLEXIBLESOLVER(n)               \
-/* Variants using Dune::FieldMatrix blocks. */      \
-template class Dune::FlexibleSolver<BM<n>, BV<n>>;  \
-/* Variants using Opm::MatrixBlock blocks. */       \
-template class Dune::FlexibleSolver<OBM<n>, BV<n>>; \
-INSTANTIATE_FLEXIBLESOLVER_CONSTRUCTOR(n)
 
+#define INSTANTIATE_FLEXIBLESOLVER(N)                  \
+template class Dune::FlexibleSolver<BM<N>, BV<N>>;     \
+template class Dune::FlexibleSolver<OBM<N>, BV<N>>;
+
+#endif // HAVE_MPI
 
 
 #endif // OPM_FLEXIBLE_SOLVER_IMPL_HEADER_INCLUDED

opm/simulators/linalg/ISTLSolverEbos.hpp

@@ -113,6 +113,8 @@ DenseMatrix transposeDenseMatrix(const DenseMatrix& M)
         typedef typename GridView::template Codim<0>::Entity Element;
         typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
         using FlexibleSolverType = Dune::FlexibleSolver<Matrix, Vector>;
+        using AbstractOperatorType = Dune::AssembledLinearOperator<Matrix, Vector, Vector>;
+        using WellModelOperator = WellModelAsLinearOperator<WellModel, Vector, Vector>;
         // Due to miscibility oil <-> gas the water eqn is the one we can replace with a pressure equation.
         static const bool waterEnabled = Indices::waterEnabled;
         static const int pindex = (waterEnabled) ? BlackOilDefaultIndexTraits::waterCompIdx : BlackOilDefaultIndexTraits::oilCompIdx;
@@ -176,17 +178,10 @@ DenseMatrix transposeDenseMatrix(const DenseMatrix& M)
 #endif
             extractParallelGridInformationToISTL(simulator_.vanguard().grid(), parallelInformation_);
             useWellConn_ = EWOMS_GET_PARAM(TypeTag, bool, MatrixAddWellContributions);
-
-            if (!useWellConn_ && useFlexible_)
-            {
-                OPM_THROW(std::logic_error, "Flexible solvers and CPR need the well contribution in the matrix. Please run with"
-                          " --matrix-add-well-contributions=true");
-            }
-
             ownersFirst_ = EWOMS_GET_PARAM(TypeTag, bool, OwnerCellsFirst);
             interiorCellNum_ = detail::numMatrixRowsToUseInSolver(simulator_.vanguard().grid(), ownersFirst_);
 
-            if ( isParallel() && (!ownersFirst_ || parameters_.linear_solver_use_amg_  || useFlexible_ ) ) {
+            if ( isParallel() && (!ownersFirst_ || parameters_.linear_solver_use_amg_) ) {
                 detail::setWellConnections(gridForConn, simulator_.vanguard().schedule().getWellsatEnd(), useWellConn_, wellConnectionsGraph_);
                 // For some reason simulator_.model().elementMapper() is not initialized at this stage
                 // Hence const auto& elemMapper = simulator_.model().elementMapper(); does not work.
@@ -741,11 +736,26 @@ DenseMatrix transposeDenseMatrix(const DenseMatrix& M)
             if (recreate_solver || !flexibleSolver_) {
                 if (isParallel()) {
 #if HAVE_MPI
-                    assert(noGhostMat_);
+                    if (useWellConn_) {
-                    flexibleSolver_.reset(new FlexibleSolverType(getMatrix(), *comm_, prm_, weightsCalculator));
+                        assert(noGhostMat_);
+                        using ParOperatorType = Dune::OverlappingSchwarzOperator<Matrix, Vector, Vector, Comm>;
+                        linearOperatorForFlexibleSolver_ = std::make_unique<ParOperatorType>(getMatrix(), *comm_);
+                        flexibleSolver_ = std::make_unique<FlexibleSolverType>(*linearOperatorForFlexibleSolver_, *comm_, prm_, weightsCalculator);
+                    } else {
+                        if (!ownersFirst_) {
+                            OPM_THROW(std::runtime_error, "In parallel, the flexible solver requires "
+                                      "--owner-cells-first=true when --matrix-add-well-contributions=false is used.");
+                        }
+                        using ParOperatorType = WellModelGhostLastMatrixAdapter<Matrix, Vector, Vector, true>;
+                        wellOperator_ = std::make_unique<WellModelOperator>(simulator_.problem().wellModel());
+                        linearOperatorForFlexibleSolver_ = std::make_unique<ParOperatorType>(getMatrix(), *wellOperator_, interiorCellNum_);
+                        flexibleSolver_ = std::make_unique<FlexibleSolverType>(*linearOperatorForFlexibleSolver_, *comm_, prm_, weightsCalculator);
+                    }
 #endif
                 } else {
-                    flexibleSolver_.reset(new FlexibleSolverType(getMatrix(), prm_, weightsCalculator));
+                    using SeqLinearOperator = Dune::MatrixAdapter<Matrix, Vector, Vector>;
+                    linearOperatorForFlexibleSolver_ = std::make_unique<SeqLinearOperator>(getMatrix());
+                    flexibleSolver_ = std::make_unique<FlexibleSolverType>(*linearOperatorForFlexibleSolver_, prm_, weightsCalculator);
                 }
             }
             else
@@ -1004,6 +1014,8 @@ DenseMatrix transposeDenseMatrix(const DenseMatrix& M)
         Vector *rhs_;
 
         std::unique_ptr<FlexibleSolverType> flexibleSolver_;
+        std::unique_ptr<AbstractOperatorType> linearOperatorForFlexibleSolver_;
+        std::unique_ptr<WellModelAsLinearOperator<WellModel, Vector, Vector>> wellOperator_;
         std::vector<int> overlapRows_;
         std::vector<int> interiorRows_;
         std::vector<std::set<int>> wellConnectionsGraph_;

opm/simulators/linalg/ISTLSolverEbosFlexible.hpp

@@ -61,18 +61,19 @@ class ISTLSolverEbosFlexible
     using Simulator = typename GET_PROP_TYPE(TypeTag, Simulator);
     using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
     using MatrixType = typename SparseMatrixAdapter::IstlMatrix;
+    using WellModel = typename GET_PROP_TYPE(TypeTag, EclWellModel);
 #if HAVE_MPI
     using Communication = Dune::OwnerOverlapCopyCommunication<int, int>;
 #else
     using Communication = int; // Dummy type.
 #endif
+    using AbstractOperatorType = Dune::AssembledLinearOperator<MatrixType, VectorType, VectorType>;
+    using WellModelOpType = WellModelAsLinearOperator<WellModel, VectorType, VectorType>;
     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;
@@ -90,6 +91,9 @@ public:
 
     explicit ISTLSolverEbosFlexible(const Simulator& simulator)
         : simulator_(simulator)
+        , ownersFirst_(EWOMS_GET_PARAM(TypeTag, bool, OwnerCellsFirst))
+        , matrixAddWellContributions_(EWOMS_GET_PARAM(TypeTag, bool, MatrixAddWellContributions))
+        , interiorCellNum_(detail::numMatrixRowsToUseInSolver(simulator_.vanguard().grid(), ownersFirst_))
     {
         parameters_.template init<TypeTag>();
         prm_ = setupPropertyTree<TypeTag>(parameters_);
@@ -134,76 +138,53 @@ public:
             parinfo->copyValuesTo(comm_->indexSet(), comm_->remoteIndices(), size, 1);
             firstcall = false;
         }
-        makeOverlapRowsInvalid(mat.istlMatrix());
+        if (isParallel() && matrixAddWellContributions_) {
+            makeOverlapRowsInvalid(mat.istlMatrix());
+        }
 #endif
-        // Decide if we should recreate the solver or just do
-        // a minimal preconditioner update.
-        const int newton_iteration = this->simulator_.model().newtonMethod().numIterations();
-        bool recreate_solver = false;
-        if (this->parameters_.cpr_reuse_setup_ == 0) {
-            // Always recreate solver.
-            recreate_solver = true;
-        } else if (this->parameters_.cpr_reuse_setup_ == 1) {
-            // Recreate solver on the first iteration of every timestep.
-            if (newton_iteration == 0) {
-                recreate_solver = true;
-            }
-        } else if (this->parameters_.cpr_reuse_setup_ == 2) {
-            // Recreate solver if the last solve used more than 10 iterations.
-            if (this->iterations() > 10) {
-                recreate_solver = true;
-            }
-        } else {
-            assert(this->parameters_.cpr_reuse_setup_ == 3);
-            assert(recreate_solver == false);
-            // Never recreate solver.
-        }
 
-        std::function<VectorType()> weightsCalculator;
+        matrix_ = &mat.istlMatrix(); // Store pointer for output if needed.
+        std::function<VectorType()> weightsCalculator = getWeightsCalculator(mat.istlMatrix(), b);
 
-        auto preconditionerType = prm_.get("preconditioner.type", "cpr");
+        if (shouldCreateSolver()) {
-        if( preconditionerType  == "cpr" ||
-            preconditionerType == "cprt"
-            )
-        {
-            bool transpose = false;
-            if(preconditionerType == "cprt"){
-                transpose = true;
-            }
-
-            auto weightsType = prm_.get("preconditioner.weight_type", "quasiimpes");
-            auto pressureIndex = this->prm_.get("preconditioner.pressure_var_index", 1);
-            if(weightsType == "quasiimpes") {
-                // weighs will be created as default in the solver
-                weightsCalculator =
-                    [&mat, transpose, pressureIndex](){
-                        return Opm::Amg::getQuasiImpesWeights<MatrixType,
-                                                              VectorType>(
-                                                                          mat.istlMatrix(),
-                                                                          pressureIndex,
-                                                                          transpose);
-                    };
-
-            }else if(weightsType == "trueimpes"  ){
-                weightsCalculator =
-                    [this, &b, pressureIndex](){
-                        return this->getTrueImpesWeights(b, pressureIndex);
-                    };
-            }else{
-                OPM_THROW(std::invalid_argument, "Weights type " << weightsType << "not implemented for cpr."
-                          << " Please use quasiimpes or trueimpes.");
-            }
-        }
-
-        if (recreate_solver || !solver_) {
             if (isParallel()) {
 #if HAVE_MPI
-                matrix_ = &mat.istlMatrix();
+                if (matrixAddWellContributions_) {
-                solver_.reset(new SolverType(mat.istlMatrix(), *comm_, prm_, weightsCalculator));
+                    using ParOperatorType = Dune::OverlappingSchwarzOperator<MatrixType, VectorType, VectorType, Communication>;
-#endif
+                    auto op = std::make_unique<ParOperatorType>(mat.istlMatrix(), *comm_);
+                    auto sol = std::make_unique<SolverType>(*op, *comm_, prm_, weightsCalculator);
+                    solver_ = std::move(sol);
+                    linear_operator_ = std::move(op);
+                } else {
+                    if (!ownersFirst_) {
+                        OPM_THROW(std::runtime_error, "In parallel, the flexible solver requires "
+                                  "--owner-cells-first=true when --matrix-add-well-contributions=false is used.");
+                    }
+                    using ParOperatorType = WellModelGhostLastMatrixAdapter<MatrixType, VectorType, VectorType, true>;
+                    auto well_op = std::make_unique<WellModelOpType>(simulator_.problem().wellModel());
+                    auto op = std::make_unique<ParOperatorType>(mat.istlMatrix(), *well_op, interiorCellNum_);
+                    auto sol = std::make_unique<SolverType>(*op, *comm_, prm_, weightsCalculator);
+                    solver_ = std::move(sol);
+                    linear_operator_ = std::move(op);
+                    well_operator_ = std::move(well_op);
+                }
+#endif // HAVE_MPI
             } else {
-                matrix_ = &mat.istlMatrix();
+                if (matrixAddWellContributions_) {
-                solver_.reset(new SolverType(mat.istlMatrix(), prm_, weightsCalculator));
+                    using SeqOperatorType = Dune::MatrixAdapter<MatrixType, VectorType, VectorType>;
+                    auto op = std::make_unique<SeqOperatorType>(mat.istlMatrix());
+                    auto sol = std::make_unique<SolverType>(*op, prm_, weightsCalculator);
+                    solver_ = std::move(sol);
+                    linear_operator_ = std::move(op);
+                } else {
+                    using SeqOperatorType = WellModelMatrixAdapter<MatrixType, VectorType, VectorType, false>;
+                    auto well_op = std::make_unique<WellModelOpType>(simulator_.problem().wellModel());
+                    auto op = std::make_unique<SeqOperatorType>(mat.istlMatrix(), *well_op);
+                    auto sol = std::make_unique<SolverType>(*op, prm_, weightsCalculator);
+                    solver_ = std::move(sol);
+                    linear_operator_ = std::move(op);
+                    well_operator_ = std::move(well_op);
+                }
             }
             rhs_ = b;
         } else {
@@ -245,6 +226,63 @@ public:
 
 protected:
 
+    bool shouldCreateSolver() const
+    {
+        // Decide if we should recreate the solver or just do
+        // a minimal preconditioner update.
+        if (!solver_) {
+            return true;
+        }
+        const int newton_iteration = this->simulator_.model().newtonMethod().numIterations();
+        bool recreate_solver = false;
+        if (this->parameters_.cpr_reuse_setup_ == 0) {
+            // Always recreate solver.
+            recreate_solver = true;
+        } else if (this->parameters_.cpr_reuse_setup_ == 1) {
+            // Recreate solver on the first iteration of every timestep.
+            if (newton_iteration == 0) {
+                recreate_solver = true;
+            }
+        } else if (this->parameters_.cpr_reuse_setup_ == 2) {
+            // Recreate solver if the last solve used more than 10 iterations.
+            if (this->iterations() > 10) {
+                recreate_solver = true;
+            }
+        } else {
+            assert(this->parameters_.cpr_reuse_setup_ == 3);
+            assert(recreate_solver == false);
+            // Never recreate solver.
+        }
+        return recreate_solver;
+    }
+
+    std::function<VectorType()> getWeightsCalculator(const MatrixType& mat, const VectorType& b) const
+    {
+        std::function<VectorType()> weightsCalculator;
+
+        auto preconditionerType = prm_.get("preconditioner.type", "cpr");
+        if (preconditionerType == "cpr" || preconditionerType == "cprt") {
+            const bool transpose = preconditionerType == "cprt";
+            const auto weightsType = prm_.get("preconditioner.weight_type", "quasiimpes");
+            const auto pressureIndex = this->prm_.get("preconditioner.pressure_var_index", 1);
+            if (weightsType == "quasiimpes") {
+                // weighs will be created as default in the solver
+                weightsCalculator = [&mat, transpose, pressureIndex]() {
+                    return Opm::Amg::getQuasiImpesWeights<MatrixType, VectorType>(mat, pressureIndex, transpose);
+                };
+            } else if (weightsType == "trueimpes") {
+                weightsCalculator = [this, &b, pressureIndex]() {
+                    return this->getTrueImpesWeights(b, pressureIndex);
+                };
+            } else {
+                OPM_THROW(std::invalid_argument,
+                          "Weights type " << weightsType << "not implemented for cpr."
+                                          << " Please use quasiimpes or trueimpes.");
+            }
+        }
+        return weightsCalculator;
+    }
+
     /// Zero out off-diagonal blocks on rows corresponding to overlap cells
     /// Diagonal blocks on ovelap rows are set to diag(1.0).
     void makeOverlapRowsInvalid(MatrixType& matrix) const
@@ -267,7 +305,7 @@ protected:
         }
     }
 
-    VectorType getTrueImpesWeights(const VectorType& b,const int pressureVarIndex)
+    VectorType getTrueImpesWeights(const VectorType& b, const int pressureVarIndex) const
     {
         VectorType weights(b.size());
         ElementContext elemCtx(simulator_);
@@ -289,14 +327,20 @@ protected:
         }
     }
 
+
     const Simulator& simulator_;
     MatrixType* matrix_;
+    std::unique_ptr<WellModelOpType> well_operator_;
+    std::unique_ptr<AbstractOperatorType> linear_operator_;
     std::unique_ptr<SolverType> solver_;
     FlowLinearSolverParameters parameters_;
     boost::property_tree::ptree prm_;
     VectorType rhs_;
     Dune::InverseOperatorResult res_;
     std::any parallelInformation_;
+    bool ownersFirst_;
+    bool matrixAddWellContributions_;
+    int interiorCellNum_;
     std::unique_ptr<Communication> comm_;
     std::vector<int> overlapRows_;
     std::vector<int> interiorRows_;

opm/simulators/linalg/PressureSolverPolicy.hpp

@@ -25,6 +25,7 @@
 #include <boost/property_tree/ptree.hpp>
 
 #include <dune/istl/solver.hh>
+#include <dune/istl/owneroverlapcopy.hh>
 
 namespace Dune
 {
@@ -55,19 +56,24 @@ namespace Amg
          * The operator will use one step of AMG to approximately solve
          * the coarse level system.
          */
-        struct PressureInverseOperator : public Dune::InverseOperator<X, X> {
+        struct PressureInverseOperator : public Dune::InverseOperator<X, X>
-            template <class Comm>
+        {
-            PressureInverseOperator(Operator& op, const boost::property_tree::ptree& prm, const Comm& comm)
+            template <typename GlobalIndex, typename LocalIndex>
+            PressureInverseOperator(Operator& op,
+                                    const boost::property_tree::ptree& prm,
+                                    const Dune::OwnerOverlapCopyCommunication<GlobalIndex, LocalIndex>& comm)
                 : linsolver_()
             {
                 assert(op.category() == Dune::SolverCategory::overlapping);
-                linsolver_ = std::make_unique<Solver>(op.getmat(), comm, prm, std::function<X()>());
+                linsolver_ = std::make_unique<Solver>(op, comm, prm, std::function<X()>());
             }
-            PressureInverseOperator(Operator& op, const boost::property_tree::ptree& prm, const SequentialInformation&)
+            PressureInverseOperator(Operator& op,
+                                    const boost::property_tree::ptree& prm,
+                                    const SequentialInformation&)
                 : linsolver_()
             {
                 assert(op.category() != Dune::SolverCategory::overlapping);
-                linsolver_ = std::make_unique<Solver>(op.getmat(), prm, std::function<X()>());
+                linsolver_ = std::make_unique<Solver>(op, prm, std::function<X()>());
             }
 
 

tests/test_flexiblesolver.cpp

@@ -75,7 +75,9 @@ testSolver(const boost::property_tree::ptree& prm, const std::string& matrix_fil
                                                                 prm.get<int>("preconditioner.pressure_var_index"),
                                                                 transpose);
               };
-    Dune::FlexibleSolver<Matrix, Vector> solver(matrix, prm, wc);
+    using SeqOperatorType = Dune::MatrixAdapter<Matrix, Vector, Vector>;
+    SeqOperatorType op(matrix);
+    Dune::FlexibleSolver<Matrix, Vector> solver(op, prm, wc);
     Vector x(rhs.size());
     Dune::InverseOperatorResult res;
     solver.apply(x, rhs, res);