add: StandardWellEquations

this is a container for the standard well equation system

CMakeLists_files.cmake

@@ -96,6 +96,7 @@ list (APPEND MAIN_SOURCE_FILES
   opm/simulators/wells/PerfData.cpp
   opm/simulators/wells/SegmentState.cpp
   opm/simulators/wells/SingleWellState.cpp
+  opm/simulators/wells/StandardWellEquations.cpp
   opm/simulators/wells/StandardWellEval.cpp
   opm/simulators/wells/StandardWellGeneric.cpp
   opm/simulators/wells/TargetCalculator.cpp
@@ -380,6 +381,8 @@ list (APPEND PUBLIC_HEADER_FILES
   opm/simulators/wells/SingleWellState.hpp
   opm/simulators/wells/StandardWell.hpp
   opm/simulators/wells/StandardWell_impl.hpp
+  opm/simulators/wells/StandardWellEquations.hpp
+  opm/simulators/wells/StandardWellEval.hpp
   opm/simulators/wells/TargetCalculator.hpp
   opm/simulators/wells/VFPHelpers.hpp
   opm/simulators/wells/VFPInjProperties.hpp

opm/simulators/wells/StandardWellEquations.cpp

@@ -0,0 +1,48 @@
+/*
+  Copyright 2017 SINTEF Digital, Mathematics and Cybernetics.
+  Copyright 2017 Statoil ASA.
+  Copyright 2016 - 2017 IRIS AS.
+
+  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/>.
+*/
+
+#include <config.h>
+#include <opm/simulators/wells/StandardWellEquations.hpp>
+
+namespace Opm
+{
+
+template<class Scalar, int numEq>
+StandardWellEquations<Scalar,numEq>::
+StandardWellEquations(const ParallelWellInfo& parallel_well_info)
+    : parallelB_(duneB_, parallel_well_info)
+{
+    duneB_.setBuildMode(OffDiagMatWell::row_wise);
+    duneC_.setBuildMode(OffDiagMatWell::row_wise),
+    invDuneD_.setBuildMode(DiagMatWell::row_wise);
+}
+
+#define INSTANCE(N) \
+template class StandardWellEquations<double,N>;
+
+INSTANCE(1)
+INSTANCE(2)
+INSTANCE(3)
+INSTANCE(4)
+INSTANCE(5)
+INSTANCE(6)
+
+}

opm/simulators/wells/StandardWellEquations.hpp

@@ -0,0 +1,83 @@
+/*
+  Copyright 2017 SINTEF Digital, Mathematics and Cybernetics.
+  Copyright 2017 Statoil ASA.
+  Copyright 2016 - 2017 IRIS AS.
+
+  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/>.
+*/
+
+
+#ifndef OPM_STANDARDWELL_EQUATIONS_HEADER_INCLUDED
+#define OPM_STANDARDWELL_EQUATIONS_HEADER_INCLUDED
+
+#include <opm/simulators/wells/WellHelpers.hpp>
+
+#include <dune/common/dynmatrix.hh>
+#include <dune/common/dynvector.hh>
+#include <dune/istl/bcrsmatrix.hh>
+#include <dune/istl/bvector.hh>
+
+namespace Opm
+{
+
+class ParallelWellInfo;
+
+template<class Scalar, int numEq>
+class StandardWellEquations
+{
+public:
+    // sparsity pattern for the matrices
+    //[A C^T    [x       =  [ res
+    // B  D ]   x_well]      res_well]
+
+    // the vector type for the res_well and x_well
+    using VectorBlockWellType = Dune::DynamicVector<Scalar>;
+    using BVectorWell = Dune::BlockVector<VectorBlockWellType>;
+
+    // the matrix type for the diagonal matrix D
+    using DiagMatrixBlockWellType = Dune::DynamicMatrix<Scalar>;
+    using DiagMatWell = Dune::BCRSMatrix<DiagMatrixBlockWellType>;
+
+    // the matrix type for the non-diagonal matrix B and C^T
+    using OffDiagMatrixBlockWellType = Dune::DynamicMatrix<Scalar>;
+    using OffDiagMatWell = Dune::BCRSMatrix<OffDiagMatrixBlockWellType>;
+
+    // block vector type
+    using BVector = Dune::BlockVector<Dune::FieldVector<Scalar,numEq>>;
+
+    StandardWellEquations(const ParallelWellInfo& parallel_well_info);
+
+    // two off-diagonal matrices
+    OffDiagMatWell duneB_;
+    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_;
+
+    // residuals of the well equations
+    BVectorWell resWell_;
+
+    // several vector used in the matrix calculation
+    mutable BVectorWell Bx_;
+    mutable BVectorWell invDrw_;
+};
+
+}
+
+#endif // OPM_STANDARDWELL_EQUATIONS_HEADER_INCLUDED

opm/simulators/wells/StandardWellEval.cpp

@@ -53,6 +53,7 @@ StandardWellEval(const WellInterfaceIndices<FluidSystem,Indices,Scalar>& baseif)
     : StandardWellGeneric<Scalar>(baseif)
     , baseif_(baseif)
     , F0_(numWellConservationEq)
+    , linSys_(baseif_.parallelWellInfo())
 {
 }
 
@@ -444,9 +445,9 @@ assembleControlEq(const WellState& well_state,
 
     // using control_eq to update the matrix and residuals
     // TODO: we should use a different index system for the well equations
-    this->resWell_[0][Bhp] = control_eq.value();
+    this->linSys_.resWell_[0][Bhp] = control_eq.value();
     for (int pv_idx = 0; pv_idx < numWellEq_; ++pv_idx) {
-        this->duneD_[0][0][Bhp][pv_idx] = control_eq.derivative(pv_idx + Indices::numEq);
+        this->linSys_.duneD_[0][0][Bhp][pv_idx] = control_eq.derivative(pv_idx + Indices::numEq);
     }
 }
 
@@ -773,7 +774,7 @@ getWellConvergence(const WellState& well_state,
     res.resize(numWellEq_);
     for (int eq_idx = 0; eq_idx < numWellEq_; ++eq_idx) {
         // magnitude of the residual matters
-        res[eq_idx] = std::abs(this->resWell_[0][eq_idx]);
+        res[eq_idx] = std::abs(this->linSys_.resWell_[0][eq_idx]);
     }
 
     std::vector<double> well_flux_residual(baseif_.numComponents());
@@ -810,7 +811,7 @@ getWellConvergence(const WellState& well_state,
     WellConvergence(baseif_).
         checkConvergenceControlEq(well_state,
                                   {1.e3, 1.e4, 1.e-4, 1.e-6, maxResidualAllowed},
-                                  std::abs(this->resWell_[0][Bhp]),
+                                  std::abs(this->linSys_.resWell_[0][Bhp]),
                                   report,
                                   deferred_logger);
 
@@ -1042,18 +1043,20 @@ init(std::vector<double>& perf_depth,
     //[A C^T    [x    =  [ res
     // B D] x_well]      res_well]
     // set the size of the matrices
-    this->duneD_.setSize(1, 1, 1);
+    this->linSys_.duneD_.setSize(1, 1, 1);
-    this->duneB_.setSize(1, num_cells, baseif_.numPerfs());
+    this->linSys_.duneB_.setSize(1, num_cells, baseif_.numPerfs());
-    this->duneC_.setSize(1, num_cells, baseif_.numPerfs());
+    this->linSys_.duneC_.setSize(1, num_cells, baseif_.numPerfs());
 
-    for (auto row=this->duneD_.createbegin(), end = this->duneD_.createend(); row!=end; ++row) {
+    for (auto row = this->linSys_.duneD_.createbegin(),
+              end = this->linSys_.duneD_.createend(); row != end; ++row) {
         // Add nonzeros for diagonal
         row.insert(row.index());
     }
     // the block size is run-time determined now
-    this->duneD_[0][0].resize(numWellEq_, numWellEq_);
+    this->linSys_.duneD_[0][0].resize(numWellEq_, numWellEq_);
 
-    for (auto row = this->duneB_.createbegin(), end = this->duneB_.createend(); row!=end; ++row) {
+    for (auto row = this->linSys_.duneB_.createbegin(),
+              end = this->linSys_.duneB_.createend(); row != end; ++row) {
         for (int perf = 0 ; perf < baseif_.numPerfs(); ++perf) {
             const int cell_idx = baseif_.cells()[perf];
             row.insert(cell_idx);
@@ -1063,11 +1066,12 @@ init(std::vector<double>& perf_depth,
     for (int perf = 0 ; perf < baseif_.numPerfs(); ++perf) {
         const int cell_idx = baseif_.cells()[perf];
          // the block size is run-time determined now
-         this->duneB_[0][cell_idx].resize(numWellEq_, Indices::numEq);
+         this->linSys_.duneB_[0][cell_idx].resize(numWellEq_, Indices::numEq);
     }
 
     // make the C^T matrix
-    for (auto row = this->duneC_.createbegin(), end = this->duneC_.createend(); row!=end; ++row) {
+    for (auto row = this->linSys_.duneC_.createbegin(),
+              end = this->linSys_.duneC_.createend(); row != end; ++row) {
         for (int perf = 0; perf < baseif_.numPerfs(); ++perf) {
             const int cell_idx = baseif_.cells()[perf];
             row.insert(cell_idx);
@@ -1076,22 +1080,22 @@ init(std::vector<double>& perf_depth,
 
     for (int perf = 0; perf < baseif_.numPerfs(); ++perf) {
         const int cell_idx = baseif_.cells()[perf];
-        this->duneC_[0][cell_idx].resize(numWellEq_, Indices::numEq);
+        this->linSys_.duneC_[0][cell_idx].resize(numWellEq_, Indices::numEq);
     }
 
-    this->resWell_.resize(1);
+    this->linSys_.resWell_.resize(1);
     // the block size of resWell_ is also run-time determined now
-    this->resWell_[0].resize(numWellEq_);
+    this->linSys_.resWell_[0].resize(numWellEq_);
 
     // resize temporary class variables
-    this->Bx_.resize( this->duneB_.N() );
+    this->linSys_.Bx_.resize(this->linSys_.duneB_.N());
-    for (unsigned i = 0; i < this->duneB_.N(); ++i) {
+    for (unsigned i = 0; i < this->linSys_.duneB_.N(); ++i) {
-        this->Bx_[i].resize(numWellEq_);
+        this->linSys_.Bx_[i].resize(numWellEq_);
     }
 
-    this->invDrw_.resize( this->duneD_.N() );
+    this->linSys_.invDrw_.resize(this->linSys_.duneD_.N());
-    for (unsigned i = 0; i < this->duneD_.N(); ++i) {
+    for (unsigned i = 0; i < this->linSys_.duneD_.N(); ++i) {
-        this->invDrw_[i].resize(numWellEq_);
+        this->linSys_.invDrw_[i].resize(numWellEq_);
     }
 }
 
@@ -1102,11 +1106,12 @@ addWellContribution(WellContributions& wellContribs) const
 {
     std::vector<int> colIndices;
     std::vector<double> nnzValues;
-    colIndices.reserve(this->duneB_.nonzeroes());
+    colIndices.reserve(this->linSys_.duneB_.nonzeroes());
-    nnzValues.reserve(this->duneB_.nonzeroes()*numStaticWellEq * Indices::numEq);
+    nnzValues.reserve(this->linSys_.duneB_.nonzeroes()*numStaticWellEq * Indices::numEq);
 
     // duneC
-    for ( auto colC = this->duneC_[0].begin(), endC = this->duneC_[0].end(); colC != endC; ++colC )
+    for (auto colC = this->linSys_.duneC_[0].begin(),
+              endC = this->linSys_.duneC_[0].end(); colC != endC; ++colC )
     {
         colIndices.emplace_back(colC.index());
         for (int i = 0; i < numStaticWellEq; ++i) {
@@ -1115,7 +1120,7 @@ addWellContribution(WellContributions& wellContribs) const
             }
         }
     }
-    wellContribs.addMatrix(WellContributions::MatrixType::C, colIndices.data(), nnzValues.data(), this->duneC_.nonzeroes());
+    wellContribs.addMatrix(WellContributions::MatrixType::C, colIndices.data(), nnzValues.data(), this->linSys_.duneC_.nonzeroes());
 
     // invDuneD
     colIndices.clear();
@@ -1124,7 +1129,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->linSys_.invDuneD_[0][0][i][j]);
         }
     }
     wellContribs.addMatrix(WellContributions::MatrixType::D, colIndices.data(), nnzValues.data(), 1);
@@ -1132,7 +1137,8 @@ addWellContribution(WellContributions& wellContribs) const
     // duneB
     colIndices.clear();
     nnzValues.clear();
-    for ( auto colB = this->duneB_[0].begin(), endB = this->duneB_[0].end(); colB != endB; ++colB )
+    for (auto colB = this->linSys_.duneB_[0].begin(),
+              endB = this->linSys_.duneB_[0].end(); colB != endB; ++colB )
     {
         colIndices.emplace_back(colB.index());
         for (int i = 0; i < numStaticWellEq; ++i) {
@@ -1141,7 +1147,14 @@ addWellContribution(WellContributions& wellContribs) const
             }
         }
     }
-    wellContribs.addMatrix(WellContributions::MatrixType::B, colIndices.data(), nnzValues.data(), this->duneB_.nonzeroes());
+    wellContribs.addMatrix(WellContributions::MatrixType::B, colIndices.data(), nnzValues.data(), this->linSys_.duneB_.nonzeroes());
+}
+
+template<class FluidSystem, class Indices, class Scalar>
+unsigned int StandardWellEval<FluidSystem,Indices,Scalar>::
+getNumBlocks() const
+{
+    return linSys_.duneB_.nonzeroes();
 }
 
 #define INSTANCE(...) \

opm/simulators/wells/StandardWellEval.hpp

@@ -23,6 +23,7 @@
 #ifndef OPM_STANDARDWELL_EVAL_HEADER_INCLUDED
 #define OPM_STANDARDWELL_EVAL_HEADER_INCLUDED
 
+#include <opm/simulators/wells/StandardWellEquations.hpp>
 #include <opm/simulators/wells/StandardWellGeneric.hpp>
 
 #include <opm/material/densead/DynamicEvaluation.hpp>
@@ -96,8 +97,12 @@ public:
     using Eval = DenseAd::Evaluation<Scalar, Indices::numEq>;
     using BVectorWell = typename StandardWellGeneric<Scalar>::BVectorWell;
 
-        /// add the contribution (C, D^-1, B matrices) of this Well to the WellContributions object
+    /// add the contribution (C, D^-1, B matrices) of this Well to the WellContributions object
-        void addWellContribution(WellContributions& wellContribs) const;
+    void addWellContribution(WellContributions& wellContribs) const;
+
+    //! \brief Returns a const reference to equation system.
+    const StandardWellEquations<Scalar,Indices::numEq>& linSys() const
+    { return linSys_; }
 
 protected:
     StandardWellEval(const WellInterfaceIndices<FluidSystem,Indices,Scalar>& baseif);
@@ -190,6 +195,8 @@ protected:
 
     // the saturations in the well bore under surface conditions at the beginning of the time step
     std::vector<double> F0_;
+
+    StandardWellEquations<Scalar,Indices::numEq> linSys_; //!< Linear equation system
 };
 
 }

opm/simulators/wells/StandardWellGeneric.cpp

@@ -52,11 +52,7 @@ StandardWellGeneric(const WellInterfaceGeneric& baseif)
     : baseif_(baseif)
     , perf_densities_(baseif_.numPerfs())
     , perf_pressure_diffs_(baseif_.numPerfs())
-    , parallelB_(duneB_, baseif_.parallelWellInfo())
 {
-    duneB_.setBuildMode(OffDiagMatWell::row_wise);
-    duneC_.setBuildMode(OffDiagMatWell::row_wise);
-    invDuneD_.setBuildMode(DiagMatWell::row_wise);
 }
 
 
@@ -149,14 +145,6 @@ computeConnectionPressureDelta()
     baseif_.parallelWellInfo().partialSumPerfValues(beg, end);
 }
 
-template<class Scalar>
-unsigned int
-StandardWellGeneric<Scalar>::
-getNumBlocks() const
-{
-    return duneB_.nonzeroes();
-}
-
 template class StandardWellGeneric<double>;
 
 }

opm/simulators/wells/StandardWellGeneric.hpp

@@ -65,11 +65,6 @@ protected:
     using OffDiagMatrixBlockWellType = Dune::DynamicMatrix<Scalar>;
     using OffDiagMatWell = Dune::BCRSMatrix<OffDiagMatrixBlockWellType>;
 
-public:
-    /// get the number of blocks of the C and B matrices, used to allocate memory in a WellContributions object
-    unsigned int getNumBlocks() const;
-
-protected:
     StandardWellGeneric(const WellInterfaceGeneric& baseif);
 
     // calculate a relaxation factor to avoid overshoot of total rates
@@ -85,28 +80,11 @@ protected:
     // Base interface reference
     const WellInterfaceGeneric& baseif_;
 
-    // residuals of the well equations
-    BVectorWell resWell_;
-
     // densities of the fluid in each perforation
     std::vector<double> perf_densities_;
     // pressure drop between different perforations
     std::vector<double> perf_pressure_diffs_;
 
-    // two off-diagonal matrices
-    OffDiagMatWell duneB_;
-    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_;
-
-    // several vector used in the matrix calculation
-    mutable BVectorWell Bx_;
-    mutable BVectorWell invDrw_;
-
     double getRho() const
     {
         return this->perf_densities_.empty() ? 0.0 : perf_densities_[0];

opm/simulators/wells/StandardWell_impl.hpp

@@ -432,10 +432,10 @@ namespace Opm
         if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return;
 
         // clear all entries
-        this->duneB_ = 0.0;
+        this->linSys_.duneB_ = 0.0;
-        this->duneC_ = 0.0;
+        this->linSys_.duneC_ = 0.0;
-        this->duneD_ = 0.0;
+        this->linSys_.duneD_ = 0.0;
-        this->resWell_ = 0.0;
+        this->linSys_.resWell_ = 0.0;
 
         assembleWellEqWithoutIterationImpl(ebosSimulator, dt, well_state, group_state, deferred_logger);
     }
@@ -488,17 +488,17 @@ namespace Opm
                 connectionRates[perf][componentIdx] = Base::restrictEval(cq_s_effective);
 
                 // subtract sum of phase fluxes in the well equations.
-                this->resWell_[0][componentIdx] += cq_s_effective.value();
+                this->linSys_.resWell_[0][componentIdx] += cq_s_effective.value();
 
                 // assemble the jacobians
                 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->linSys_.duneC_[0][cell_idx][pvIdx][componentIdx] -= cq_s_effective.derivative(pvIdx+Indices::numEq); // intput in transformed matrix
-                    this->duneD_[0][0][componentIdx][pvIdx] += cq_s_effective.derivative(pvIdx+Indices::numEq);
+                    this->linSys_.duneD_[0][0][componentIdx][pvIdx] += cq_s_effective.derivative(pvIdx+Indices::numEq);
                 }
 
                 for (int pvIdx = 0; pvIdx < Indices::numEq; ++pvIdx) {
-                    this->duneB_[0][cell_idx][componentIdx][pvIdx] += cq_s_effective.derivative(pvIdx);
+                    this->linSys_.duneB_[0][cell_idx][componentIdx][pvIdx] += cq_s_effective.derivative(pvIdx);
                 }
 
                 // Store the perforation phase flux for later usage.
@@ -512,7 +512,7 @@ namespace Opm
 
             if constexpr (has_zFraction) {
                 for (int pvIdx = 0; pvIdx < this->numWellEq_; ++pvIdx) {
-                    this->duneC_[0][cell_idx][pvIdx][Indices::contiZfracEqIdx] -= cq_s_zfrac_effective.derivative(pvIdx+Indices::numEq);
+                    this->linSys_.duneC_[0][cell_idx][pvIdx][Indices::contiZfracEqIdx] -= cq_s_zfrac_effective.derivative(pvIdx+Indices::numEq);
                 }
             }
         }
@@ -529,7 +529,7 @@ namespace Opm
         }
 
         // accumulate resWell_ and duneD_ in parallel to get effects of all perforations (might be distributed)
-        wellhelpers::sumDistributedWellEntries(this->duneD_[0][0], this->resWell_[0],
+        wellhelpers::sumDistributedWellEntries(this->linSys_.duneD_[0][0], this->linSys_.resWell_[0],
                                                this->parallel_well_info_.communication());
         // add vol * dF/dt + Q to the well equations;
         for (int componentIdx = 0; componentIdx < numWellConservationEq; ++componentIdx) {
@@ -542,9 +542,9 @@ namespace Opm
             }
             resWell_loc -= this->getQs(componentIdx) * this->well_efficiency_factor_;
             for (int pvIdx = 0; pvIdx < this->numWellEq_; ++pvIdx) {
-                this->duneD_[0][0][componentIdx][pvIdx] += resWell_loc.derivative(pvIdx+Indices::numEq);
+                this->linSys_.duneD_[0][0][componentIdx][pvIdx] += resWell_loc.derivative(pvIdx+Indices::numEq);
             }
-            this->resWell_[0][componentIdx] += resWell_loc.value();
+            this->linSys_.resWell_[0][componentIdx] += resWell_loc.value();
         }
 
         const auto& summaryState = ebosSimulator.vanguard().summaryState();
@@ -554,13 +554,13 @@ namespace Opm
 
         // do the local inversion of D.
         try {
-            this->invDuneD_ = this->duneD_; // Not strictly need if not cpr with well contributions is used
+            this->linSys_.invDuneD_ = this->linSys_.duneD_; // Not strictly need if not cpr with well contributions is used
-            detail::invertMatrix(this->invDuneD_[0][0]);
+            detail::invertMatrix(this->linSys_.invDuneD_[0][0]);
         } catch (NumericalProblem&) {
             // for singular matrices, use identity as the inverse
-            this->invDuneD_[0][0] = 0.0;
+            this->linSys_.invDuneD_[0][0] = 0.0;
-            for (size_t i = 0; i < this->invDuneD_[0][0].rows(); ++i) {
+            for (size_t i = 0; i < this->linSys_.invDuneD_[0][0].rows(); ++i) {
-                this->invDuneD_[0][0][i][i] = 1.0;
+                this->linSys_.invDuneD_[0][0][i][i] = 1.0;
             }
         } catch( ... ) {
             OPM_DEFLOG_THROW(NumericalIssue,"Error when inverting local well equations for well " + name(), deferred_logger);
@@ -1689,7 +1689,7 @@ namespace Opm
         // which is why we do not put the assembleWellEq here.
         BVectorWell dx_well(1);
         dx_well[0].resize(this->numWellEq_);
-        this->invDuneD_.mv(this->resWell_, dx_well);
+        this->linSys_.invDuneD_.mv(this->linSys_.resWell_, dx_well);
 
         updateWellState(dx_well, well_state, deferred_logger);
     }
@@ -1725,20 +1725,20 @@ namespace Opm
             // Contributions are already in the matrix itself
             return;
         }
-        assert( this->Bx_.size() == this->duneB_.N() );
+        assert(this->linSys_.Bx_.size() == this->linSys_.duneB_.N());
-        assert( this->invDrw_.size() == this->invDuneD_.N() );
+        assert(this->linSys_.invDrw_.size() == this->linSys_.invDuneD_.N());
 
         // Bx_ = duneB_ * x
-        this->parallelB_.mv(x, this->Bx_);
+        this->linSys_.parallelB_.mv(x, this->linSys_.Bx_);
 
         // invDBx = invDuneD_ * Bx_
         // TODO: with this, we modified the content of the invDrw_.
         // Is it necessary to do this to save some memory?
-        BVectorWell& invDBx = this->invDrw_;
+        BVectorWell& invDBx = this->linSys_.invDrw_;
-        this->invDuneD_.mv(this->Bx_, invDBx);
+        this->linSys_.invDuneD_.mv(this->linSys_.Bx_, invDBx);
 
         // Ax = Ax - duneC_^T * invDBx
-        this->duneC_.mmtv(invDBx,Ax);
+        this->linSys_.duneC_.mmtv(invDBx,Ax);
     }
 
 
@@ -1751,12 +1751,12 @@ namespace Opm
     {
         if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return;
 
-        assert( this->invDrw_.size() == this->invDuneD_.N() );
+        assert(this->linSys_.invDrw_.size() == this->linSys_.invDuneD_.N());
 
         // invDrw_ = invDuneD_ * resWell_
-        this->invDuneD_.mv(this->resWell_, this->invDrw_);
+        this->linSys_.invDuneD_.mv(this->linSys_.resWell_, this->linSys_.invDrw_);
         // r = r - duneC_^T * invDrw_
-        this->duneC_.mmtv(this->invDrw_, r);
+        this->linSys_.duneC_.mmtv(this->linSys_.invDrw_, r);
     }
 
     template<typename TypeTag>
@@ -1766,11 +1766,11 @@ namespace Opm
     {
         if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return;
 
-        BVectorWell resWell = this->resWell_;
+        BVectorWell resWell = this->linSys_.resWell_;
         // resWell = resWell - B * x
-        this->parallelB_.mmv(x, resWell);
+        this->linSys_.parallelB_.mmv(x, resWell);
         // xw = D^-1 * resWell
-        this->invDuneD_.mv(resWell, xw);
+        this->linSys_.invDuneD_.mv(resWell, xw);
     }
 
 
@@ -2178,13 +2178,15 @@ namespace Opm
         // at (0,j) only if this well has a perforation at cell j.
         typename SparseMatrixAdapter::MatrixBlock tmpMat;
         Dune::DynamicMatrix<Scalar> tmp;
-        for ( auto colC = this->duneC_[0].begin(), endC = this->duneC_[0].end(); colC != endC; ++colC )
+        for (auto colC = this->linSys_.duneC_[0].begin(),
+                  endC = this->linSys_.duneC_[0].end(); colC != endC; ++colC)
         {
             const auto row_index = colC.index();
 
-            for ( auto colB = this->duneB_[0].begin(), endB = this->duneB_[0].end(); colB != endB; ++colB )
+            for (auto colB = this->linSys_.duneB_[0].begin(),
+                      endB = this->linSys_.duneB_[0].end(); colB != endB; ++colB)
             {
-                detail::multMatrix(this->invDuneD_[0][0],  (*colB), tmp);
+                detail::multMatrix(this->linSys_.invDuneD_[0][0],  (*colB), tmp);
                 detail::negativeMultMatrixTransposed((*colC), tmp, tmpMat);
                 jacobian.addToBlock( row_index, colB.index(), tmpMat );
             }
@@ -2217,12 +2219,13 @@ namespace Opm
         int nperf = 0;
         auto cell_weights = weights[0];// not need for not(use_well_weights)
         cell_weights = 0.0;
-        assert(this->duneC_.M() == weights.size());
+        assert(this->linSys_.duneC_.M() == weights.size());
-        const int welldof_ind = this->duneC_.M() + this->index_of_well_;
+        const int welldof_ind = this->linSys_.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 ){
+        if (!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) {
+            for (auto colC = this->linSys_.duneC_[0].begin(),
+                      endC = this->linSys_.duneC_[0].end(); colC != endC; ++colC) {
                 const auto row_ind = colC.index();
                 const auto& bw = weights[row_ind];
                 double matel = 0;        
@@ -2244,7 +2247,7 @@ namespace Opm
         bweights[0] = 0.0;
         double diagElem = 0;
         {            
-            if ( use_well_weights ){
+            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
@@ -2252,7 +2255,7 @@ namespace Opm
                 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 = this->linSys_.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;
@@ -2277,7 +2280,7 @@ namespace Opm
                     }
                     bweights[0][blockSz-1] = 0.0;
                     diagElem = 0.0;
-                    const auto& locmat = this->duneD_[0][0]; 
+                    const auto& locmat = this->linSys_.duneD_[0][0];
                     for (size_t i = 0; i < cell_weights.size(); ++i) {
                         diagElem += locmat[i][bhp_var_index]*cell_weights[i];
                     }
@@ -2289,7 +2292,8 @@ namespace Opm
         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) {
+            for (auto colB = this->linSys_.duneB_[0].begin(),
+                      endB = this->linSys_.duneB_[0].end(); colB != endB; ++colB) {
                 const auto col_index = colB.index();
                 const auto& bw = bweights[0];
                 double matel = 0;
@@ -2466,7 +2470,7 @@ namespace Opm
 
         // equation for the water velocity
         const EvalWell eq_wat_vel = this->primary_variables_evaluation_[wat_vel_index] - water_velocity;
-        this->resWell_[0][wat_vel_index] = eq_wat_vel.value();
+        this->linSys_.resWell_[0][wat_vel_index] = eq_wat_vel.value();
 
         const auto& ws = well_state.well(this->index_of_well_);
         const auto& perf_data = ws.perf_data;
@@ -2481,15 +2485,15 @@ namespace Opm
         const EvalWell eq_pskin = this->primary_variables_evaluation_[pskin_index]
                                   - pskin(throughput, this->primary_variables_evaluation_[wat_vel_index], poly_conc, deferred_logger);
 
-        this->resWell_[0][pskin_index] = eq_pskin.value();
+        this->linSys_.resWell_[0][pskin_index] = eq_pskin.value();
         for (int pvIdx = 0; pvIdx < this->numWellEq_; ++pvIdx) {
-            this->duneD_[0][0][wat_vel_index][pvIdx] = eq_wat_vel.derivative(pvIdx+Indices::numEq);
+            this->linSys_.duneD_[0][0][wat_vel_index][pvIdx] = eq_wat_vel.derivative(pvIdx+Indices::numEq);
-            this->duneD_[0][0][pskin_index][pvIdx] = eq_pskin.derivative(pvIdx+Indices::numEq);
+            this->linSys_.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
         for (int pvIdx = 0; pvIdx < Indices::numEq; ++pvIdx) {
-            this->duneB_[0][cell_idx][wat_vel_index][pvIdx] = eq_wat_vel.derivative(pvIdx);
+            this->linSys_.duneB_[0][cell_idx][wat_vel_index][pvIdx] = eq_wat_vel.derivative(pvIdx);
         }
     }