Averaging multiple layers of NTG values

This commit is contained in:
Tor Harald Sandve 2015-03-19 12:24:18 +01:00
parent 65199735e2
commit 388ce7548e

View File

@ -222,23 +222,25 @@ namespace Opm
const int nx = grid.cartdims[0]; const int nx = grid.cartdims[0];
const int ny = grid.cartdims[1]; const int ny = grid.cartdims[1];
const int cartesianCellIdx = global_cell[cellIdx]; const int cartesianCellIdx = global_cell[cellIdx];
const int cartesianCellIdxAbove = cartesianCellIdx - nx*ny;
// Average properties when the cell above exist const double cellVolume = eclgrid->getCellVolume(cartesianCellIdx);
// and has porv less than the MINPV threshold ntg[cartesianCellIdx] *= cellVolume;
if ( cartesianCellIdxAbove > 0 && 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] > 0 &&
porv[cartesianCellIdxAbove] < eclgrid->getMinpvValue() ) { porv[cartesianCellIdxAbove] < eclgrid->getMinpvValue() ) {
// Volume weighted arithmetic average of NTG // Volume weighted arithmetic average of NTG
const double cellVolume = eclgrid->getCellVolume(cartesianCellIdx);
const double cellAboveVolume = eclgrid->getCellVolume(cartesianCellIdxAbove); const double cellAboveVolume = eclgrid->getCellVolume(cartesianCellIdxAbove);
const double inv_totalCellVolume = 1/(cellVolume + cellAboveVolume); totalCellVolume += cellAboveVolume;
ntg[cartesianCellIdx] *= cellAboveVolume;
ntg[cartesianCellIdx] += ntg[cartesianCellIdxAbove]*cellAboveVolume; ntg[cartesianCellIdx] += ntg[cartesianCellIdxAbove]*cellAboveVolume;
ntg[cartesianCellIdx] *= inv_totalCellVolume; cartesianCellIdxAbove -= nx*ny;
} }
ntg[cartesianCellIdx] /= totalCellVolume;
} }
} }