diff --git a/opm/autodiff/GeoProps.hpp b/opm/autodiff/GeoProps.hpp index fa37751e3..7a82c0d08 100644 --- a/opm/autodiff/GeoProps.hpp +++ b/opm/autodiff/GeoProps.hpp @@ -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 + void minPvFillProps_(const Grid &grid, + Opm::EclipseStateConstPtr eclState, + std::vector &ntg); + Vector pvol_ ; Vector trans_; Vector gpot_ ; @@ -197,6 +208,41 @@ namespace Opm }; + template + inline void DerivedGeology::minPvFillProps_(const GridType &grid, + Opm::EclipseStateConstPtr eclState, + std::vector &ntg) + { + + int numCells = Opm::AutoDiffGrid::numCells(grid); + const int* global_cell = Opm::UgGridHelpers::globalCell(grid); + EclipseGridConstPtr eclgrid = eclState->getEclipseGrid(); + std::vector 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 inline void DerivedGeology::multiplyHalfIntersections_(const GridType &grid,