00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020 #ifndef OPM_GEOPROPS_HEADER_INCLUDED
00021 #define OPM_GEOPROPS_HEADER_INCLUDED
00022
00023 #include <opm/core/grid.h>
00024 #include <opm/autodiff/GridHelpers.hpp>
00025 #include <opm/common/ErrorMacros.hpp>
00026
00027 #include <opm/core/pressure/tpfa/TransTpfa.hpp>
00028
00029 #include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
00030 #include <opm/parser/eclipse/EclipseState/Grid/EclipseGrid.hpp>
00031 #include <opm/parser/eclipse/EclipseState/Grid/GridProperty.hpp>
00032 #include <opm/parser/eclipse/EclipseState/Grid/NNC.hpp>
00033 #include <opm/parser/eclipse/EclipseState/Grid/TransMult.hpp>
00034 #include <opm/core/grid/PinchProcessor.hpp>
00035 #include <opm/common/utility/platform_dependent/disable_warnings.h>
00036 #include <opm/output/data/Cells.hpp>
00037 #include <opm/output/data/Solution.hpp>
00038
00039 #include <Eigen/Eigen>
00040
00041 #ifdef HAVE_OPM_GRID
00042 #include <dune/common/version.hh>
00043 #include <dune/grid/CpGrid.hpp>
00044 #include <dune/grid/common/mcmgmapper.hh>
00045 #endif
00046
00047 #include <opm/common/utility/platform_dependent/reenable_warnings.h>
00048
00049 #include <cstddef>
00050
00051 namespace Opm
00052 {
00053
00059 class DerivedGeology
00060 {
00061 public:
00062 typedef Eigen::ArrayXd Vector;
00063
00066 template <class Props, class Grid>
00067 DerivedGeology(const Grid& grid,
00068 const Props& props ,
00069 const EclipseState& eclState,
00070 const bool use_local_perm,
00071 const double* grav = 0
00072
00073 )
00074 : pvol_ (Opm::AutoDiffGrid::numCells(grid))
00075 , trans_(Opm::AutoDiffGrid::numFaces(grid))
00076 , gpot_ (Vector::Zero(Opm::AutoDiffGrid::cell2Faces(grid).noEntries(), 1))
00077 , z_(Opm::AutoDiffGrid::numCells(grid))
00078 , use_local_perm_(use_local_perm)
00079 {
00080 update(grid, props, eclState, grav);
00081 }
00082
00084 template <class Props, class Grid>
00085 void update(const Grid& grid,
00086 const Props& props ,
00087 const EclipseState& eclState,
00088 const double* grav)
00089
00090 {
00091 int numCells = AutoDiffGrid::numCells(grid);
00092 int numFaces = AutoDiffGrid::numFaces(grid);
00093 const int *cartDims = AutoDiffGrid::cartDims(grid);
00094 int numCartesianCells =
00095 cartDims[0]
00096 * cartDims[1]
00097 * cartDims[2];
00098
00099
00100 std::vector<double> multpv(numCartesianCells, 1.0);
00101 const auto& eclProps = eclState.get3DProperties();
00102 if (eclProps.hasDeckDoubleGridProperty("MULTPV")) {
00103 multpv = eclProps.getDoubleGridProperty("MULTPV").getData();
00104 }
00105
00106
00107 std::vector<double> ntg(numCartesianCells, 1.0);
00108 if (eclProps.hasDeckDoubleGridProperty("NTG")) {
00109 ntg = eclProps.getDoubleGridProperty("NTG").getData();
00110 }
00111
00112
00113 const auto& eclgrid = eclState.getInputGrid();
00114
00115
00116 computePoreVolume_(grid, eclState);
00117
00118
00119 nnc_ = eclState.getInputNNC();
00120
00121
00122 Vector htrans(AutoDiffGrid::numCellFaces(grid));
00123 Grid* ug = const_cast<Grid*>(& grid);
00124
00125 if (! use_local_perm_) {
00126 tpfa_htrans_compute(ug, props.permeability(), htrans.data());
00127 }
00128 else {
00129 tpfa_loc_trans_compute_(grid,eclgrid, props.permeability(),htrans);
00130 }
00131
00132
00133
00134
00135
00136
00137 bool opmfil = eclgrid.getMinpvMode() == MinpvMode::ModeEnum::OpmFIL;
00138
00139 opmfil = true;
00140 if (opmfil) {
00141 minPvFillProps_(grid, eclState, ntg);
00142 }
00143
00144 std::vector<double> mult;
00145 multiplyHalfIntersections_(grid, eclState, ntg, htrans, mult);
00146
00147 if (!opmfil && eclgrid.isPinchActive()) {
00148
00149 pinchProcess_(grid, eclState, htrans, numCells);
00150 }
00151
00152
00153
00154 tpfa_trans_compute(ug, htrans.data(), trans_.data());
00155
00156
00157
00158 for (int faceIdx = 0; faceIdx < numFaces; faceIdx++) {
00159 trans_[faceIdx] *= mult[faceIdx];
00160 }
00161
00162
00163 noncartesian_ = nnc_;
00164 exportNncStructure(grid);
00165
00166
00167 for (int c = 0; c<numCells; ++c){
00168 z_[c] = Opm::UgGridHelpers::cellCenterDepth(grid, c);
00169 }
00170
00171
00172
00173 std::fill(gravity_, gravity_ + 3, 0.0);
00174 if (grav != 0) {
00175 const typename Vector::Index nd = AutoDiffGrid::dimensions(grid);
00176 typedef typename AutoDiffGrid::ADCell2FacesTraits<Grid>::Type Cell2Faces;
00177 Cell2Faces c2f=AutoDiffGrid::cell2Faces(grid);
00178
00179 std::size_t i = 0;
00180 for (typename Vector::Index c = 0; c < numCells; ++c) {
00181 const double* const cc = AutoDiffGrid::cellCentroid(grid, c);
00182
00183 typename Cell2Faces::row_type faces=c2f[c];
00184 typedef typename Cell2Faces::row_type::iterator Iter;
00185
00186 for (Iter f=faces.begin(), end=faces.end(); f!=end; ++f, ++i) {
00187 auto fc = AutoDiffGrid::faceCentroid(grid, *f);
00188
00189 for (typename Vector::Index d = 0; d < nd; ++d) {
00190 gpot_[i] += grav[d] * (fc[d] - cc[d]);
00191 }
00192 }
00193 }
00194 std::copy(grav, grav + nd, gravity_);
00195 }
00196 }
00197
00198
00199
00200
00201 const Vector& poreVolume() const { return pvol_ ;}
00202 const Vector& transmissibility() const { return trans_ ;}
00203 const Vector& gravityPotential() const { return gpot_ ;}
00204 const Vector& z() const { return z_ ;}
00205 const double* gravity() const { return gravity_;}
00206 Vector& poreVolume() { return pvol_ ;}
00207 Vector& transmissibility() { return trans_ ;}
00208 const NNC& nnc() const { return nnc_;}
00209 const NNC& nonCartesianConnections() const { return noncartesian_;}
00210
00211
00230
00231 template <class Grid>
00232 data::Solution simProps( const Grid& grid ) const {
00233 using namespace UgGridHelpers;
00234 const int* dims = cartDims( grid );
00235 const int globalSize = dims[0] * dims[1] * dims[2];
00236 const auto& trans = this->transmissibility( );
00237
00238 data::CellData tranx = {UnitSystem::measure::transmissibility, std::vector<double>( globalSize ), data::TargetType::INIT};
00239 data::CellData trany = {UnitSystem::measure::transmissibility, std::vector<double>( globalSize ), data::TargetType::INIT};
00240 data::CellData tranz = {UnitSystem::measure::transmissibility, std::vector<double>( globalSize ), data::TargetType::INIT};
00241
00242 size_t num_faces = numFaces(grid);
00243 auto fc = faceCells(grid);
00244 for (size_t i = 0; i < num_faces; ++i) {
00245 auto c1 = std::min( fc(i,0) , fc(i,1));
00246 auto c2 = std::max( fc(i,0) , fc(i,1));
00247
00248 if (c1 == -1 || c2 == -1)
00249 continue;
00250
00251 c1 = globalCell(grid) ? globalCell(grid)[c1] : c1;
00252 c2 = globalCell(grid) ? globalCell(grid)[c2] : c2;
00253
00254 if ((c2 - c1) == 1) {
00255 tranx.data[c1] = trans[i];
00256 }
00257
00258 if ((c2 - c1) == dims[0]) {
00259 trany.data[c1] = trans[i];
00260 }
00261
00262 if ((c2 - c1) == dims[0]*dims[1]) {
00263 tranz.data[c1] = trans[i];
00264 }
00265 }
00266
00267 return { {"TRANX" , tranx},
00268 {"TRANY" , trany} ,
00269 {"TRANZ" , tranz } };
00270 }
00271
00272
00273 private:
00274 template <class Grid>
00275 void multiplyHalfIntersections_(const Grid &grid,
00276 const EclipseState& eclState,
00277 const std::vector<double> &ntg,
00278 Vector &halfIntersectTransmissibility,
00279 std::vector<double> &intersectionTransMult);
00280
00281 template <class Grid>
00282 void tpfa_loc_trans_compute_(const Grid &grid,
00283 const EclipseGrid& eclGrid,
00284 const double* perm,
00285 Vector &hTrans);
00286
00287 template <class Grid>
00288 void minPvFillProps_(const Grid &grid,
00289 const EclipseState& eclState,
00290 std::vector<double> &ntg);
00291
00292 template <class GridType>
00293 void computePoreVolume_(const GridType &grid,
00294 const EclipseState& eclState)
00295 {
00296 int numCells = Opm::AutoDiffGrid::numCells(grid);
00297 const int* globalCell = Opm::UgGridHelpers::globalCell(grid);
00298 const auto& eclGrid = eclState.getInputGrid();
00299 const int nx = eclGrid.getNX();
00300 const int ny = eclGrid.getNY();
00301
00302
00303 const std::vector<double>& porvData =
00304 eclState.get3DProperties().getDoubleGridProperty("PORV").getData();
00305 pvol_.resize(numCells);
00306
00307
00308 const std::vector<int>& actnumData =
00309 eclState.get3DProperties().getIntGridProperty("ACTNUM").getData();
00310
00311
00312 for (int cellIdx = 0; cellIdx < numCells; ++cellIdx) {
00313 const int cellCartIdx = globalCell[cellIdx];
00314
00315 double cellPoreVolume = porvData[cellCartIdx];
00316
00317 if (eclGrid.getMinpvMode() == MinpvMode::ModeEnum::OpmFIL) {
00318
00319
00320 for (int aboveCellCartIdx = cellCartIdx - nx*ny;
00321 aboveCellCartIdx >= 0;
00322 aboveCellCartIdx -= nx*ny)
00323 {
00324 if (porvData[aboveCellCartIdx] >= eclGrid.getMinpvValue()) {
00325
00326
00327 break;
00328 }
00329
00330 const double aboveCellVolume = eclGrid.getCellVolume(aboveCellCartIdx);
00331 if (actnumData[aboveCellCartIdx] == 0 && aboveCellVolume > 1e-6) {
00332
00333
00334 break;
00335 }
00336
00337 cellPoreVolume += porvData[aboveCellCartIdx];
00338 }
00339 }
00340
00341 pvol_[cellIdx] = cellPoreVolume;
00342 }
00343 }
00344
00345 template <class Grid>
00346 void pinchProcess_(const Grid& grid,
00347 const Opm::EclipseState& eclState,
00348 const Vector& htrans,
00349 int numCells);
00350
00351
00353 template <typename Grid>
00354 bool cartesianAdjacent(const Grid& grid, int g1, int g2) {
00355 int diff = std::abs(g1 - g2);
00356
00357 const int * dimens = UgGridHelpers::cartDims(grid);
00358 if (diff == 1)
00359 return true;
00360 if (diff == dimens[0])
00361 return true;
00362 if (diff == dimens[0] * dimens[1])
00363 return true;
00364
00365 return false;
00366 }
00367
00368
00374 template <typename Grid>
00375 void exportNncStructure(const Grid& grid) {
00376
00377 using namespace UgGridHelpers;
00378
00379 size_t num_faces = numFaces(grid);
00380
00381 auto fc = faceCells(grid);
00382 for (size_t i = 0; i < num_faces; ++i) {
00383 auto c1 = fc(i, 0);
00384 auto c2 = fc(i, 1);
00385
00386 if (c1 == -1 || c2 == -1)
00387 continue;
00388
00389 c1 = globalCell(grid) ? globalCell(grid)[c1] : c1;
00390 c2 = globalCell(grid) ? globalCell(grid)[c2] : c2;
00391 if (!cartesianAdjacent(grid, c1, c2)) {
00392
00393
00394 noncartesian_.addNNC(c1, c2, trans_[i]);
00395 }
00396 }
00397 }
00398
00399 Vector pvol_ ;
00400 Vector trans_;
00401 Vector gpot_ ;
00402 Vector z_;
00403 double gravity_[3];
00404 bool use_local_perm_;
00405
00406
00407 NNC nnc_;
00408
00409 NNC noncartesian_;
00410 };
00411
00412 template <class GridType>
00413 inline void DerivedGeology::minPvFillProps_(const GridType &grid,
00414 const EclipseState& eclState,
00415 std::vector<double> &ntg)
00416 {
00417
00418 int numCells = Opm::AutoDiffGrid::numCells(grid);
00419 const int* global_cell = Opm::UgGridHelpers::globalCell(grid);
00420 const int* cartdims = Opm::UgGridHelpers::cartDims(grid);
00421 const auto& eclgrid = eclState.getInputGrid();
00422 const auto& porv = eclState.get3DProperties().getDoubleGridProperty("PORV").getData();
00423 const auto& actnum = eclState.get3DProperties().getIntGridProperty("ACTNUM").getData();
00424 for (int cellIdx = 0; cellIdx < numCells; ++cellIdx) {
00425 const int nx = cartdims[0];
00426 const int ny = cartdims[1];
00427 const int cartesianCellIdx = global_cell[cellIdx];
00428
00429 const double cellVolume = eclgrid.getCellVolume(cartesianCellIdx);
00430 ntg[cartesianCellIdx] *= cellVolume;
00431 double totalCellVolume = cellVolume;
00432
00433
00434
00435 int cartesianCellIdxAbove = cartesianCellIdx - nx*ny;
00436 while ( cartesianCellIdxAbove >= 0 &&
00437 actnum[cartesianCellIdxAbove] > 0 &&
00438 porv[cartesianCellIdxAbove] < eclgrid.getMinpvValue() ) {
00439
00440
00441 const double cellAboveVolume = eclgrid.getCellVolume(cartesianCellIdxAbove);
00442 totalCellVolume += cellAboveVolume;
00443 ntg[cartesianCellIdx] += ntg[cartesianCellIdxAbove]*cellAboveVolume;
00444 cartesianCellIdxAbove -= nx*ny;
00445 }
00446 ntg[cartesianCellIdx] /= totalCellVolume;
00447 }
00448 }
00449
00450
00451
00452
00453 template <class GridType>
00454 inline void DerivedGeology::pinchProcess_(const GridType& grid,
00455 const Opm::EclipseState& eclState,
00456 const Vector& htrans,
00457 int numCells)
00458 {
00459
00460
00461 auto eclgrid = eclState.getInputGrid();
00462 auto& eclProps = eclState.get3DProperties();
00463 const double minpv = eclgrid.getMinpvValue();
00464 const double thickness = eclgrid.getPinchThresholdThickness();
00465 auto transMode = eclgrid.getPinchOption();
00466 auto multzMode = eclgrid.getMultzOption();
00467 PinchProcessor<GridType> pinch(minpv, thickness, transMode, multzMode);
00468
00469 std::vector<double> htrans_copy(htrans.size());
00470 std::copy_n(htrans.data(), htrans.size(), htrans_copy.begin());
00471
00472 std::vector<int> actnum;
00473 eclgrid.exportACTNUM(actnum);
00474
00475 const auto& transMult = eclState.getTransMult();
00476 std::vector<double> multz(numCells, 0.0);
00477 const int* global_cell = Opm::UgGridHelpers::globalCell(grid);
00478
00479 for (int i = 0; i < numCells; ++i) {
00480 multz[i] = transMult.getMultiplier(global_cell[i], Opm::FaceDir::ZPlus);
00481 }
00482
00483
00484 const auto& porv = eclProps.getDoubleGridProperty("PORV").getData();
00485 pinch.process(grid, htrans_copy, actnum, multz, porv, nnc_);
00486 }
00487
00488
00489
00490
00491 template <class GridType>
00492 inline void DerivedGeology::multiplyHalfIntersections_(const GridType &grid,
00493 const EclipseState& eclState,
00494 const std::vector<double> &ntg,
00495 Vector &halfIntersectTransmissibility,
00496 std::vector<double> &intersectionTransMult)
00497 {
00498 int numCells = Opm::AutoDiffGrid::numCells(grid);
00499
00500 int numIntersections = Opm::AutoDiffGrid::numFaces(grid);
00501 intersectionTransMult.resize(numIntersections);
00502 std::fill(intersectionTransMult.begin(), intersectionTransMult.end(), 1.0);
00503
00504 const TransMult& multipliers = eclState.getTransMult();
00505 auto cell2Faces = Opm::UgGridHelpers::cell2Faces(grid);
00506 auto faceCells = Opm::AutoDiffGrid::faceCells(grid);
00507 const int* global_cell = Opm::UgGridHelpers::globalCell(grid);
00508 int cellFaceIdx = 0;
00509
00510 for (int cellIdx = 0; cellIdx < numCells; ++cellIdx) {
00511
00512 auto cellFacesRange = cell2Faces[cellIdx];
00513
00514 for(auto cellFaceIter = cellFacesRange.begin(), cellFaceEnd = cellFacesRange.end();
00515 cellFaceIter != cellFaceEnd; ++cellFaceIter, ++cellFaceIdx)
00516 {
00517
00518 int cartesianCellIdx = global_cell[cellIdx];
00519
00520
00521 int faceIdx = *cellFaceIter;
00522
00523
00524 int faceTag = Opm::UgGridHelpers::faceTag(grid, cellFaceIter);
00525
00526
00527 Opm::FaceDir::DirEnum faceDirection;
00528 if (faceTag == 0)
00529 faceDirection = Opm::FaceDir::XMinus;
00530 else if (faceTag == 1)
00531 faceDirection = Opm::FaceDir::XPlus;
00532 else if (faceTag == 2)
00533 faceDirection = Opm::FaceDir::YMinus;
00534 else if (faceTag == 3)
00535 faceDirection = Opm::FaceDir::YPlus;
00536 else if (faceTag == 4)
00537 faceDirection = Opm::FaceDir::ZMinus;
00538 else if (faceTag == 5)
00539 faceDirection = Opm::FaceDir::ZPlus;
00540 else
00541 OPM_THROW(std::logic_error, "Unhandled face direction: " << faceTag);
00542
00543
00544 switch (faceDirection) {
00545 case Opm::FaceDir::XMinus:
00546 case Opm::FaceDir::XPlus:
00547 case Opm::FaceDir::YMinus:
00548 case Opm::FaceDir::YPlus:
00549 halfIntersectTransmissibility[cellFaceIdx] *= ntg[cartesianCellIdx];
00550 break;
00551 default:
00552
00553 break;
00554 }
00555
00556
00557 intersectionTransMult[faceIdx] *=
00558 multipliers.getMultiplier(cartesianCellIdx, faceDirection);
00559
00560
00561 const int cellIdxInside = faceCells(faceIdx, 0);
00562 const int cellIdxOutside = faceCells(faceIdx, 1);
00563
00564
00565 if (cellIdxInside < 0 || cellIdxOutside < 0) {
00566 continue;
00567 }
00568 const int cartesianCellIdxInside = global_cell[cellIdxInside];
00569 const int cartesianCellIdxOutside = global_cell[cellIdxOutside];
00570
00571 if (cartesianCellIdx == cartesianCellIdxInside) {
00572 intersectionTransMult[faceIdx] *= multipliers.getRegionMultiplier(cartesianCellIdxInside,cartesianCellIdxOutside,faceDirection);
00573 }
00574
00575
00576 }
00577 }
00578 }
00579
00580 template <class GridType>
00581 inline void DerivedGeology::tpfa_loc_trans_compute_(const GridType& grid,
00582 const EclipseGrid& eclGrid,
00583 const double* perm,
00584 Vector& hTrans){
00585
00586
00587
00588
00589
00590
00591 int numCells = AutoDiffGrid::numCells(grid);
00592 int cellFaceIdx = 0;
00593 auto cell2Faces = Opm::UgGridHelpers::cell2Faces(grid);
00594 auto faceCells = Opm::UgGridHelpers::faceCells(grid);
00595
00596 for (int cellIdx = 0; cellIdx < numCells; ++cellIdx) {
00597
00598 auto cellFacesRange = cell2Faces[cellIdx];
00599
00600 for(auto cellFaceIter = cellFacesRange.begin(), cellFaceEnd = cellFacesRange.end();
00601 cellFaceIter != cellFaceEnd; ++cellFaceIter, ++cellFaceIdx)
00602 {
00603
00604 const int faceIdx = *cellFaceIter;
00605
00606
00607 const int faceTag = Opm::UgGridHelpers::faceTag(grid, cellFaceIter);
00608
00609
00610 const int d = std::floor(faceTag/2) * 4;
00611
00612
00613 double dist = 0.0;
00614 double cn = 0.0;
00615 double sgn = 2.0 * (faceCells(faceIdx, 0) == cellIdx) - 1;
00616 const int dim = Opm::UgGridHelpers::dimensions(grid);
00617
00618 int cartesianCellIdx = AutoDiffGrid::globalCell(grid)[cellIdx];
00619 auto cellCenter = eclGrid.getCellCenter(cartesianCellIdx);
00620 const auto& faceCenter = Opm::UgGridHelpers::faceCenterEcl(grid, cellIdx, faceTag);
00621 const auto& faceAreaNormalEcl = Opm::UgGridHelpers::faceAreaNormalEcl(grid, faceIdx);
00622
00623 for (int indx = 0; indx < dim; ++indx) {
00624 const double Ci = faceCenter[indx] - cellCenter[indx];
00625 dist += Ci*Ci;
00626 cn += sgn * Ci * faceAreaNormalEcl[ indx ];
00627 }
00628
00629 if (cn < 0){
00630 switch (d) {
00631 case 0:
00632 OPM_MESSAGE("Warning: negative X-transmissibility value in cell: " << cellIdx << " replace by absolute value") ;
00633 break;
00634 case 4:
00635 OPM_MESSAGE("Warning: negative Y-transmissibility value in cell: " << cellIdx << " replace by absolute value") ;
00636 break;
00637 case 8:
00638 OPM_MESSAGE("Warning: negative Z-transmissibility value in cell: " << cellIdx << " replace by absolute value") ;
00639 break;
00640 default:
00641 OPM_THROW(std::logic_error, "Inconsistency in the faceTag in cell: " << cellIdx);
00642
00643 }
00644 cn = -cn;
00645 }
00646 hTrans[cellFaceIdx] = perm[cellIdx*dim*dim + d] * cn / dist;
00647
00648 }
00649 }
00650
00651 }
00652
00653 }
00654
00655
00656
00657 #endif // OPM_GEOPROPS_HEADER_INCLUDED