From 21c897fbca321b4990571c314a2cc9fc026dd21c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Tue, 29 Aug 2023 10:43:59 +0200 Subject: [PATCH] Treat Numerical Aquifer Connections as NNCs for Output Purposes Connections between reservoir cells and numerical aquifer cells, or between numerical aquifer cells when multiple such cells define a single numerical aquifer, should always be treated as NNCs for output purposes and end up in the (NNC1,NNC2,TRANNNC) output arrays. To this end, make a special purpose predicate to identify numerical aquifer connections when forming the output NNC and transmissibility arrays and act accordingly in member functions 'computeTrans_()' and 'exportNncStructure_()'. While here, also pick up the NNC transmissibility value from 'globalTrans()' since multiplier operations like MULTREGT might have affected the explicit values entered in the NNC, EDITNNC, and EDITNNCR keywords. This is in preparation of properly incorporating such multipliers in follow-up work. Finally, fix a subtle problem caused by using 'std::abs()' to check for non-zero connections. When accounting for explicit NNCs, it might happen that the final transmissibility would become negative with a sufficiently large absolute value that 'abs(t) > threshold' would be true. This would result in outputting a negative transmissibility value to the NNC arrays which would confuse result processors. --- ebos/eclgenericwriter_impl.hh | 176 +++++++++++++++++++++++++--------- 1 file changed, 132 insertions(+), 44 deletions(-) diff --git a/ebos/eclgenericwriter_impl.hh b/ebos/eclgenericwriter_impl.hh index 9955df1f7..93251b384 100644 --- a/ebos/eclgenericwriter_impl.hh +++ b/ebos/eclgenericwriter_impl.hh @@ -55,6 +55,7 @@ #include #include #include +#include #include #include #include @@ -263,17 +264,17 @@ writeInit(const std::function& map) } template -data::Solution EclGenericWriter:: +data::Solution +EclGenericWriter:: computeTrans_(const std::unordered_map& cartesianToActive, const std::function& map) const { const auto& cartMapper = *equilCartMapper_; const auto& cartDims = cartMapper.cartesianDimensions(); - const int globalSize = cartDims[0] * cartDims[1] * cartDims[2]; auto tranx = data::CellData { UnitSystem::measure::transmissibility, - std::vector(globalSize, 0.0), + std::vector(cartDims[0] * cartDims[1] * cartDims[2], 0.0), data::TargetType::INIT }; @@ -281,15 +282,24 @@ computeTrans_(const std::unordered_map& cartesianToActive, auto tranz = tranx; using GlobalGridView = typename EquilGrid::LeafGridView; - const GlobalGridView& globalGridView = equilGrid_->leafGridView(); using GlobElementMapper = Dune::MultipleCodimMultipleGeomTypeMapper; - GlobElementMapper globalElemMapper(globalGridView, Dune::mcmgElementLayout()); + const GlobalGridView& globalGridView = this->equilGrid_->leafGridView(); + const GlobElementMapper globalElemMapper { globalGridView, Dune::mcmgElementLayout() }; + + auto isNumAquCell = [numAquCell = this->eclState_.aquifer().hasNumericalAquifer() + ? this->eclState_.aquifer().numericalAquifers().allAquiferCellIds() + : std::vector{}] + (const std::size_t cellIdx) + { + return std::binary_search(numAquCell.begin(), numAquCell.end(), cellIdx); + }; for (const auto& elem : elements(globalGridView)) { for (const auto& is : intersections(globalGridView, elem)) { if (!is.neighbor()) continue; // intersection is on the domain boundary + // Not 'const' because remapped if 'map' is non-null. unsigned c1 = globalElemMapper.index(is.inside()); unsigned c2 = globalElemMapper.index(is.outside()); @@ -299,6 +309,15 @@ computeTrans_(const std::unordered_map& cartesianToActive, // Ordering of compressed and uncompressed index should be the same const int cartIdx1 = cartMapper.cartesianIndex( c1 ); const int cartIdx2 = cartMapper.cartesianIndex( c2 ); + + if (isNumAquCell(cartIdx1) || isNumAquCell(cartIdx2)) { + // Connections involving numerical aquifers are always NNCs + // for the purpose of file output. This holds even for + // connections between cells like (I,J,K) and (I+1,J,K) + // which are nominally neighbours in the Cartesian grid. + continue; + } + // Ordering of compressed and uncompressed index should be the same assert(cartIdx1 <= cartIdx2); int gc1 = std::min(cartIdx1, cartIdx2); @@ -334,41 +353,105 @@ computeTrans_(const std::unordered_map& cartesianToActive, } template -std::vector EclGenericWriter:: -exportNncStructure_(const std::unordered_map& cartesianToActive, const std::function& map) const +std::vector +EclGenericWriter:: +exportNncStructure_(const std::unordered_map& cartesianToActive, + const std::function& map) const { - std::size_t nx = eclState_.getInputGrid().getNX(); - std::size_t ny = eclState_.getInputGrid().getNY(); - auto nncData = eclState_.getInputNNC().input(); - const auto& unitSystem = eclState_.getDeckUnitSystem(); + auto isNumAquCell = [numAquCell = this->eclState_.aquifer().hasNumericalAquifer() + ? this->eclState_.aquifer().numericalAquifers().allAquiferCellIds() + : std::vector{}] + (const std::size_t cellIdx) + { + return std::binary_search(numAquCell.begin(), numAquCell.end(), cellIdx); + }; - for( const auto& entry : nncData ) { - // test whether NNC is not a neighboring connection - // cell2>=cell1 holds due to sortNncAndApplyEditnnc - assert( entry.cell2 >= entry.cell1 ); - auto cellDiff = entry.cell2 - entry.cell1; + auto isNumAquConn = [&isNumAquCell](const std::size_t cellIdx1, + const std::size_t cellIdx2) + { + return isNumAquCell(cellIdx1) || isNumAquCell(cellIdx2); + }; - if (cellDiff != 1 && cellDiff != nx && cellDiff != nx*ny) { - auto tt = unitSystem.from_si(UnitSystem::measure::transmissibility, entry.trans); - // Eclipse ignores NNCs (with EDITNNC applied) that are small. Seems like the threshold is 1.0e-6 - if ( tt >= 1.0e-6 ) - outputNnc_.emplace_back(entry.cell1, entry.cell2, entry.trans); + auto isCartesianNeighbour = [nx = this->eclState_.getInputGrid().getNX(), + ny = this->eclState_.getInputGrid().getNY()] + (const std::size_t cellIdx1, const std::size_t cellIdx2) + { + const auto cellDiff = cellIdx2 - cellIdx1; + + return (cellDiff == 1) + || (cellDiff == nx) + || (cellDiff == nx * ny); + }; + + auto activeCell = [&cartesianToActive](const std::size_t cellIdx) + { + auto pos = cartesianToActive.find(cellIdx); + return (pos == cartesianToActive.end()) ? -1 : pos->second; + }; + + const auto& nncData = this->eclState_.getInputNNC().input(); + const auto& unitSystem = this->eclState_.getDeckUnitSystem(); + + for (const auto& entry : nncData) { + // Ignore most explicit NNCs between otherwise neighbouring cells. + // We keep NNCs that involve cells with numerical aquifers even if + // these might be between neighbouring cells in the Cartesian + // grid--e.g., between cells (I,J,K) and (I+1,J,K). All such + // connections should be written to NNC output arrays provided the + // transmissibility value is sufficiently large. + // + // The condition cell2 >= cell1 holds by construction of nncData. + assert (entry.cell2 >= entry.cell1); + + if (! isCartesianNeighbour(entry.cell1, entry.cell2) || + isNumAquConn(entry.cell1, entry.cell2)) + { + // Pick up transmissibility value from 'globalTrans()' since + // multiplier keywords like MULTREGT might have impacted the + // values entered in primary sources like NNC/EDITNNC/EDITNNCR. + const auto c1 = activeCell(entry.cell1); + const auto c2 = activeCell(entry.cell2); + + if ((c1 < 0) || (c2 < 0)) { + // Connection between inactive cells? Unexpected at this + // level. Might consider 'throw'ing if this happens... + continue; + } + + const auto trans = this->globalTrans().transmissibility(c1, c2); + const auto tt = unitSystem + .from_si(UnitSystem::measure::transmissibility, trans); + + // ECLIPSE ignores NNCs (with EDITNNC/EDITNNCR applied) with + // small transmissibility values. Seems like the threshold is + // 1.0e-6 in output units. + if (std::isnormal(tt) && ! (tt < 1.0e-6)) { + this->outputNnc_.emplace_back(entry.cell1, entry.cell2, trans); + } } } + auto isDirectNeighbours = [&isCartesianNeighbour, &cartesianToActive, + cartDims = &this->cartMapper_.cartesianDimensions()] + (const std::size_t cellIdx1, const std::size_t cellIdx2) + { + return isCartesianNeighbour(cellIdx1, cellIdx2) + || directVerticalNeighbors(*cartDims, cartesianToActive, cellIdx1, cellIdx2); + }; + using GlobalGridView = typename EquilGrid::LeafGridView; - const GlobalGridView& globalGridView = equilGrid_->leafGridView(); using GlobElementMapper = Dune::MultipleCodimMultipleGeomTypeMapper; - GlobElementMapper globalElemMapper(globalGridView, Dune::mcmgElementLayout()); + const GlobalGridView& globalGridView = this->equilGrid_->leafGridView(); + const GlobElementMapper globalElemMapper { globalGridView, Dune::mcmgElementLayout() }; // Cartesian index mapper for the serial I/O grid - const auto& equilCartMapper = *equilCartMapper_; - const auto& cartDims = cartMapper_.cartesianDimensions(); + const auto& equilCartMapper = *equilCartMapper_; for (const auto& elem : elements(globalGridView)) { for (const auto& is : intersections(globalGridView, elem)) { if (!is.neighbor()) continue; // intersection is on the domain boundary + // Not 'const' because remapped if 'map' is non-null. unsigned c1 = globalElemMapper.index(is.inside()); unsigned c2 = globalElemMapper.index(is.outside()); @@ -381,38 +464,43 @@ exportNncStructure_(const std::unordered_map& cartesianToActive, const if ( cc2 < cc1 ) std::swap(cc1, cc2); - auto cellDiff = cc2 - cc1; - // Re-ordering in case of non-empty mapping between equilGrid to grid if (map) { c1 = map(c1); // equilGridToGrid map c2 = map(c2); } - if (cellDiff != 1 && - cellDiff != nx && - 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); - auto candidate = std::lower_bound(nncData.begin(), nncData.end(), NNCdata(cc1, cc2, 0.0)); + if (isNumAquConn(cc1, cc2) || ! isDirectNeighbours(cc1, cc2)) { + // We need to check whether an NNC for this face was also + // specified via the NNC keyword in the deck. + auto t = this->globalTrans().transmissibility(c1, c2); + auto candidate = std::lower_bound(nncData.begin(), nncData.end(), + NNCdata { cc1, cc2, 0.0 }); - while ( candidate != nncData.end() && candidate->cell1 == cc1 - && candidate->cell2 == cc2) { + while ((candidate != nncData.end()) && + (candidate->cell1 == cc1) && + (candidate->cell2 == cc2)) + { t -= candidate->trans; ++candidate; } - // eclipse ignores NNCs with zero transmissibility (different threshold than for NNC - // with corresponding EDITNNC above). In addition we do set small transmissibilties - // to zero when setting up the simulator. These will be ignored here, too. - auto tt = unitSystem.from_si(UnitSystem::measure::transmissibility, std::abs(t)); - if ( tt > 1e-12 ) - outputNnc_.push_back({cc1, cc2, t}); + + // ECLIPSE ignores NNCs with zero transmissibility + // (different threshold than for NNC with corresponding + // EDITNNC above). In addition we do set small + // transmissibilities to zero when setting up the simulator. + // These will be ignored here, too. + const auto tt = unitSystem + .from_si(UnitSystem::measure::transmissibility, t); + + if (std::isnormal(tt) && (tt > 1.0e-12)) { + this->outputNnc_.emplace_back(cc1, cc2, t); + } } } } - return outputNnc_; + + return this->outputNnc_; } template