Merge pull request #2361 from akva2/removeNTGfill

remove minpvFillNtg a temporary solution while not supporting standard minpv
This commit is contained in:
Atgeirr Flø Rasmussen 2020-02-24 13:34:43 +01:00 committed by GitHub
commit 8cd9da0d83
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
2 changed files with 1 additions and 77 deletions

View File

@ -2366,7 +2366,6 @@ private:
const auto& simulator = this->simulator();
const auto& vanguard = simulator.vanguard();
const auto& eclState = vanguard.eclState();
const auto& eclGrid = eclState.getInputGrid();
size_t numDof = this->model().numGridDof();
@ -2375,37 +2374,10 @@ private:
const auto& fp = eclState.fieldProps();
const std::vector<double> porvData = fp.porv(true);
const std::vector<int> actnumData = fp.actnum();
int nx = eclGrid.getNX();
int ny = eclGrid.getNY();
for (size_t dofIdx = 0; dofIdx < numDof; ++ dofIdx) {
unsigned cartElemIdx = vanguard.cartesianIndex(dofIdx);
Scalar poreVolume = porvData[cartElemIdx];
// sum up the pore volume of the active cell and all inactive ones above it
// which were disabled due to their pore volume being too small. If energy is
// conserved, cells are not disabled due to a too small pore volume because
// such cells still store and conduct energy.
if (!enableEnergy && eclGrid.getMinpvMode() == Opm::MinpvMode::ModeEnum::OpmFIL) {
const std::vector<Scalar>& minPvVector = eclGrid.getMinpvVector();
for (int aboveElemCartIdx = static_cast<int>(cartElemIdx) - nx*ny;
aboveElemCartIdx >= 0;
aboveElemCartIdx -= nx*ny)
{
if (porvData[aboveElemCartIdx] >= minPvVector[aboveElemCartIdx])
// the cartesian element above exhibits a pore volume which larger or
// equal to the minimum one
break;
Scalar aboveElemVolume = eclGrid.getCellVolume(aboveElemCartIdx);
if (actnumData[aboveElemCartIdx] == 0 && aboveElemVolume > 1e-3)
// stop at explicitly disabled elements, but only if their volume is
// greater than 10^-3 m^3
break;
poreVolume += porvData[aboveElemCartIdx];
}
}
// we define the porosity as the accumulated pore volume divided by the
// geometric volume of the element. Note that -- in pathetic cases -- it can
// be larger than 1.0!

View File

@ -136,8 +136,7 @@ public:
ElementMapper elemMapper(gridView, Dune::mcmgElementLayout());
// get the ntg values, the ntg values are modified for the cells merged with minpv
std::vector<double> ntg;
minPvFillNtg_(ntg);
const std::vector<double>& ntg = eclState.fieldProps().get_global_double("NTG");
unsigned numElements = elemMapper.size();
@ -844,53 +843,6 @@ private:
}
}
void minPvFillNtg_(std::vector<double>& averageNtg) const
{
// compute volume weighted arithmetic average of NTG for
// cells merged as an result of minpv.
const auto& eclState = vanguard_.eclState();
const auto& eclGrid = eclState.getInputGrid();
bool opmfil = eclGrid.getMinpvMode() == Opm::MinpvMode::OpmFIL;
std::vector<double> ntg;
ntg = eclState.fieldProps().get_global_double("NTG");
// just return the unmodified ntg if opmfil is not used
averageNtg = ntg;
if (!opmfil)
return;
const auto& cartMapper = vanguard_.cartesianIndexMapper();
const auto& cartDims = cartMapper.cartesianDimensions();
const auto& actnum = eclState.fieldProps().actnum();
const auto& porv = eclState.fieldProps().porv(true);
assert(dimWorld > 1);
const size_t nxny = cartDims[0] * cartDims[1];
for (size_t cartesianCellIdx = 0; cartesianCellIdx < ntg.size(); ++cartesianCellIdx) {
// use the original ntg values for the inactive cells
if (!actnum[cartesianCellIdx])
continue;
// Average properties as long as there exist cells above
// that has pore volume less than the MINPV threshold
const double cellVolume = eclGrid.getCellVolume(cartesianCellIdx);
double ntgCellVolume = ntg[cartesianCellIdx] * cellVolume;
double totalCellVolume = cellVolume;
int cartesianCellIdxAbove = cartesianCellIdx - nxny;
while (cartesianCellIdxAbove >= 0 &&
actnum[cartesianCellIdxAbove] > 0 &&
porv[cartesianCellIdxAbove] < eclGrid.getMinpvVector()[cartesianCellIdxAbove])
{
// Volume weighted arithmetic average of NTG
const double cellAboveVolume = eclGrid.getCellVolume(cartesianCellIdxAbove);
totalCellVolume += cellAboveVolume;
ntgCellVolume += ntg[cartesianCellIdxAbove]*cellAboveVolume;
cartesianCellIdxAbove -= nxny;
}
averageNtg[cartesianCellIdx] = ntgCellVolume / totalCellVolume;
}
}
const Vanguard& vanguard_;
Scalar transmissibilityThreshold_;
std::vector<DimMatrix> permeability_;