20 #ifndef OPM_GEOPROPS_HEADER_INCLUDED
21 #define OPM_GEOPROPS_HEADER_INCLUDED
23 #include <opm/core/grid.h>
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/common/utility/platform_dependent/disable_warnings.h>
33 #include <Eigen/Eigen>
35 #ifdef HAVE_DUNE_CORNERPOINT
36 #include <dune/common/version.hh>
37 #include <dune/grid/CpGrid.hpp>
38 #include <dune/grid/common/mcmgmapper.hh>
41 #include <opm/common/utility/platform_dependent/reenable_warnings.h>
60 template <
class Props,
class Gr
id>
63 Opm::EclipseStateConstPtr eclState,
64 const bool use_local_perm,
65 const double* grav = 0
68 : pvol_ (
Opm::AutoDiffGrid::numCells(grid))
69 , trans_(
Opm::AutoDiffGrid::numFaces(grid))
70 , gpot_ (Vector::Zero(
Opm::AutoDiffGrid::cell2Faces(grid).noEntries(), 1))
71 , z_(
Opm::AutoDiffGrid::numCells(grid))
73 int numCells = AutoDiffGrid::numCells(grid);
74 int numFaces = AutoDiffGrid::numFaces(grid);
75 const int *cartDims = AutoDiffGrid::cartDims(grid);
76 int numCartesianCells =
82 std::vector<double> multpv(numCartesianCells, 1.0);
83 if (eclState->hasDoubleGridProperty(
"MULTPV")) {
84 multpv = eclState->getDoubleGridProperty(
"MULTPV")->getData();
88 std::vector<double> ntg(numCartesianCells, 1.0);
89 if (eclState->hasDoubleGridProperty(
"NTG")) {
90 ntg = eclState->getDoubleGridProperty(
"NTG")->getData();
96 EclipseGridConstPtr eclgrid = eclState->getEclipseGrid();
100 for (
int cellIdx = 0; cellIdx < numCells; ++cellIdx) {
101 int cartesianCellIdx = AutoDiffGrid::globalCell(grid)[cellIdx];
103 props.porosity()[cellIdx]
104 * multpv[cartesianCellIdx]
105 * ntg[cartesianCellIdx];
106 if (eclgrid->getMinpvMode() == MinpvMode::ModeEnum::OpmFIL) {
107 pvol_[cellIdx] *= AutoDiffGrid::cellVolume(grid, cellIdx);
109 pvol_[cellIdx] *= eclgrid->getCellVolume(cartesianCellIdx);
117 minPvFillProps_(grid, eclState,ntg);
121 Vector htrans(AutoDiffGrid::numCellFaces(grid));
122 Grid* ug =
const_cast<Grid*
>(& grid);
124 if (! use_local_perm) {
125 tpfa_htrans_compute(ug, props.permeability(), htrans.data());
128 tpfa_loc_trans_compute_(grid,props.permeability(),htrans);
131 std::vector<double> mult;
132 multiplyHalfIntersections_(grid, eclState, ntg, htrans, mult);
136 tpfa_trans_compute(ug, htrans.data(), trans_.data());
140 for (
int faceIdx = 0; faceIdx < numFaces; faceIdx++) {
141 trans_[faceIdx] *= mult[faceIdx];
145 for (
int c = 0; c<numCells; ++c){
146 z_[c] = Opm::UgGridHelpers::cellCentroidCoordinate(grid, c, 2);
151 std::fill(gravity_, gravity_ + 3, 0.0);
153 const typename Vector::Index nd = AutoDiffGrid::dimensions(grid);
155 Cell2Faces c2f=AutoDiffGrid::cell2Faces(grid);
158 for (
typename Vector::Index c = 0; c < numCells; ++c) {
159 const double*
const cc = AutoDiffGrid::cellCentroid(grid, c);
161 typename Cell2Faces::row_type faces=c2f[c];
162 typedef typename Cell2Faces::row_type::iterator Iter;
164 for (Iter f=faces.begin(), end=faces.end(); f!=end; ++f, ++i) {
165 auto fc = AutoDiffGrid::faceCentroid(grid, *f);
167 for (
typename Vector::Index d = 0; d < nd; ++d) {
168 gpot_[i] += grav[d] * (fc[d] - cc[d]);
172 std::copy(grav, grav + nd, gravity_);
179 const Vector&
z()
const {
return z_ ;}
180 const double*
gravity()
const {
return gravity_;}
185 template <
class Gr
id>
186 void multiplyHalfIntersections_(
const Grid &grid,
187 Opm::EclipseStateConstPtr eclState,
188 const std::vector<double> &ntg,
189 Vector &halfIntersectTransmissibility,
190 std::vector<double> &intersectionTransMult);
192 template <
class Gr
id>
193 void tpfa_loc_trans_compute_(
const Grid &grid,
197 template <
class Gr
id>
198 void minPvFillProps_(
const Grid &grid,
199 Opm::EclipseStateConstPtr eclState,
200 std::vector<double> &ntg);
213 template <
class Gr
idType>
214 inline void DerivedGeology::minPvFillProps_(
const GridType &grid,
215 Opm::EclipseStateConstPtr eclState,
216 std::vector<double> &ntg)
219 int numCells = Opm::AutoDiffGrid::numCells(grid);
220 const int* global_cell = Opm::UgGridHelpers::globalCell(grid);
221 const int* cartdims = Opm::UgGridHelpers::cartDims(grid);
222 EclipseGridConstPtr eclgrid = eclState->getEclipseGrid();
223 std::vector<double> porv = eclState->getDoubleGridProperty(
"PORV")->getData();
224 for (
int cellIdx = 0; cellIdx < numCells; ++cellIdx) {
225 const int nx = cartdims[0];
226 const int ny = cartdims[1];
227 const int cartesianCellIdx = global_cell[cellIdx];
229 const double cellVolume = eclgrid->getCellVolume(cartesianCellIdx);
230 ntg[cartesianCellIdx] *= cellVolume;
231 double totalCellVolume = cellVolume;
235 int cartesianCellIdxAbove = cartesianCellIdx - nx*ny;
236 while ( cartesianCellIdxAbove >= 0 &&
237 porv[cartesianCellIdxAbove] > 0 &&
238 porv[cartesianCellIdxAbove] < eclgrid->getMinpvValue() ) {
241 const double cellAboveVolume = eclgrid->getCellVolume(cartesianCellIdxAbove);
242 totalCellVolume += cellAboveVolume;
243 ntg[cartesianCellIdx] += ntg[cartesianCellIdxAbove]*cellAboveVolume;
244 cartesianCellIdxAbove -= nx*ny;
246 ntg[cartesianCellIdx] /= totalCellVolume;
250 template <
class Gr
idType>
251 inline void DerivedGeology::multiplyHalfIntersections_(
const GridType &grid,
252 Opm::EclipseStateConstPtr eclState,
253 const std::vector<double> &ntg,
254 Vector &halfIntersectTransmissibility,
255 std::vector<double> &intersectionTransMult)
257 int numCells = Opm::AutoDiffGrid::numCells(grid);
259 int numIntersections = Opm::AutoDiffGrid::numFaces(grid);
260 intersectionTransMult.resize(numIntersections);
261 std::fill(intersectionTransMult.begin(), intersectionTransMult.end(), 1.0);
263 std::shared_ptr<const Opm::TransMult> multipliers = eclState->getTransMult();
264 auto cell2Faces = Opm::UgGridHelpers::cell2Faces(grid);
265 auto faceCells = Opm::AutoDiffGrid::faceCells(grid);
266 const int* global_cell = Opm::UgGridHelpers::globalCell(grid);
269 for (
int cellIdx = 0; cellIdx < numCells; ++cellIdx) {
271 auto cellFacesRange = cell2Faces[cellIdx];
273 for(
auto cellFaceIter = cellFacesRange.begin(), cellFaceEnd = cellFacesRange.end();
274 cellFaceIter != cellFaceEnd; ++cellFaceIter, ++cellFaceIdx)
277 int cartesianCellIdx = global_cell[cellIdx];
280 int faceIdx = *cellFaceIter;
283 int faceTag = Opm::UgGridHelpers::faceTag(grid, cellFaceIter);
286 Opm::FaceDir::DirEnum faceDirection;
288 faceDirection = Opm::FaceDir::XMinus;
289 else if (faceTag == 1)
290 faceDirection = Opm::FaceDir::XPlus;
291 else if (faceTag == 2)
292 faceDirection = Opm::FaceDir::YMinus;
293 else if (faceTag == 3)
294 faceDirection = Opm::FaceDir::YPlus;
295 else if (faceTag == 4)
296 faceDirection = Opm::FaceDir::ZMinus;
297 else if (faceTag == 5)
298 faceDirection = Opm::FaceDir::ZPlus;
300 OPM_THROW(std::logic_error,
"Unhandled face direction: " << faceTag);
303 switch (faceDirection) {
304 case Opm::FaceDir::XMinus:
305 case Opm::FaceDir::XPlus:
306 case Opm::FaceDir::YMinus:
307 case Opm::FaceDir::YPlus:
308 halfIntersectTransmissibility[cellFaceIdx] *= ntg[cartesianCellIdx];
316 intersectionTransMult[faceIdx] *=
317 multipliers->getMultiplier(cartesianCellIdx, faceDirection);
320 const int cellIdxInside = faceCells(faceIdx, 0);
321 const int cellIdxOutside = faceCells(faceIdx, 1);
324 if (cellIdxInside < 0 || cellIdxOutside < 0) {
327 const int cartesianCellIdxInside = global_cell[cellIdxInside];
328 const int cartesianCellIdxOutside = global_cell[cellIdxOutside];
330 if (cartesianCellIdx == cartesianCellIdxInside) {
331 intersectionTransMult[faceIdx] *= multipliers->getRegionMultiplier(cartesianCellIdxInside,cartesianCellIdxOutside,faceDirection);
339 template <
class Gr
idType>
340 inline void DerivedGeology::tpfa_loc_trans_compute_(
const GridType& grid,
349 int numCells = AutoDiffGrid::numCells(grid);
351 auto cell2Faces = Opm::UgGridHelpers::cell2Faces(grid);
352 auto faceCells = Opm::UgGridHelpers::faceCells(grid);
354 for (
int cellIdx = 0; cellIdx < numCells; ++cellIdx) {
356 auto cellFacesRange = cell2Faces[cellIdx];
358 for(
auto cellFaceIter = cellFacesRange.begin(), cellFaceEnd = cellFacesRange.end();
359 cellFaceIter != cellFaceEnd; ++cellFaceIter, ++cellFaceIdx)
362 const int faceIdx = *cellFaceIter;
365 const int faceTag = Opm::UgGridHelpers::faceTag(grid, cellFaceIter);
368 const int d = std::floor(faceTag/2) * 4;
373 double sgn = 2.0 * (faceCells(faceIdx, 0) == cellIdx) - 1;
374 const int dim = Opm::UgGridHelpers::dimensions(grid);
376 const double* faceNormal = Opm::UgGridHelpers::faceNormal(grid, faceIdx);
377 #if HAVE_DUNE_CORNERPOINT
379 Dune::FieldVector< double, 3 > scaledFaceNormal( 0 );
380 for (
int indx = 0; indx < dim; ++indx) {
381 scaledFaceNormal[ indx ] = faceNormal[ indx ];
384 scaledFaceNormal /= scaledFaceNormal.two_norm();
386 scaledFaceNormal *= Opm::UgGridHelpers::faceArea(grid, faceIdx);
388 const double* scaledFaceNormal = faceNormal;
391 for (
int indx = 0; indx < dim; ++indx) {
392 const double Ci = Opm::UgGridHelpers::faceCentroid(grid, faceIdx)[indx] -
393 Opm::UgGridHelpers::cellCentroidCoordinate(grid, cellIdx, indx);
395 cn += sgn * Ci * scaledFaceNormal[ indx ];
401 OPM_MESSAGE(
"Warning: negative X-transmissibility value in cell: " << cellIdx <<
" replace by absolute value") ;
404 OPM_MESSAGE(
"Warning: negative Y-transmissibility value in cell: " << cellIdx <<
" replace by absolute value") ;
407 OPM_MESSAGE(
"Warning: negative Z-transmissibility value in cell: " << cellIdx <<
" replace by absolute value") ;
410 OPM_THROW(std::logic_error,
"Inconsistency in the faceTag in cell: " << cellIdx);
415 hTrans[cellFaceIdx] = perm[cellIdx*dim*dim + d] * cn / dist;
425 #endif // OPM_GEOPROPS_HEADER_INCLUDED
Definition: GeoProps.hpp:53
const Vector & transmissibility() const
Definition: GeoProps.hpp:177
Definition: AdditionalObjectDeleter.hpp:22
DerivedGeology(const Grid &grid, const Props &props, Opm::EclipseStateConstPtr eclState, const bool use_local_perm, const double *grav=0)
Definition: GeoProps.hpp:61
Mapping of the grid type to the type of the cell to faces mapping.
Definition: GridHelpers.hpp:64
Vector & poreVolume()
Definition: GeoProps.hpp:181
const Vector & poreVolume() const
Definition: GeoProps.hpp:176
const Vector & gravityPotential() const
Definition: GeoProps.hpp:178
Eigen::ArrayXd Vector
Definition: GeoProps.hpp:56
const double * gravity() const
Definition: GeoProps.hpp:180
const Vector & z() const
Definition: GeoProps.hpp:179
Vector & transmissibility()
Definition: GeoProps.hpp:182