20 #ifndef OPM_GEOPROPS_HEADER_INCLUDED 21 #define OPM_GEOPROPS_HEADER_INCLUDED 23 #include <opm/core/grid.h> 24 #include <opm/autodiff/GridHelpers.hpp> 25 #include <opm/common/ErrorMacros.hpp> 27 #include <opm/core/pressure/tpfa/TransTpfa.hpp> 29 #include <opm/parser/eclipse/EclipseState/EclipseState.hpp> 30 #include <opm/parser/eclipse/EclipseState/Grid/EclipseGrid.hpp> 31 #include <opm/parser/eclipse/EclipseState/Grid/GridProperty.hpp> 32 #include <opm/parser/eclipse/EclipseState/Grid/NNC.hpp> 33 #include <opm/parser/eclipse/EclipseState/Grid/TransMult.hpp> 34 #include <opm/core/grid/PinchProcessor.hpp> 35 #include <opm/common/utility/platform_dependent/disable_warnings.h> 36 #include <opm/output/data/Cells.hpp> 37 #include <opm/output/data/Solution.hpp> 39 #include <Eigen/Eigen> 42 #include <dune/common/version.hh> 43 #include <dune/grid/CpGrid.hpp> 44 #include <dune/grid/common/mcmgmapper.hh> 47 #include <opm/common/utility/platform_dependent/reenable_warnings.h> 62 typedef Eigen::ArrayXd Vector;
66 template <
class Props,
class Gr
id>
69 const EclipseState& eclState,
70 const bool use_local_perm,
71 const double* grav = 0
74 : pvol_ (
Opm::AutoDiffGrid::numCells(grid))
75 , trans_(
Opm::AutoDiffGrid::numFaces(grid))
76 , gpot_ (Vector::Zero(
Opm::AutoDiffGrid::cell2Faces(grid).noEntries(), 1))
77 , z_(
Opm::AutoDiffGrid::numCells(grid))
78 , use_local_perm_(use_local_perm)
80 update(grid, props, eclState, grav);
84 template <
class Props,
class Gr
id>
87 const EclipseState& eclState,
91 int numCells = AutoDiffGrid::numCells(grid);
92 int numFaces = AutoDiffGrid::numFaces(grid);
93 const int *cartDims = AutoDiffGrid::cartDims(grid);
94 int numCartesianCells =
100 std::vector<double> multpv(numCartesianCells, 1.0);
101 const auto& eclProps = eclState.get3DProperties();
102 if (eclProps.hasDeckDoubleGridProperty(
"MULTPV")) {
103 multpv = eclProps.getDoubleGridProperty(
"MULTPV").getData();
107 std::vector<double> ntg(numCartesianCells, 1.0);
108 if (eclProps.hasDeckDoubleGridProperty(
"NTG")) {
109 ntg = eclProps.getDoubleGridProperty(
"NTG").getData();
113 const auto& eclgrid = eclState.getInputGrid();
116 computePoreVolume_(grid, eclState);
119 nnc_ = eclState.getInputNNC();
122 Vector htrans(AutoDiffGrid::numCellFaces(grid));
123 Grid* ug =
const_cast<Grid*
>(& grid);
125 if (! use_local_perm_) {
126 tpfa_htrans_compute(ug, props.permeability(), htrans.data());
129 tpfa_loc_trans_compute_(grid,eclgrid, props.permeability(),htrans);
137 bool opmfil = eclgrid.getMinpvMode() == MinpvMode::ModeEnum::OpmFIL;
141 minPvFillProps_(grid, eclState, ntg);
144 std::vector<double> mult;
145 multiplyHalfIntersections_(grid, eclState, ntg, htrans, mult);
147 if (!opmfil && eclgrid.isPinchActive()) {
149 pinchProcess_(grid, eclState, htrans, numCells);
154 tpfa_trans_compute(ug, htrans.data(), trans_.data());
158 for (
int faceIdx = 0; faceIdx < numFaces; faceIdx++) {
159 trans_[faceIdx] *= mult[faceIdx];
163 noncartesian_ = nnc_;
164 exportNncStructure(grid);
167 for (
int c = 0; c<numCells; ++c){
168 z_[c] = Opm::UgGridHelpers::cellCenterDepth(grid, c);
173 std::fill(gravity_, gravity_ + 3, 0.0);
175 const typename Vector::Index nd = AutoDiffGrid::dimensions(grid);
177 Cell2Faces c2f=AutoDiffGrid::cell2Faces(grid);
180 for (
typename Vector::Index c = 0; c < numCells; ++c) {
181 const double*
const cc = AutoDiffGrid::cellCentroid(grid, c);
183 typename Cell2Faces::row_type faces=c2f[c];
184 typedef typename Cell2Faces::row_type::iterator Iter;
186 for (Iter f=faces.begin(), end=faces.end(); f!=end; ++f, ++i) {
187 auto fc = AutoDiffGrid::faceCentroid(grid, *f);
189 for (
typename Vector::Index d = 0; d < nd; ++d) {
190 gpot_[i] += grav[d] * (fc[d] - cc[d]);
194 std::copy(grav, grav + nd, gravity_);
201 const Vector& poreVolume()
const {
return pvol_ ;}
202 const Vector& transmissibility()
const {
return trans_ ;}
203 const Vector& gravityPotential()
const {
return gpot_ ;}
204 const Vector& z()
const {
return z_ ;}
205 const double* gravity()
const {
return gravity_;}
206 Vector& poreVolume() {
return pvol_ ;}
207 Vector& transmissibility() {
return trans_ ;}
208 const NNC& nnc()
const {
return nnc_;}
209 const NNC& nonCartesianConnections()
const {
return noncartesian_;}
231 template <
class Gr
id>
232 data::Solution
simProps(
const Grid& grid )
const {
233 using namespace UgGridHelpers;
234 const int* dims = cartDims( grid );
235 const int globalSize = dims[0] * dims[1] * dims[2];
236 const auto& trans = this->transmissibility( );
238 data::CellData tranx = {UnitSystem::measure::transmissibility, std::vector<double>( globalSize ), data::TargetType::INIT};
239 data::CellData trany = {UnitSystem::measure::transmissibility, std::vector<double>( globalSize ), data::TargetType::INIT};
240 data::CellData tranz = {UnitSystem::measure::transmissibility, std::vector<double>( globalSize ), data::TargetType::INIT};
242 size_t num_faces = numFaces(grid);
243 auto fc = faceCells(grid);
244 for (
size_t i = 0; i < num_faces; ++i) {
245 auto c1 = std::min( fc(i,0) , fc(i,1));
246 auto c2 = std::max( fc(i,0) , fc(i,1));
248 if (c1 == -1 || c2 == -1)
251 c1 = globalCell(grid) ? globalCell(grid)[c1] : c1;
252 c2 = globalCell(grid) ? globalCell(grid)[c2] : c2;
254 if ((c2 - c1) == 1) {
255 tranx.data[c1] = trans[i];
258 if ((c2 - c1) == dims[0]) {
259 trany.data[c1] = trans[i];
262 if ((c2 - c1) == dims[0]*dims[1]) {
263 tranz.data[c1] = trans[i];
267 return { {
"TRANX" , tranx},
269 {
"TRANZ" , tranz } };
274 template <
class Gr
id>
275 void multiplyHalfIntersections_(
const Grid &grid,
276 const EclipseState& eclState,
277 const std::vector<double> &ntg,
278 Vector &halfIntersectTransmissibility,
279 std::vector<double> &intersectionTransMult);
281 template <
class Gr
id>
282 void tpfa_loc_trans_compute_(
const Grid &grid,
283 const EclipseGrid& eclGrid,
287 template <
class Gr
id>
288 void minPvFillProps_(
const Grid &grid,
289 const EclipseState& eclState,
290 std::vector<double> &ntg);
292 template <
class Gr
idType>
293 void computePoreVolume_(
const GridType &grid,
294 const EclipseState& eclState)
296 int numCells = Opm::AutoDiffGrid::numCells(grid);
297 const int* globalCell = Opm::UgGridHelpers::globalCell(grid);
298 const auto& eclGrid = eclState.getInputGrid();
299 const int nx = eclGrid.getNX();
300 const int ny = eclGrid.getNY();
303 const std::vector<double>& porvData =
304 eclState.get3DProperties().getDoubleGridProperty(
"PORV").getData();
305 pvol_.resize(numCells);
308 const std::vector<int>& actnumData =
309 eclState.get3DProperties().getIntGridProperty(
"ACTNUM").getData();
312 for (
int cellIdx = 0; cellIdx < numCells; ++cellIdx) {
313 const int cellCartIdx = globalCell[cellIdx];
315 double cellPoreVolume = porvData[cellCartIdx];
317 if (eclGrid.getMinpvMode() == MinpvMode::ModeEnum::OpmFIL) {
320 for (
int aboveCellCartIdx = cellCartIdx - nx*ny;
321 aboveCellCartIdx >= 0;
322 aboveCellCartIdx -= nx*ny)
324 if (porvData[aboveCellCartIdx] >= eclGrid.getMinpvValue()) {
330 const double aboveCellVolume = eclGrid.getCellVolume(aboveCellCartIdx);
331 if (actnumData[aboveCellCartIdx] == 0 && aboveCellVolume > 1e-6) {
337 cellPoreVolume += porvData[aboveCellCartIdx];
341 pvol_[cellIdx] = cellPoreVolume;
345 template <
class Gr
id>
346 void pinchProcess_(
const Grid& grid,
347 const Opm::EclipseState& eclState,
348 const Vector& htrans,
353 template <
typename Gr
id>
354 bool cartesianAdjacent(
const Grid& grid,
int g1,
int g2) {
355 int diff = std::abs(g1 - g2);
357 const int * dimens = UgGridHelpers::cartDims(grid);
360 if (diff == dimens[0])
362 if (diff == dimens[0] * dimens[1])
374 template <
typename Gr
id>
375 void exportNncStructure(
const Grid& grid) {
377 using namespace UgGridHelpers;
379 size_t num_faces = numFaces(grid);
381 auto fc = faceCells(grid);
382 for (
size_t i = 0; i < num_faces; ++i) {
386 if (c1 == -1 || c2 == -1)
389 c1 = globalCell(grid) ? globalCell(grid)[c1] : c1;
390 c2 = globalCell(grid) ? globalCell(grid)[c2] : c2;
391 if (!cartesianAdjacent(grid, c1, c2)) {
394 noncartesian_.addNNC(c1, c2, trans_[i]);
404 bool use_local_perm_;
412 template <
class Gr
idType>
413 inline void DerivedGeology::minPvFillProps_(
const GridType &grid,
414 const EclipseState& eclState,
415 std::vector<double> &ntg)
418 int numCells = Opm::AutoDiffGrid::numCells(grid);
419 const int* global_cell = Opm::UgGridHelpers::globalCell(grid);
420 const int* cartdims = Opm::UgGridHelpers::cartDims(grid);
421 const auto& eclgrid = eclState.getInputGrid();
422 const auto& porv = eclState.get3DProperties().getDoubleGridProperty(
"PORV").getData();
423 const auto& actnum = eclState.get3DProperties().getIntGridProperty(
"ACTNUM").getData();
424 for (
int cellIdx = 0; cellIdx < numCells; ++cellIdx) {
425 const int nx = cartdims[0];
426 const int ny = cartdims[1];
427 const int cartesianCellIdx = global_cell[cellIdx];
429 const double cellVolume = eclgrid.getCellVolume(cartesianCellIdx);
430 ntg[cartesianCellIdx] *= cellVolume;
431 double totalCellVolume = cellVolume;
435 int cartesianCellIdxAbove = cartesianCellIdx - nx*ny;
436 while ( cartesianCellIdxAbove >= 0 &&
437 actnum[cartesianCellIdxAbove] > 0 &&
438 porv[cartesianCellIdxAbove] < eclgrid.getMinpvValue() ) {
441 const double cellAboveVolume = eclgrid.getCellVolume(cartesianCellIdxAbove);
442 totalCellVolume += cellAboveVolume;
443 ntg[cartesianCellIdx] += ntg[cartesianCellIdxAbove]*cellAboveVolume;
444 cartesianCellIdxAbove -= nx*ny;
446 ntg[cartesianCellIdx] /= totalCellVolume;
453 template <
class Gr
idType>
454 inline void DerivedGeology::pinchProcess_(
const GridType& grid,
455 const Opm::EclipseState& eclState,
456 const Vector& htrans,
461 auto eclgrid = eclState.getInputGrid();
462 auto& eclProps = eclState.get3DProperties();
463 const double minpv = eclgrid.getMinpvValue();
464 const double thickness = eclgrid.getPinchThresholdThickness();
465 auto transMode = eclgrid.getPinchOption();
466 auto multzMode = eclgrid.getMultzOption();
467 PinchProcessor<GridType> pinch(minpv, thickness, transMode, multzMode);
469 std::vector<double> htrans_copy(htrans.size());
470 std::copy_n(htrans.data(), htrans.size(), htrans_copy.begin());
472 std::vector<int> actnum;
473 eclgrid.exportACTNUM(actnum);
475 const auto& transMult = eclState.getTransMult();
476 std::vector<double> multz(numCells, 0.0);
477 const int* global_cell = Opm::UgGridHelpers::globalCell(grid);
479 for (
int i = 0; i < numCells; ++i) {
480 multz[i] = transMult.getMultiplier(global_cell[i], Opm::FaceDir::ZPlus);
484 const auto& porv = eclProps.getDoubleGridProperty(
"PORV").getData();
485 pinch.process(grid, htrans_copy, actnum, multz, porv, nnc_);
491 template <
class Gr
idType>
492 inline void DerivedGeology::multiplyHalfIntersections_(
const GridType &grid,
493 const EclipseState& eclState,
494 const std::vector<double> &ntg,
495 Vector &halfIntersectTransmissibility,
496 std::vector<double> &intersectionTransMult)
498 int numCells = Opm::AutoDiffGrid::numCells(grid);
500 int numIntersections = Opm::AutoDiffGrid::numFaces(grid);
501 intersectionTransMult.resize(numIntersections);
502 std::fill(intersectionTransMult.begin(), intersectionTransMult.end(), 1.0);
504 const TransMult& multipliers = eclState.getTransMult();
505 auto cell2Faces = Opm::UgGridHelpers::cell2Faces(grid);
506 auto faceCells = Opm::AutoDiffGrid::faceCells(grid);
507 const int* global_cell = Opm::UgGridHelpers::globalCell(grid);
510 for (
int cellIdx = 0; cellIdx < numCells; ++cellIdx) {
512 auto cellFacesRange = cell2Faces[cellIdx];
514 for(
auto cellFaceIter = cellFacesRange.begin(), cellFaceEnd = cellFacesRange.end();
515 cellFaceIter != cellFaceEnd; ++cellFaceIter, ++cellFaceIdx)
518 int cartesianCellIdx = global_cell[cellIdx];
521 int faceIdx = *cellFaceIter;
524 int faceTag = Opm::UgGridHelpers::faceTag(grid, cellFaceIter);
527 Opm::FaceDir::DirEnum faceDirection;
529 faceDirection = Opm::FaceDir::XMinus;
530 else if (faceTag == 1)
531 faceDirection = Opm::FaceDir::XPlus;
532 else if (faceTag == 2)
533 faceDirection = Opm::FaceDir::YMinus;
534 else if (faceTag == 3)
535 faceDirection = Opm::FaceDir::YPlus;
536 else if (faceTag == 4)
537 faceDirection = Opm::FaceDir::ZMinus;
538 else if (faceTag == 5)
539 faceDirection = Opm::FaceDir::ZPlus;
541 OPM_THROW(std::logic_error,
"Unhandled face direction: " << faceTag);
544 switch (faceDirection) {
545 case Opm::FaceDir::XMinus:
546 case Opm::FaceDir::XPlus:
547 case Opm::FaceDir::YMinus:
548 case Opm::FaceDir::YPlus:
549 halfIntersectTransmissibility[cellFaceIdx] *= ntg[cartesianCellIdx];
557 intersectionTransMult[faceIdx] *=
558 multipliers.getMultiplier(cartesianCellIdx, faceDirection);
561 const int cellIdxInside = faceCells(faceIdx, 0);
562 const int cellIdxOutside = faceCells(faceIdx, 1);
565 if (cellIdxInside < 0 || cellIdxOutside < 0) {
568 const int cartesianCellIdxInside = global_cell[cellIdxInside];
569 const int cartesianCellIdxOutside = global_cell[cellIdxOutside];
571 if (cartesianCellIdx == cartesianCellIdxInside) {
572 intersectionTransMult[faceIdx] *= multipliers.getRegionMultiplier(cartesianCellIdxInside,cartesianCellIdxOutside,faceDirection);
580 template <
class Gr
idType>
581 inline void DerivedGeology::tpfa_loc_trans_compute_(
const GridType& grid,
582 const EclipseGrid& eclGrid,
591 int numCells = AutoDiffGrid::numCells(grid);
593 auto cell2Faces = Opm::UgGridHelpers::cell2Faces(grid);
594 auto faceCells = Opm::UgGridHelpers::faceCells(grid);
596 for (
int cellIdx = 0; cellIdx < numCells; ++cellIdx) {
598 auto cellFacesRange = cell2Faces[cellIdx];
600 for(
auto cellFaceIter = cellFacesRange.begin(), cellFaceEnd = cellFacesRange.end();
601 cellFaceIter != cellFaceEnd; ++cellFaceIter, ++cellFaceIdx)
604 const int faceIdx = *cellFaceIter;
607 const int faceTag = Opm::UgGridHelpers::faceTag(grid, cellFaceIter);
610 const int d = std::floor(faceTag/2) * 4;
615 double sgn = 2.0 * (faceCells(faceIdx, 0) == cellIdx) - 1;
616 const int dim = Opm::UgGridHelpers::dimensions(grid);
618 int cartesianCellIdx = AutoDiffGrid::globalCell(grid)[cellIdx];
619 auto cellCenter = eclGrid.getCellCenter(cartesianCellIdx);
620 const auto& faceCenter = Opm::UgGridHelpers::faceCenterEcl(grid, cellIdx, faceTag);
621 const auto& faceAreaNormalEcl = Opm::UgGridHelpers::faceAreaNormalEcl(grid, faceIdx);
623 for (
int indx = 0; indx < dim; ++indx) {
624 const double Ci = faceCenter[indx] - cellCenter[indx];
626 cn += sgn * Ci * faceAreaNormalEcl[ indx ];
632 OPM_MESSAGE(
"Warning: negative X-transmissibility value in cell: " << cellIdx <<
" replace by absolute value") ;
635 OPM_MESSAGE(
"Warning: negative Y-transmissibility value in cell: " << cellIdx <<
" replace by absolute value") ;
638 OPM_MESSAGE(
"Warning: negative Z-transmissibility value in cell: " << cellIdx <<
" replace by absolute value") ;
641 OPM_THROW(std::logic_error,
"Inconsistency in the faceTag in cell: " << cellIdx);
646 hTrans[cellFaceIdx] = perm[cellIdx*dim*dim + d] * cn / dist;
657 #endif // OPM_GEOPROPS_HEADER_INCLUDED data::Solution simProps(const Grid &grid) const
Most properties are loaded by the parser, and managed by the EclipseState class in the opm-parser...
Definition: GeoProps.hpp:232
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition: AdditionalObjectDeleter.hpp:22
void update(const Grid &grid, const Props &props, const EclipseState &eclState, const double *grav)
compute the all geological properties at a given report step
Definition: GeoProps.hpp:85
DerivedGeology(const Grid &grid, const Props &props, const EclipseState &eclState, const bool use_local_perm, const double *grav=0)
Construct contained derived geological properties from grid and property information.
Definition: GeoProps.hpp:67
Class containing static geological properties that are derived from grid and petrophysical properties...
Definition: GeoProps.hpp:59
Mapping of the grid type to the type of the cell to faces mapping.
Definition: GridHelpers.hpp:65