All Classes Namespaces Files Functions Variables Typedefs Enumerator Pages
GeoProps.hpp
1 /*
2  Copyright 2013 SINTEF ICT, Applied Mathematics.
3 
4  This file is part of the Open Porous Media project (OPM).
5 
6  OPM is free software: you can redistribute it and/or modify
7  it under the terms of the GNU General Public License as published by
8  the Free Software Foundation, either version 3 of the License, or
9  (at your option) any later version.
10 
11  OPM is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  GNU General Public License for more details.
15 
16  You should have received a copy of the GNU General Public License
17  along with OPM. If not, see <http://www.gnu.org/licenses/>.
18 */
19 
20 #ifndef OPM_GEOPROPS_HEADER_INCLUDED
21 #define OPM_GEOPROPS_HEADER_INCLUDED
22 
23 #include <opm/core/grid.h>
24 #include <opm/autodiff/GridHelpers.hpp>
25 #include <opm/common/ErrorMacros.hpp>
26 //#include <opm/core/pressure/tpfa/trans_tpfa.h>
27 #include <opm/core/pressure/tpfa/TransTpfa.hpp>
28 
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>
38 
39 #include <Eigen/Eigen>
40 
41 #ifdef HAVE_OPM_GRID
42 #include <dune/common/version.hh>
43 #include <dune/grid/CpGrid.hpp>
44 #include <dune/grid/common/mcmgmapper.hh>
45 #endif
46 
47 #include <opm/common/utility/platform_dependent/reenable_warnings.h>
48 
49 #include <cstddef>
50 
51 namespace Opm
52 {
53 
60  {
61  public:
62  typedef Eigen::ArrayXd Vector;
63 
66  template <class Props, class Grid>
67  DerivedGeology(const Grid& grid,
68  const Props& props ,
69  const EclipseState& eclState,
70  const bool use_local_perm,
71  const double* grav = 0
72 
73  )
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)
79  {
80  update(grid, props, eclState, grav);
81  }
82 
84  template <class Props, class Grid>
85  void update(const Grid& grid,
86  const Props& props ,
87  const EclipseState& eclState,
88  const double* grav)
89 
90  {
91  int numCells = AutoDiffGrid::numCells(grid);
92  int numFaces = AutoDiffGrid::numFaces(grid);
93  const int *cartDims = AutoDiffGrid::cartDims(grid);
94  int numCartesianCells =
95  cartDims[0]
96  * cartDims[1]
97  * cartDims[2];
98 
99  // get the pore volume multipliers from the EclipseState
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();
104  }
105 
106  // get the net-to-gross cell thickness from the EclipseState
107  std::vector<double> ntg(numCartesianCells, 1.0);
108  if (eclProps.hasDeckDoubleGridProperty("NTG")) {
109  ntg = eclProps.getDoubleGridProperty("NTG").getData();
110  }
111 
112  // Get grid from parser.
113  const auto& eclgrid = eclState.getInputGrid();
114 
115  // update the pore volume of all active cells in the grid
116  computePoreVolume_(grid, eclState);
117 
118  // Non-neighbour connections.
119  nnc_ = eclState.getInputNNC();
120 
121  // Transmissibility
122  Vector htrans(AutoDiffGrid::numCellFaces(grid));
123  Grid* ug = const_cast<Grid*>(& grid);
124 
125  if (! use_local_perm_) {
126  tpfa_htrans_compute(ug, props.permeability(), htrans.data());
127  }
128  else {
129  tpfa_loc_trans_compute_(grid,eclgrid, props.permeability(),htrans);
130  }
131 
132  // Use volume weighted arithmetic average of the NTG values for
133  // the cells effected by the current OPM cpgrid process algorithm
134  // for MINPV. Note that the change does not effect the pore volume calculations
135  // as the pore volume is currently defaulted to be comparable to ECLIPSE, but
136  // only the transmissibility calculations.
137  bool opmfil = eclgrid.getMinpvMode() == MinpvMode::ModeEnum::OpmFIL;
138  // opmfil is hardcoded to be true. i.e the volume weighting is always used
139  opmfil = true;
140  if (opmfil) {
141  minPvFillProps_(grid, eclState, ntg);
142  }
143 
144  std::vector<double> mult;
145  multiplyHalfIntersections_(grid, eclState, ntg, htrans, mult);
146 
147  if (!opmfil && eclgrid.isPinchActive()) {
148  // opmfil is hardcoded to be true. i.e the pinch processor is never used
149  pinchProcess_(grid, eclState, htrans, numCells);
150  }
151 
152  // combine the half-face transmissibilites into the final face
153  // transmissibilites.
154  tpfa_trans_compute(ug, htrans.data(), trans_.data());
155 
156  // multiply the face transmissibilities with their appropriate
157  // transmissibility multipliers
158  for (int faceIdx = 0; faceIdx < numFaces; faceIdx++) {
159  trans_[faceIdx] *= mult[faceIdx];
160  }
161 
162  // Create the set of noncartesian connections.
163  noncartesian_ = nnc_;
164  exportNncStructure(grid);
165 
166  // Compute z coordinates
167  for (int c = 0; c<numCells; ++c){
168  z_[c] = Opm::UgGridHelpers::cellCenterDepth(grid, c);
169  }
170 
171 
172  // Gravity potential
173  std::fill(gravity_, gravity_ + 3, 0.0);
174  if (grav != 0) {
175  const typename Vector::Index nd = AutoDiffGrid::dimensions(grid);
176  typedef typename AutoDiffGrid::ADCell2FacesTraits<Grid>::Type Cell2Faces;
177  Cell2Faces c2f=AutoDiffGrid::cell2Faces(grid);
178 
179  std::size_t i = 0;
180  for (typename Vector::Index c = 0; c < numCells; ++c) {
181  const double* const cc = AutoDiffGrid::cellCentroid(grid, c);
182 
183  typename Cell2Faces::row_type faces=c2f[c];
184  typedef typename Cell2Faces::row_type::iterator Iter;
185 
186  for (Iter f=faces.begin(), end=faces.end(); f!=end; ++f, ++i) {
187  auto fc = AutoDiffGrid::faceCentroid(grid, *f);
188 
189  for (typename Vector::Index d = 0; d < nd; ++d) {
190  gpot_[i] += grav[d] * (fc[d] - cc[d]);
191  }
192  }
193  }
194  std::copy(grav, grav + nd, gravity_);
195  }
196  }
197 
198 
199 
200 
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_;}
210 
211 
230 
231  template <class Grid>
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( );
237 
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};
241 
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));
247 
248  if (c1 == -1 || c2 == -1)
249  continue;
250 
251  c1 = globalCell(grid) ? globalCell(grid)[c1] : c1;
252  c2 = globalCell(grid) ? globalCell(grid)[c2] : c2;
253 
254  if ((c2 - c1) == 1) {
255  tranx.data[c1] = trans[i];
256  }
257 
258  if ((c2 - c1) == dims[0]) {
259  trany.data[c1] = trans[i];
260  }
261 
262  if ((c2 - c1) == dims[0]*dims[1]) {
263  tranz.data[c1] = trans[i];
264  }
265  }
266 
267  return { {"TRANX" , tranx},
268  {"TRANY" , trany} ,
269  {"TRANZ" , tranz } };
270  }
271 
272 
273  private:
274  template <class Grid>
275  void multiplyHalfIntersections_(const Grid &grid,
276  const EclipseState& eclState,
277  const std::vector<double> &ntg,
278  Vector &halfIntersectTransmissibility,
279  std::vector<double> &intersectionTransMult);
280 
281  template <class Grid>
282  void tpfa_loc_trans_compute_(const Grid &grid,
283  const EclipseGrid& eclGrid,
284  const double* perm,
285  Vector &hTrans);
286 
287  template <class Grid>
288  void minPvFillProps_(const Grid &grid,
289  const EclipseState& eclState,
290  std::vector<double> &ntg);
291 
292  template <class GridType>
293  void computePoreVolume_(const GridType &grid,
294  const EclipseState& eclState)
295  {
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();
301 
302  // the "raw" pore volume.
303  const std::vector<double>& porvData =
304  eclState.get3DProperties().getDoubleGridProperty("PORV").getData();
305  pvol_.resize(numCells);
306 
307  // the "activation number" grid property
308  const std::vector<int>& actnumData =
309  eclState.get3DProperties().getIntGridProperty("ACTNUM").getData();
310 
311 
312  for (int cellIdx = 0; cellIdx < numCells; ++cellIdx) {
313  const int cellCartIdx = globalCell[cellIdx];
314 
315  double cellPoreVolume = porvData[cellCartIdx];
316 
317  if (eclGrid.getMinpvMode() == MinpvMode::ModeEnum::OpmFIL) {
318  // Sum the pore volumes of the cells above which have been deactivated
319  // because their volume less is less than the MINPV threshold
320  for (int aboveCellCartIdx = cellCartIdx - nx*ny;
321  aboveCellCartIdx >= 0;
322  aboveCellCartIdx -= nx*ny)
323  {
324  if (porvData[aboveCellCartIdx] >= eclGrid.getMinpvValue()) {
325  // stop if we encounter a cell which has a pore volume which is
326  // at least as large as the minimum one
327  break;
328  }
329 
330  const double aboveCellVolume = eclGrid.getCellVolume(aboveCellCartIdx);
331  if (actnumData[aboveCellCartIdx] == 0 && aboveCellVolume > 1e-6) {
332  // stop at explicitly disabled cells, but only if their volume is
333  // greater than 10^-6 m^3
334  break;
335  }
336 
337  cellPoreVolume += porvData[aboveCellCartIdx];
338  }
339  }
340 
341  pvol_[cellIdx] = cellPoreVolume;
342  }
343  }
344 
345  template <class Grid>
346  void pinchProcess_(const Grid& grid,
347  const Opm::EclipseState& eclState,
348  const Vector& htrans,
349  int numCells);
350 
351 
353  template <typename Grid>
354  bool cartesianAdjacent(const Grid& grid, int g1, int g2) {
355  int diff = std::abs(g1 - g2);
356 
357  const int * dimens = UgGridHelpers::cartDims(grid);
358  if (diff == 1)
359  return true;
360  if (diff == dimens[0])
361  return true;
362  if (diff == dimens[0] * dimens[1])
363  return true;
364 
365  return false;
366  }
367 
368 
374  template <typename Grid>
375  void exportNncStructure(const Grid& grid) {
376  // we use numFaces, numCells, cell2Faces, globalCell from UgGridHelpers
377  using namespace UgGridHelpers;
378 
379  size_t num_faces = numFaces(grid);
380 
381  auto fc = faceCells(grid);
382  for (size_t i = 0; i < num_faces; ++i) {
383  auto c1 = fc(i, 0);
384  auto c2 = fc(i, 1);
385 
386  if (c1 == -1 || c2 == -1)
387  continue; // face on grid boundary
388  // translate from active cell idx (ac1,ac2) to global cell idx
389  c1 = globalCell(grid) ? globalCell(grid)[c1] : c1;
390  c2 = globalCell(grid) ? globalCell(grid)[c2] : c2;
391  if (!cartesianAdjacent(grid, c1, c2)) {
392  // suppose c1,c2 is specified in ECLIPSE input
393  // we here overwrite its trans by grid's
394  noncartesian_.addNNC(c1, c2, trans_[i]);
395  }
396  }
397  }
398 
399  Vector pvol_ ;
400  Vector trans_;
401  Vector gpot_ ;
402  Vector z_;
403  double gravity_[3]; // Size 3 even if grid is 2-dim.
404  bool use_local_perm_;
405 
406  // Non-neighboring connections
407  NNC nnc_;
408  // Non-cartesian connections
409  NNC noncartesian_;
410  };
411 
412  template <class GridType>
413  inline void DerivedGeology::minPvFillProps_(const GridType &grid,
414  const EclipseState& eclState,
415  std::vector<double> &ntg)
416  {
417 
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];
428 
429  const double cellVolume = eclgrid.getCellVolume(cartesianCellIdx);
430  ntg[cartesianCellIdx] *= cellVolume;
431  double totalCellVolume = cellVolume;
432 
433  // Average properties as long as there exist cells above
434  // that has pore volume less than the MINPV threshold
435  int cartesianCellIdxAbove = cartesianCellIdx - nx*ny;
436  while ( cartesianCellIdxAbove >= 0 &&
437  actnum[cartesianCellIdxAbove] > 0 &&
438  porv[cartesianCellIdxAbove] < eclgrid.getMinpvValue() ) {
439 
440  // Volume weighted arithmetic average of NTG
441  const double cellAboveVolume = eclgrid.getCellVolume(cartesianCellIdxAbove);
442  totalCellVolume += cellAboveVolume;
443  ntg[cartesianCellIdx] += ntg[cartesianCellIdxAbove]*cellAboveVolume;
444  cartesianCellIdxAbove -= nx*ny;
445  }
446  ntg[cartesianCellIdx] /= totalCellVolume;
447  }
448  }
449 
450 
451 
452 
453  template <class GridType>
454  inline void DerivedGeology::pinchProcess_(const GridType& grid,
455  const Opm::EclipseState& eclState,
456  const Vector& htrans,
457  int numCells)
458  {
459  // NOTE that this function is currently never invoked due to
460  // opmfil being hardcoded to be true.
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);
468 
469  std::vector<double> htrans_copy(htrans.size());
470  std::copy_n(htrans.data(), htrans.size(), htrans_copy.begin());
471 
472  std::vector<int> actnum;
473  eclgrid.exportACTNUM(actnum);
474 
475  const auto& transMult = eclState.getTransMult();
476  std::vector<double> multz(numCells, 0.0);
477  const int* global_cell = Opm::UgGridHelpers::globalCell(grid);
478 
479  for (int i = 0; i < numCells; ++i) {
480  multz[i] = transMult.getMultiplier(global_cell[i], Opm::FaceDir::ZPlus);
481  }
482 
483  // Note the pore volume from eclState is used and not the pvol_ calculated above
484  const auto& porv = eclProps.getDoubleGridProperty("PORV").getData();
485  pinch.process(grid, htrans_copy, actnum, multz, porv, nnc_);
486  }
487 
488 
489 
490 
491  template <class GridType>
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)
497  {
498  int numCells = Opm::AutoDiffGrid::numCells(grid);
499 
500  int numIntersections = Opm::AutoDiffGrid::numFaces(grid);
501  intersectionTransMult.resize(numIntersections);
502  std::fill(intersectionTransMult.begin(), intersectionTransMult.end(), 1.0);
503 
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);
508  int cellFaceIdx = 0;
509 
510  for (int cellIdx = 0; cellIdx < numCells; ++cellIdx) {
511  // loop over all logically-Cartesian faces of the current cell
512  auto cellFacesRange = cell2Faces[cellIdx];
513 
514  for(auto cellFaceIter = cellFacesRange.begin(), cellFaceEnd = cellFacesRange.end();
515  cellFaceIter != cellFaceEnd; ++cellFaceIter, ++cellFaceIdx)
516  {
517  // the index of the current cell in arrays for the logically-Cartesian grid
518  int cartesianCellIdx = global_cell[cellIdx];
519 
520  // The index of the face in the compressed grid
521  int faceIdx = *cellFaceIter;
522 
523  // the logically-Cartesian direction of the face
524  int faceTag = Opm::UgGridHelpers::faceTag(grid, cellFaceIter);
525 
526  // Translate the C face tag into the enum used by opm-parser's TransMult class
527  Opm::FaceDir::DirEnum faceDirection;
528  if (faceTag == 0) // left
529  faceDirection = Opm::FaceDir::XMinus;
530  else if (faceTag == 1) // right
531  faceDirection = Opm::FaceDir::XPlus;
532  else if (faceTag == 2) // back
533  faceDirection = Opm::FaceDir::YMinus;
534  else if (faceTag == 3) // front
535  faceDirection = Opm::FaceDir::YPlus;
536  else if (faceTag == 4) // bottom
537  faceDirection = Opm::FaceDir::ZMinus;
538  else if (faceTag == 5) // top
539  faceDirection = Opm::FaceDir::ZPlus;
540  else
541  OPM_THROW(std::logic_error, "Unhandled face direction: " << faceTag);
542 
543  // Account for NTG in horizontal one-sided transmissibilities
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];
550  break;
551  default:
552  // do nothing for the top and bottom faces
553  break;
554  }
555 
556  // Multiplier contribution on this face for MULT[XYZ] logical cartesian multipliers
557  intersectionTransMult[faceIdx] *=
558  multipliers.getMultiplier(cartesianCellIdx, faceDirection);
559 
560  // Multiplier contribution on this fase for region multipliers
561  const int cellIdxInside = faceCells(faceIdx, 0);
562  const int cellIdxOutside = faceCells(faceIdx, 1);
563 
564  // Do not apply region multipliers in the case of boundary connections
565  if (cellIdxInside < 0 || cellIdxOutside < 0) {
566  continue;
567  }
568  const int cartesianCellIdxInside = global_cell[cellIdxInside];
569  const int cartesianCellIdxOutside = global_cell[cellIdxOutside];
570  // Only apply the region multipliers from the inside
571  if (cartesianCellIdx == cartesianCellIdxInside) {
572  intersectionTransMult[faceIdx] *= multipliers.getRegionMultiplier(cartesianCellIdxInside,cartesianCellIdxOutside,faceDirection);
573  }
574 
575 
576  }
577  }
578  }
579 
580  template <class GridType>
581  inline void DerivedGeology::tpfa_loc_trans_compute_(const GridType& grid,
582  const EclipseGrid& eclGrid,
583  const double* perm,
584  Vector& hTrans){
585 
586  // Using Local coordinate system for the transmissibility calculations
587  // hTrans(cellFaceIdx) = K(cellNo,j) * sum( C(:,i) .* N(:,j), 2) / sum(C.*C, 2)
588  // where K is a diagonal permeability tensor, C is the distance from cell centroid
589  // to face centroid and N is the normal vector pointing outwards with norm equal to the face area.
590  // Off-diagonal permeability values are ignored without warning
591  int numCells = AutoDiffGrid::numCells(grid);
592  int cellFaceIdx = 0;
593  auto cell2Faces = Opm::UgGridHelpers::cell2Faces(grid);
594  auto faceCells = Opm::UgGridHelpers::faceCells(grid);
595 
596  for (int cellIdx = 0; cellIdx < numCells; ++cellIdx) {
597  // loop over all logically-Cartesian faces of the current cell
598  auto cellFacesRange = cell2Faces[cellIdx];
599 
600  for(auto cellFaceIter = cellFacesRange.begin(), cellFaceEnd = cellFacesRange.end();
601  cellFaceIter != cellFaceEnd; ++cellFaceIter, ++cellFaceIdx)
602  {
603  // The index of the face in the compressed grid
604  const int faceIdx = *cellFaceIter;
605 
606  // the logically-Cartesian direction of the face
607  const int faceTag = Opm::UgGridHelpers::faceTag(grid, cellFaceIter);
608 
609  // d = 0: XPERM d = 4: YPERM d = 8: ZPERM ignores off-diagonal permeability values.
610  const int d = std::floor(faceTag/2) * 4;
611 
612  // compute the half transmissibility
613  double dist = 0.0;
614  double cn = 0.0;
615  double sgn = 2.0 * (faceCells(faceIdx, 0) == cellIdx) - 1;
616  const int dim = Opm::UgGridHelpers::dimensions(grid);
617 
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);
622 
623  for (int indx = 0; indx < dim; ++indx) {
624  const double Ci = faceCenter[indx] - cellCenter[indx];
625  dist += Ci*Ci;
626  cn += sgn * Ci * faceAreaNormalEcl[ indx ];
627  }
628 
629  if (cn < 0){
630  switch (d) {
631  case 0:
632  OPM_MESSAGE("Warning: negative X-transmissibility value in cell: " << cellIdx << " replace by absolute value") ;
633  break;
634  case 4:
635  OPM_MESSAGE("Warning: negative Y-transmissibility value in cell: " << cellIdx << " replace by absolute value") ;
636  break;
637  case 8:
638  OPM_MESSAGE("Warning: negative Z-transmissibility value in cell: " << cellIdx << " replace by absolute value") ;
639  break;
640  default:
641  OPM_THROW(std::logic_error, "Inconsistency in the faceTag in cell: " << cellIdx);
642 
643  }
644  cn = -cn;
645  }
646  hTrans[cellFaceIdx] = perm[cellIdx*dim*dim + d] * cn / dist;
647 
648  }
649  }
650 
651  }
652 
653 }
654 
655 
656 
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
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