diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake
index d082b95a0..ef0f23eec 100644
--- a/CMakeLists_files.cmake
+++ b/CMakeLists_files.cmake
@@ -42,12 +42,12 @@ list (APPEND MAIN_SOURCE_FILES
   opm/simulators/flow/SimulatorFullyImplicitBlackoilEbos.cpp
   opm/simulators/flow/ValidationFunctions.cpp
   opm/simulators/linalg/ExtractParallelGridInformationToISTL.cpp
-  opm/simulators/linalg/FlexibleSolver1.cpp
-  opm/simulators/linalg/FlexibleSolver2.cpp
-  opm/simulators/linalg/FlexibleSolver3.cpp
-  opm/simulators/linalg/FlexibleSolver4.cpp
-  opm/simulators/linalg/FlexibleSolver5.cpp
-  opm/simulators/linalg/FlexibleSolver6.cpp
+#  opm/simulators/linalg/FlexibleSolver1.cpp
+#  opm/simulators/linalg/FlexibleSolver2.cpp
+#  opm/simulators/linalg/FlexibleSolver3.cpp
+#  opm/simulators/linalg/FlexibleSolver4.cpp
+#  opm/simulators/linalg/FlexibleSolver5.cpp
+#  opm/simulators/linalg/FlexibleSolver6.cpp
   opm/simulators/linalg/PropertyTree.cpp
   opm/simulators/linalg/setupPropertyTree.cpp
   opm/simulators/utils/PartiallySupportedFlowKeywords.cpp
diff --git a/opm/simulators/linalg/FlexibleSolver.hpp b/opm/simulators/linalg/FlexibleSolver.hpp
index 925e01412..a0bcebd75 100644
--- a/opm/simulators/linalg/FlexibleSolver.hpp
+++ b/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>
-class FlexibleSolver : public Dune::InverseOperator<VectorTypeT, VectorTypeT>
+template <class Operator>
+class FlexibleSolver : public Dune::InverseOperator<typename Operator::domain_type,
+                                                    typename Operator::range_type>
 {
 public:
-    using MatrixType = MatrixTypeT;
-    using VectorType = VectorTypeT;
+    using VectorType = typename Operator::domain_type; // Assuming symmetry: domain == range
 
-    /// 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_;
-    std::shared_ptr<AbstractOperatorType> linearoperator_for_precond_;
+    Operator* linearoperator_for_solver_;
     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
diff --git a/opm/simulators/linalg/FlexibleSolver_impl.hpp b/opm/simulators/linalg/FlexibleSolver_impl.hpp
index f8f5bd8fa..e2d927d05 100644
--- a/opm/simulators/linalg/FlexibleSolver_impl.hpp
+++ b/opm/simulators/linalg/FlexibleSolver_impl.hpp
@@ -36,9 +36,9 @@
 namespace Dune
 {
     /// Create a sequential solver.
-    template <class MatrixType, class VectorType>
-    FlexibleSolver<MatrixType, VectorType>::
-    FlexibleSolver(AbstractOperatorType& op,
+    template <class Operator>
+    FlexibleSolver<Operator>::
+    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(AbstractOperatorType& op,
+    FlexibleSolver<Operator>::
+    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>::
-    initOpPrecSp(AbstractOperatorType& op,
+    FlexibleSolver<Operator>::
+    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,
-                                                                                    child ? *child : Opm::PropertyTree(),
-                                                                                    weightsCalculator,
-                                                                                    comm,
-                                                                                    pressureIndex);
+        preconditioner_ = Opm::PreconditionerFactory<Operator, Comm>::create(op,
+                                                                             child ? *child : Opm::PropertyTree(),
+                                                                             weightsCalculator,
+                                                                             comm,
+                                                                             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>::
-    initOpPrecSp(AbstractOperatorType& op,
+    FlexibleSolver<Operator>::
+    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,
-                                                                              child ? *child : Opm::PropertyTree(),
-                                                                              weightsCalculator,
-                                                                              pressureIndex);
+        preconditioner_ = Opm::PreconditionerFactory<Operator>::create(op,
+                                                                       child ? *child : Opm::PropertyTree(),
+                                                                       weightsCalculator,
+                                                                       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>::
-    init(AbstractOperatorType& op,
+    FlexibleSolver<Operator>::
+    init(Operator& op,
          const Comm& comm,
          const Opm::PropertyTree& prm,
          const std::function<VectorType()> weightsCalculator,
@@ -203,7 +198,7 @@ namespace Dune
 
 
 // Macros to simplify explicit instantiation of FlexibleSolver for various block sizes.
-
+/*
 template <int N>
 using BV = Dune::BlockVector<Dune::FieldVector<double, N>>;
 template <int N>
@@ -240,6 +235,8 @@ template class Dune::FlexibleSolver<BM<N>, BV<N>>;     \
 template class Dune::FlexibleSolver<OBM<N>, BV<N>>;
 
 #endif // HAVE_MPI
+*/
 
+#define INSTANTIATE_FLEXIBLESOLVER(N) do(){}
 
 #endif // OPM_FLEXIBLE_SOLVER_IMPL_HEADER_INCLUDED
diff --git a/opm/simulators/linalg/ISTLSolverEbos.hpp b/opm/simulators/linalg/ISTLSolverEbos.hpp
index 890fcc4ae..78960d797 100644
--- a/opm/simulators/linalg/ISTLSolverEbos.hpp
+++ b/opm/simulators/linalg/ISTLSolverEbos.hpp
@@ -87,9 +87,11 @@ 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 ParOperatorType = WellModelGhostLastMatrixAdapter<Matrix, Vector, Vector, true>;
         using ElementMapper = GetPropType<TypeTag, Properties::ElementMapper>;
         constexpr static std::size_t pressureIndex = GetPropType<TypeTag, Properties::Indices>::pressureSwitchIdx;
 
@@ -381,35 +383,47 @@ namespace Opm
 #if HAVE_MPI
                     if (useWellConn_) {
                         using ParOperatorType = Dune::OverlappingSchwarzOperator<Matrix, Vector, Vector, Comm>;
-                        linearOperatorForFlexibleSolver_ = std::make_unique<ParOperatorType>(getMatrix(), *comm_);
-                        flexibleSolver_ = std::make_unique<FlexibleSolverType>(*linearOperatorForFlexibleSolver_, *comm_, prm_, weightsCalculator,
-                                                                               pressureIndex);
+                        auto op = std::make_unique<ParOperatorType>(getMatrix(), *comm_);
+                        using FlexibleSolverType = Dune::FlexibleSolver<ParOperatorType>;
+                        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_);
-                        flexibleSolver_ = std::make_unique<FlexibleSolverType>(*linearOperatorForFlexibleSolver_, *comm_, prm_, weightsCalculator,
-                                                                               pressureIndex);
+                        auto op = std::make_unique<ParOperatorType>(getMatrix(), *wellOperator_, interiorCellNum_);
+                        using FlexibleSolverType = Dune::FlexibleSolver<ParOperatorType>;
+                        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());
-                        flexibleSolver_ = std::make_unique<FlexibleSolverType>(*linearOperatorForFlexibleSolver_, prm_, weightsCalculator,
-                                                                               pressureIndex);
+                        auto op = std::make_unique<SeqOperatorType>(getMatrix());
+                        using FlexibleSolverType = Dune::FlexibleSolver<SeqOperatorType>;
+                        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_);
-                        flexibleSolver_ = std::make_unique<FlexibleSolverType>(*linearOperatorForFlexibleSolver_, prm_, weightsCalculator,
-                                                                               pressureIndex);
+                        auto op = std::make_unique<SeqOperatorType>(getMatrix(), *wellOperator_);
+                        using FlexibleSolverType = Dune::FlexibleSolver<SeqOperatorType>;
+                        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();
             }
         }
 
@@ -452,9 +466,11 @@ namespace Opm
             using namespace std::string_literals;
 
             auto preconditionerType = prm_.get("preconditioner.type"s, "cpr"s);
-            if (preconditionerType == "cpr" || preconditionerType == "cprt") {
-                const bool transpose = preconditionerType == "cprt";
+            if (preconditionerType == "cpr" || preconditionerType == "cprt"
+                || preconditionerType == "cprw" || preconditionerType == "cprwt") {
+                const bool transpose = preconditionerType == "cprt" || preconditionerType == "cprwt";
                 const auto weightsType = prm_.get("preconditioner.weight_type"s, "quasiimpes"s);
+                const auto pressureIndex = this->prm_.get("preconditioner.pressure_var_index", 1);
                 if (weightsType == "quasiimpes") {
                     // weights will be created as default in the solver
                     // assignment p = pressureIndex prevent compiler warning about
@@ -598,8 +614,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_;
diff --git a/opm/simulators/linalg/OwningTwoLevelPreconditioner.hpp b/opm/simulators/linalg/OwningTwoLevelPreconditioner.hpp
index dd9eddcf8..f3c0e3c5f 100644
--- a/opm/simulators/linalg/OwningTwoLevelPreconditioner.hpp
+++ b/opm/simulators/linalg/OwningTwoLevelPreconditioner.hpp
@@ -50,7 +50,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>
@@ -183,7 +183,7 @@ private:
                                                   ParCoarseOperatorType>;
     using LevelTransferPolicy = Opm::PressureTransferPolicy<OperatorType, CoarseOperatorType, Communication, transpose>;
     using CoarseSolverPolicy = Dune::Amg::PressureSolverPolicy<CoarseOperatorType,
-                                                               FlexibleSolver<PressureMatrixType, PressureVectorType>,
+                                                               FlexibleSolver<CoarseOperatorType>,
                                                                LevelTransferPolicy>;
     using TwoLevelMethod
         = Dune::Amg::TwoLevelMethodCpr<OperatorType, CoarseSolverPolicy, Dune::Preconditioner<VectorType, VectorType>>;
diff --git a/opm/simulators/linalg/OwningTwoLevelPreconditionerWell.hpp b/opm/simulators/linalg/OwningTwoLevelPreconditionerWell.hpp
new file mode 100644
index 000000000..a338ebfb3
--- /dev/null
+++ b/opm/simulators/linalg/OwningTwoLevelPreconditionerWell.hpp
@@ -0,0 +1,231 @@
+/*
+  Copyright 2019 SINTEF Digital, Mathematics and Cybernetics.
+
+  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/PreconditionerWithUpdate.hpp>
+#include <opm/simulators/linalg/PressureSolverPolicy.hpp>
+#include <opm/simulators/linalg/PressureBhpTransferPolicy.hpp>
+#include <opm/simulators/linalg/getQuasiImpesWeights.hpp>
+#include <opm/simulators/linalg/twolevelmethodcpr.hh>
+
+#include <opm/common/ErrorMacros.hpp>
+
+#include <dune/common/fmatrix.hh>
+#include <dune/istl/bcrsmatrix.hh>
+#include <dune/istl/paamg/amg.hh>
+
+#include <boost/property_tree/ptree.hpp>
+
+#include <fstream>
+#include <type_traits>
+
+
+namespace Opm
+{
+// Circular dependency between PreconditionerFactory [which can make an OwningTwoLevelPreconditioner]
+// and OwningTwoLevelPreconditioner [which uses PreconditionerFactory to choose the fine-level smoother]
+// must be broken, accomplished by forward-declaration here.
+//template <class Operator, class Comm = Dune::Amg::SequentialInformation>
+//class PreconditionerFactory;
+}
+
+namespace Dune
+{
+
+
+// Must forward-declare FlexibleSolver as we want to use it as solver for the pressure system.
+template <class Operator>
+class FlexibleSolver;
+
+// template <typename T, typename A, int i>
+// std::ostream& operator<<(std::ostream& out,
+//                          const BlockVector< FieldVector< T, i >, A >&  vector)
+// {
+//    Dune::writeMatrixMarket(vector, out);
+//    return out;
+// }
+
+//     template <typename T, typename A, int i>
+//     std::istream& operator>>(std::istream& input,
+//                              BlockVector< FieldVector< T, i >, A >&  vector)
+// {
+//    Dune::readMatrixMarket(vector, input);
+//    return input;
+// }
+
+/// A version of the two-level preconditioner that is:
+/// - Self-contained, because it owns its policy components.
+/// - Flexible, because it uses the runtime-flexible solver
+///   and preconditioner factory.
+template <class OperatorType,
+          class VectorType,
+          bool transpose = false,
+          class Communication = Dune::Amg::SequentialInformation>
+class OwningTwoLevelPreconditionerWell : public Dune::PreconditionerWithUpdate<VectorType, VectorType>
+{
+public:
+    using pt = boost::property_tree::ptree;
+    using MatrixType = typename OperatorType::matrix_type;
+    using PrecFactory = Opm::PreconditionerFactory<OperatorType, Communication>;
+    using AbstractOperatorType = Dune::AssembledLinearOperator<MatrixType, VectorType, VectorType>;
+
+    OwningTwoLevelPreconditionerWell(const OperatorType& linearoperator, const pt& prm,
+                                 const std::function<VectorType()> weightsCalculator)
+        : linear_operator_(linearoperator)
+        , finesmoother_(PrecFactory::create(linearoperator,
+                                            prm.get_child_optional("finesmoother")?
+                                            prm.get_child("finesmoother"): pt()))
+        , comm_(nullptr)
+        , weightsCalculator_(weightsCalculator)
+        , weights_(weightsCalculator())
+        , levelTransferPolicy_(dummy_comm_, weights_, prm)//.get<int>("pressure_var_index"))
+        , coarseSolverPolicy_(prm.get_child_optional("coarsesolver")? prm.get_child("coarsesolver") : pt())
+        , twolevel_method_(linearoperator,
+                           finesmoother_,
+                           levelTransferPolicy_,
+                           coarseSolverPolicy_,
+                           prm.get<int>("pre_smooth", transpose? 1 : 0),
+                           prm.get<int>("post_smooth", transpose? 0 : 1))
+        , prm_(prm)
+    {
+        if (prm.get<int>("verbosity", 0) > 10) {
+            std::string filename = prm.get<std::string>("weights_filename", "impes_weights.txt");
+            std::ofstream outfile(filename);
+            if (!outfile) {
+                OPM_THROW(std::ofstream::failure, "Could not write weights to file " << filename << ".");
+            }
+            Dune::writeMatrixMarket(weights_, outfile);
+        }
+    }
+
+    OwningTwoLevelPreconditionerWell(const OperatorType& linearoperator,
+                                     const pt& prm,
+                                     const std::function<VectorType()> weightsCalculator, const Communication& comm)
+        : linear_operator_(linearoperator)
+        , finesmoother_(PrecFactory::create(linearoperator,
+                                            prm.get_child_optional("finesmoother")?
+                                            prm.get_child("finesmoother"): pt(), comm))
+        , comm_(&comm)
+        , weightsCalculator_(weightsCalculator)
+        , weights_(weightsCalculator())
+        , levelTransferPolicy_(*comm_, weights_, prm)//.get<int>("pressure_var_index", 1))
+        , coarseSolverPolicy_(prm.get_child_optional("coarsesolver")? prm.get_child("coarsesolver") : pt())
+        , twolevel_method_(linearoperator,
+                           finesmoother_,
+                           levelTransferPolicy_,
+                           coarseSolverPolicy_,
+                           prm.get<int>("pre_smooth", transpose? 1 : 0),
+                           prm.get<int>("post_smooth", transpose? 0 : 1))
+        , prm_(prm)
+    {
+        if (prm.get<int>("verbosity", 0) > 10 && comm.communicator().rank() == 0) {
+            auto filename = prm.get<std::string>("weights_filename", "impes_weights.txt");
+            std::ofstream outfile(filename);
+            if (!outfile) {
+                OPM_THROW(std::ofstream::failure, "Could not write weights to file " << filename << ".");
+            }
+            Dune::writeMatrixMarket(weights_, outfile);
+        }
+    }
+
+    virtual void pre(VectorType& x, VectorType& b) override
+    {
+        twolevel_method_.pre(x, b);
+    }
+
+    virtual void apply(VectorType& v, const VectorType& d) override
+    {
+        twolevel_method_.apply(v, d);
+    }
+
+    virtual void post(VectorType& x) override
+    {
+        twolevel_method_.post(x);
+    }
+
+    virtual void update() override
+    {
+        weights_ = weightsCalculator_();
+        updateImpl(comm_);
+    }
+
+    virtual Dune::SolverCategory::Category category() const override
+    {
+        return linear_operator_.category();
+    }
+
+private:
+    using PressureMatrixType = Dune::BCRSMatrix<Dune::FieldMatrix<double, 1, 1>>;
+    using PressureVectorType = Dune::BlockVector<Dune::FieldVector<double, 1>>;
+    using SeqCoarseOperatorType = Dune::MatrixAdapter<PressureMatrixType, PressureVectorType, PressureVectorType>;
+    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::PressureBhpTransferPolicy<OperatorType, CoarseOperatorType, Communication, transpose>;
+    using CoarseSolverPolicy = Dune::Amg::PressureSolverPolicy<CoarseOperatorType,
+                                                               FlexibleSolver<CoarseOperatorType>,
+                                                               LevelTransferPolicy>;
+    using TwoLevelMethod
+        = Dune::Amg::TwoLevelMethodCpr<OperatorType, CoarseSolverPolicy, Dune::Preconditioner<VectorType, VectorType>>;
+
+    // Handling parallel vs serial instantiation of preconditioner factory.
+    template <class Comm>
+    void updateImpl(const Comm*)
+    {
+        // Parallel case.
+        // using ParOperatorType = Dune::OverlappingSchwarzOperator<MatrixType, VectorType, VectorType, Comm>;
+        // auto op_prec = std::make_shared<ParOperatorType>(linear_operator_.getmat(), *comm_);
+        auto child = prm_.get_child_optional("finesmoother");
+        finesmoother_ = PrecFactory::create(linear_operator_, child ? *child : pt(), *comm_);
+        twolevel_method_.updatePreconditioner(finesmoother_, coarseSolverPolicy_);
+        // linearoperator_for_precond_ = op_prec;
+    }
+
+    void updateImpl(const Dune::Amg::SequentialInformation*)
+    {
+        // Serial case.
+        // using SeqOperatorType = Dune::MatrixAdapter<MatrixType, VectorType, VectorType>;
+        // auto op_prec = std::make_shared<SeqOperatorType>(linear_operator_.getmat());
+        auto child = prm_.get_child_optional("finesmoother");
+        finesmoother_ = PrecFactory::create(linear_operator_, child ? *child : pt());
+        twolevel_method_.updatePreconditioner(finesmoother_, coarseSolverPolicy_);
+        // linearoperator_for_precond_ = op_prec;
+    }
+
+    const OperatorType& linear_operator_;
+    std::shared_ptr<Dune::Preconditioner<VectorType, VectorType>> finesmoother_;
+    const Communication* comm_;
+    std::function<VectorType()> weightsCalculator_;
+    VectorType weights_;
+    LevelTransferPolicy levelTransferPolicy_;
+    CoarseSolverPolicy coarseSolverPolicy_;
+    TwoLevelMethod twolevel_method_;
+    boost::property_tree::ptree prm_;
+    Communication dummy_comm_;
+    //std::shared_ptr<AbstractOperatorType> linearoperator_for_precond_;
+};
+
+} // namespace Dune
+
+
+
+
+
diff --git a/opm/simulators/linalg/PreconditionerFactory.hpp b/opm/simulators/linalg/PreconditionerFactory.hpp
index d90f180dc..eb4d0b48b 100644
--- a/opm/simulators/linalg/PreconditionerFactory.hpp
+++ b/opm/simulators/linalg/PreconditionerFactory.hpp
@@ -24,10 +24,12 @@
 
 #include <opm/simulators/linalg/OwningBlockPreconditioner.hpp>
 #include <opm/simulators/linalg/OwningTwoLevelPreconditioner.hpp>
+#include <opm/simulators/linalg/OwningTwoLevelPreconditionerWell.hpp>
 #include <opm/simulators/linalg/ParallelOverlappingILU0.hpp>
 #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 <dune/istl/paamg/amg.hh>
 #include <dune/istl/paamg/kamg.hh>
@@ -328,6 +330,27 @@ private:
             }
             return std::make_shared<OwningTwoLevelPreconditioner<O, V, true, 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");
+                             }
+                             return std::make_shared<OwningTwoLevelPreconditionerWell<O, V, false, Comm>>(
+                                 op, prm, weightsCalculator, comm);
+                         });
+            doAddCreator("cprwt",
+                         [](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");
+                             }
+                             return std::make_shared<OwningTwoLevelPreconditionerWell<O, V, true, Comm>>(
+                                 op, prm, weightsCalculator, comm);
+                         });
+        }
     }
 
     // Add a useful default set of preconditioners to the factory.
@@ -449,6 +472,21 @@ 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) {
+                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<OwningTwoLevelPreconditionerWell<O, V, false>>(op, prm, weightsCalculator);
+            });
+            doAddCreator("cprwt", [](const O& op, const P& prm, const std::function<Vector()>& 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");
+                }
+                return std::make_shared<OwningTwoLevelPreconditionerWell<O, V, true>>(op, prm, weightsCalculator);
+            });
+            }
+
         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())
                                 {
diff --git a/opm/simulators/linalg/PressureBhpTransferPolicy.hpp b/opm/simulators/linalg/PressureBhpTransferPolicy.hpp
new file mode 100644
index 000000000..6279ec8ea
--- /dev/null
+++ b/opm/simulators/linalg/PressureBhpTransferPolicy.hpp
@@ -0,0 +1,254 @@
+/*
+  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 num_proc = comm.communicator().size();
+        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 * num_proc + 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>();
+    }
+
+
+
+    template <class FineOperator, class CoarseOperator, class Communication, bool transpose = false>
+    class PressureBhpTransferPolicy : public Dune::Amg::LevelTransferPolicyCpr<FineOperator, CoarseOperator>
+    {
+    public:
+        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 boost::property_tree::ptree& prm)
+            : communication_(&const_cast<Communication&>(comm))
+            , weights_(weights)
+            , prm_(prm)
+            , pressure_var_index_(prm_.get<int>("pressure_var_index"))
+        {
+        }
+
+        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
+            fineOperator.addWellPressureEquations(*coarseLevelMatrix_, 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_;
+    boost::property_tree::ptree prm_;
+    const int pressure_var_index_;
+    std::shared_ptr<Communication> coarseLevelCommunication_;
+    std::shared_ptr<typename CoarseOperator::matrix_type> coarseLevelMatrix_;
+    
+};
+
+} // namespace Opm
+
diff --git a/opm/simulators/linalg/WellOperators.hpp b/opm/simulators/linalg/WellOperators.hpp
index f161f0dc2..e8565e623 100644
--- a/opm/simulators/linalg/WellOperators.hpp
+++ b/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>
-class WellModelAsLinearOperator : public Dune::LinearOperator<X, Y>
+template <class X, class Y>
+class LinearOperatorExtra : public Dune::LinearOperator<X, Y>
 {
 public:
-    using Base = Dune::LinearOperator<X, Y>;
-    using field_type = typename Base::field_type;
+    using PressureMatrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, 1, 1>>;
+    virtual void addWellPressureEquations(PressureMatrix& jacobian, const X& 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 = Dune::BCRSMatrix<Dune::FieldMatrix<double, 1, 1>>;
     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 override
+    {
+        wellMod_.addWellPressureEquations(jacobian, weights);
+    }
+    void addWellPressureEquationsStruct(PressureMatrix& jacobian) const override
+    {
+        wellMod_.addWellPressureEquationsStruct(jacobian);
+    }
+    int getNumberOfExtraEquations() const override
+    {
+        return wellMod_.numLocalWells();
+    }
+
 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<Dune::FieldMatrix<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
+    {
+        wellOper_.addWellPressureEquations(jacobian, 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<Dune::FieldMatrix<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
+    {
+        wellOper_.addWellPressureEquations(jacobian, 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_;
 };
 
diff --git a/opm/simulators/wells/BlackoilWellModel.hpp b/opm/simulators/wells/BlackoilWellModel.hpp
index 2d6f37946..fa9d4af49 100644
--- a/opm/simulators/wells/BlackoilWellModel.hpp
+++ b/opm/simulators/wells/BlackoilWellModel.hpp
@@ -294,6 +294,29 @@ namespace Opm {
             WellInterfacePtr getWell(const std::string& well_name) const;
             bool hasWell(const std::string& well_name) const;
 
+            // The number of components in the model.
+            int numComponents() const;
+
+            int numLocalWells() const;
+
+            int numPhases() const;
+
+            using PressureMatrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, 1, 1>>;
+
+            void addWellPressureEquations(PressureMatrix& jacobian, const BVector& weights) const
+            {
+                for (const auto& well : well_container_) {
+                    well->addWellPressureEquations(jacobian, weights);
+                }
+            }
+
+            void addWellPressureEquationsStruct(PressureMatrix& jacobian) const
+            {
+                for (const auto& well : well_container_) {
+                    well->addWellPressureEquationsStruct(jacobian);
+                }
+            }
+
             void initGliftEclWellMap(GLiftEclWells &ecl_well_map);
 
             /// \brief Get list of local nonshut wells
@@ -386,12 +409,7 @@ namespace Opm {
             void calculateProductivityIndexValues(const WellInterface<TypeTag>* wellPtr,
                                                   DeferredLogger& deferred_logger);
 
-            // The number of components in the model.
-            int numComponents() const;
-
-            int reportStepIndex() const;
-
-            void assembleWellEq(const double dt, DeferredLogger& deferred_logger);
+            void assembleWellEq(const std::vector<Scalar>& B_avg, const double dt, Opm::DeferredLogger& deferred_logger);
 
             bool maybeDoGasLiftOptimize(DeferredLogger& deferred_logger);
 
diff --git a/opm/simulators/wells/StandardWell.hpp b/opm/simulators/wells/StandardWell.hpp
index d3d68db95..1507eb322 100644
--- a/opm/simulators/wells/StandardWell.hpp
+++ b/opm/simulators/wells/StandardWell.hpp
@@ -94,6 +94,7 @@ namespace Opm
         using PolymerModule =  BlackOilPolymerModule<TypeTag>;
         using FoamModule = BlackOilFoamModule<TypeTag>;
         using BrineModule = BlackOilBrineModule<TypeTag>;
+        using PressureMatrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, 1, 1>>;
 
         // number of the conservation equations
         static constexpr int numWellConservationEq = Indices::numPhases + Indices::numSolvents;
@@ -178,7 +179,11 @@ 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 addWellPressureEquationsStruct(PressureMatrix& mat) const override;
+
+        virtual void addWellPressureEquations(PressureMatrix& mat, const BVector& x) const override;
 
         // iterate well equations with the specified control until converged
         bool iterateWellEqWithControl(const Simulator& ebosSimulator,
diff --git a/opm/simulators/wells/StandardWell_impl.hpp b/opm/simulators/wells/StandardWell_impl.hpp
index 64b725987..446c36172 100644
--- a/opm/simulators/wells/StandardWell_impl.hpp
+++ b/opm/simulators/wells/StandardWell_impl.hpp
@@ -2170,6 +2170,96 @@ namespace Opm
 
 
 
+    template <typename TypeTag>
+    void
+    StandardWell<TypeTag>::addWellPressureEquationsStruct(PressureMatrix& jacobian) const
+    {
+        // sustem is the pressur variant of
+        // We need to change matrx A as follows
+        // A CT
+        // B D
+        // we need to add the elemenst of CT
+        // then we need to ad the quasiimpes type well equation for B D if the well is not
+        // BHP contolled
+        const int welldof_ind = duneC_.M() + index_of_well_;
+        for (auto colC = duneC_[0].begin(), endC = duneC_[0].end(); colC != endC; ++colC) {
+            const auto row_index = colC.index();
+            double matel = 0;
+            jacobian.entry(row_index, welldof_ind) = matel;
+        }
+
+        jacobian.entry(welldof_ind, welldof_ind) = 0.0;
+
+        // set the matrix elements for well reservoir coupling
+        for (auto colB = duneB_[0].begin(), endB = duneB_[0].end(); colB != endB; ++colB) {
+            const auto col_index = colB.index();
+            double matel = 0;
+            jacobian.entry(welldof_ind, col_index) = matel;
+        }
+    }
+
+
+
+
+    template <typename TypeTag>
+    void
+    StandardWell<TypeTag>::addWellPressureEquations(PressureMatrix& jacobian, const BVector& weights) const
+    {
+        // sustem is the pressur variant of
+        // We need to change matrx A as follows
+        // A CT
+        // B D
+        // we need to add the elemenst of CT
+        // then we need to ad the quasiimpes type well equation for B D if the well is not
+        // BHP contolled
+        int bhp_var_index = Bhp;
+        assert(duneC_.M() == weights.size());
+        const int welldof_ind = duneC_.M() + index_of_well_;
+        for (auto colC = duneC_[0].begin(), endC = 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;
+        }
+        // make quasipes weights for bhp it should be trival
+        using VectorBlockType = typename BVector::block_type;
+        using MatrixBlockType = DiagMatrixBlockWellType;
+        VectorBlockType bweights(0.0);
+        double diagElem = 0;
+        {
+            // const DiagMatrixBlockWellType& invA = invDuneD_[0][0];
+            VectorBlockType rhs(0.0);
+            rhs[bhp_var_index] = 1.0;
+            MatrixBlockType inv_diag_block = invDuneD_[0][0];
+            auto inv_diag_block_transpose = Opm::wellhelpers::transposeDenseDynMatrix(inv_diag_block);
+            // diag_block_transpose.solve(bweights, rhs);
+            inv_diag_block_transpose.mv(rhs, bweights);
+            // NB how to scale to make it most symmetric
+            // double abs_max = *std::max_element(
+            //    bweights.begin(), bweights.end(), [](double a, double b) { return std::fabs(a) < std::fabs(b); });
+            // bweights /= std::fabs(abs_max);
+        }
+        //
+        jacobian[welldof_ind][welldof_ind] = 1.0;
+        // set the matrix elements for well reservoir coupling
+        for (auto colB = duneB_[0].begin(), endB = duneB_[0].end(); colB != endB; ++colB) {
+            const auto col_index = colB.index();
+            const auto& bw = bweights;
+            double matel = 0;
+            for (size_t i = 0; i < bw.size(); ++i) {
+                matel += (*colB)[bhp_var_index][i] * bw[i];
+            }
+            jacobian[welldof_ind][col_index] = matel;
+        }
+    }
+
+
+
+
     template<typename TypeTag>
     typename StandardWell<TypeTag>::EvalWell
     StandardWell<TypeTag>::
diff --git a/opm/simulators/wells/WellHelpers.hpp b/opm/simulators/wells/WellHelpers.hpp
index b566f92b1..9240ec540 100644
--- a/opm/simulators/wells/WellHelpers.hpp
+++ b/opm/simulators/wells/WellHelpers.hpp
@@ -233,8 +233,19 @@ namespace Opm {
         }
 
 
-    } // namespace wellhelpers
+        template <class DenseMatrix>
+        DenseMatrix transposeDenseDynMatrix(const DenseMatrix& M)
+        {
+            DenseMatrix tmp{M.cols(), M.rows()};
+            for (int i = 0; i < M.rows(); ++i) {
+                for (int j = 0; j < M.cols(); ++j) {
+                    tmp[j][i] = M[i][j];
+                }
+            }
+            return tmp;
+        }
 
+    } // namespace wellhelpers
 }
 
 #endif
diff --git a/opm/simulators/wells/WellInterface.hpp b/opm/simulators/wells/WellInterface.hpp
index 0c14522c2..1db58c0bd 100644
--- a/opm/simulators/wells/WellInterface.hpp
+++ b/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<Dune::FieldMatrix<double, 1, 1>>;
 
     using RateConverterType =
     typename WellInterfaceFluidSystem<FluidSystem>::RateConverterType;
@@ -224,6 +225,15 @@ public:
     // Add well contributions to matrix
     virtual void addWellContributions(SparseMatrixAdapter&) const = 0;
 
+    virtual void addWellPressureEquationsStruct(PressureMatrix&) const
+    {
+    }
+
+    virtual void addWellPressureEquations(PressureMatrix&,
+                                          const BVector& x) const
+    {
+    }
+
     void addCellRates(RateVector& rates, int cellIdx) const;
 
     Scalar volumetricSurfaceRateForConnection(int cellIdx, int phaseIdx) const;
diff --git a/tests/test_flexiblesolver.cpp b/tests/test_flexiblesolver.cpp
index 53a5973b6..d4a3974d9 100644
--- a/tests/test_flexiblesolver.cpp
+++ b/tests/test_flexiblesolver.cpp
@@ -74,7 +74,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);