Merge pull request #2458 from blattms/use-indices-for-modifying-matrices

Do not use local id set to index matrices.

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)
                 {