Do not use local id set to index matrices.

That the local ids were consecutive and starting from 0 was just
a coincidence and they should never be used to access linear systems
or vectors. This commit fixes this by using the correct mappers instead.

Note the we removed some computations from the constructor of
ISTLSolverEbosCpr as it inherits from ISTLSolverEbos and the operations
already happnen in constructor of the base class.

opm/simulators/linalg/ISTLSolverEbos.hpp

@@ -364,7 +364,13 @@ protected:
 
             if (!ownersFirst_ || parameters_.linear_solver_use_amg_  || parameters_.use_cpr_ ) {
                 detail::setWellConnections(gridForConn, wellsForConn, useWellConn, wellConnectionsGraph_);
-                detail::findOverlapAndInterior(gridForConn, overlapRows_, interiorRows_);
+                // For some reason simulator_.model().elementMapper() is not initialized at this stage
+                // Hence const auto& elemMapper = simulator_.model().elementMapper(); does not work.
+                // Set it up manually
+                using ElementMapper =
+                    Dune::MultipleCodimMultipleGeomTypeMapper<GridView>;
+                ElementMapper elemMapper(simulator_.vanguard().grid().leafGridView(), Dune::mcmgElementLayout());
+                detail::findOverlapAndInterior(gridForConn, elemMapper, overlapRows_, interiorRows_);
                 if (gridForConn.comm().size() > 1) {
 
                     noGhostAdjacency();
@@ -787,6 +793,12 @@ protected:
         void noGhostAdjacency()
         {
             const auto& grid = simulator_.vanguard().grid();
+            // For some reason simulator_.model().elementMapper() is not initialized at this stage.
+            // Hence const auto& elemMapper = simulator_.model().elementMapper(); does not work.
+            // Set it up manually
+            using ElementMapper =
+                Dune::MultipleCodimMultipleGeomTypeMapper<GridView>;
+            ElementMapper elemMapper(simulator_.vanguard().grid().leafGridView(), Dune::mcmgElementLayout());
             typedef typename Matrix::size_type size_type;
             size_type numCells = grid.size( 0 );
             noGhostMat_.reset(new Matrix(numCells, numCells, Matrix::random));
@@ -794,7 +806,6 @@ protected:
             std::vector<std::set<size_type>> pattern;
             pattern.resize(numCells);
 
-            const auto& lid = grid.localIdSet();
             const auto& gridView = grid.leafGridView();
             auto elemIt = gridView.template begin<0>();
             const auto& elemEndIt = gridView.template end<0>();
@@ -803,7 +814,7 @@ protected:
             for (; elemIt != elemEndIt; ++elemIt)
             {
                 const auto& elem = *elemIt;
-                size_type idx = lid.id(elem);
+                size_type idx = elemMapper.index(elem);
                 pattern[idx].insert(idx);
 
                 // Add well non-zero connections
@@ -822,7 +833,7 @@ protected:
                         //check if face has neighbor
                         if (is->neighbor())
                         {
-                            size_type nid = lid.id(is->outside());
+                            size_type nid = elemMapper.index(is->outside());
                             pattern[idx].insert(nid);
                         }
                     }

opm/simulators/linalg/ISTLSolverEbosCpr.hpp

@@ -97,10 +97,7 @@ namespace Opm
         ///                                with dune-istl the information about the parallelization.
         explicit ISTLSolverEbosCpr(const Simulator& simulator)
             : SuperClass(simulator), oldMat()
-        {
+        {}
-            extractParallelGridInformationToISTL(this->simulator_.vanguard().grid(), this->parallelInformation_);
-            detail::findOverlapAndInterior(this->simulator_.vanguard().grid(), this->overlapRows_, this->interiorRows_);
-        }
 
         void prepare(const SparseMatrixAdapter& M, Vector& b)
         {

opm/simulators/linalg/ISTLSolverEbosFlexible.hpp

@@ -54,6 +54,7 @@ namespace Opm
 template <class TypeTag>
 class ISTLSolverEbosFlexible
 {
+    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
     using SparseMatrixAdapter = typename GET_PROP_TYPE(TypeTag, SparseMatrixAdapter);
     using VectorType = typename GET_PROP_TYPE(TypeTag, GlobalEqVector);
     using Simulator = typename GET_PROP_TYPE(TypeTag, Simulator);
@@ -77,7 +78,13 @@ public:
         parameters_.template init<TypeTag>();
         prm_ = setupPropertyTree(parameters_);
         extractParallelGridInformationToISTL(simulator_.vanguard().grid(), parallelInformation_);
-        detail::findOverlapAndInterior(simulator_.vanguard().grid(), overlapRows_, interiorRows_);
+        // For some reason simulator_.model().elementMapper() is not initialized at this stage
+        // Hence const auto& elemMapper = simulator_.model().elementMapper(); does not work.
+        // Set it up manually
+        using ElementMapper =
+            Dune::MultipleCodimMultipleGeomTypeMapper<GridView>;
+        ElementMapper elemMapper(simulator_.vanguard().grid().leafGridView(), Dune::mcmgElementLayout());
+        detail::findOverlapAndInterior(simulator_.vanguard().grid(), elemMapper, overlapRows_, interiorRows_);
 #if HAVE_MPI
         if (parallelInformation_.type() == typeid(ParallelISTLInformation)) {
             // Parallel case.

opm/simulators/linalg/findOverlapRowsAndColumns.hpp

@@ -82,16 +82,14 @@ namespace detail
     /// \param grid The grid where we look for overlap cells.
     /// \param overlapRows List where overlap rows are stored.
     /// \param interiorRows List where overlap rows are stored.
-    template<class Grid>
+    template<class Grid, class Mapper>
-    void findOverlapAndInterior(const Grid& grid, std::vector<int>& overlapRows,
+    void findOverlapAndInterior(const Grid& grid, const Mapper& mapper, std::vector<int>& overlapRows,
                                 std::vector<int>& interiorRows)
     {
         //only relevant in parallel case.
         if ( grid.comm().size() > 1)
         {
             //Numbering of cells
-            const auto& lid = grid.localIdSet();
-
             const auto& gridView = grid.leafGridView();
             auto elemIt = gridView.template begin<0>();
             const auto& elemEndIt = gridView.template end<0>();
@@ -100,7 +98,7 @@ namespace detail
             for (; elemIt != elemEndIt; ++elemIt)
             {
                 const auto& elem = *elemIt;
-                int lcell = lid.id(elem);
+                int lcell = mapper.index(elem);
 
                 if (elem.partitionType() != Dune::InteriorEntity)
                 {