OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

Treat Numerical Aquifer Connections as NNCs for Output Purposes

Browse Source 
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.
This commit is contained in:
Bård Skaflestad
2023-08-29 10:43:59 +02:00
parent c1a6fba113
commit 21c897fbca
1 changed files with 132 additions and 44 deletions
Download Patch File Download Diff File Expand all files Collapse all files

176
ebos/eclgenericwriter_impl.hh
View File

@@ -55,6 +55,7 @@
#include <algorithm>
#include <array>
#include <cassert>
#include <cmath>
#include <functional>
#include <map>
#include <memory>
@@ -263,17 +264,17 @@ writeInit(const std::function<unsigned int(unsigned int)>& map)
}
template<class Grid, class EquilGrid, class GridView, class ElementMapper, class Scalar>
data::Solution EclGenericWriter<Grid,EquilGrid,GridView,ElementMapper,Scalar>::
data::Solution
EclGenericWriter<Grid,EquilGrid,GridView,ElementMapper,Scalar>::
computeTrans_(const std::unordered_map<int,int>& cartesianToActive,
const std::function<unsigned int(unsigned int)>& 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<double>(globalSize, 0.0),
std::vector<double>(cartDims[0] * cartDims[1] * cartDims[2], 0.0),
data::TargetType::INIT
};
@@ -281,15 +282,24 @@ computeTrans_(const std::unordered_map<int,int>& cartesianToActive,
auto tranz = tranx;
using GlobalGridView = typename EquilGrid::LeafGridView;
const GlobalGridView& globalGridView = equilGrid_->leafGridView();
using GlobElementMapper = Dune::MultipleCodimMultipleGeomTypeMapper<GlobalGridView>;
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<std::size_t>{}]
(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<int,int>& 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<int,int>& cartesianToActive,
}
template<class Grid, class EquilGrid, class GridView, class ElementMapper, class Scalar>
std::vector<NNCdata> EclGenericWriter<Grid,EquilGrid,GridView,ElementMapper,Scalar>::
exportNncStructure_(const std::unordered_map<int,int>& cartesianToActive, const std::function<unsigned int(unsigned int)>& map) const
std::vector<NNCdata>
EclGenericWriter<Grid,EquilGrid,GridView,ElementMapper,Scalar>::
exportNncStructure_(const std::unordered_map<int,int>& cartesianToActive,
const std::function<unsigned int(unsigned int)>& 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<std::size_t>{}]
(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<GlobalGridView>;
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<int,int>& 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<class Grid, class EquilGrid, class GridView, class ElementMapper, class Scalar>