From 1c2d3fbcc748a3f21ca892bf4158de507bf58433 Mon Sep 17 00:00:00 2001 From: Markus Blatt Date: Thu, 12 Mar 2020 13:59:46 +0100 Subject: [PATCH] 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 | 19 +++++++++++++++---- opm/simulators/linalg/ISTLSolverEbosCpr.hpp | 5 +---- .../linalg/ISTLSolverEbosFlexible.hpp | 9 ++++++++- .../linalg/findOverlapRowsAndColumns.hpp | 8 +++----- 4 files changed, 27 insertions(+), 14 deletions(-) diff --git a/opm/simulators/linalg/ISTLSolverEbos.hpp b/opm/simulators/linalg/ISTLSolverEbos.hpp index b05fe6f90..fa15b37cd 100644 --- a/opm/simulators/linalg/ISTLSolverEbos.hpp +++ b/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; + 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; + 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> 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); } } diff --git a/opm/simulators/linalg/ISTLSolverEbosCpr.hpp b/opm/simulators/linalg/ISTLSolverEbosCpr.hpp index 40f2bc959..7fb7e2c9f 100644 --- a/opm/simulators/linalg/ISTLSolverEbosCpr.hpp +++ b/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) { diff --git a/opm/simulators/linalg/ISTLSolverEbosFlexible.hpp b/opm/simulators/linalg/ISTLSolverEbosFlexible.hpp index 48445e868..27fde5cd7 100644 --- a/opm/simulators/linalg/ISTLSolverEbosFlexible.hpp +++ b/opm/simulators/linalg/ISTLSolverEbosFlexible.hpp @@ -54,6 +54,7 @@ namespace Opm template 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(); 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; + 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. diff --git a/opm/simulators/linalg/findOverlapRowsAndColumns.hpp b/opm/simulators/linalg/findOverlapRowsAndColumns.hpp index 28eea46d6..c177ec5f7 100644 --- a/opm/simulators/linalg/findOverlapRowsAndColumns.hpp +++ b/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 - void findOverlapAndInterior(const Grid& grid, std::vector& overlapRows, + template + void findOverlapAndInterior(const Grid& grid, const Mapper& mapper, std::vector& overlapRows, std::vector& 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) {