Merge pull request #3937 from atgeirr/pressure-bhp-cpr

Pressure bhp cpr

compareECLFiles.cmake

@@ -1482,7 +1482,7 @@ if(MPI_FOUND)
                                        ABS_TOL 0.17
                                        REL_TOL ${coarse_rel_tol_parallel}
                                        DIR aquifer-num
-                                       TEST_ARGS --tolerance-cnv=0.00003 --time-step-control=pid --linsolver=cpr)
+                                       TEST_ARGS --tolerance-cnv=0.000003 --time-step-control=pid --linsolver=cpr)
 
   add_test_compare_parallel_simulation(CASENAME numerical_aquifer_3d_1aqu
                                        FILENAME 3D_1AQU_3CELLS

ebos/eclgenerictracermodel.cc

@@ -63,7 +63,7 @@ struct TracerSolverSelector
 {
     using Comm = Dune::OwnerOverlapCopyCommunication<int, int>;
     using TracerOperator = Dune::OverlappingSchwarzOperator<M, V, V, Comm>;
-    using type = Dune::FlexibleSolver<M, V>;
+    using type = Dune::FlexibleSolver<TracerOperator>;
 };
 template<class Vector, class Grid, class Matrix>
 std::tuple<std::unique_ptr<Dune::OverlappingSchwarzOperator<Matrix,Vector,Vector,
@@ -83,7 +83,7 @@ createParallelFlexibleSolver(const Dune::CpGrid& grid, const Matrix& M, const Pr
 {
         using TracerOperator = Dune::OverlappingSchwarzOperator<Matrix,Vector,Vector,
                                                                 Dune::OwnerOverlapCopyCommunication<int,int>>;
-        using TracerSolver = Dune::FlexibleSolver<Matrix, Vector>;
+        using TracerSolver = Dune::FlexibleSolver<TracerOperator>;
         const auto& cellComm = grid.cellCommunication();
         auto op = std::make_unique<TracerOperator>(M, cellComm);
         auto dummyWeights = [](){ return Vector();};

ebos/eclgenerictracermodel.hh

@@ -48,7 +48,7 @@ class EclipseState;
 template<class Grid, class GridView, class DofMapper, class Stencil, class Scalar>
 class EclGenericTracerModel {
 public:
-    using TracerMatrix = Dune::BCRSMatrix<Dune::FieldMatrix<Scalar, 1, 1>>;
+    using TracerMatrix = Dune::BCRSMatrix<Opm::MatrixBlock<Scalar, 1, 1>>;
     using TracerVector = Dune::BlockVector<Dune::FieldVector<Scalar,1>>;
     using CartesianIndexMapper = Dune::CartesianIndexMapper<Grid>;
     static const int dimWorld = Grid::dimensionworld;

opm-simulators-prereqs.cmake

@@ -13,6 +13,7 @@ set (opm-simulators_CONFIG_VAR
   HAVE_VEXCL
   HAVE_SUITESPARSE_UMFPACK_H
   HAVE_DUNE_ISTL
+  DUNE_ISTL_WITH_CHECKING
   DUNE_ISTL_VERSION_MAJOR
   DUNE_ISTL_VERSION_MINOR
   DUNE_ISTL_VERSION_REVISION

opm/simulators/linalg/FlexibleSolver.hpp

@@ -35,27 +35,25 @@ namespace Dune
 /// (operator, scalar product, iterative solver and preconditioner) and sets
 /// them up based on runtime parameters, using the PreconditionerFactory for
 /// setting up preconditioners.
-template <class MatrixTypeT, class VectorTypeT>
+template <class Operator>
-class FlexibleSolver : public Dune::InverseOperator<VectorTypeT, VectorTypeT>
+class FlexibleSolver : public Dune::InverseOperator<typename Operator::domain_type,
+                                                    typename Operator::range_type>
 {
 public:
-    using MatrixType = MatrixTypeT;
+    using VectorType = typename Operator::domain_type; // Assuming symmetry: domain == range
-    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(AbstractOperatorType& op,
+    FlexibleSolver(Operator& op,
                    const Opm::PropertyTree& prm,
                    const std::function<VectorType()>& weightsCalculator,
                    std::size_t pressureIndex);
 
     /// Create a parallel solver (if Comm is e.g. OwnerOverlapCommunication).
     template <class Comm>
-    FlexibleSolver(AbstractOperatorType& op,
+    FlexibleSolver(Operator& op,
                    const Comm& comm,
                    const Opm::PropertyTree& prm,
                    const std::function<VectorType()>& weightsCalculator,
@@ -76,11 +74,11 @@ private:
 
     // Machinery for making sequential or parallel operators/preconditioners/scalar products.
     template <class Comm>
-    void initOpPrecSp(AbstractOperatorType& op, const Opm::PropertyTree& prm,
+    void initOpPrecSp(Operator& op, const Opm::PropertyTree& prm,
                       const std::function<VectorType()> weightsCalculator, const Comm& comm,
                       std::size_t pressureIndex);
 
-    void initOpPrecSp(AbstractOperatorType& op, const Opm::PropertyTree& prm,
+    void initOpPrecSp(Operator& op, const Opm::PropertyTree& prm,
                       const std::function<VectorType()> weightsCalculator, const Dune::Amg::SequentialInformation&,
                       std::size_t pressureIndex);
 
@@ -89,14 +87,13 @@ private:
     // Main initialization routine.
     // Call with Comm == Dune::Amg::SequentialInformation to get a serial solver.
     template <class Comm>
-    void init(AbstractOperatorType& op,
+    void init(Operator& op,
               const Comm& comm,
               const Opm::PropertyTree& prm,
               const std::function<VectorType()> weightsCalculator,
               std::size_t pressureIndex);
 
-    AbstractOperatorType* linearoperator_for_solver_;
+    Operator* 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_;
@@ -104,6 +101,7 @@ private:
 
 } // namespace Dune
 
+//#include <opm/simulators/linalg/FlexibleSolver_impl.hpp>
 
 
 #endif // OPM_FLEXIBLE_SOLVER_HEADER_INCLUDED

opm/simulators/linalg/FlexibleSolver_impl.hpp

@@ -36,9 +36,9 @@
 namespace Dune
 {
     /// Create a sequential solver.
-    template <class MatrixType, class VectorType>
+    template <class Operator>
-    FlexibleSolver<MatrixType, VectorType>::
+    FlexibleSolver<Operator>::
-    FlexibleSolver(AbstractOperatorType& op,
+    FlexibleSolver(Operator& op,
                    const Opm::PropertyTree& prm,
                    const std::function<VectorType()>& weightsCalculator,
                    std::size_t pressureIndex)
@@ -48,10 +48,10 @@ namespace Dune
     }
 
     /// Create a parallel solver (if Comm is e.g. OwnerOverlapCommunication).
-    template <class MatrixType, class VectorType>
+    template <class Operator>
     template <class Comm>
-    FlexibleSolver<MatrixType, VectorType>::
+    FlexibleSolver<Operator>::
-    FlexibleSolver(AbstractOperatorType& op,
+    FlexibleSolver(Operator& op,
                    const Comm& comm,
                    const Opm::PropertyTree& prm,
                    const std::function<VectorType()>& weightsCalculator,
@@ -60,89 +60,83 @@ namespace Dune
         init(op, comm, prm, weightsCalculator, pressureIndex);
     }
 
-    template <class MatrixType, class VectorType>
+    template <class Operator>
     void
-    FlexibleSolver<MatrixType, VectorType>::
+    FlexibleSolver<Operator>::
     apply(VectorType& x, VectorType& rhs, Dune::InverseOperatorResult& res)
     {
         linsolver_->apply(x, rhs, res);
     }
 
-    template <class MatrixType, class VectorType>
+    template <class Operator>
     void
-    FlexibleSolver<MatrixType, VectorType>::
+    FlexibleSolver<Operator>::
     apply(VectorType& x, VectorType& rhs, double reduction, Dune::InverseOperatorResult& res)
     {
         linsolver_->apply(x, rhs, reduction, res);
     }
 
     /// Access the contained preconditioner.
-    template <class MatrixType, class VectorType>
+    template <class Operator>
     auto
-    FlexibleSolver<MatrixType, VectorType>::
+    FlexibleSolver<Operator>::
     preconditioner() -> AbstractPrecondType&
     {
         return *preconditioner_;
     }
 
-    template <class MatrixType, class VectorType>
+    template <class Operator>
     Dune::SolverCategory::Category
-    FlexibleSolver<MatrixType, VectorType>::
+    FlexibleSolver<Operator>::
     category() const
     {
         return linearoperator_for_solver_->category();
     }
 
     // Machinery for making sequential or parallel operators/preconditioners/scalar products.
-    template <class MatrixType, class VectorType>
+    template <class Operator>
     template <class Comm>
     void
-    FlexibleSolver<MatrixType, VectorType>::
+    FlexibleSolver<Operator>::
-    initOpPrecSp(AbstractOperatorType& op,
+    initOpPrecSp(Operator& op,
                  const Opm::PropertyTree& prm,
                  const std::function<VectorType()> weightsCalculator,
                  const Comm& comm,
                  std::size_t pressureIndex)
     {
         // Parallel case.
-        using ParOperatorType = Dune::OverlappingSchwarzOperator<MatrixType, VectorType, VectorType, Comm>;
         linearoperator_for_solver_ = &op;
-        auto op_prec = std::make_shared<ParOperatorType>(op.getmat(), comm);
         auto child = prm.get_child_optional("preconditioner");
-        preconditioner_ = Opm::PreconditionerFactory<ParOperatorType, Comm>::create(*op_prec,
+        preconditioner_ = Opm::PreconditionerFactory<Operator, Comm>::create(op,
-                                                                                    child ? *child : Opm::PropertyTree(),
+                                                                             child ? *child : Opm::PropertyTree(),
-                                                                                    weightsCalculator,
+                                                                             weightsCalculator,
-                                                                                    comm,
+                                                                             comm,
-                                                                                    pressureIndex);
+                                                                             pressureIndex);
         scalarproduct_ = Dune::createScalarProduct<VectorType, Comm>(comm, op.category());
-        linearoperator_for_precond_ = op_prec;
     }
 
-    template <class MatrixType, class VectorType>
+    template <class Operator>
     void
-    FlexibleSolver<MatrixType, VectorType>::
+    FlexibleSolver<Operator>::
-    initOpPrecSp(AbstractOperatorType& op,
+    initOpPrecSp(Operator& op,
                  const Opm::PropertyTree& prm,
                  const std::function<VectorType()> weightsCalculator,
                  const Dune::Amg::SequentialInformation&,
                  std::size_t pressureIndex)
     {
         // Sequential case.
-        using SeqOperatorType = Dune::MatrixAdapter<MatrixType, VectorType, VectorType>;
         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(*op_prec,
+        preconditioner_ = Opm::PreconditionerFactory<Operator>::create(op,
-                                                                              child ? *child : Opm::PropertyTree(),
+                                                                       child ? *child : Opm::PropertyTree(),
-                                                                              weightsCalculator,
+                                                                       weightsCalculator,
-                                                                              pressureIndex);
+                                                                       pressureIndex);
         scalarproduct_ = std::make_shared<Dune::SeqScalarProduct<VectorType>>();
-        linearoperator_for_precond_ = op_prec;
     }
 
-    template <class MatrixType, class VectorType>
+    template <class Operator>
     void
-    FlexibleSolver<MatrixType, VectorType>::
+    FlexibleSolver<Operator>::
     initSolver(const Opm::PropertyTree& prm, const bool is_iorank)
     {
         const double tol = prm.get<double>("tol", 1e-2);
@@ -175,6 +169,7 @@ namespace Dune
 #if HAVE_SUITESPARSE_UMFPACK
         } else if (solver_type == "umfpack") {
             bool dummy = false;
+            using MatrixType = std::remove_const_t<std::remove_reference_t<decltype(linearoperator_for_solver_->getmat())>>;
             linsolver_.reset(new Dune::UMFPack<MatrixType>(linearoperator_for_solver_->getmat(), verbosity, dummy));
 #endif
         } else {
@@ -185,11 +180,11 @@ namespace Dune
 
     // Main initialization routine.
     // Call with Comm == Dune::Amg::SequentialInformation to get a serial solver.
-    template <class MatrixType, class VectorType>
+    template <class Operator>
     template <class Comm>
     void
-    FlexibleSolver<MatrixType, VectorType>::
+    FlexibleSolver<Operator>::
-    init(AbstractOperatorType& op,
+    init(Operator& op,
          const Comm& comm,
          const Opm::PropertyTree& prm,
          const std::function<VectorType()> weightsCalculator,
@@ -204,42 +199,55 @@ namespace Dune
 
 // Macros to simplify explicit instantiation of FlexibleSolver for various block sizes.
 
+// Vectors and matrices.
 template <int N>
 using BV = Dune::BlockVector<Dune::FieldVector<double, N>>;
 template <int N>
-using BM = Dune::BCRSMatrix<Dune::FieldMatrix<double, N, N>>;
-template <int N>
 using OBM = Dune::BCRSMatrix<Opm::MatrixBlock<double, N, N>>;
 
+// Sequential operators.
+template <int N>
+using SeqOpM = Dune::MatrixAdapter<OBM<N>, BV<N>, BV<N>>;
+template <int N>
+using SeqOpW = Opm::WellModelMatrixAdapter<OBM<N>, BV<N>, BV<N>, false>;
+
 #if HAVE_MPI
 
+// Parallel communicator and operators.
 using Comm = Dune::OwnerOverlapCopyCommunication<int, int>;
+template <int N>
+using ParOpM = Dune::OverlappingSchwarzOperator<OBM<N>, BV<N>, BV<N>, Comm>;
+template <int N>
+using ParOpW = Opm::WellModelGhostLastMatrixAdapter<OBM<N>, BV<N>, BV<N>, true>;
 
 // 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)                  \
+#define INSTANTIATE_FLEXIBLESOLVER_OP(Operator)                                                          \
-template class Dune::FlexibleSolver<BM<N>, BV<N>>;     \
+template class Dune::FlexibleSolver<Operator>;                                                           \
-template class Dune::FlexibleSolver<OBM<N>, BV<N>>;    \
+template Dune::FlexibleSolver<Operator>::FlexibleSolver(Operator& op,                                    \
-template Dune::FlexibleSolver<BM<N>, BV<N>>::FlexibleSolver(AbstractOperatorType& op,                         \
+                                                        const Comm& comm,                                \
-                                                            const Comm& comm,                                 \
+                                                        const Opm::PropertyTree& prm,                    \
-                                                            const Opm::PropertyTree& prm,                     \
+                                                        const std::function<typename Operator::domain_type()>& weightsCalculator, \
-                                                            const std::function<BV<N>()>& weightsCalculator, \
+                                                        std::size_t pressureIndex);
-std::size_t);                                                                 \
+#define INSTANTIATE_FLEXIBLESOLVER(N)     \
-template Dune::FlexibleSolver<OBM<N>, BV<N>>::FlexibleSolver(AbstractOperatorType& op,                        \
+INSTANTIATE_FLEXIBLESOLVER_OP(SeqOpM<N>); \
-                                                             const Comm& comm,                                \
+INSTANTIATE_FLEXIBLESOLVER_OP(SeqOpW<N>); \
-                                                             const Opm::PropertyTree& prm,                    \
+INSTANTIATE_FLEXIBLESOLVER_OP(ParOpM<N>); \
-                                                             const std::function<BV<N>()>& weightsCalculator, \
+INSTANTIATE_FLEXIBLESOLVER_OP(ParOpW<N>);
-std::size_t);
 
 #else // HAVE_MPI
 
-#define INSTANTIATE_FLEXIBLESOLVER(N)                  \
+#define INSTANTIATE_FLEXIBLESOLVER_OP(Operator) \
-template class Dune::FlexibleSolver<BM<N>, BV<N>>;     \
+template class Dune::FlexibleSolver<Operator>;
-template class Dune::FlexibleSolver<OBM<N>, BV<N>>;
+
+#define INSTANTIATE_FLEXIBLESOLVER(N)     \
+INSTANTIATE_FLEXIBLESOLVER_OP(SeqOpM<N>); \
+INSTANTIATE_FLEXIBLESOLVER_OP(SeqOpW<N>);
 
 #endif // HAVE_MPI
 
 
+
 #endif // OPM_FLEXIBLE_SOLVER_IMPL_HEADER_INCLUDED

opm/simulators/linalg/FlowLinearSolverParameters.hpp

@@ -114,6 +114,10 @@ struct CprReuseSetup {
     using type = UndefinedProperty;
 };
 template<class TypeTag, class MyTypeTag>
+struct CprReuseInterval {
+    using type = UndefinedProperty;
+};
+template<class TypeTag, class MyTypeTag>
 struct Linsolver {
     using type = UndefinedProperty;
 };
@@ -137,7 +141,6 @@ template<class TypeTag, class MyTypeTag>
 struct FpgaBitstream {
     using type = UndefinedProperty;
 };
-
 template<class TypeTag>
 struct LinearSolverReduction<TypeTag, TTag::FlowIstlSolverParams> {
     using type = GetPropType<TypeTag, Scalar>;
@@ -213,6 +216,10 @@ struct CprReuseSetup<TypeTag, TTag::FlowIstlSolverParams> {
     static constexpr int value = 3;
 };
 template<class TypeTag>
+struct CprReuseInterval<TypeTag, TTag::FlowIstlSolverParams> {
+    static constexpr int value = 10;
+};
+template<class TypeTag>
 struct Linsolver<TypeTag, TTag::FlowIstlSolverParams> {
     static constexpr auto value = "ilu0";
 };
@@ -263,8 +270,9 @@ namespace Opm
         std::string accelerator_mode_;
         int bda_device_id_;
         int opencl_platform_id_;
-        int cpr_max_ell_iter_ = 20;
+        int cpr_max_ell_iter_;
-        int cpr_reuse_setup_ = 0;
+        int cpr_reuse_setup_;
+        int cpr_reuse_interval_;
         std::string opencl_ilu_reorder_;
         std::string fpga_bitstream_;
 
@@ -287,6 +295,7 @@ namespace Opm
             scale_linear_system_ = EWOMS_GET_PARAM(TypeTag, bool, ScaleLinearSystem);
             cpr_max_ell_iter_  =  EWOMS_GET_PARAM(TypeTag, int, CprMaxEllIter);
             cpr_reuse_setup_  =  EWOMS_GET_PARAM(TypeTag, int, CprReuseSetup);
+            cpr_reuse_interval_  =  EWOMS_GET_PARAM(TypeTag, int, CprReuseInterval);
             linsolver_ = EWOMS_GET_PARAM(TypeTag, std::string, Linsolver);
             accelerator_mode_ = EWOMS_GET_PARAM(TypeTag, std::string, AcceleratorMode);
             bda_device_id_ = EWOMS_GET_PARAM(TypeTag, int, BdaDeviceId);
@@ -312,7 +321,8 @@ namespace Opm
             EWOMS_REGISTER_PARAM(TypeTag, bool, LinearSolverIgnoreConvergenceFailure, "Continue with the simulation like nothing happened after the linear solver did not converge");
             EWOMS_REGISTER_PARAM(TypeTag, bool, ScaleLinearSystem, "Scale linear system according to equation scale and primary variable types");
             EWOMS_REGISTER_PARAM(TypeTag, int, CprMaxEllIter, "MaxIterations of the elliptic pressure part of the cpr solver");
-            EWOMS_REGISTER_PARAM(TypeTag, int, CprReuseSetup, "Reuse preconditioner setup. Valid options are 0: recreate the preconditioner for every linear solve, 1: recreate once every timestep, 2: recreate if last linear solve took more than 10 iterations, 3: never recreate");
+            EWOMS_REGISTER_PARAM(TypeTag, int, CprReuseSetup, "Reuse preconditioner setup. Valid options are 0: recreate the preconditioner for every linear solve, 1: recreate once every timestep, 2: recreate if last linear solve took more than 10 iterations, 3: never recreate, 4: recreated every CprReuseInterval");
+            EWOMS_REGISTER_PARAM(TypeTag, int, CprReuseInterval, "Reuse preconditioner interval. Used when CprReuseSetup is set to 4, then the preconditioner will be fully recreated instead of reused every N linear solve, where N is this parameter.");
             EWOMS_REGISTER_PARAM(TypeTag, std::string, Linsolver, "Configuration of solver. Valid options are: ilu0 (default), cpr (an alias for cpr_trueimpes), cpr_quasiimpes, cpr_trueimpes or amg. Alternatively, you can request a configuration to be read from a JSON file by giving the filename here, ending with '.json.'");
             EWOMS_REGISTER_PARAM(TypeTag, std::string, AcceleratorMode, "Use GPU (cusparseSolver or openclSolver) or FPGA (fpgaSolver) as the linear solver, usage: '--accelerator-mode=[none|cusparse|opencl|fpga|amgcl]'");
             EWOMS_REGISTER_PARAM(TypeTag, int, BdaDeviceId, "Choose device ID for cusparseSolver or openclSolver, use 'nvidia-smi' or 'clinfo' to determine valid IDs");

opm/simulators/linalg/ISTLSolverEbos.hpp

@@ -87,8 +87,9 @@ namespace Opm
         using Matrix = typename SparseMatrixAdapter::IstlMatrix;
         using ThreadManager = GetPropType<TypeTag, Properties::ThreadManager>;
         using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
-        using FlexibleSolverType = Dune::FlexibleSolver<Matrix, Vector>;
+        using AbstractSolverType = Dune::InverseOperator<Vector, Vector>;
         using AbstractOperatorType = Dune::AssembledLinearOperator<Matrix, Vector, Vector>;
+        using AbstractPreconditionerType = Dune::PreconditionerWithUpdate<Vector, Vector>;
         using WellModelOperator = WellModelAsLinearOperator<WellModel, Vector, Vector>;
         using ElementMapper = GetPropType<TypeTag, Properties::ElementMapper>;
         constexpr static std::size_t pressureIndex = GetPropType<TypeTag, Properties::Indices>::pressureSwitchIdx;
@@ -118,6 +119,7 @@ namespace Opm
         explicit ISTLSolverEbos(const Simulator& simulator)
             : simulator_(simulator),
               iterations_( 0 ),
+              calls_( 0 ),
               converged_(false),
               matrix_()
         {
@@ -257,6 +259,7 @@ namespace Opm
         }
 
         bool solve(Vector& x) {
+            calls_ += 1;
             // Write linear system if asked for.
             const int verbosity = prm_.get<int>("verbosity", 0);
             const bool write_matrix = verbosity > 10;
@@ -386,35 +389,47 @@ namespace Opm
 #if HAVE_MPI
                     if (useWellConn_) {
                         using ParOperatorType = Dune::OverlappingSchwarzOperator<Matrix, Vector, Vector, Comm>;
-                        linearOperatorForFlexibleSolver_ = std::make_unique<ParOperatorType>(getMatrix(), *comm_);
+                        auto op = std::make_unique<ParOperatorType>(getMatrix(), *comm_);
-                        flexibleSolver_ = std::make_unique<FlexibleSolverType>(*linearOperatorForFlexibleSolver_, *comm_, prm_, weightsCalculator,
+                        using FlexibleSolverType = Dune::FlexibleSolver<ParOperatorType>;
-                                                                               pressureIndex);
+                        auto sol = std::make_unique<FlexibleSolverType>(*op, *comm_, prm_, weightsCalculator, pressureIndex);
+                        preconditionerForFlexibleSolver_ = &(sol->preconditioner());
+                        linearOperatorForFlexibleSolver_ = std::move(op);
+                        flexibleSolver_ = std::move(sol);
                     } else {
                         using ParOperatorType = WellModelGhostLastMatrixAdapter<Matrix, Vector, Vector, true>;
                         wellOperator_ = std::make_unique<WellModelOperator>(simulator_.problem().wellModel());
-                        linearOperatorForFlexibleSolver_ = std::make_unique<ParOperatorType>(getMatrix(), *wellOperator_, interiorCellNum_);
+                        auto op = std::make_unique<ParOperatorType>(getMatrix(), *wellOperator_, interiorCellNum_);
-                        flexibleSolver_ = std::make_unique<FlexibleSolverType>(*linearOperatorForFlexibleSolver_, *comm_, prm_, weightsCalculator,
+                        using FlexibleSolverType = Dune::FlexibleSolver<ParOperatorType>;
-                                                                               pressureIndex);
+                        auto sol = std::make_unique<FlexibleSolverType>(*op, *comm_, prm_, weightsCalculator, pressureIndex);
+                        preconditionerForFlexibleSolver_ = &(sol->preconditioner());
+                        linearOperatorForFlexibleSolver_ = std::move(op);
+                        flexibleSolver_ = std::move(sol);
                     }
 #endif
                 } else {
                     if (useWellConn_) {
                         using SeqOperatorType = Dune::MatrixAdapter<Matrix, Vector, Vector>;
-                        linearOperatorForFlexibleSolver_ = std::make_unique<SeqOperatorType>(getMatrix());
+                        auto op = std::make_unique<SeqOperatorType>(getMatrix());
-                        flexibleSolver_ = std::make_unique<FlexibleSolverType>(*linearOperatorForFlexibleSolver_, prm_, weightsCalculator,
+                        using FlexibleSolverType = Dune::FlexibleSolver<SeqOperatorType>;
-                                                                               pressureIndex);
+                        auto sol = std::make_unique<FlexibleSolverType>(*op, prm_, weightsCalculator, pressureIndex);
+                        preconditionerForFlexibleSolver_ = &(sol->preconditioner());
+                        linearOperatorForFlexibleSolver_ = std::move(op);
+                        flexibleSolver_ = std::move(sol);
                     } else {
                         using SeqOperatorType = WellModelMatrixAdapter<Matrix, Vector, Vector, false>;
                         wellOperator_ = std::make_unique<WellModelOperator>(simulator_.problem().wellModel());
-                        linearOperatorForFlexibleSolver_ = std::make_unique<SeqOperatorType>(getMatrix(), *wellOperator_);
+                        auto op = std::make_unique<SeqOperatorType>(getMatrix(), *wellOperator_);
-                        flexibleSolver_ = std::make_unique<FlexibleSolverType>(*linearOperatorForFlexibleSolver_, prm_, weightsCalculator,
+                        using FlexibleSolverType = Dune::FlexibleSolver<SeqOperatorType>;
-                                                                               pressureIndex);
+                        auto sol = std::make_unique<FlexibleSolverType>(*op, prm_, weightsCalculator, pressureIndex);
+                        preconditionerForFlexibleSolver_ = &(sol->preconditioner());
+                        linearOperatorForFlexibleSolver_ = std::move(op);
+                        flexibleSolver_ = std::move(sol);
                     }
                 }
             }
             else
             {
-                flexibleSolver_->preconditioner().update();
+                preconditionerForFlexibleSolver_->update();
             }
         }
 
@@ -441,10 +456,27 @@ namespace Opm
                 // Recreate solver if the last solve used more than 10 iterations.
                 return this->iterations() > 10;
             }
+            if (this->parameters_.cpr_reuse_setup_ == 3) {
+                // Recreate solver if the last solve used more than 10 iterations.
+                return false;
+            }
+            if (this->parameters_.cpr_reuse_setup_ == 4) {
+                // Recreate solver every 'step' solve calls.
+                const int step = this->parameters_.cpr_reuse_interval_;
+                const bool create = ((calls_ % step) == 0);
+                return create;
+            }
 
-            // Otherwise, do not recreate solver.
+            // If here, we have an invalid parameter.
-            assert(this->parameters_.cpr_reuse_setup_ == 3);
+            const bool on_io_rank = (simulator_.gridView().comm().rank() == 0);
+            std::string msg = "Invalid value: " + std::to_string(this->parameters_.cpr_reuse_setup_)
+                + " for --cpr-reuse-setup parameter, run with --help to see allowed values.";
+            if (on_io_rank) {
+                OpmLog::error(msg);
+            }
+            throw std::runtime_error(msg);
 
+            // Never reached.
             return false;
         }
 
@@ -457,8 +489,9 @@ namespace Opm
             using namespace std::string_literals;
 
             auto preconditionerType = prm_.get("preconditioner.type"s, "cpr"s);
-            if (preconditionerType == "cpr" || preconditionerType == "cprt") {
+            if (preconditionerType == "cpr" || preconditionerType == "cprt"
-                const bool transpose = preconditionerType == "cprt";
+                || preconditionerType == "cprw" || preconditionerType == "cprwt") {
+                const bool transpose = preconditionerType == "cprt" || preconditionerType == "cprwt";
                 const auto weightsType = prm_.get("preconditioner.weight_type"s, "quasiimpes"s);
                 if (weightsType == "quasiimpes") {
                     // weights will be created as default in the solver
@@ -594,6 +627,7 @@ namespace Opm
 
         const Simulator& simulator_;
         mutable int iterations_;
+        mutable int calls_;
         mutable bool converged_;
         std::any parallelInformation_;
 
@@ -603,8 +637,9 @@ namespace Opm
 
         std::unique_ptr<Matrix> blockJacobiForGPUILU0_;
 
-        std::unique_ptr<FlexibleSolverType> flexibleSolver_;
+        std::unique_ptr<AbstractSolverType> flexibleSolver_;
         std::unique_ptr<AbstractOperatorType> linearOperatorForFlexibleSolver_;
+        AbstractPreconditionerType* preconditionerForFlexibleSolver_;
         std::unique_ptr<WellModelAsLinearOperator<WellModel, Vector, Vector>> wellOperator_;
         std::vector<int> overlapRows_;
         std::vector<int> interiorRows_;

opm/simulators/linalg/MatrixMarketSpecializations.hpp

@@ -44,19 +44,20 @@ namespace MatrixMarketImpl
         }
     };
 
+    template<class M>
+    struct mm_multipliers;
+
+    template <typename T, int i, int j, typename A>
+    struct mm_multipliers<BCRSMatrix<Opm::MatrixBlock<T,i,j>, A>>
+    {
+        enum {
+            rows = i,
+            cols = j
+        };
+    };
+
 } // namespace MatrixMarketImpl
 
-template<class M>
-struct mm_multipliers;
-
-template <typename T, int i, int j, typename A>
-struct mm_multipliers<BCRSMatrix<Opm::MatrixBlock<T,i,j>, A>>
-{
-    enum {
-        rows = i,
-        cols = j
-    };
-};
 
 } // namespace Dune
 

opm/simulators/linalg/OwningTwoLevelPreconditioner.hpp

@@ -22,8 +22,6 @@
 
 #include <opm/simulators/linalg/PreconditionerWithUpdate.hpp>
 #include <opm/simulators/linalg/PressureSolverPolicy.hpp>
-#include <opm/simulators/linalg/PressureTransferPolicy.hpp>
-#include <opm/simulators/linalg/getQuasiImpesWeights.hpp>
 #include <opm/simulators/linalg/twolevelmethodcpr.hh>
 
 #include <opm/common/ErrorMacros.hpp>
@@ -50,7 +48,7 @@ namespace Dune
 
 
 // Must forward-declare FlexibleSolver as we want to use it as solver for the pressure system.
-template <class MatrixTypeT, class VectorTypeT>
+template <class Operator>
 class FlexibleSolver;
 
 template <typename T, typename A, int i>
@@ -75,13 +73,14 @@ std::ostream& operator<<(std::ostream& out,
 ///   and preconditioner factory.
 template <class OperatorType,
           class VectorType,
-          bool transpose = false,
+          class LevelTransferPolicy,
           class Communication = Dune::Amg::SequentialInformation>
 class OwningTwoLevelPreconditioner : public Dune::PreconditionerWithUpdate<VectorType, VectorType>
 {
 public:
     using MatrixType = typename OperatorType::matrix_type;
     using PrecFactory = Opm::PreconditionerFactory<OperatorType, Communication>;
+    using AbstractOperatorType = Dune::AssembledLinearOperator<MatrixType, VectorType, VectorType>;
 
     OwningTwoLevelPreconditioner(const OperatorType& linearoperator, const Opm::PropertyTree& prm,
                                  const std::function<VectorType()> weightsCalculator,
@@ -94,14 +93,14 @@ public:
         , comm_(nullptr)
         , weightsCalculator_(weightsCalculator)
         , weights_(weightsCalculator())
-        , levelTransferPolicy_(dummy_comm_, weights_, pressureIndex)
+        , levelTransferPolicy_(dummy_comm_, weights_, prm, pressureIndex)
-        , coarseSolverPolicy_(prm.get_child_optional("coarsesolver")? prm.get_child("coarsesolver") : Opm::PropertyTree())
+        , coarseSolverPolicy_(prm.get_child_optional("coarsesolver") ? prm.get_child("coarsesolver") : Opm::PropertyTree())
         , twolevel_method_(linearoperator,
                            finesmoother_,
                            levelTransferPolicy_,
                            coarseSolverPolicy_,
-                           prm.get<int>("pre_smooth", transpose? 1 : 0),
+                           prm.get<int>("pre_smooth", 0),
-                           prm.get<int>("post_smooth", transpose? 0 : 1))
+                           prm.get<int>("post_smooth", 1))
         , prm_(prm)
     {
         if (prm.get<int>("verbosity", 0) > 10) {
@@ -126,14 +125,14 @@ public:
         , comm_(&comm)
         , weightsCalculator_(weightsCalculator)
         , weights_(weightsCalculator())
-        , levelTransferPolicy_(*comm_, weights_, pressureIndex)
+        , levelTransferPolicy_(*comm_, weights_, prm, pressureIndex)
-        , coarseSolverPolicy_(prm.get_child_optional("coarsesolver")? prm.get_child("coarsesolver") : Opm::PropertyTree())
+        , coarseSolverPolicy_(prm.get_child_optional("coarsesolver") ? prm.get_child("coarsesolver") : Opm::PropertyTree())
         , twolevel_method_(linearoperator,
                            finesmoother_,
                            levelTransferPolicy_,
                            coarseSolverPolicy_,
-                           prm.get<int>("pre_smooth", transpose? 1 : 0),
+                           prm.get<int>("pre_smooth", 0),
-                           prm.get<int>("post_smooth", transpose? 0 : 1))
+                           prm.get<int>("post_smooth", 1))
         , prm_(prm)
     {
         if (prm.get<int>("verbosity", 0) > 10 && comm.communicator().rank() == 0) {
@@ -173,17 +172,9 @@ public:
     }
 
 private:
-    using PressureMatrixType = Dune::BCRSMatrix<Dune::FieldMatrix<double, 1, 1>>;
+    using CoarseOperator = typename LevelTransferPolicy::CoarseOperator;
-    using PressureVectorType = Dune::BlockVector<Dune::FieldVector<double, 1>>;
+    using CoarseSolverPolicy = Dune::Amg::PressureSolverPolicy<CoarseOperator,
-    using SeqCoarseOperatorType = Dune::MatrixAdapter<PressureMatrixType, PressureVectorType, PressureVectorType>;
+                                                               FlexibleSolver<CoarseOperator>,
-    using ParCoarseOperatorType
-        = Dune::OverlappingSchwarzOperator<PressureMatrixType, PressureVectorType, PressureVectorType, Communication>;
-    using CoarseOperatorType = std::conditional_t<std::is_same<Communication, Dune::Amg::SequentialInformation>::value,
-                                                  SeqCoarseOperatorType,
-                                                  ParCoarseOperatorType>;
-    using LevelTransferPolicy = Opm::PressureTransferPolicy<OperatorType, CoarseOperatorType, Communication, transpose>;
-    using CoarseSolverPolicy = Dune::Amg::PressureSolverPolicy<CoarseOperatorType,
-                                                               FlexibleSolver<PressureMatrixType, PressureVectorType>,
                                                                LevelTransferPolicy>;
     using TwoLevelMethod
         = Dune::Amg::TwoLevelMethodCpr<OperatorType, CoarseSolverPolicy, Dune::Preconditioner<VectorType, VectorType>>;

opm/simulators/linalg/PreconditionerFactory.hpp

@@ -28,6 +28,11 @@
 #include <opm/simulators/linalg/PreconditionerWithUpdate.hpp>
 #include <opm/simulators/linalg/PropertyTree.hpp>
 #include <opm/simulators/linalg/amgcpr.hh>
+#include <opm/simulators/linalg/WellOperators.hpp>
+
+#include <opm/simulators/linalg/PressureTransferPolicy.hpp>
+#include <opm/simulators/linalg/PressureBhpTransferPolicy.hpp>
+
 
 #include <dune/istl/paamg/amg.hh>
 #include <dune/istl/paamg/kamg.hh>
@@ -318,7 +323,8 @@ private:
             {
                 OPM_THROW(std::logic_error, "Pressure index out of bounds. It needs to specified for CPR");
             }
-            return std::make_shared<OwningTwoLevelPreconditioner<O, V, false, Comm>>(op, prm, weightsCalculator, pressureIndex, comm);
+            using LevelTransferPolicy = Opm::PressureTransferPolicy<O, Comm, false>;
+            return std::make_shared<OwningTwoLevelPreconditioner<O, V, LevelTransferPolicy, Comm>>(op, prm, weightsCalculator, pressureIndex, comm);
         });
         doAddCreator("cprt", [](const O& op, const P& prm, const std::function<Vector()> weightsCalculator, std::size_t pressureIndex, const C& comm) {
             assert(weightsCalculator);
@@ -326,8 +332,22 @@ private:
             {
                 OPM_THROW(std::logic_error, "Pressure index out of bounds. It needs to specified for CPR");
             }
-            return std::make_shared<OwningTwoLevelPreconditioner<O, V, true, Comm>>(op, prm, weightsCalculator, pressureIndex, comm);
+            using LevelTransferPolicy = Opm::PressureTransferPolicy<O, Comm, true>;
+            return std::make_shared<OwningTwoLevelPreconditioner<O, V, LevelTransferPolicy, Comm>>(op, prm, weightsCalculator, pressureIndex, comm);
         });
+
+        if constexpr (std::is_same_v<O, WellModelGhostLastMatrixAdapter<M, V, V, true>>) {
+            doAddCreator("cprw",
+                         [](const O& op, const P& prm, const std::function<Vector()> weightsCalculator, std::size_t pressureIndex, const C& comm) {
+                             assert(weightsCalculator);
+                             if (pressureIndex == std::numeric_limits<std::size_t>::max()) {
+                                 OPM_THROW(std::logic_error, "Pressure index out of bounds. It needs to specified for CPR");
+                             }
+                             using LevelTransferPolicy = Opm::PressureBhpTransferPolicy<O, Comm, false>;
+                             return std::make_shared<OwningTwoLevelPreconditioner<O, V, LevelTransferPolicy, Comm>>(
+                                 op, prm, weightsCalculator, pressureIndex, comm);
+                         });
+        }
     }
 
     // Add a useful default set of preconditioners to the factory.
@@ -449,19 +469,31 @@ private:
                 return wrapPreconditioner<Dune::Amg::FastAMG<O, V>>(op, crit, parms);
             });
         }
+        if constexpr (std::is_same_v<O, WellModelMatrixAdapter<M, V, V, false>>) {
+            doAddCreator("cprw", [](const O& op, const P& prm, const std::function<Vector()>& weightsCalculator, std::size_t pressureIndex) {
+                if (pressureIndex == std::numeric_limits<std::size_t>::max()) {
+                    OPM_THROW(std::logic_error, "Pressure index out of bounds. It needs to specified for CPR");
+                }
+                using LevelTransferPolicy = Opm::PressureBhpTransferPolicy<O, Dune::Amg::SequentialInformation, false>;
+                return std::make_shared<OwningTwoLevelPreconditioner<O, V, LevelTransferPolicy>>(op, prm, weightsCalculator, pressureIndex);
+            });
+            }
+
         doAddCreator("cpr", [](const O& op, const P& prm, const std::function<Vector()>& weightsCalculator, std::size_t pressureIndex) {
                                 if (pressureIndex == std::numeric_limits<std::size_t>::max())
                                 {
                                     OPM_THROW(std::logic_error, "Pressure index out of bounds. It needs to specified for CPR");
                                 }
-                                return std::make_shared<OwningTwoLevelPreconditioner<O, V, false>>(op, prm, weightsCalculator, pressureIndex);
+                                using LevelTransferPolicy = Opm::PressureTransferPolicy<O, Dune::Amg::SequentialInformation, false>;
+                                return std::make_shared<OwningTwoLevelPreconditioner<O, V, LevelTransferPolicy>>(op, prm, weightsCalculator, pressureIndex);
         });
         doAddCreator("cprt", [](const O& op, const P& prm, const std::function<Vector()>& weightsCalculator, std::size_t pressureIndex) {
                                 if (pressureIndex == std::numeric_limits<std::size_t>::max())
                                 {
                                     OPM_THROW(std::logic_error, "Pressure index out of bounds. It needs to specified for CPR");
                                 }
-                                 return std::make_shared<OwningTwoLevelPreconditioner<O, V, true>>(op, prm, weightsCalculator, pressureIndex);
+                                using LevelTransferPolicy = Opm::PressureTransferPolicy<O, Dune::Amg::SequentialInformation, true>;
+                                return std::make_shared<OwningTwoLevelPreconditioner<O, V, LevelTransferPolicy>>(op, prm, weightsCalculator, pressureIndex);
         });
     }
 

opm/simulators/linalg/PressureBhpTransferPolicy.hpp

@@ -0,0 +1,270 @@
+/*
+  Copyright 2019 SINTEF Digital, Mathematics and Cybernetics.
+  Copyright 2019 Dr. Blatt - HPC-Simulation-Software & Services.
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#pragma once
+
+#include <opm/simulators/linalg/twolevelmethodcpr.hh>
+#include <dune/istl/paamg/pinfo.hh>
+
+namespace Opm
+{
+    template <class Communication>
+    void extendCommunicatorWithWells(const Communication& comm,
+                                     std::shared_ptr<Communication>& commRW,
+                                     const int nw)
+    {
+        // used for extending the coarse communicator pattern
+        using IndexSet = typename Communication::ParallelIndexSet;
+        using LocalIndex = typename IndexSet::LocalIndex;
+        const IndexSet& indset = comm.indexSet();
+        IndexSet& indset_rw = commRW->indexSet();
+        const int max_nw = comm.communicator().max(nw);
+        const int rank = comm.communicator().rank();
+        int glo_max = 0;
+        size_t loc_max = 0;
+        indset_rw.beginResize();
+        for (auto ind = indset.begin(), indend = indset.end(); ind != indend; ++ind) {
+            indset_rw.add(ind->global(), LocalIndex(ind->local(), ind->local().attribute(), true));
+            const int glo = ind->global();
+            const size_t loc = ind->local().local();
+            glo_max = std::max(glo_max, glo);
+            loc_max = std::max(loc_max, loc);
+        }
+        const int global_max = comm.communicator().max(glo_max);
+        // used append the welldofs at the end
+        assert(loc_max + 1 == indset.size());
+        size_t local_ind = loc_max + 1;
+        for (int i = 0; i < nw; ++i) {
+            // need to set unique global number
+            const size_t v = global_max + max_nw * rank + i + 1;
+            // set to true to not have problems with higher levels if growing of domains is used
+            indset_rw.add(v, LocalIndex(local_ind, Dune::OwnerOverlapCopyAttributeSet::owner, true));
+            ++local_ind;
+        }
+        indset_rw.endResize();
+        assert(indset_rw.size() == indset.size() + nw);
+        // assume same communication pattern
+        commRW->remoteIndices().setNeighbours(comm.remoteIndices().getNeighbours());
+        commRW->remoteIndices().template rebuild<true>();
+        //commRW->clearInterfaces(); may need this for correct rebuild
+    }
+
+    namespace Details
+    {
+        using PressureMatrixType = Dune::BCRSMatrix<Opm::MatrixBlock<double, 1, 1>>;
+        using PressureVectorType = Dune::BlockVector<Dune::FieldVector<double, 1>>;
+        using SeqCoarseOperatorType = Dune::MatrixAdapter<PressureMatrixType, PressureVectorType, PressureVectorType>;
+        template <class Comm>
+        using ParCoarseOperatorType
+            = Dune::OverlappingSchwarzOperator<PressureMatrixType, PressureVectorType, PressureVectorType, Comm>;
+        template <class Comm>
+        using CoarseOperatorType = std::conditional_t<std::is_same<Comm, Dune::Amg::SequentialInformation>::value,
+                                                      SeqCoarseOperatorType,
+                                                      ParCoarseOperatorType<Comm>>;
+    } // namespace Details
+
+    template <class FineOperator, class Communication, bool transpose = false>
+    class PressureBhpTransferPolicy : public Dune::Amg::LevelTransferPolicyCpr<FineOperator, Details::CoarseOperatorType<Communication>>
+    {
+    public:
+        typedef typename Details::CoarseOperatorType<Communication> CoarseOperator;
+        typedef Dune::Amg::LevelTransferPolicyCpr<FineOperator, CoarseOperator> ParentType;
+        typedef Communication ParallelInformation;
+        typedef typename FineOperator::domain_type FineVectorType;
+
+    public:
+        PressureBhpTransferPolicy(const Communication& comm,
+                                  const FineVectorType& weights,
+                                  const Opm::PropertyTree& prm,
+                                  const std::size_t pressureIndex)
+            : communication_(&const_cast<Communication&>(comm))
+            , weights_(weights)
+            , prm_(prm)
+            , pressure_var_index_(pressureIndex)
+        {
+        }
+
+        virtual void createCoarseLevelSystem(const FineOperator& fineOperator) override
+        {
+            using CoarseMatrix = typename CoarseOperator::matrix_type;
+            const auto& fineLevelMatrix = fineOperator.getmat();
+            const auto& nw = fineOperator.getNumberOfExtraEquations();
+            if (prm_.get<bool>("add_wells")) {
+                const size_t average_elements_per_row
+                    = static_cast<size_t>(std::ceil(fineLevelMatrix.nonzeroes() / fineLevelMatrix.N()));
+                const double overflow_fraction = 1.2;
+                coarseLevelMatrix_.reset(new CoarseMatrix(fineLevelMatrix.N() + nw,
+                                                          fineLevelMatrix.M() + nw,
+                                                          average_elements_per_row,
+                                                          overflow_fraction,
+                                                          CoarseMatrix::implicit));
+                int rownum = 0;
+                for (const auto& row : fineLevelMatrix) {
+                    for (auto col = row.begin(), cend = row.end(); col != cend; ++col) {
+                        coarseLevelMatrix_->entry(rownum, col.index()) = 0.0;
+                    }
+                    ++rownum;
+                }
+            } else {
+                coarseLevelMatrix_.reset(
+                    new CoarseMatrix(fineLevelMatrix.N(), fineLevelMatrix.M(), CoarseMatrix::row_wise));
+                auto createIter = coarseLevelMatrix_->createbegin();
+                for (const auto& row : fineLevelMatrix) {
+                    for (auto col = row.begin(), cend = row.end(); col != cend; ++col) {
+                        createIter.insert(col.index());
+                    }
+                    ++createIter;
+                }
+            }
+#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 7)
+        if constexpr (std::is_same_v<Communication, Dune::Amg::SequentialInformation>) {
+            coarseLevelCommunication_ = std::make_shared<Communication>();
+        } else {
+            coarseLevelCommunication_ = std::make_shared<Communication>(
+                communication_->communicator(), communication_->category(), false);
+        }
+#else
+        if constexpr (std::is_same_v<Communication, Dune::Amg::SequentialInformation>) {
+            coarseLevelCommunication_ = std::make_shared<Communication>();
+        } else {
+            coarseLevelCommunication_ = std::make_shared<Communication>(
+                communication_->communicator(), communication_->getSolverCategory(), false);
+        }
+#endif
+        if (prm_.get<bool>("add_wells")) {
+            fineOperator.addWellPressureEquationsStruct(*coarseLevelMatrix_);
+            coarseLevelMatrix_->compress(); // all elemenst should be set
+            if constexpr (!std::is_same_v<Communication, Dune::Amg::SequentialInformation>) {
+                extendCommunicatorWithWells(*communication_, coarseLevelCommunication_, nw);
+            }
+        }
+        calculateCoarseEntries(fineOperator);
+
+        this->lhs_.resize(this->coarseLevelMatrix_->M());
+        this->rhs_.resize(this->coarseLevelMatrix_->N());
+        using OperatorArgs = typename Dune::Amg::ConstructionTraits<CoarseOperator>::Arguments;
+#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 7)
+        OperatorArgs oargs(coarseLevelMatrix_, *coarseLevelCommunication_);
+        this->operator_ = Dune::Amg::ConstructionTraits<CoarseOperator>::construct(oargs);
+#else
+        OperatorArgs oargs(*coarseLevelMatrix_, *coarseLevelCommunication_);
+        this->operator_.reset(Dune::Amg::ConstructionTraits<CoarseOperator>::construct(oargs));
+#endif
+    }
+
+    virtual void calculateCoarseEntries(const FineOperator& fineOperator) override
+    {
+        const auto& fineMatrix = fineOperator.getmat();
+        *coarseLevelMatrix_ = 0;
+        auto rowCoarse = coarseLevelMatrix_->begin();
+        for (auto row = fineMatrix.begin(), rowEnd = fineMatrix.end(); row != rowEnd; ++row, ++rowCoarse) {
+            assert(row.index() == rowCoarse.index());
+            auto entryCoarse = rowCoarse->begin();
+            for (auto entry = row->begin(), entryEnd = row->end(); entry != entryEnd; ++entry, ++entryCoarse) {
+                assert(entry.index() == entryCoarse.index());
+                double matrix_el = 0;
+                if (transpose) {
+                    const auto& bw = weights_[entry.index()];
+                    for (size_t i = 0; i < bw.size(); ++i) {
+                        matrix_el += (*entry)[pressure_var_index_][i] * bw[i];
+                    }
+                } else {
+                    const auto& bw = weights_[row.index()];
+                    for (size_t i = 0; i < bw.size(); ++i) {
+                        matrix_el += (*entry)[i][pressure_var_index_] * bw[i];
+                    }
+                }
+                (*entryCoarse) = matrix_el;
+            }
+        }
+        if (prm_.get<bool>("add_wells")) {
+            assert(transpose == false); // not implemented
+            bool use_well_weights = prm_.get<bool>("use_well_weights");
+            fineOperator.addWellPressureEquations(*coarseLevelMatrix_, weights_, use_well_weights);
+#ifndef NDEBUG
+            std::advance(rowCoarse, fineOperator.getNumberOfExtraEquations());
+            assert(rowCoarse == coarseLevelMatrix_->end());
+#endif
+
+        }
+    }
+
+    virtual void moveToCoarseLevel(const typename ParentType::FineRangeType& fine) override
+    {
+        //NB we iterate over fine assumming welldofs is at the end
+        // Set coarse vector to zero
+        this->rhs_ = 0;
+
+        auto end = fine.end(), begin = fine.begin();
+
+        for (auto block = begin; block != end; ++block) {
+            const auto& bw = weights_[block.index()];
+            double rhs_el = 0.0;
+            if (transpose) {
+                rhs_el = (*block)[pressure_var_index_];
+            } else {
+                for (size_t i = 0; i < block->size(); ++i) {
+                    rhs_el += (*block)[i] * bw[i];
+                }
+            }
+            this->rhs_[block - begin] = rhs_el;
+        }
+
+        this->lhs_ = 0;
+    }
+
+    virtual void moveToFineLevel(typename ParentType::FineDomainType& fine) override
+    {
+        //NB we iterate over fine assumming welldofs is at the end
+        auto end = fine.end(), begin = fine.begin();
+
+        for (auto block = begin; block != end; ++block) {
+            if (transpose) {
+                const auto& bw = weights_[block.index()];
+                for (size_t i = 0; i < block->size(); ++i) {
+                    (*block)[i] = this->lhs_[block - begin] * bw[i];
+                }
+            } else {
+                (*block)[pressure_var_index_] = this->lhs_[block - begin];
+            }
+        }
+    }
+
+    virtual PressureBhpTransferPolicy* clone() const override
+    {
+        return new PressureBhpTransferPolicy(*this);
+    }
+
+    const Communication& getCoarseLevelCommunication() const
+    {
+        return *coarseLevelCommunication_;
+    }
+
+private:
+    Communication* communication_;
+    const FineVectorType& weights_;
+    PropertyTree prm_;
+    const int pressure_var_index_;
+    std::shared_ptr<Communication> coarseLevelCommunication_;
+    std::shared_ptr<typename CoarseOperator::matrix_type> coarseLevelMatrix_;
+};
+
+} // namespace Opm
+

opm/simulators/linalg/PressureTransferPolicy.hpp

@@ -23,21 +23,44 @@
 
 
 #include <opm/simulators/linalg/twolevelmethodcpr.hh>
+#include <opm/simulators/linalg/PropertyTree.hpp>
+#include <opm/simulators/linalg/matrixblock.hh>
 
 
 namespace Opm
 {
 
-template <class FineOperator, class CoarseOperator, class Communication, bool transpose = false>
+namespace Details
-class PressureTransferPolicy : public Dune::Amg::LevelTransferPolicyCpr<FineOperator, CoarseOperator>
+{
+    using PressureMatrixType = Dune::BCRSMatrix<Opm::MatrixBlock<double, 1, 1>>;
+    using PressureVectorType = Dune::BlockVector<Dune::FieldVector<double, 1>>;
+    using SeqCoarseOperatorType = Dune::MatrixAdapter<PressureMatrixType, PressureVectorType, PressureVectorType>;
+    template <class Comm>
+    using ParCoarseOperatorType
+        = Dune::OverlappingSchwarzOperator<PressureMatrixType, PressureVectorType, PressureVectorType, Comm>;
+    template <class Comm>
+    using CoarseOperatorType = std::conditional_t<std::is_same<Comm, Dune::Amg::SequentialInformation>::value,
+                                                  SeqCoarseOperatorType,
+                                                  ParCoarseOperatorType<Comm>>;
+} // namespace Details
+
+
+
+template <class FineOperator, class Communication, bool transpose = false>
+class PressureTransferPolicy
+    : public Dune::Amg::LevelTransferPolicyCpr<FineOperator, Details::CoarseOperatorType<Communication>>
 {
 public:
+    typedef typename Details::CoarseOperatorType<Communication> CoarseOperator;
     typedef Dune::Amg::LevelTransferPolicyCpr<FineOperator, CoarseOperator> ParentType;
     typedef Communication ParallelInformation;
     typedef typename FineOperator::domain_type FineVectorType;
 
 public:
-    PressureTransferPolicy(const Communication& comm, const FineVectorType& weights, int pressure_var_index)
+    PressureTransferPolicy(const Communication& comm,
+                           const FineVectorType& weights,
+                           const Opm::PropertyTree& /*prm*/,
+                           int pressure_var_index)
         : communication_(&const_cast<Communication&>(comm))
         , weights_(weights)
         , pressure_var_index_(pressure_var_index)

opm/simulators/linalg/WellOperators.hpp

@@ -45,30 +45,40 @@ namespace Opm
 /// WellModelMatrixAdapter separately for each TypeTag, as it will only
 /// depend on the matrix and vector types involved, which typically are
 /// just one for each block size with block sizes 1-4.
-template <class WellModel, class X, class Y>
+template <class X, class Y>
-class WellModelAsLinearOperator : public Dune::LinearOperator<X, Y>
+class LinearOperatorExtra : public Dune::LinearOperator<X, Y>
 {
 public:
-    using Base = Dune::LinearOperator<X, Y>;
+    using PressureMatrix = Dune::BCRSMatrix<Opm::MatrixBlock<double, 1, 1>>;
-    using field_type = typename Base::field_type;
+    virtual void addWellPressureEquations(PressureMatrix& jacobian, const X& weights,const bool use_well_weights) const = 0;
+    virtual void addWellPressureEquationsStruct(PressureMatrix& jacobian) const = 0;
+    virtual int getNumberOfExtraEquations() const = 0;
+};
 
+template <class WellModel, class X, class Y>
+class WellModelAsLinearOperator : public Opm::LinearOperatorExtra<X, Y>
+{
+public:
+    using Base = Opm::LinearOperatorExtra<X, Y>;
+    using field_type = typename Base::field_type;
+    using PressureMatrix = typename Base::PressureMatrix;
     explicit WellModelAsLinearOperator(const WellModel& wm)
         : wellMod_(wm)
     {
     }
     /*! \brief apply operator to x:  \f$ y = A(x) \f$
-          The input vector is consistent and the output must also be
+       The input vector is consistent and the output must also be
        consistent on the interior+border partition.
      */
-    void apply (const X& x, Y& y) const override
+    void apply(const X& x, Y& y) const override
     {
-        wellMod_.apply(x, y );
+        wellMod_.apply(x, y);
     }
 
     //! apply operator to x, scale and add:  \f$ y = y + \alpha A(x) \f$
-    virtual void applyscaleadd (field_type alpha, const X& x, Y& y) const override
+    virtual void applyscaleadd(field_type alpha, const X& x, Y& y) const override
     {
-        wellMod_.applyScaleAdd( alpha, x, y );
+        wellMod_.applyScaleAdd(alpha, x, y);
     }
 
     /// Category for operator.
@@ -80,6 +90,19 @@ public:
     {
         return Dune::SolverCategory::sequential;
     }
+    void addWellPressureEquations(PressureMatrix& jacobian, const X& weights,const bool use_well_weights) const override
+    {
+        wellMod_.addWellPressureEquations(jacobian, weights, use_well_weights);
+    }
+    void addWellPressureEquationsStruct(PressureMatrix& jacobian) const override
+    {
+        wellMod_.addWellPressureEquationsStruct(jacobian);
+    }
+    int getNumberOfExtraEquations() const override
+    {
+        return wellMod_.numLocalWellsEnd();
+    }
+
 private:
     const WellModel& wellMod_;
 };
@@ -103,7 +126,7 @@ public:
   typedef X domain_type;
   typedef Y range_type;
   typedef typename X::field_type field_type;
-
+  using PressureMatrix = Dune::BCRSMatrix<Opm::MatrixBlock<double, 1, 1>>;
 #if HAVE_MPI
   typedef Dune::OwnerOverlapCopyCommunication<int,int> communication_type;
 #else
@@ -118,7 +141,7 @@ public:
 
   //! constructor: just store a reference to a matrix
   WellModelMatrixAdapter (const M& A,
-                          const Dune::LinearOperator<X, Y>& wellOper,
+                          const Opm::LinearOperatorExtra<X, Y>& wellOper,
                           const std::shared_ptr< communication_type >& comm = std::shared_ptr< communication_type >())
       : A_( A ), wellOper_( wellOper ), comm_(comm)
   {}
@@ -153,9 +176,22 @@ public:
 
   virtual const matrix_type& getmat() const override { return A_; }
 
+    void addWellPressureEquations(PressureMatrix& jacobian, const X& weights,const bool use_well_weights) const
+    {
+        wellOper_.addWellPressureEquations(jacobian, weights, use_well_weights);
+    }
+    void addWellPressureEquationsStruct(PressureMatrix& jacobian) const
+    {
+        wellOper_.addWellPressureEquationsStruct(jacobian);
+    }
+    int getNumberOfExtraEquations() const
+    {
+        return wellOper_.getNumberOfExtraEquations();
+    }
+
 protected:
   const matrix_type& A_ ;
-  const Dune::LinearOperator<X, Y>& wellOper_;
+  const Opm::LinearOperatorExtra<X, Y>& wellOper_;
   std::shared_ptr< communication_type > comm_;
 };
 
@@ -176,7 +212,7 @@ public:
     typedef X domain_type;
     typedef Y range_type;
     typedef typename X::field_type field_type;
-
+    using PressureMatrix = Dune::BCRSMatrix<Opm::MatrixBlock<double, 1, 1>>;
 #if HAVE_MPI
     typedef Dune::OwnerOverlapCopyCommunication<int,int> communication_type;
 #else
@@ -192,7 +228,7 @@ public:
 
     //! constructor: just store a reference to a matrix
     WellModelGhostLastMatrixAdapter (const M& A,
-                                     const Dune::LinearOperator<X, Y>& wellOper,
+                                     const Opm::LinearOperatorExtra<X, Y>& wellOper,
                                      const size_t interiorSize )
         : A_( A ), wellOper_( wellOper ), interiorSize_(interiorSize)
     {}
@@ -230,6 +266,19 @@ public:
 
     virtual const matrix_type& getmat() const override { return A_; }
 
+    void addWellPressureEquations(PressureMatrix& jacobian, const X& weights,const bool use_well_weights) const
+    {
+        wellOper_.addWellPressureEquations(jacobian, weights, use_well_weights);
+    }
+    void addWellPressureEquationsStruct(PressureMatrix& jacobian) const
+    {
+        wellOper_.addWellPressureEquationsStruct(jacobian);
+    }
+    int getNumberOfExtraEquations() const
+    {
+        return wellOper_.getNumberOfExtraEquations();
+    }
+
 protected:
     void ghostLastProject(Y& y) const
     {
@@ -239,7 +288,7 @@ protected:
     }
 
     const matrix_type& A_ ;
-    const Dune::LinearOperator<X, Y>& wellOper_;
+    const Opm::LinearOperatorExtra< X, Y>& wellOper_;
     size_t interiorSize_;
 };
 

opm/simulators/linalg/setupPropertyTree.cpp

@@ -75,6 +75,14 @@ setupPropertyTree(FlowLinearSolverParameters p, // Note: copying the parameters
         return setupCPR(conf, p);
     }
 
+    if ((conf == "cprw")) {      
+        if (!LinearSolverMaxIterSet) {
+            // Use our own default unless it was explicitly overridden by user.
+            p.linear_solver_maxiter_ = 20;
+        }      
+        return setupCPRW(conf, p);
+    }
+
     if (conf == "amg") {
         return setupAMG(conf, p);
     }
@@ -95,6 +103,50 @@ setupPropertyTree(FlowLinearSolverParameters p, // Note: copying the parameters
               << " Please use ilu0, cpr, cpr_trueimpes, cpr_quasiimpes or isai");
 }
 
+PropertyTree
+setupCPRW(const std::string& /*conf*/, const FlowLinearSolverParameters& p)
+{
+    using namespace std::string_literals;
+    PropertyTree prm;
+    prm.put("maxiter", p.linear_solver_maxiter_);
+    prm.put("tol", p.linear_solver_reduction_);
+    prm.put("verbosity", p.linear_solver_verbosity_);
+    prm.put("solver", "bicgstab"s);
+    prm.put("preconditioner.type", "cprw"s);
+    prm.put("preconditioner.use_well_weights", "false"s);
+    prm.put("preconditioner.add_wells", "true"s);
+    prm.put("preconditioner.weight_type", "trueimpes"s);
+    prm.put("preconditioner.finesmoother.type", "ParOverILU0"s);
+    prm.put("preconditioner.finesmoother.relaxation", 1.0);
+    prm.put("preconditioner.verbosity", 0);
+    prm.put("preconditioner.coarsesolver.maxiter", 1);
+    prm.put("preconditioner.coarsesolver.tol", 1e-1);
+    prm.put("preconditioner.coarsesolver.solver", "loopsolver"s);
+    prm.put("preconditioner.coarsesolver.verbosity", 0);
+    prm.put("preconditioner.coarsesolver.preconditioner.type", "amg"s);
+    prm.put("preconditioner.coarsesolver.preconditioner.alpha", 0.333333333333);
+    prm.put("preconditioner.coarsesolver.preconditioner.relaxation", 1.0);
+    prm.put("preconditioner.coarsesolver.preconditioner.iterations", p.cpr_max_ell_iter_);
+    prm.put("preconditioner.coarsesolver.preconditioner.coarsenTarget", 1200);
+    prm.put("preconditioner.coarsesolver.preconditioner.pre_smooth", 1);
+    prm.put("preconditioner.coarsesolver.preconditioner.post_smooth", 1);
+    prm.put("preconditioner.coarsesolver.preconditioner.beta", 0.0);
+    prm.put("preconditioner.coarsesolver.preconditioner.smoother", "ILU0"s);
+    prm.put("preconditioner.coarsesolver.preconditioner.verbosity", 0);
+    prm.put("preconditioner.coarsesolver.preconditioner.maxlevel", 15);
+    prm.put("preconditioner.coarsesolver.preconditioner.skip_isolated", 0);
+    // We request to accumulate data to 1 process always as our matrix
+    // graph might be unsymmetric and hence not supported by the PTScotch/ParMetis
+    // calls in DUNE. Accumulating to 1 skips PTScotch/ParMetis
+    prm.put("preconditioner.coarsesolver.preconditioner.accumulate", 1);
+    prm.put("preconditioner.coarsesolver.preconditioner.prolongationdamping", 1.0);
+    prm.put("preconditioner.coarsesolver.preconditioner.maxdistance", 2);
+    prm.put("preconditioner.coarsesolver.preconditioner.maxconnectivity", 15);
+    prm.put("preconditioner.coarsesolver.preconditioner.maxaggsize", 6);
+    prm.put("preconditioner.coarsesolver.preconditioner.minaggsize", 4);
+    return prm;
+}
+
 PropertyTree
 setupCPR(const std::string& conf, const FlowLinearSolverParameters& p)
 {

opm/simulators/linalg/setupPropertyTree.hpp

@@ -32,6 +32,7 @@ PropertyTree setupPropertyTree(FlowLinearSolverParameters p,
                                bool LinearSolverMaxIterSet,
                                bool CprMaxEllIterSet);
 
+PropertyTree setupCPRW(const std::string& conf, const FlowLinearSolverParameters& p);    
 PropertyTree setupCPR(const std::string& conf, const FlowLinearSolverParameters& p);
 PropertyTree setupAMG(const std::string& conf, const FlowLinearSolverParameters& p);
 PropertyTree setupILU(const std::string& conf, const FlowLinearSolverParameters& p);

opm/simulators/wells/BlackoilWellModel.hpp

@@ -112,7 +112,7 @@ namespace Opm {
                 typename BlackoilWellModelGeneric::GLiftWellStateMap;
             using GLiftEclWells = typename GasLiftGroupInfo::GLiftEclWells;
             using GLiftSyncGroups = typename GasLiftSingleWellGeneric::GLiftSyncGroups;
-
+            constexpr static std::size_t pressureVarIndex = GetPropType<TypeTag, Properties::Indices>::pressureSwitchIdx;
             typedef typename BaseAuxiliaryModule<TypeTag>::NeighborSet NeighborSet;
 
             static const int numEq = Indices::numEq;
@@ -127,7 +127,7 @@ namespace Opm {
             typedef Dune::FieldVector<Scalar, numEq    > VectorBlockType;
             typedef Dune::BlockVector<VectorBlockType> BVector;
 
-            typedef Dune::FieldMatrix<Scalar, numEq, numEq > MatrixBlockType;
+            // typedef Dune::FieldMatrix<Scalar, numEq, numEq > MatrixBlockType;
 
             typedef BlackOilPolymerModule<TypeTag> PolymerModule;
             typedef BlackOilMICPModule<TypeTag> MICPModule;
@@ -264,12 +264,7 @@ namespace Opm {
 
             const SimulatorReportSingle& lastReport() const;
 
-            void addWellContributions(SparseMatrixAdapter& jacobian) const
+            void addWellContributions(SparseMatrixAdapter& jacobian) const;
-            {
-                for ( const auto& well: well_container_ ) {
-                    well->addWellContributions(jacobian);
-                }
-            }
 
             // called at the beginning of a report step
             void beginReportStep(const int time_step);
@@ -294,14 +289,26 @@ namespace Opm {
             WellInterfacePtr getWell(const std::string& well_name) const;
             bool hasWell(const std::string& well_name) const;
 
+            using PressureMatrix = Dune::BCRSMatrix<Opm::MatrixBlock<double, 1, 1>>;
+
+            int numLocalWellsEnd() const;
+
+            void addWellPressureEquations(PressureMatrix& jacobian, const BVector& weights,const bool use_well_weights) const;
+
+            std::vector<std::vector<int>> getMaxWellConnections() const;
+
+            void addWellPressureEquationsStruct(PressureMatrix& jacobian) const;
+
             void initGliftEclWellMap(GLiftEclWells &ecl_well_map);
 
             /// \brief Get list of local nonshut wells
-            const std::vector<WellInterfacePtr>& localNonshutWells()
+            const std::vector<WellInterfacePtr>& localNonshutWells() const
             {
                 return well_container_;
             }
 
+            int numLocalNonshutWells() const;
+
         protected:
             Simulator& ebosSimulator_;
 

opm/simulators/wells/BlackoilWellModel_impl.hpp

@@ -1173,10 +1173,130 @@ namespace Opm {
         Ax.axpy( alpha, scaleAddRes_ );
     }
 
+    template<typename TypeTag>
+    void
+    BlackoilWellModel<TypeTag>::
+    addWellContributions(SparseMatrixAdapter& jacobian) const
+    {
+        for ( const auto& well: well_container_ ) {
+            well->addWellContributions(jacobian);
+        }
+    }
 
+    template<typename TypeTag>
+    void
+    BlackoilWellModel<TypeTag>::
+    addWellPressureEquations(PressureMatrix& jacobian, const BVector& weights,const bool use_well_weights) const
+    {
+        int nw =  this->numLocalWellsEnd();
+        int rdofs = local_num_cells_;
+        for ( int i = 0; i < nw; i++ ){
+            int wdof = rdofs + i; 
+            jacobian[wdof][wdof] = 1.0;// better scaling ?
+        }
 
+        for ( const auto& well : well_container_ ) {
+            well->addWellPressureEquations(jacobian, weights, pressureVarIndex, use_well_weights, this->wellState());
+        }
+    }
 
+    template<typename TypeTag>
+    int
+    BlackoilWellModel<TypeTag>::
+    numLocalWellsEnd() const
+    {
+        auto w = schedule().getWellsatEnd();
+        w.erase(std::remove_if(w.begin(), w.end(), not_on_process_), w.end());
+        return w.size();
+    }
 
+    template<typename TypeTag>
+    std::vector<std::vector<int>>
+    BlackoilWellModel<TypeTag>::
+    getMaxWellConnections() const
+    {
+        std::vector<std::vector<int>> wells;
+        // Create cartesian to compressed mapping
+        const auto& globalCell = grid().globalCell();
+        const auto& cartesianSize = grid().logicalCartesianSize();
+
+        auto size = cartesianSize[0]*cartesianSize[1]*cartesianSize[2];
+
+        std::vector<int> cartesianToCompressed(size, -1);
+        auto begin = globalCell.begin();
+
+        for ( auto cell = begin, end = globalCell.end(); cell != end; ++cell )
+        {
+            cartesianToCompressed[ *cell ] = cell - begin;
+        }
+
+        auto schedule_wells = schedule().getWellsatEnd();
+        schedule_wells.erase(std::remove_if(schedule_wells.begin(), schedule_wells.end(), not_on_process_), schedule_wells.end());
+        wells.reserve(schedule_wells.size());
+
+        // initialize the additional cell connections introduced by wells.
+        for ( const auto& well : schedule_wells )
+        {
+            std::vector<int> compressed_well_perforations;
+            // All possible completions of the well
+            const auto& completionSet = well.getConnections();
+            compressed_well_perforations.reserve(completionSet.size());
+
+            for ( size_t c=0; c < completionSet.size(); c++ )
+            {
+                const auto& completion = completionSet.get(c);
+                int i = completion.getI();
+                int j = completion.getJ();
+                int k = completion.getK();
+                int cart_grid_idx = i + cartesianSize[0]*(j + cartesianSize[1]*k);
+                int compressed_idx = cartesianToCompressed[cart_grid_idx];
+
+                if ( compressed_idx >= 0 ) // Ignore completions in inactive/remote cells.
+                {
+                    compressed_well_perforations.push_back(compressed_idx);
+                }
+            }
+
+            // also include wells with no perforations in case
+            std::sort(compressed_well_perforations.begin(),
+                      compressed_well_perforations.end());
+
+            wells.push_back(compressed_well_perforations);
+            
+        }
+        return wells;
+    }
+
+    template<typename TypeTag>
+    void
+    BlackoilWellModel<TypeTag>::
+    addWellPressureEquationsStruct(PressureMatrix& jacobian) const
+    {
+        int nw =  this->numLocalWellsEnd();
+        int rdofs = local_num_cells_;
+        for(int i=0; i < nw; i++){
+            int wdof = rdofs + i; 
+            jacobian.entry(wdof,wdof) = 1.0;// better scaling ?
+        }
+        std::vector<std::vector<int>> wellconnections = getMaxWellConnections();
+        for(int i=0; i < nw; i++){
+            const auto& perfcells = wellconnections[i];
+            for(int perfcell : perfcells){
+                int wdof = rdofs + i; 
+                jacobian.entry(wdof,perfcell) = 0.0;
+                jacobian.entry(perfcell, wdof) = 0.0;
+            }
+        }
+    }
+
+    template<typename TypeTag>
+    int
+    BlackoilWellModel<TypeTag>::
+    numLocalNonshutWells() const
+    {
+        return well_container_.size();
+    }
+    
     template<typename TypeTag>
     void
     BlackoilWellModel<TypeTag>::

opm/simulators/wells/MultisegmentWell.hpp

@@ -74,6 +74,7 @@ namespace Opm
         using MSWEval::GTotal;
         using MSWEval::SPres;
         using MSWEval::numWellEq;
+        using typename Base::PressureMatrix;
 
         MultisegmentWell(const Well& well,
                          const ParallelWellInfo& pw_info,
@@ -137,6 +138,12 @@ namespace Opm
                                              DeferredLogger& deferred_logger) const override;
 
         virtual void  addWellContributions(SparseMatrixAdapter& jacobian) const override;
+        
+        virtual void addWellPressureEquations(PressureMatrix& mat,
+                                              const BVector& x,
+                                              const int pressureVarIndex,
+                                              const bool use_well_weights,
+                                              const WellState& well_state) const override;
 
         virtual std::vector<double> computeCurrentWellRates(const Simulator& ebosSimulator,
                                                             DeferredLogger& deferred_logger) const override;

opm/simulators/wells/MultisegmentWell_impl.hpp

@@ -752,6 +752,73 @@ namespace Opm
     }
 
 
+   template<typename TypeTag>
+    void
+    MultisegmentWell<TypeTag>::
+    addWellPressureEquations(PressureMatrix& jacobian,
+                             const BVector& weights,
+                             const int pressureVarIndex,
+                             const bool /*use_well_weights*/,
+                             const WellState& well_state) const
+    {
+        // Add the pressure contribution to the cpr system for the well
+
+        // Add for coupling from well to reservoir
+        const auto seg_pressure_var_ind  = this->SPres; 
+        const int welldof_ind = this->duneC_.M() + this->index_of_well_;
+        if(not(this->isPressureControlled(well_state))){
+            for (size_t rowC = 0; rowC < this->duneC_.N(); ++rowC) {
+                for (auto colC = this->duneC_[rowC].begin(), endC = this->duneC_[rowC].end(); colC != endC; ++colC) {
+                    const auto row_index = colC.index();
+                    const auto& bw = weights[row_index];
+                    double matel = 0.0;
+                
+                    for(size_t i = 0; i< bw.size(); ++i){
+                        matel += bw[i]*(*colC)[seg_pressure_var_ind][i];
+                    }
+                    jacobian[row_index][welldof_ind] += matel;
+                
+                }
+            }
+        }
+        // make cpr weights for well by pure avarage of reservoir weights of the perforations
+        if(not(this->isPressureControlled(well_state))){        
+            auto well_weight = weights[0]*0.0;
+            int num_perfs = 0;
+            for (size_t rowB = 0; rowB < this->duneB_.N(); ++rowB) {
+                for (auto colB = this->duneB_[rowB].begin(), endB = this->duneB_[rowB].end(); colB != endB; ++colB) {
+                    const auto col_index = colB.index();
+                    const auto& bw = weights[col_index];
+                    well_weight += bw;
+                    num_perfs += 1;
+                }
+            }
+        
+            well_weight /= num_perfs;
+            assert(num_perfs>0);
+        
+            // Add for coupling from reservoir to well and caclulate diag elelement corresping to incompressible standard well
+            double diag_ell = 0.0;
+            for (size_t rowB = 0; rowB < this->duneB_.N(); ++rowB) {
+                const auto& bw = well_weight;
+                for (auto colB = this->duneB_[rowB].begin(), endB = this->duneB_[rowB].end(); colB != endB; ++colB) {
+                    const auto col_index = colB.index();               
+                    double matel = 0.0;
+                    for(size_t i = 0; i< bw.size(); ++i){
+                        matel += bw[i] *(*colB)[i][pressureVarIndex];
+                    }
+                    jacobian[welldof_ind][col_index] += matel;
+                    diag_ell -= matel;
+                }
+            }
+        
+            assert(diag_ell > 0.0);
+            jacobian[welldof_ind][welldof_ind] = diag_ell;
+        }else{
+            jacobian[welldof_ind][welldof_ind] = 1.0; // maybe we could have used diag_ell if calculated
+        }
+    }
+    
 
     template<typename TypeTag>
     template<class Value>

opm/simulators/wells/StandardWell.hpp

@@ -94,6 +94,7 @@ namespace Opm
         using PolymerModule =  BlackOilPolymerModule<TypeTag>;
         using FoamModule = BlackOilFoamModule<TypeTag>;
         using BrineModule = BlackOilBrineModule<TypeTag>;
+        using typename Base::PressureMatrix;
 
         // number of the conservation equations
         static constexpr int numWellConservationEq = Indices::numPhases + Indices::numSolvents;
@@ -121,6 +122,7 @@ namespace Opm
         using Eval = typename StdWellEval::Eval;
         using EvalWell = typename StdWellEval::EvalWell;
         using BVectorWell = typename StdWellEval::BVectorWell;
+        using DiagMatrixBlockWellType = typename StdWellEval::DiagMatrixBlockWellType;
 
         StandardWell(const Well& well,
                      const ParallelWellInfo& pw_info,
@@ -178,7 +180,13 @@ namespace Opm
                                              WellState& well_state,
                                              DeferredLogger& deferred_logger) const override;
 
-        virtual void  addWellContributions(SparseMatrixAdapter& mat) const override;
+        virtual void addWellContributions(SparseMatrixAdapter& mat) const override;
+
+        virtual void addWellPressureEquations(PressureMatrix& mat,
+                                              const BVector& x,
+                                              const int pressureVarIndex,
+                                              const bool use_well_weights,
+                                              const WellState& well_state) const override;
 
         // iterate well equations with the specified control until converged
         bool iterateWellEqWithControl(const Simulator& ebosSimulator,

opm/simulators/wells/StandardWellEval.cpp

@@ -433,7 +433,7 @@ assembleControlEq(const WellState& well_state,
     // TODO: we should use a different index system for the well equations
     this->resWell_[0][Bhp] = control_eq.value();
     for (int pv_idx = 0; pv_idx < numWellEq_; ++pv_idx) {
-        this->invDuneD_[0][0][Bhp][pv_idx] = control_eq.derivative(pv_idx + Indices::numEq);
+        this->duneD_[0][0][Bhp][pv_idx] = control_eq.derivative(pv_idx + Indices::numEq);
     }
 }
 
@@ -1056,16 +1056,16 @@ init(std::vector<double>& perf_depth,
     //[A C^T    [x    =  [ res
     // B D] x_well]      res_well]
     // set the size of the matrices
-    this->invDuneD_.setSize(1, 1, 1);
+    this->duneD_.setSize(1, 1, 1);
     this->duneB_.setSize(1, num_cells, baseif_.numPerfs());
     this->duneC_.setSize(1, num_cells, baseif_.numPerfs());
 
-    for (auto row=this->invDuneD_.createbegin(), end = this->invDuneD_.createend(); row!=end; ++row) {
+    for (auto row=this->duneD_.createbegin(), end = this->duneD_.createend(); row!=end; ++row) {
         // Add nonzeros for diagonal
         row.insert(row.index());
     }
     // the block size is run-time determined now
-    this->invDuneD_[0][0].resize(numWellEq_, numWellEq_);
+    this->duneD_[0][0].resize(numWellEq_, numWellEq_);
 
     for (auto row = this->duneB_.createbegin(), end = this->duneB_.createend(); row!=end; ++row) {
         for (int perf = 0 ; perf < baseif_.numPerfs(); ++perf) {
@@ -1103,8 +1103,8 @@ init(std::vector<double>& perf_depth,
         this->Bx_[i].resize(numWellEq_);
     }
 
-    this->invDrw_.resize( this->invDuneD_.N() );
+    this->invDrw_.resize( this->duneD_.N() );
-    for (unsigned i = 0; i < this->invDuneD_.N(); ++i) {
+    for (unsigned i = 0; i < this->duneD_.N(); ++i) {
         this->invDrw_[i].resize(numWellEq_);
     }
 }
@@ -1139,7 +1139,7 @@ addWellContribution(WellContributions& wellContribs) const
     for (int i = 0; i < numStaticWellEq; ++i)
     {
         for (int j = 0; j < numStaticWellEq; ++j) {
-            nnzValues.emplace_back(this->invDuneD_[0][0][i][j]);
+            nnzValues.emplace_back(this->duneD_[0][0][i][j]);
         }
     }
     wellContribs.addMatrix(WellContributions::MatrixType::D, colIndices.data(), nnzValues.data(), 1);

opm/simulators/wells/StandardWellGeneric.hpp

@@ -124,6 +124,7 @@ protected:
     OffDiagMatWell duneC_;
     // diagonal matrix for the well
     DiagMatWell invDuneD_;
+    DiagMatWell duneD_;
 
     // Wrapper for the parallel application of B for distributed wells
     wellhelpers::ParallelStandardWellB<Scalar> parallelB_;

opm/simulators/wells/StandardWell_impl.hpp

@@ -431,7 +431,7 @@ namespace Opm
         // clear all entries
         this->duneB_ = 0.0;
         this->duneC_ = 0.0;
-        this->invDuneD_ = 0.0;
+        this->duneD_ = 0.0;
         this->resWell_ = 0.0;
 
         assembleWellEqWithoutIterationImpl(ebosSimulator, dt, well_state, group_state, deferred_logger);
@@ -490,7 +490,7 @@ namespace Opm
                 for (int pvIdx = 0; pvIdx < this->numWellEq_; ++pvIdx) {
                     // also need to consider the efficiency factor when manipulating the jacobians.
                     this->duneC_[0][cell_idx][pvIdx][componentIdx] -= cq_s_effective.derivative(pvIdx+Indices::numEq); // intput in transformed matrix
-                    this->invDuneD_[0][0][componentIdx][pvIdx] += cq_s_effective.derivative(pvIdx+Indices::numEq);
+                    this->duneD_[0][0][componentIdx][pvIdx] += cq_s_effective.derivative(pvIdx+Indices::numEq);
                 }
 
                 for (int pvIdx = 0; pvIdx < Indices::numEq; ++pvIdx) {
@@ -524,8 +524,8 @@ namespace Opm
             ws.vaporized_wat_rate = comm.sum(ws.vaporized_wat_rate);
         }
 
-        // accumulate resWell_ and invDuneD_ in parallel to get effects of all perforations (might be distributed)
+        // accumulate resWell_ and duneD_ in parallel to get effects of all perforations (might be distributed)
-        wellhelpers::sumDistributedWellEntries(this->invDuneD_[0][0], this->resWell_[0],
+        wellhelpers::sumDistributedWellEntries(this->duneD_[0][0], this->resWell_[0],
                                                this->parallel_well_info_.communication());
         // add vol * dF/dt + Q to the well equations;
         for (int componentIdx = 0; componentIdx < numWellConservationEq; ++componentIdx) {
@@ -538,7 +538,7 @@ namespace Opm
             }
             resWell_loc -= this->getQs(componentIdx) * this->well_efficiency_factor_;
             for (int pvIdx = 0; pvIdx < this->numWellEq_; ++pvIdx) {
-                this->invDuneD_[0][0][componentIdx][pvIdx] += resWell_loc.derivative(pvIdx+Indices::numEq);
+                this->duneD_[0][0][componentIdx][pvIdx] += resWell_loc.derivative(pvIdx+Indices::numEq);
             }
             this->resWell_[0][componentIdx] += resWell_loc.value();
         }
@@ -550,6 +550,7 @@ namespace Opm
 
         // do the local inversion of D.
         try {
+            this->invDuneD_ = this->duneD_; // Not strictly need if not cpr with well contributions is used
             Dune::ISTLUtility::invertMatrix(this->invDuneD_[0][0]);
         } catch( ... ) {
             OPM_DEFLOG_THROW(NumericalIssue,"Error when inverting local well equations for well " + name(), deferred_logger);
@@ -2166,9 +2167,116 @@ namespace Opm
         }
     }
 
+    
+    template <typename TypeTag>
+    void
+    StandardWell<TypeTag>::addWellPressureEquations(PressureMatrix& jacobian,
+                                                    const BVector& weights,
+                                                    const int pressureVarIndex,
+                                                    const bool use_well_weights,
+                                                    const WellState& well_state) const
+    {
+        // This adds pressure quation for cpr
+        // For use_well_weights=true
+        //    weights lamda = inv(D)'v  v = 0 v(bhpInd) = 1
+        //    the well equations are summed i lambda' B(:,pressureVarINd) -> B  lambda'*D(:,bhpInd) -> D
+        // For use_well_weights = false
+        //    weights lambda = \sum_i w /n where ths sum is over weights of all perforation cells
+        //    in the case of pressure controlled trivial equations are used and bhp  C=B=0
+        //    then the flow part of the well equations are summed lambda'*B(1:n,pressureVarInd) -> B lambda'*D(1:n,bhpInd) -> D
+        // For bouth
+        //    C -> w'C(:,bhpInd) where w is weights of the perforation cell
+        
+        // Add the well contributions in cpr
+        // use_well_weights is a quasiimpes formulation which is not implemented in multisegment
+        int bhp_var_index = Bhp;
+        int nperf = 0;
+        auto cell_weights = weights[0]*0.0;// not need for not(use_well_weights)
+        assert(this->duneC_.M() == weights.size());
+        const int welldof_ind = this->duneC_.M() + this->index_of_well_;
+        // do not assume anything about pressure controlled with use_well_weights (work fine with the assumtion also)
+        if( not( this->isPressureControlled(well_state) ) || use_well_weights ){
+            // make coupling for reservoir to well            
+            for (auto colC = this->duneC_[0].begin(), endC = this->duneC_[0].end(); colC != endC; ++colC) {
+                const auto row_ind = colC.index();
+                const auto& bw = weights[row_ind];
+                double matel = 0;        
+                assert((*colC).M() == bw.size());
+                for (size_t i = 0; i < bw.size(); ++i) {
+                    matel += (*colC)[bhp_var_index][i] * bw[i];
+                }
+                
+                jacobian[row_ind][welldof_ind] = matel;
+                cell_weights += bw;
+                nperf += 1;
+            }
+        }
+        cell_weights /= nperf;
+        
+        BVectorWell  bweights(1);
+        size_t blockSz = this->numWellEq_;
+        bweights[0].resize(blockSz);
+        bweights[0] = 0.0;
+        double diagElem = 0;
+        {            
+            if ( use_well_weights ){
+                // calculate weighs and set diagonal element
+                //NB! use this options without treating pressure controlled separated
+                //NB! calculate quasiimpes well weights NB do not work well with trueimpes reservoir weights
+                double abs_max = 0;
+                BVectorWell rhs(1);           
+                rhs[0].resize(blockSz);
+                rhs[0][bhp_var_index] = 1.0;
+                DiagMatrixBlockWellType inv_diag_block = this->invDuneD_[0][0];
+                DiagMatrixBlockWellType inv_diag_block_transpose = Opm::wellhelpers::transposeDenseDynMatrix(inv_diag_block);
+                for (size_t i = 0; i < blockSz; ++i) {
+                    bweights[0][i] = 0;
+                    for (size_t j = 0; j < blockSz; ++j) {
+                        bweights[0][i] += inv_diag_block_transpose[i][j]*rhs[0][j];
+                    }
+                    abs_max = std::max(abs_max, std::fabs(bweights[0][i]));
+                }
+                assert( abs_max > 0.0 );
+                for (size_t i = 0; i < blockSz; ++i) {
+                    bweights[0][i] /= abs_max;
+                }
+                diagElem = 1.0/abs_max;
+            }else{
+                // set diagonal element
+                if( this->isPressureControlled(well_state) ){
+                    bweights[0][blockSz-1] = 1.0;
+                    diagElem = 1.0;// better scaling could have used the calculation below if weights were calculated
+                }else{
+                    for (size_t i = 0; i < cell_weights.size(); ++i) {
+                        bweights[0][i] = cell_weights[i];
+                    }
+                    bweights[0][blockSz-1] = 0.0;
+                    diagElem = 0.0;
+                    const auto& locmat = this->duneD_[0][0]; 
+                    for (size_t i = 0; i < cell_weights.size(); ++i) {
+                        diagElem += locmat[i][bhp_var_index]*cell_weights[i];
+                    }
+                    
+                }                 
+            }
+        }
+        //
+        jacobian[welldof_ind][welldof_ind] = diagElem;
+        // set the matrix elements for well reservoir coupling
+        if( not( this->isPressureControlled(well_state) ) || use_well_weights ){
+            for (auto colB = this->duneB_[0].begin(), endB = this->duneB_[0].end(); colB != endB; ++colB) {
+                const auto col_index = colB.index();
+                const auto& bw = bweights[0];
+                double matel = 0;
+                for (size_t i = 0; i < bw.size(); ++i) {
+                     matel += (*colB)[i][pressureVarIndex] * bw[i];
+                }
+                jacobian[welldof_ind][col_index] = matel;
+            }
+        }
+    }
 
-
+    
-
 
     template<typename TypeTag>
     typename StandardWell<TypeTag>::EvalWell
@@ -2350,8 +2458,8 @@ namespace Opm
 
         this->resWell_[0][pskin_index] = eq_pskin.value();
         for (int pvIdx = 0; pvIdx < this->numWellEq_; ++pvIdx) {
-            this->invDuneD_[0][0][wat_vel_index][pvIdx] = eq_wat_vel.derivative(pvIdx+Indices::numEq);
+            this->duneD_[0][0][wat_vel_index][pvIdx] = eq_wat_vel.derivative(pvIdx+Indices::numEq);
-            this->invDuneD_[0][0][pskin_index][pvIdx] = eq_pskin.derivative(pvIdx+Indices::numEq);
+            this->duneD_[0][0][pskin_index][pvIdx] = eq_pskin.derivative(pvIdx+Indices::numEq);
         }
 
         // the water velocity is impacted by the reservoir primary varaibles. It needs to enter matrix B

opm/simulators/wells/WellHelpers.hpp

@@ -232,9 +232,21 @@ namespace Opm {
             return cube;
         }
 
+        // explicite transpose of  dense matrix due to compilation problems
+        // used for caclulating quasiimpes well weights
+        template <class DenseMatrix>
+        DenseMatrix transposeDenseDynMatrix(const DenseMatrix& M)
+        {
+            DenseMatrix tmp{M.cols(), M.rows()};
+            for (size_t i = 0; i < M.rows(); ++i) {
+                for (size_t j = 0; j < M.cols(); ++j) {
+                    tmp[j][i] = M[i][j];
+                }
+            }
+            return tmp;
+        }
 
     } // namespace wellhelpers
-
 }
 
 #endif

opm/simulators/wells/WellInterface.hpp

@@ -95,6 +95,7 @@ public:
     using MatrixBlockType = Dune::FieldMatrix<Scalar, Indices::numEq, Indices::numEq>;
     using BVector = Dune::BlockVector<VectorBlockType>;
     using Eval = DenseAd::Evaluation<Scalar, /*size=*/Indices::numEq>;
+    using PressureMatrix = Dune::BCRSMatrix<Opm::MatrixBlock<double, 1, 1>>;
 
     using RateConverterType =
     typename WellInterfaceFluidSystem<FluidSystem>::RateConverterType;
@@ -224,6 +225,14 @@ public:
     // Add well contributions to matrix
     virtual void addWellContributions(SparseMatrixAdapter&) const = 0;
 
+    virtual bool isPressureControlled(const WellState& well_state) const;
+    
+    virtual void addWellPressureEquations(PressureMatrix& mat,
+                                          const BVector& x,
+                                          const int pressureVarIndex,
+                                          const bool use_well_weights,
+                                          const WellState& well_state) const = 0;
+
     void addCellRates(RateVector& rates, int cellIdx) const;
 
     Scalar volumetricSurfaceRateForConnection(int cellIdx, int phaseIdx) const;

opm/simulators/wells/WellInterface_impl.hpp

@@ -534,7 +534,33 @@ namespace Opm
         assembleWellEqWithoutIteration(ebosSimulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
     }
 
-
+    template<typename TypeTag>
+    bool
+    WellInterface<TypeTag>::isPressureControlled(const WellState& well_state) const
+    {
+         bool thp_controlled_well = false;
+        bool bhp_controlled_well = false;
+        const auto& ws = well_state.well(this->index_of_well_);
+        if (this->isInjector()) {
+            const Well::InjectorCMode& current = ws.injection_cmode;
+            if (current == Well::InjectorCMode::THP) {
+                thp_controlled_well = true;
+            }
+            if (current == Well::InjectorCMode::BHP) {
+                bhp_controlled_well = true;
+            }
+        } else {
+            const Well::ProducerCMode& current = ws.production_cmode;
+            if (current == Well::ProducerCMode::THP) {
+                thp_controlled_well = true;
+            }
+            if (current == Well::ProducerCMode::BHP) {
+                bhp_controlled_well = true;
+            }
+        }
+        bool ispressureControlled =  (bhp_controlled_well || thp_controlled_well);
+        return ispressureControlled;
+    }
 
     template<typename TypeTag>
     void

tests/test_blackoil_amg.cpp

@@ -23,6 +23,7 @@
 #include <boost/test/unit_test.hpp>
 #include <opm/simulators/linalg/PreconditionerFactory.hpp>
 #include <opm/simulators/linalg/FlexibleSolver.hpp>
+#include <opm/simulators/linalg/getQuasiImpesWeights.hpp>
 
 #include <dune/common/parallel/mpihelper.hh>
 #include <dune/common/fmatrix.hh>

tests/test_flexiblesolver.cpp

@@ -30,6 +30,7 @@
 
 #include <opm/simulators/linalg/FlexibleSolver.hpp>
 #include <opm/simulators/linalg/getQuasiImpesWeights.hpp>
+#include <opm/simulators/linalg/matrixblock.hh>
 #include <opm/simulators/linalg/PropertyTree.hpp>
 
 #include <dune/common/fmatrix.hh>
@@ -44,7 +45,7 @@ template <int bz>
 Dune::BlockVector<Dune::FieldVector<double, bz>>
 testSolver(const Opm::PropertyTree& prm, const std::string& matrix_filename, const std::string& rhs_filename)
 {
-    using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, bz, bz>>;
+    using Matrix = Dune::BCRSMatrix<Opm::MatrixBlock<double, bz, bz>>;
     using Vector = Dune::BlockVector<Dune::FieldVector<double, bz>>;
     Matrix matrix;
     {
@@ -52,7 +53,8 @@ testSolver(const Opm::PropertyTree& prm, const std::string& matrix_filename, con
         if (!mfile) {
             throw std::runtime_error("Could not read matrix file");
         }
-        readMatrixMarket(matrix, mfile);
+        using M = Dune::BCRSMatrix<Dune::FieldMatrix<double, bz, bz>>;
+        readMatrixMarket(reinterpret_cast<M&>(matrix), mfile); // Hack to avoid hassle
     }
     Vector rhs;
     {
@@ -74,7 +76,7 @@ testSolver(const Opm::PropertyTree& prm, const std::string& matrix_filename, con
 
     using SeqOperatorType = Dune::MatrixAdapter<Matrix, Vector, Vector>;
     SeqOperatorType op(matrix);
-    Dune::FlexibleSolver<Matrix, Vector> solver(op, prm, wc, 1);
+    Dune::FlexibleSolver<SeqOperatorType> solver(op, prm, wc, 1);
     Vector x(rhs.size());
     Dune::InverseOperatorResult res;
     solver.apply(x, rhs, res);

tests/test_preconditionerfactory.cpp

@@ -28,9 +28,13 @@
 #if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6) && \
     BOOST_VERSION / 100 % 1000 > 48
 
+#include <opm/simulators/linalg/matrixblock.hh>
+#include <opm/simulators/linalg/ilufirstelement.hh>
+
 #include <opm/simulators/linalg/PreconditionerFactory.hpp>
 #include <opm/simulators/linalg/PropertyTree.hpp>
 #include <opm/simulators/linalg/FlexibleSolver.hpp>
+#include <opm/simulators/linalg/getQuasiImpesWeights.hpp>
 
 #include <dune/common/fvector.hh>
 #include <dune/istl/bvector.hh>
@@ -70,7 +74,7 @@ template <int bz>
 Dune::BlockVector<Dune::FieldVector<double, bz>>
 testPrec(const Opm::PropertyTree& prm, const std::string& matrix_filename, const std::string& rhs_filename)
 {
-    using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, bz, bz>>;
+    using Matrix = Dune::BCRSMatrix<Opm::MatrixBlock<double, bz, bz>>;
     using Vector = Dune::BlockVector<Dune::FieldVector<double, bz>>;
     Matrix matrix;
     {
@@ -78,7 +82,8 @@ testPrec(const Opm::PropertyTree& prm, const std::string& matrix_filename, const
         if (!mfile) {
             throw std::runtime_error("Could not read matrix file");
         }
-        readMatrixMarket(matrix, mfile);
+        using M = Dune::BCRSMatrix<Dune::FieldMatrix<double, bz, bz>>;
+        readMatrixMarket(reinterpret_cast<M&>(matrix), mfile); // Hack to avoid hassle
     }
     Vector rhs;
     {
@@ -161,7 +166,7 @@ BOOST_AUTO_TEST_CASE(TestDefaultPreconditionerFactory)
 
 
 template <int bz>
-using M = Dune::BCRSMatrix<Dune::FieldMatrix<double, bz, bz>>;
+using M = Dune::BCRSMatrix<Opm::MatrixBlock<double, bz, bz>>;
 template <int bz>
 using V = Dune::BlockVector<Dune::FieldVector<double, bz>>;
 template <int bz>
@@ -275,15 +280,16 @@ template <int bz>
 Dune::BlockVector<Dune::FieldVector<double, bz>>
 testPrecRepeating(const Opm::PropertyTree& prm, const std::string& matrix_filename, const std::string& rhs_filename)
 {
-    using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, bz, bz>>;
+    using Matrix = M<bz>;
-    using Vector = Dune::BlockVector<Dune::FieldVector<double, bz>>;
+    using Vector = V<bz>;
     Matrix matrix;
     {
         std::ifstream mfile(matrix_filename);
         if (!mfile) {
             throw std::runtime_error("Could not read matrix file");
         }
-        readMatrixMarket(matrix, mfile);
+        using M = Dune::BCRSMatrix<Dune::FieldMatrix<double, bz, bz>>;
+        readMatrixMarket(reinterpret_cast<M&>(matrix), mfile); // Hack to avoid hassle
     }
     Vector rhs;
     {