From 337c637fa01f330a42d5655be928e4ef88bb60db Mon Sep 17 00:00:00 2001 From: Markus Blatt Date: Thu, 16 May 2019 14:58:49 +0200 Subject: [PATCH] Write transmissibility between direct vertical neighbors into TRANZ. In the case where two were direct vertical neighbors in the grid but not in the underlying cartesian grid (e.g. because of MINPV or pinch outs), we treated them as NNCs and wrote the transmissibilty to TRANNC. With this patch we detect this situation (two neighbor cells with identical i and i and no active cells between them) and do not create an NNC in the eclipse output files but write the transmissibility to TRANZ. --- ebos/eclwriter.hh | 77 ++++++++++++++++++++++++++++++++++++++++++----- 1 file changed, 69 insertions(+), 8 deletions(-) diff --git a/ebos/eclwriter.hh b/ebos/eclwriter.hh index 8f7ddb266..b2cd47be2 100644 --- a/ebos/eclwriter.hh +++ b/ebos/eclwriter.hh @@ -43,6 +43,7 @@ #include #include +#include #include #include @@ -69,6 +70,55 @@ class EclWriter; template class EclOutputBlackOilModule; +/*! + * \brief Detect whether two cells are direct vertical neighbours. + * + * I.e. have the same i and j index and all cartesian cells between them + * along the vertical column are inactive. + * + * \tparam CM The type of the cartesian index mapper. + * \param cartMapper The mapper onto cartesian indices. + * \param cartesianToActive The mapping of cartesian indices to active indices. + * \param smallGlobalIndex The cartesian cell index of the cell with smaller index + * \param largeGlobalIndex The cartesian cell index of the cell with larger index + * \return True if the cells have the same i and j indices and all cartesian cells + * between them are inactive. + */ +inline +bool directVerticalNeighbors(const std::array& cartDims, + const std::unordered_map& cartesianToActive, + int smallGlobalIndex, int largeGlobalIndex) +{ + assert(smallGlobalIndex <= largeGlobalIndex); + std::array ijk1, ijk2; + auto globalToIjk = [cartDims](int gc) { + std::array ijk; + ijk[0] = gc % cartDims[0]; + gc /= cartDims[0]; + ijk[1] = gc % cartDims[1]; + ijk[2] = gc / cartDims[1]; + return ijk; + }; + ijk1 = globalToIjk(smallGlobalIndex); + ijk2 = globalToIjk(largeGlobalIndex); + assert(ijk2[2]>=ijk1[2]); + + if ( ijk1[0] == ijk2[0] && ijk1[1] == ijk2[1] && (ijk2[2] - ijk1[2]) > 1) + { + assert((largeGlobalIndex-smallGlobalIndex)%(cartDims[0]*cartDims[1])==0); + for ( int gi = smallGlobalIndex + cartDims[0] * cartDims[1]; gi < largeGlobalIndex; + gi += cartDims[0] * cartDims[1] ) + { + if ( cartesianToActive.find( gi ) != cartesianToActive.end() ) + { + return false; + } + } + return true; + } else + return false; +} + /*! * \ingroup EclBlackOilSimulator * @@ -145,7 +195,9 @@ public: std::map > integerVectors; if (collectToIORank_.isParallel()) integerVectors.emplace("MPI_RANK", collectToIORank_.globalRanks()); - eclIO_->writeInitial(computeTrans_(), integerVectors, exportNncStructure_()); + auto cartMap = Opm::cartesianToCompressed(globalGrid_.size(0), + Opm::UgGridHelpers::globalCell(globalGrid_)); + eclIO_->writeInitial(computeTrans_(cartMap), integerVectors, exportNncStructure_(cartMap)); } } @@ -305,7 +357,7 @@ private: static bool enableEclOutput_() { return EWOMS_GET_PARAM(TypeTag, bool, EnableEclOutput); } - Opm::data::Solution computeTrans_() const + Opm::data::Solution computeTrans_(const std::unordered_map& cartesianToActive) const { const auto& cartMapper = simulator_.vanguard().cartesianIndexMapper(); const auto& cartDims = cartMapper.cartesianDimensions(); @@ -361,16 +413,23 @@ private: if (c1 > c2) continue; // we only need to handle each connection once, thank you. - + // Ordering of compressed and uncompressed index should be the same + assert(cartesianCellIdx[c1] <= cartesianCellIdx[c2]); int gc1 = std::min(cartesianCellIdx[c1], cartesianCellIdx[c2]); int gc2 = std::max(cartesianCellIdx[c1], cartesianCellIdx[c2]); - if (gc2 - gc1 == 1) + + if (gc2 - gc1 == 1) { tranx.data[gc1] = globalTrans->transmissibility(c1, c2); + continue; // skip other if clauses as they are false, last one needs some computation + } - if (gc2 - gc1 == cartDims[0]) + if (gc2 - gc1 == cartDims[0]) { trany.data[gc1] = globalTrans->transmissibility(c1, c2); + continue; // skipt next if clause as it needs some computation + } - if (gc2 - gc1 == cartDims[0]*cartDims[1]) + if ( gc2 - gc1 == cartDims[0]*cartDims[1] || + directVerticalNeighbors(cartDims, cartesianToActive, gc1, gc2)) tranz.data[gc1] = globalTrans->transmissibility(c1, c2); } } @@ -380,7 +439,7 @@ private: {"TRANZ", tranz}}; } - Opm::NNC exportNncStructure_() const + Opm::NNC exportNncStructure_(const std::unordered_map& cartesianToActive) const { std::size_t nx = eclState().getInputGrid().getNX(); std::size_t ny = eclState().getInputGrid().getNY(); @@ -433,6 +492,7 @@ private: globalTrans = &simulator_.problem().eclTransmissibilities(); } + auto cartDims = simulator_.vanguard().cartesianIndexMapper().cartesianDimensions(); auto elemIt = globalGridView.template begin(); const auto& elemEndIt = globalGridView.template end(); for (; elemIt != elemEndIt; ++ elemIt) { @@ -465,7 +525,8 @@ private: if (cellDiff != 1 && cellDiff != nx && - cellDiff != nx*ny) { + cellDiff != nx*ny && + ! directVerticalNeighbors(cartDims, cartesianToActive, cc1, cc2)) { // We need to check whether an NNC for this face was also specified // via the NNC keyword in the deck (i.e. in the first origNncSize entries. auto t = globalTrans->transmissibility(c1, c2);