Merge pull request #335 from totto82/minpv-fill

Use average NTG values for merged cells due to MINPV algorithm
This commit is contained in:
Atgeirr Flø Rasmussen 2015-04-13 10:20:36 +02:00
commit e2a920b1fe

View File

@ -109,6 +109,12 @@ namespace Opm
pvol_[cellIdx] *= eclgrid->getCellVolume(cartesianCellIdx);
}
}
// Use volume weighted arithmetic average of the NTG values for
// the cells effected by the current OPM cpgrid process algorithm
// for MINPV. Note that the change does not effect the pore volume calculations
// as the pore volume is currently defaulted to be comparable to ECLIPSE, but
// only the transmissibility calculations.
minPvFillProps_(grid, eclState,ntg);
// Transmissibility
@ -186,6 +192,11 @@ namespace Opm
const double* perm,
Vector &hTrans);
template <class Grid>
void minPvFillProps_(const Grid &grid,
Opm::EclipseStateConstPtr eclState,
std::vector<double> &ntg);
Vector pvol_ ;
Vector trans_;
Vector gpot_ ;
@ -197,6 +208,41 @@ namespace Opm
};
template <class GridType>
inline void DerivedGeology::minPvFillProps_(const GridType &grid,
Opm::EclipseStateConstPtr eclState,
std::vector<double> &ntg)
{
int numCells = Opm::AutoDiffGrid::numCells(grid);
const int* global_cell = Opm::UgGridHelpers::globalCell(grid);
EclipseGridConstPtr eclgrid = eclState->getEclipseGrid();
std::vector<double> porv = eclState->getDoubleGridProperty("PORV")->getData();
for (int cellIdx = 0; cellIdx < numCells; ++cellIdx) {
const int nx = grid.cartdims[0];
const int ny = grid.cartdims[1];
const int cartesianCellIdx = global_cell[cellIdx];
const double cellVolume = eclgrid->getCellVolume(cartesianCellIdx);
ntg[cartesianCellIdx] *= cellVolume;
double totalCellVolume = cellVolume;
// Average properties as long as there exist cells above
// that has pore volume less than the MINPV threshold
int cartesianCellIdxAbove = cartesianCellIdx - nx*ny;
while ( cartesianCellIdxAbove >= 0 &&
porv[cartesianCellIdxAbove] > 0 &&
porv[cartesianCellIdxAbove] < eclgrid->getMinpvValue() ) {
// Volume weighted arithmetic average of NTG
const double cellAboveVolume = eclgrid->getCellVolume(cartesianCellIdxAbove);
totalCellVolume += cellAboveVolume;
ntg[cartesianCellIdx] += ntg[cartesianCellIdxAbove]*cellAboveVolume;
cartesianCellIdxAbove -= nx*ny;
}
ntg[cartesianCellIdx] /= totalCellVolume;
}
}
template <class GridType>
inline void DerivedGeology::multiplyHalfIntersections_(const GridType &grid,