-Fixed buges related to cpr with wells.

- change interfaces to have access to pressureVarIndex
- added option in cmake files to propagate checking in dune-istl

opm-simulators-prereqs.cmake

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

opm/simulators/linalg/PressureBhpTransferPolicy.hpp

@@ -183,11 +183,12 @@ namespace Opm
         if (prm_.get<bool>("add_wells")) {
             assert(transpose == false); // not implemented
             fineOperator.addWellPressureEquations(*coarseLevelMatrix_, weights_);
-        }
 #ifndef NDEBUG
-        std::advance(rowCoarse, fineOperator.getNumberOfExtraEquations());
+            std::advance(rowCoarse, fineOperator.getNumberOfExtraEquations());
-        assert(rowCoarse == coarseLevelMatrix_->end());
+            assert(rowCoarse == coarseLevelMatrix_->end());
 #endif
+
+        }
     }
 
     virtual void moveToCoarseLevel(const typename ParentType::FineRangeType& fine) override

opm/simulators/wells/BlackoilWellModel.hpp

@@ -112,7 +112,7 @@ namespace Opm {
                 typename BlackoilWellModelGeneric::GLiftWellStateMap;
             using GLiftEclWells = typename GasLiftGroupInfo::GLiftEclWells;
             using GLiftSyncGroups = typename GasLiftSingleWellGeneric::GLiftSyncGroups;
-
+            constexpr static std::size_t pressureVarIndex = GetPropType<TypeTag, Properties::Indices>::pressureSwitchIdx;
             typedef typename BaseAuxiliaryModule<TypeTag>::NeighborSet NeighborSet;
 
             static const int numEq = Indices::numEq;
@@ -299,7 +299,7 @@ namespace Opm {
             void addWellPressureEquations(PressureMatrix& jacobian, const BVector& weights) const
             {
                 for (const auto& well : well_container_) {
-                    well->addWellPressureEquations(jacobian, weights);
+                    well->addWellPressureEquations(jacobian, weights, pressureVarIndex);
                 }
             }
 

opm/simulators/wells/StandardWell.hpp

@@ -122,6 +122,7 @@ namespace Opm
         using Eval = typename StdWellEval::Eval;
         using EvalWell = typename StdWellEval::EvalWell;
         using BVectorWell = typename StdWellEval::BVectorWell;
+        using DiagMatrixBlockWellType = typename StdWellEval::DiagMatrixBlockWellType;
 
         StandardWell(const Well& well,
                      const ParallelWellInfo& pw_info,
@@ -183,7 +184,7 @@ namespace Opm
 
         virtual void addWellPressureEquationsStruct(PressureMatrix& mat) const override;
 
-        virtual void addWellPressureEquations(PressureMatrix& mat, const BVector& x) const override;
+        virtual void addWellPressureEquations(PressureMatrix& mat, const BVector& x,const int pressureVarIndex) const override;
 
         // iterate well equations with the specified control until converged
         bool iterateWellEqWithControl(const Simulator& ebosSimulator,

opm/simulators/wells/StandardWell_impl.hpp

@@ -2203,7 +2203,7 @@ namespace Opm
 
     template <typename TypeTag>
     void
-    StandardWell<TypeTag>::addWellPressureEquations(PressureMatrix& jacobian, const BVector& weights) const
+    StandardWell<TypeTag>::addWellPressureEquations(PressureMatrix& jacobian, const BVector& weights,const int pressureVarIndex) const
     {
         // sustem is the pressur variant of
         // We need to change matrx A as follows
@@ -2226,17 +2226,28 @@ namespace Opm
             jacobian[row_ind][welldof_ind] = matel;
         }
         // make quasipes weights for bhp it should be trival
-        using VectorBlockType = typename BVector::block_type;
+        //using VectorBlockType = BVectorWell;
-        VectorBlockType bweights(0.0);
+        //VectorBlockType
+        BVectorWell  bweights(1);
+        size_t blockSz = this->numWellEq_;
+        bweights[0].resize(blockSz);
         double diagElem = 0;
         {
             // const DiagMatrixBlockWellType& invA = invDuneD_[0][0];
-            VectorBlockType rhs(0.0);
+            BVectorWell rhs(1);           
-            rhs[bhp_var_index] = 1.0;
+            rhs[0].resize(blockSz);
-            auto inv_diag_block = this->invDuneD_[0][0];
+            rhs[0][bhp_var_index] = 1.0;
-            auto inv_diag_block_transpose = Opm::wellhelpers::transposeDenseDynMatrix(inv_diag_block);
+            DiagMatrixBlockWellType inv_diag_block = this->invDuneD_[0][0];
+            DiagMatrixBlockWellType inv_diag_block_transpose = Opm::wellhelpers::transposeDenseDynMatrix(inv_diag_block);
             // diag_block_transpose.solve(bweights, rhs);
-            inv_diag_block_transpose.mv(rhs, bweights);
+            //HACK due to template errors
+            for (size_t i = 0; i < blockSz; ++i) {
+                bweights[0][i] = 0;
+                for (size_t j = 0; j < blockSz; ++j) {
+                    bweights[0][i] += inv_diag_block_transpose[i][j]*rhs[0][j];
+                }
+            }
+            //inv_diag_block_transpose.mv(rhs[0], bweights[0]);
             // 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); });
@@ -2247,10 +2258,11 @@ namespace Opm
         // set the matrix elements for well reservoir coupling
         for (auto colB = this->duneB_[0].begin(), endB = this->duneB_[0].end(); colB != endB; ++colB) {
             const auto col_index = colB.index();
-            const auto& bw = bweights;
+            const auto& bw = bweights[0];
             double matel = 0;
             for (size_t i = 0; i < bw.size(); ++i) {
-                matel += (*colB)[bhp_var_index][i] * bw[i];
+                const double w = bw[i];
+                matel += (*colB)[i][pressureVarIndex] * bw[i];
             }
             jacobian[welldof_ind][col_index] = matel;
         }

opm/simulators/wells/WellInterface.hpp

@@ -230,7 +230,7 @@ public:
     }
 
     virtual void addWellPressureEquations(PressureMatrix&,
-                                          const BVector& x) const
+                                          const BVector& x,const int pressureVarIndex) const
     {
     }