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

fix style issues in EclTransmissibility and make the editNnc handling a bit simpler

Browse Source 
the new editNNC code might be slightly slower, but I doubt that the
difference is even measureable for real decks.
This commit is contained in:
Andreas Lauser 2019-02-01 17:33:30 +01:00
parent 0aa9a2c6b8
commit 80e04aa8b4
1 changed files with 113 additions and 122 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

235
ebos/ecltransmissibility.hh
View File

@ -91,15 +91,13 @@ class EclTransmissibility
static const unsigned elemIdxShift = 32; // bits static const unsigned elemIdxShift = 32; // bits
using NncData = Opm::NNCdata;
public: public:
EclTransmissibility(const Vanguard& vanguard) EclTransmissibility(const Vanguard& vanguard)
: vanguard_(vanguard) : vanguard_(vanguard)
{ {
const Opm::UnitSystem& unitSystem = vanguard_.deck().getActiveUnitSystem(); const Opm::UnitSystem& unitSystem = vanguard_.deck().getActiveUnitSystem();
transmissibility_threshold_ = unitSystem.parse("Transmissibility").getSIScaling() * 1e-6; transmissibilityThreshold_ = unitSystem.parse("Transmissibility").getSIScaling() * 1e-6;
} }
/*! /*!
@ -331,11 +329,10 @@ public:
// Then the smallest multiplier is applied. // Then the smallest multiplier is applied.
// Default is to apply the top and bottom multiplier // Default is to apply the top and bottom multiplier
bool useSmallestMultiplier = eclGrid.getMultzOption() == Opm::PinchMode::ModeEnum::ALL; bool useSmallestMultiplier = eclGrid.getMultzOption() == Opm::PinchMode::ModeEnum::ALL;
if (useSmallestMultiplier) { if (useSmallestMultiplier)
applyAllZMultipliers_(trans, insideFaceIdx, insideCartElemIdx, outsideCartElemIdx, transMult, cartDims); applyAllZMultipliers_(trans, insideFaceIdx, insideCartElemIdx, outsideCartElemIdx, transMult, cartDims);
} else { else
applyMultipliers_(trans, insideFaceIdx, insideCartElemIdx, transMult); applyMultipliers_(trans, insideFaceIdx, insideCartElemIdx, transMult);
}
// ... and outside elements // ... and outside elements
applyMultipliers_(trans, outsideFaceIdx, outsideCartElemIdx, transMult); applyMultipliers_(trans, outsideFaceIdx, outsideCartElemIdx, transMult);
@ -385,7 +382,7 @@ public:
int cartElemIdx = vanguard_.cartesianIndexMapper().cartesianIndex(elemIdx); int cartElemIdx = vanguard_.cartesianIndexMapper().cartesianIndex(elemIdx);
globalToLocal[cartElemIdx] = elemIdx; globalToLocal[cartElemIdx] = elemIdx;
} }
applyEditNNCToGridTrans_(elemMapper, globalToLocal); applyEditNncToGridTrans_(elemMapper, globalToLocal);
applyNncToGridTrans_(globalToLocal); applyNncToGridTrans_(globalToLocal);
//remove very small non-neighbouring transmissibilities //remove very small non-neighbouring transmissibilities
@ -432,18 +429,19 @@ public:
private: private:
void removeSmallNonCartesianTransmissibilities_() { void removeSmallNonCartesianTransmissibilities_()
{
const auto& cartMapper = vanguard_.cartesianIndexMapper(); const auto& cartMapper = vanguard_.cartesianIndexMapper();
const auto& cartDims = cartMapper.cartesianDimensions(); const auto& cartDims = cartMapper.cartesianDimensions();
for ( auto&& trans: trans_ ){ for (auto&& trans: trans_) {
if (trans.second < transmissibility_threshold_) { if (trans.second < transmissibilityThreshold_) {
const auto& id = trans.first; const auto& id = trans.first;
const auto& elements = isIdReverse_(id); const auto& elements = isIdReverse_(id);
int gc1 = std::min(cartMapper.cartesianIndex(elements.first), cartMapper.cartesianIndex(elements.second)); int gc1 = std::min(cartMapper.cartesianIndex(elements.first), cartMapper.cartesianIndex(elements.second));
int gc2 = std::max(cartMapper.cartesianIndex(elements.first), cartMapper.cartesianIndex(elements.second)); int gc2 = std::max(cartMapper.cartesianIndex(elements.first), cartMapper.cartesianIndex(elements.second));
// only adjust the NNCs // only adjust the NNCs
if (gc2 - gc1 == 1 || gc2 - gc1 == cartDims[0] || gc2 - gc1 == cartDims[0]*cartDims[1] ) if (gc2 - gc1 == 1 || gc2 - gc1 == cartDims[0] || gc2 - gc1 == cartDims[0]*cartDims[1])
continue; continue;
//remove transmissibilities less than the threshold (by default 1e-6 in the deck's unit system) //remove transmissibilities less than the threshold (by default 1e-6 in the deck's unit system)
@ -452,16 +450,20 @@ private:
} }
} }
void applyAllZMultipliers_(Scalar& trans, unsigned insideFaceIdx, unsigned insideCartElemIdx, unsigned outsideCartElemIdx, void applyAllZMultipliers_(Scalar& trans,
const Opm::TransMult& transMult, const std::array<int, dimWorld>& cartDims){ unsigned insideFaceIdx,
unsigned insideCartElemIdx,
if (insideFaceIdx > 3) { //TOP or BOTTOM unsigned outsideCartElemIdx,
const Opm::TransMult& transMult,
const std::array<int, dimWorld>& cartDims)
{
if (insideFaceIdx > 3) { // top or or bottom
Scalar mult = 1e20; Scalar mult = 1e20;
unsigned cartElemIdx = insideCartElemIdx; unsigned cartElemIdx = insideCartElemIdx;
// pick the smallest multiplier while looking down the pillar untill reaching the other end of the connection // pick the smallest multiplier while looking down the pillar untill reaching the other end of the connection
// for the inbetween cells we apply it from both sides // for the inbetween cells we apply it from both sides
while (cartElemIdx != outsideCartElemIdx) { while (cartElemIdx != outsideCartElemIdx) {
if (insideFaceIdx == 4 || cartElemIdx !=insideCartElemIdx ) if (insideFaceIdx == 4 || cartElemIdx !=insideCartElemIdx )
mult = std::min(mult, transMult.getMultiplier(cartElemIdx, Opm::FaceDir::ZMinus)); mult = std::min(mult, transMult.getMultiplier(cartElemIdx, Opm::FaceDir::ZMinus));
if (insideFaceIdx == 5 || cartElemIdx !=insideCartElemIdx) if (insideFaceIdx == 5 || cartElemIdx !=insideCartElemIdx)
mult = std::min(mult, transMult.getMultiplier(cartElemIdx, Opm::FaceDir::ZPlus)); mult = std::min(mult, transMult.getMultiplier(cartElemIdx, Opm::FaceDir::ZPlus));
@ -469,13 +471,13 @@ private:
cartElemIdx += cartDims[0]*cartDims[1]; cartElemIdx += cartDims[0]*cartDims[1];
} }
trans *= mult; trans *= mult;
} else {
applyMultipliers_(trans, insideFaceIdx, insideCartElemIdx, transMult);
} }
else
applyMultipliers_(trans, insideFaceIdx, insideCartElemIdx, transMult);
} }
void updateFromEclState_(){ void updateFromEclState_()
{
const auto& gridView = vanguard_.gridView(); const auto& gridView = vanguard_.gridView();
const auto& cartMapper = vanguard_.cartesianIndexMapper(); const auto& cartMapper = vanguard_.cartesianIndexMapper();
const auto& cartDims = cartMapper.cartesianDimensions(); const auto& cartDims = cartMapper.cartesianDimensions();
@ -517,31 +519,28 @@ private:
int gc2 = std::max(cartMapper.cartesianIndex(c1), cartMapper.cartesianIndex(c2)); int gc2 = std::max(cartMapper.cartesianIndex(c1), cartMapper.cartesianIndex(c2));
if (gc2 - gc1 == 1) { if (gc2 - gc1 == 1) {
if (inputTranx.deckAssigned()) { if (inputTranx.deckAssigned())
// set simulator internal transmissibilities to values from inputTranx // set simulator internal transmissibilities to values from inputTranx
trans_[isId] = inputTranx.iget(gc1); trans_[isId] = inputTranx.iget(gc1);
} else { else
// Scale transmissibilities with scale factor from inputTranx // Scale transmissibilities with scale factor from inputTranx
trans_[isId] *= inputTranx.iget(gc1); trans_[isId] *= inputTranx.iget(gc1);
}
} }
else if (gc2 - gc1 == cartDims[0]) { else if (gc2 - gc1 == cartDims[0]) {
if (inputTrany.deckAssigned()) { if (inputTrany.deckAssigned())
// set simulator internal transmissibilities to values from inputTrany // set simulator internal transmissibilities to values from inputTrany
trans_[isId] = inputTrany.iget(gc1); trans_[isId] = inputTrany.iget(gc1);
} else { else
// Scale transmissibilities with scale factor from inputTrany // Scale transmissibilities with scale factor from inputTrany
trans_[isId] *= inputTrany.iget(gc1); trans_[isId] *= inputTrany.iget(gc1);
}
} }
else if (gc2 - gc1 == cartDims[0]*cartDims[1]) { else if (gc2 - gc1 == cartDims[0]*cartDims[1]) {
if (inputTranz.deckAssigned()) { if (inputTranz.deckAssigned())
// set simulator internal transmissibilities to values from inputTranz // set simulator internal transmissibilities to values from inputTranz
trans_[isId] = inputTranz.iget(gc1); trans_[isId] = inputTranz.iget(gc1);
} else { else
// Scale transmissibilities with scale factor from inputTranz // Scale transmissibilities with scale factor from inputTranz
trans_[isId] *= inputTranz.iget(gc1); trans_[isId] *= inputTranz.iget(gc1);
}
} }
//else.. We don't support modification of NNC at the moment. //else.. We don't support modification of NNC at the moment.
} }
@ -550,15 +549,15 @@ private:
template <class Intersection> template <class Intersection>
void computeFaceProperties( const Intersection& intersection, void computeFaceProperties(const Intersection& intersection,
const int insideElemIdx, const int insideElemIdx,
const int insideFaceIdx, const int insideFaceIdx,
const int outsideElemIdx, const int outsideElemIdx,
const int outsideFaceIdx, const int outsideFaceIdx,
DimVector& faceCenterInside, DimVector& faceCenterInside,
DimVector& faceCenterOutside, DimVector& faceCenterOutside,
DimVector& faceAreaNormal, DimVector& faceAreaNormal,
/*isCpGrid=*/std::false_type) const /*isCpGrid=*/std::false_type) const
{ {
// default implementation for DUNE grids // default implementation for DUNE grids
const auto& geometry = intersection.geometry(); const auto& geometry = intersection.geometry();
@ -570,19 +569,19 @@ private:
} }
template <class Intersection> template <class Intersection>
void computeFaceProperties( const Intersection& intersection, void computeFaceProperties(const Intersection& intersection,
const int insideElemIdx, const int insideElemIdx,
const int insideFaceIdx, const int insideFaceIdx,
const int outsideElemIdx, const int outsideElemIdx,
const int outsideFaceIdx, const int outsideFaceIdx,
DimVector& faceCenterInside, DimVector& faceCenterInside,
DimVector& faceCenterOutside, DimVector& faceCenterOutside,
DimVector& faceAreaNormal, DimVector& faceAreaNormal,
/*isCpGrid=*/std::true_type) const /*isCpGrid=*/std::true_type) const
{ {
int faceIdx = intersection.id(); int faceIdx = intersection.id();
faceCenterInside = vanguard_.grid().faceCenterEcl(insideElemIdx,insideFaceIdx); faceCenterInside = vanguard_.grid().faceCenterEcl(insideElemIdx, insideFaceIdx);
faceCenterOutside = vanguard_.grid().faceCenterEcl(outsideElemIdx,outsideFaceIdx); faceCenterOutside = vanguard_.grid().faceCenterEcl(outsideElemIdx, outsideFaceIdx);
faceAreaNormal = vanguard_.grid().faceAreaNormalEcl(faceIdx); faceAreaNormal = vanguard_.grid().faceAreaNormalEcl(faceIdx);
} }
@ -600,50 +599,53 @@ private:
* and the second the NNCs not resembled by faces of the grid. NNCs specified for * and the second the NNCs not resembled by faces of the grid. NNCs specified for
* inactive cells are omitted in these vectors. * inactive cells are omitted in these vectors.
*/ */
std::tuple<std::vector<NncData>, std::vector<NncData> > std::tuple<std::vector<Opm::NNCdata>, std::vector<Opm::NNCdata> >
applyNncToGridTrans_(const std::vector<int>& cartesianToCompressed) applyNncToGridTrans_(const std::vector<int>& cartesianToCompressed)
{ {
// First scale NNCs with EDITNNC. // First scale NNCs with EDITNNC.
std::vector<NncData> processedNnc, unprocessedNnc; std::vector<Opm::NNCdata> unprocessedNnc;
std::vector<Opm::NNCdata> processedNnc;
const auto& nnc = vanguard_.eclState().getInputNNC(); const auto& nnc = vanguard_.eclState().getInputNNC();
if ( ! nnc.hasNNC() ) if (!nnc.hasNNC())
{
return make_tuple(processedNnc, unprocessedNnc); return make_tuple(processedNnc, unprocessedNnc);
}
auto nncData = nnc.nncdata(); auto nncData = nnc.nncdata();
auto editnncData = vanguard_.eclState().getInputEDITNNC().data(); auto editnncData = vanguard_.eclState().getInputEDITNNC().data();
auto nncLess = [](const NncData& d1, const NncData& d2){ auto nncLess =
return ( d1.cell1 < d2.cell1 ) || [](const Opm::NNCdata& d1, const Opm::NNCdata& d2)
( d1.cell1 == d2.cell1 && d1.cell2 < d2.cell2 ); {
}; return
(d1.cell1 < d2.cell1)
|| (d1.cell1 == d2.cell1 && d1.cell2 < d2.cell2);
};
std::sort(nncData.begin(), nncData.end(), nncLess); std::sort(nncData.begin(), nncData.end(), nncLess);
auto candidate = nncData.begin(); auto candidate = nncData.begin();
for ( const auto& edit: editnncData ) for (const auto& edit: editnncData) {
{ auto printNncWarning =
auto printNncWarning = [](int c1, int c2){ [](int c1, int c2)
std::ostringstream sstr; {
sstr << "Cannot edit NNC from " << c1 << " to " << c2 std::ostringstream sstr;
<< " as it does not exist"; sstr << "Cannot edit NNC from " << c1 << " to " << c2
Opm::OpmLog::warning(sstr.str()); << " as it does not exist";
}; Opm::OpmLog::warning(sstr.str());
if ( candidate == nncData.end() ) };
{ if (candidate == nncData.end()) {
// no more NNCs left // no more NNCs left
printNncWarning(edit.cell1, edit.cell2); printNncWarning(edit.cell1, edit.cell2);
continue; continue;
} }
if ( candidate->cell1 != edit.cell1 || candidate->cell2 != edit.cell2 ) if (candidate->cell1 != edit.cell1 || candidate->cell2 != edit.cell2) {
{ candidate = std::lower_bound(candidate, nncData.end(), Opm::NNCdata(edit.cell1, edit.cell2, 0), nncLess);
candidate = std::lower_bound(candidate, nncData.end(), NncData(edit.cell1, edit.cell2, 0), nncLess); if (candidate == nncData.end()) {
if ( candidate == nncData.end() )
{
// no more NNCs left // no more NNCs left
printNncWarning(edit.cell1, edit.cell2); printNncWarning(edit.cell1, edit.cell2);
continue; continue;
} }
} }
auto firstCandidate = candidate; auto firstCandidate = candidate;
while ( candidate != nncData.end() && candidate->cell1 == edit.cell1 && candidate->cell2 == edit.cell2 ) while (candidate != nncData.end()
&& candidate->cell1 == edit.cell1
&& candidate->cell2 == edit.cell2)
{ {
candidate->trans *= edit.trans; candidate->trans *= edit.trans;
++candidate; ++candidate;
@ -653,40 +655,30 @@ private:
candidate = firstCandidate; candidate = firstCandidate;
} }
for (const auto& nncEntry : nnc.nncdata()) for (const auto& nncEntry : nnc.nncdata()) {
{
auto c1 = nncEntry.cell1; auto c1 = nncEntry.cell1;
auto c2 = nncEntry.cell2; auto c2 = nncEntry.cell2;
auto low = cartesianToCompressed[c1]; auto low = cartesianToCompressed[c1];
auto high = cartesianToCompressed[c2]; auto high = cartesianToCompressed[c2];
if ( low > high) if (low > high)
{
std::swap(low, high); std::swap(low, high);
}
if ( low == -1 && high == -1 ) if (low == -1 && high == -1)
{
// Silently discard as it is not between active cells // Silently discard as it is not between active cells
continue; continue;
}
if ( low == -1 || high == -1) if (low == -1 || high == -1)
{
OPM_THROW(std::logic_error, "NNC between active and inactive cells (" << OPM_THROW(std::logic_error, "NNC between active and inactive cells (" <<
low << " -> " << high); low << " -> " << high);
}
auto candidate = trans_.find(isId_(low, high)); auto candidate = trans_.find(isId_(low, high));
if ( candidate == trans_.end() ) if (candidate == trans_.end())
{
// This NNC is not resembled by the grid. Save it for later // This NNC is not resembled by the grid. Save it for later
// processing with local cell values // processing with local cell values
unprocessedNnc.push_back({c1, c2, nncEntry.trans}); unprocessedNnc.push_back({c1, c2, nncEntry.trans});
} else {
else
{
// NNC is represented by the grid and might be a neighboring connection // NNC is represented by the grid and might be a neighboring connection
// In this case the transmissibilty is added to the value already // In this case the transmissibilty is added to the value already
// set or computed. // set or computed.
@ -698,41 +690,37 @@ private:
} }
/// \brief Multiplies the grid transmissibilities according to EDITNNC. /// \brief Multiplies the grid transmissibilities according to EDITNNC.
void applyEditNNCToGridTrans_(const ElementMapper& elementMapper, void applyEditNncToGridTrans_(const ElementMapper& elementMapper,
const std::vector<int>& globalToLocal) const std::vector<int>& globalToLocal)
{ {
const auto& editNNC = vanguard_.eclState().getInputEDITNNC(); const auto& editNnc = vanguard_.eclState().getInputEDITNNC();
if ( editNNC.empty() ) if (editNnc.empty())
{
return; return;
}
// editNNC is supposed to only reference non-neighboring connections and not // editNnc is supposed to only reference non-neighboring connections and not
// neighboring connections. Use all entries for scaling if there is an NNC. // neighboring connections. Use all entries for scaling if there is an NNC.
// variable nnc incremented in loop body. // variable nnc incremented in loop body.
for (auto nnc = editNNC.data().begin(), end = editNNC.data().end(); nnc != end; ) auto nnc = editNnc.data().begin();
{ auto end = editNnc.data().end();
auto c1 = nnc->cell1, c2 = nnc->cell2; while (nnc != end) {
auto low = globalToLocal[c1], high = globalToLocal[c2]; auto c1 = nnc->cell1;
if ( low > high) auto c2 = nnc->cell2;
{ auto low = globalToLocal[c1];
auto high = globalToLocal[c2];
if (low > high)
std::swap(low, high); std::swap(low, high);
}
auto candidate = trans_.find(isId_(low, high));
if ( candidate == trans_.end() ) auto candidate = trans_.find(isId_(low, high));
{ if (candidate == trans_.end()) {
std::ostringstream sstr; std::ostringstream sstr;
sstr << "Cannot edit NNC from " << c1 << " to " << c2 sstr << "Cannot edit NNC from " << c1 << " to " << c2
<< " as it does not exist"; << " as it does not exist";
Opm::OpmLog::warning(sstr.str()); Opm::OpmLog::warning(sstr.str());
++nnc; ++nnc;
} }
else else {
{
// NNC exists // NNC exists
while ( nnc!= end && c1==nnc->cell1 && c2==nnc->cell2 ) while (nnc!= end && c1==nnc->cell1 && c2==nnc->cell2) {
{
candidate->second *= nnc->trans; candidate->second *= nnc->trans;
++nnc; ++nnc;
} }
@ -831,7 +819,9 @@ private:
return x; return x;
} }
void applyMultipliers_(Scalar& trans, unsigned faceIdx, unsigned cartElemIdx, void applyMultipliers_(Scalar& trans,
unsigned faceIdx,
unsigned cartElemIdx,
const Opm::TransMult& transMult) const const Opm::TransMult& transMult) const
{ {
// apply multiplyer for the transmissibility of the face. (the // apply multiplyer for the transmissibility of the face. (the
@ -861,7 +851,9 @@ private:
} }
} }
void applyNtg_(Scalar& trans, unsigned faceIdx, unsigned cartElemIdx, void applyNtg_(Scalar& trans,
unsigned faceIdx,
unsigned cartElemIdx,
const std::vector<double>& ntg) const const std::vector<double>& ntg) const
{ {
// apply multiplyer for the transmissibility of the face. (the // apply multiplyer for the transmissibility of the face. (the
@ -886,9 +878,8 @@ private:
} }
} }
void minPvFillNtg_(std::vector<double>& averageNtg) const { void minPvFillNtg_(std::vector<double>& averageNtg) const
{
// compute volume weighted arithmetic average of NTG for // compute volume weighted arithmetic average of NTG for
// cells merged as an result of minpv. // cells merged as an result of minpv.
const auto& eclState = vanguard_.eclState(); const auto& eclState = vanguard_.eclState();
@ -911,8 +902,7 @@ private:
const auto& cartDims = cartMapper.cartesianDimensions(); const auto& cartDims = cartMapper.cartesianDimensions();
assert(dimWorld > 1); assert(dimWorld > 1);
const size_t nxny = cartDims[0] * cartDims[1]; const size_t nxny = cartDims[0] * cartDims[1];
for (size_t cartesianCellIdx = 0; cartesianCellIdx < ntg.size(); ++cartesianCellIdx) for (size_t cartesianCellIdx = 0; cartesianCellIdx < ntg.size(); ++cartesianCellIdx) {
{
// use the original ntg values for the inactive cells // use the original ntg values for the inactive cells
if (!actnum[cartesianCellIdx]) if (!actnum[cartesianCellIdx])
continue; continue;
@ -923,22 +913,23 @@ private:
double ntgCellVolume = ntg[cartesianCellIdx] * cellVolume; double ntgCellVolume = ntg[cartesianCellIdx] * cellVolume;
double totalCellVolume = cellVolume; double totalCellVolume = cellVolume;
int cartesianCellIdxAbove = cartesianCellIdx - nxny; int cartesianCellIdxAbove = cartesianCellIdx - nxny;
while ( cartesianCellIdxAbove >= 0 && while (cartesianCellIdxAbove >= 0 &&
actnum[cartesianCellIdxAbove] > 0 && actnum[cartesianCellIdxAbove] > 0 &&
porv[cartesianCellIdxAbove] < eclGrid.getMinpvVector()[cartesianCellIdxAbove] ) { porv[cartesianCellIdxAbove] < eclGrid.getMinpvVector()[cartesianCellIdxAbove])
{
// Volume weighted arithmetic average of NTG // Volume weighted arithmetic average of NTG
const double cellAboveVolume = eclGrid.getCellVolume(cartesianCellIdxAbove); const double cellAboveVolume = eclGrid.getCellVolume(cartesianCellIdxAbove);
totalCellVolume += cellAboveVolume; totalCellVolume += cellAboveVolume;
ntgCellVolume += ntg[cartesianCellIdxAbove]*cellAboveVolume; ntgCellVolume += ntg[cartesianCellIdxAbove]*cellAboveVolume;
cartesianCellIdxAbove -= nxny; cartesianCellIdxAbove -= nxny;
} }
averageNtg[cartesianCellIdx] = ntgCellVolume / totalCellVolume; averageNtg[cartesianCellIdx] = ntgCellVolume / totalCellVolume;
} }
} }
const Vanguard& vanguard_; const Vanguard& vanguard_;
Scalar transmissibility_threshold_; Scalar transmissibilityThreshold_;
std::vector<DimMatrix> permeability_; std::vector<DimMatrix> permeability_;
std::unordered_map<std::uint64_t, Scalar> trans_; std::unordered_map<std::uint64_t, Scalar> trans_;
std::map<std::pair<unsigned, unsigned>, Scalar> transBoundary_; std::map<std::pair<unsigned, unsigned>, Scalar> transBoundary_;