PinchProcessor.hpp
Go to the documentation of this file.
1 /*
2  Copyright 2015 Statoil ASA.
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_PINCHPROCESSOR_HEADER_INCLUDED
21 #define OPM_PINCHPROCESSOR_HEADER_INCLUDED
22 
23 #include <opm/common/ErrorMacros.hpp>
24 #include <opm/core/grid/GridHelpers.hpp>
25 #include <opm/parser/eclipse/EclipseState/Grid/NNC.hpp>
26 #include <opm/parser/eclipse/EclipseState/Grid/FaceDir.hpp>
27 #include <opm/parser/eclipse/EclipseState/Grid/PinchMode.hpp>
28 #include <opm/core/utility/Units.hpp>
29 #include <array>
30 #include <iostream>
31 #include <algorithm>
32 #include <unordered_map>
33 #include <limits>
34 
35 namespace Opm
36 {
37 
38  template <class Grid>
39  class PinchProcessor
40  {
41  public:
47  PinchProcessor(const double minpvValue,
48  const double thickness,
49  const PinchMode::ModeEnum transMode,
50  const PinchMode::ModeEnum multzMode);
62  void process(const Grid& grid,
63  const std::vector<double>& htrans,
64  const std::vector<int>& actnum,
65  const std::vector<double>& multz,
66  const std::vector<double>& pv,
67  NNC& nnc);
68 
69  private:
70  double minpvValue_;
71  double thickness_;
72  PinchMode::ModeEnum transMode_;
73  PinchMode::ModeEnum multzMode_;
74 
76  std::vector<int> getMinpvCells_(const std::vector<int>& actnum,
77  const std::vector<double>& pv);
78 
80  int interface_(const Grid& grid,
81  const int cellIdx1,
82  const int cellIdx2);
83 
85  int interface_(const Grid& grid,
86  const int cellIdx,
87  const Opm::FaceDir::DirEnum& faceDir);
88 
90  std::vector<std::vector<int> >
91  getPinchoutsColumn_(const Grid& grid,
92  const std::vector<int>& actnum,
93  const std::vector<double>& pv);
94 
96  int getGlobalIndex_(const int i, const int j, const int k, const int* dims);
97 
99  std::array<int, 3> getCartIndex_(const int idx,
100  const int* dims);
101 
103  std::vector<double> transCompute_(const Grid& grid,
104  const std::vector<double>& htrans,
105  const std::vector<int>& pinCells,
106  const std::vector<int>& pinFaces);
107 
109  std::vector<int> getHfIdxMap_(const Grid& grid);
110 
112  int getActiveCellIdx_(const Grid& grid,
113  const int globalIdx);
114 
116  void transTopbot_(const Grid& grid,
117  const std::vector<double>& htrans,
118  const std::vector<int>& actnum,
119  const std::vector<double>& multz,
120  const std::vector<double>& pv,
121  NNC& nnc);
122 
124  std::unordered_multimap<int, double> multzOptions_(const Grid& grid,
125  const std::vector<int>& pinCells,
126  const std::vector<int>& pinFaces,
127  const std::vector<double>& multz,
128  const std::vector<std::vector<int> >& seg);
129 
131  void applyMultz_(std::vector<double>& trans,
132  const std::unordered_multimap<int, double>& multzmap);
133 
134  };
135 
136 
137 
138  template <class Grid>
139  inline PinchProcessor<Grid>::PinchProcessor(const double minpv,
140  const double thickness,
141  const PinchMode::ModeEnum transMode,
142  const PinchMode::ModeEnum multzMode)
143  {
144  minpvValue_ = minpv;
145  thickness_ = thickness;
146  transMode_ = transMode;
147  multzMode_ = multzMode;
148  }
149 
150 
151 
152  template <class Grid>
153  inline int PinchProcessor<Grid>::getGlobalIndex_(const int i, const int j, const int k, const int* dims)
154  {
155  return i + dims[0] * (j + dims[1] * k);
156  }
157 
158 
159 
160  template <class Grid>
161  inline std::array<int, 3> PinchProcessor<Grid>::getCartIndex_(const int idx,
162  const int* dims)
163  {
164  std::array<int, 3> ijk;
165  ijk[0] = (idx % dims[0]);
166  ijk[1] = ((idx / dims[0]) % dims[1]);
167  ijk[2] = ((idx / dims[0]) / dims[1]);
168 
169  return ijk;
170  }
171 
172 
173 
174  template<class Grid>
175  inline int PinchProcessor<Grid>::interface_(const Grid& grid,
176  const int cellIdx1,
177  const int cellIdx2)
178  {
179  const auto cell_faces = Opm::UgGridHelpers::cell2Faces(grid);
180  int commonFace = -1;
181  const int actCellIdx1 = getActiveCellIdx_(grid, cellIdx1);
182  const int actCellIdx2 = getActiveCellIdx_(grid, cellIdx2);
183  const auto cellFacesRange1 = cell_faces[actCellIdx1];
184  const auto cellFacesRange2 = cell_faces[actCellIdx2];
185  for (const auto& f1 : cellFacesRange1) {
186  for (const auto& f2 : cellFacesRange2) {
187  if (f1 == f2) {
188  commonFace = f1;
189  break;
190  }
191  }
192  }
193 
194  if (commonFace == -1) {
195  const auto dims = Opm::UgGridHelpers::cartDims(grid);
196  const auto ijk1 = getCartIndex_(cellIdx1, dims);
197  const auto ijk2 = getCartIndex_(cellIdx2, dims);
198 
199  OPM_THROW(std::logic_error, "Couldn't find the common face for cell "
200  << cellIdx1<< "("<<ijk1[0]<<","<<ijk1[1]<<","<<ijk1[2]<<")"
201  << " and " << cellIdx2<<"("<<ijk2[0]<<","<<ijk2[1]<<","<<ijk2[2]<<")");
202  }
203 
204  return commonFace;
205  }
206 
207 
208 
209  template<class Grid>
210  inline int PinchProcessor<Grid>::interface_(const Grid& grid,
211  const int cellIdx,
212  const Opm::FaceDir::DirEnum& faceDir)
213  {
214  const auto actCellIdx = getActiveCellIdx_(grid, cellIdx);
215  const auto cell_faces = Opm::UgGridHelpers::cell2Faces(grid);
216  const auto cellFacesRange = cell_faces[actCellIdx];
217  int faceIdx = -1;
218  for (auto cellFaceIter = cellFacesRange.begin(); cellFaceIter != cellFacesRange.end(); ++cellFaceIter) {
219  int tag = Opm::UgGridHelpers::faceTag(grid, cellFaceIter);
220  if ( (faceDir == Opm::FaceDir::ZMinus && tag == 4) || (faceDir == Opm::FaceDir::ZPlus && tag == 5) ) {
221  faceIdx = *cellFaceIter;
222  }
223  }
224 
225  if (faceIdx == -1) {
226  OPM_THROW(std::logic_error, "Couldn't find the face for cell ." << cellIdx);
227  }
228 
229  return faceIdx;
230  }
231 
232 
233 
234  template<class Grid>
235  inline std::vector<int> PinchProcessor<Grid>::getMinpvCells_(const std::vector<int>& actnum,
236  const std::vector<double>& pv)
237  {
238  std::vector<int> minpvCells(pv.size(), 0);
239  for (int idx = 0; idx < static_cast<int>(pv.size()); ++idx) {
240  if (actnum[idx]) {
241  if (pv[idx] < minpvValue_) {
242  minpvCells[idx] = 1;
243  }
244  }
245  }
246  return minpvCells;
247  }
248 
249 
250 
251  template<class Grid>
252  inline std::vector<int> PinchProcessor<Grid>::getHfIdxMap_(const Grid& grid)
253  {
254  std::vector<int> hf_ix(2*Opm::UgGridHelpers::numFaces(grid), -1);
255  const auto& f2c = Opm::UgGridHelpers::faceCells(grid);
256  const auto& cf = Opm::UgGridHelpers::cell2Faces(grid);
257 
258  for (int c = 0, i = 0; c < Opm::UgGridHelpers::numCells(grid); ++c) {
259  for (const auto& f: cf[c]) {
260  const auto off = 0 + (f2c(f, 0) != c);
261  hf_ix[2*f + off] = i++;
262  }
263  }
264  return hf_ix;
265  }
266 
267 
268 
269  template<class Grid>
270  inline int PinchProcessor<Grid>::getActiveCellIdx_(const Grid& grid,
271  const int globalIdx)
272  {
273  const int nc = Opm::UgGridHelpers::numCells(grid);
274  const int* global_cell = Opm::UgGridHelpers::globalCell(grid);
275  int idx = -1;
276  for (int i = 0; i < nc; ++i) {
277  if (global_cell[i] == globalIdx) {
278  idx = i;
279  break;
280  }
281  }
282  return idx;
283  }
284 
285 
286 
287  template<class Grid>
288  inline std::vector<double> PinchProcessor<Grid>::transCompute_(const Grid& grid,
289  const std::vector<double>& htrans,
290  const std::vector<int>& pinCells,
291  const std::vector<int>& pinFaces)
292  {
293  const int nc = Opm::UgGridHelpers::numCells(grid);
294  const int nf = Opm::UgGridHelpers::numFaces(grid);
295  std::vector<double> trans(nf, 0);
296  int cellFaceIdx = 0;
297  auto cell_faces = Opm::UgGridHelpers::cell2Faces(grid);
298  const auto& hfmap = getHfIdxMap_(grid);
299  const auto& f2c = Opm::UgGridHelpers::faceCells(grid);
300  for (int cellIdx = 0; cellIdx < nc; ++cellIdx) {
301  auto cellFacesRange = cell_faces[cellIdx];
302  for (auto cellFaceIter = cellFacesRange.begin(); cellFaceIter != cellFacesRange.end(); ++cellFaceIter, ++cellFaceIdx) {
303  const int faceIdx = *cellFaceIter;
304  const auto pos = std::find(pinFaces.begin(), pinFaces.end(), faceIdx);
305  if (pos == pinFaces.end()) {
306  trans[faceIdx] += 1. / htrans[cellFaceIdx];
307  } else {
308  const int idx1 = std::distance(std::begin(pinFaces), pos);
309  int idx2;
310  if (idx1 % 2 == 0) {
311  idx2 = idx1 + 1;
312  } else {
313  idx2 = idx1 - 1;
314  }
315  const int f1 = hfmap[2*pinFaces[idx1] + (f2c(pinFaces[idx1], 0) != getActiveCellIdx_(grid, pinCells[idx1]))];
316  const int f2 = hfmap[2*pinFaces[idx2] + (f2c(pinFaces[idx2], 0) != getActiveCellIdx_(grid, pinCells[idx2]))];
317  trans[faceIdx] = (1. / htrans[f1] + 1. / htrans[f2]);
318  trans[pinFaces[idx2]] = trans[faceIdx];
319  }
320  }
321  }
322 
323  for (auto f = 0; f < nf; ++f) {
324  trans[f] = 1. / trans[f];
325  }
326 
327  return trans;
328  }
329 
330 
331 
332  template<class Grid>
333  inline std::vector<std::vector<int>> PinchProcessor<Grid>::getPinchoutsColumn_(const Grid& grid,
334  const std::vector<int>& actnum,
335  const std::vector<double>& pv)
336  {
337  const int* dims = Opm::UgGridHelpers::cartDims(grid);
338  std::vector<int> minpvCells = getMinpvCells_(actnum, pv);
339  std::vector<std::vector<int>> segment;
340  for (int z = 0; z < dims[2]; ++z) {
341  for (int y = 0; y < dims[1]; ++y) {
342  for (int x = 0; x < dims[0]; ++x) {
343  const int c = getGlobalIndex_(x, y, z, dims);
344  std::vector<int> seg;
345  if (minpvCells[c]) {
346  seg.push_back(c);
347  minpvCells[c] = 0;
348  for (int zz = z+1; zz < dims[2]; ++zz) {
349  const int cc = getGlobalIndex_(x, y, zz, dims);
350  if (minpvCells[cc]) {
351  seg.push_back(cc);
352  minpvCells[cc] = 0;
353  } else {
354  break;
355  }
356  }
357  segment.push_back(seg);
358  }
359  }
360  }
361  }
362 
363  return segment;
364  }
365 
366 
367 
368  template<class Grid>
369  inline void PinchProcessor<Grid>::transTopbot_(const Grid& grid,
370  const std::vector<double>& htrans,
371  const std::vector<int>& actnum,
372  const std::vector<double>& multz,
373  const std::vector<double>& pv,
374  NNC& nnc)
375  {
376  const int* dims = Opm::UgGridHelpers::cartDims(grid);
377  std::vector<int> pinFaces;
378  std::vector<int> pinCells;
379  std::vector<std::vector<int> > newSeg;
380  auto minpvSeg = getPinchoutsColumn_(grid, actnum, pv);
381  for (auto& seg : minpvSeg) {
382  std::array<int, 3> ijk1 = getCartIndex_(seg.front(), dims);
383  std::array<int, 3> ijk2 = getCartIndex_(seg.back(), dims);
384  auto tmp = seg;
385  if ((ijk1[2]-1) >= 0 && (ijk2[2]+1) < dims[2]) {
386  int topCell = getGlobalIndex_(ijk1[0], ijk1[1], ijk1[2]-1, dims);
387  int botCell = getGlobalIndex_(ijk2[0], ijk2[1], ijk2[2]+1, dims);
391  if (!actnum[topCell]) {
392  seg.insert(seg.begin(), topCell);
393  for (int topk = ijk1[2]-2; topk > 0; --topk) {
394  topCell = getGlobalIndex_(ijk1[0], ijk1[1], topk, dims);
395  if (actnum[topCell]) {
396  break;
397  } else {
398  auto it = seg.begin();
399  seg.insert(it, topCell);
400  }
401  }
402  pinFaces.push_back(interface_(grid, topCell, Opm::FaceDir::ZPlus));
403  } else {
404  pinFaces.push_back(interface_(grid, topCell, seg.front()));
405  }
406  tmp.insert(tmp.begin(), topCell);
407  newSeg.push_back(tmp);
408  pinCells.push_back(topCell);
409  if (!actnum[botCell]) {
410  seg.push_back(botCell);
411  for (int botk = ijk2[2]+2; botk < dims[2]; ++botk) {
412  botCell = getGlobalIndex_(ijk2[0], ijk2[1], botk, dims);
413  if (actnum[botCell]) {
414  break;
415  } else {
416  seg.push_back(botCell);
417  }
418  }
419  pinFaces.push_back(interface_(grid, botCell, Opm::FaceDir::ZMinus));
420  } else {
421  pinFaces.push_back(interface_(grid, seg.back(), botCell));
422  }
423  pinCells.push_back(botCell);
424  }
425  }
426 
427  auto faceTrans = transCompute_(grid, htrans, pinCells, pinFaces);
428  auto multzmap = multzOptions_(grid, pinCells, pinFaces, multz, newSeg);
429  applyMultz_(faceTrans, multzmap);
430  for (int i = 0; i < static_cast<int>(pinCells.size())/2; ++i) {
431  nnc.addNNC(static_cast<int>(pinCells[2*i]), static_cast<int>(pinCells[2*i+1]), faceTrans[pinFaces[2*i]]);
432  }
433  }
434 
435 
436 
437 
438  template<class Grid>
439  inline std::unordered_multimap<int, double> PinchProcessor<Grid>::multzOptions_(const Grid& grid,
440  const std::vector<int>& pinCells,
441  const std::vector<int>& pinFaces,
442  const std::vector<double>& multz,
443  const std::vector<std::vector<int> >& segs)
444  {
445  std::unordered_multimap<int, double> multzmap;
446  if (multzMode_ == PinchMode::ModeEnum::TOP) {
447  for (int i = 0; i < static_cast<int>(pinFaces.size())/2; ++i) {
448  multzmap.insert(std::make_pair(pinFaces[2*i], multz[getActiveCellIdx_(grid, pinCells[2*i])]));
449  multzmap.insert(std::make_pair(pinFaces[2*i+1],multz[getActiveCellIdx_(grid, pinCells[2*i])]));
450  }
451  } else if (multzMode_ == PinchMode::ModeEnum::ALL) {
452  for (auto& seg : segs) {
453  //find the min multz in seg cells.
454  auto multzValue = std::numeric_limits<double>::max();
455  for (auto& cellIdx : seg) {
456  auto activeIdx = getActiveCellIdx_(grid, cellIdx);
457  if (activeIdx != -1) {
458  multzValue = std::min(multzValue, multz[activeIdx]);
459  }
460  }
461  //find the right face.
462  auto index = std::distance(std::begin(pinCells), std::find(pinCells.begin(), pinCells.end(), seg.front()));
463  multzmap.insert(std::make_pair(pinFaces[index], multzValue));
464  multzmap.insert(std::make_pair(pinFaces[index+1], multzValue));
465  }
466  }
467 
468  return multzmap;
469  }
470 
471 
472 
473  template<class Grid>
474  inline void PinchProcessor<Grid>::applyMultz_(std::vector<double>& trans,
475  const std::unordered_multimap<int, double>& multzmap)
476  {
477  for (auto& x : multzmap) {
478  trans[x.first] *= x.second;
479  }
480  }
481 
482 
483 
484  template<class Grid>
485  inline void PinchProcessor<Grid>::process(const Grid& grid,
486  const std::vector<double>& htrans,
487  const std::vector<int>& actnum,
488  const std::vector<double>& multz,
489  const std::vector<double>& pv,
490  NNC& nnc)
491  {
492  transTopbot_(grid, htrans, actnum, multz, pv, nnc);
493  }
494 
495 } // namespace Opm
496 #endif // OPM_PINCHPROCESSOR_HEADER_INCLUDED
Opm::UgGridHelpers::faceCells
FaceCellTraits< UnstructuredGrid >::Type faceCells(const UnstructuredGrid &grid)
Get the face to cell mapping of a grid.
Opm
Definition: AnisotropicEikonal.hpp:43
Opm::UgGridHelpers::faceTag
int faceTag(const UnstructuredGrid &grid, boost::iterator_range< const int * >::const_iterator cell_face)
Get Eclipse Cartesian tag of a face.
Opm::UgGridHelpers::numCells
int numCells(const UnstructuredGrid &grid)
Get the number of cells of a grid.
Opm::UgGridHelpers::numFaces
int numFaces(const UnstructuredGrid &grid)
Get the number of faces of a grid.
Opm::PinchProcessor::process
void process(const Grid &grid, const std::vector< double > &htrans, const std::vector< int > &actnum, const std::vector< double > &multz, const std::vector< double > &pv, NNC &nnc)
Definition: PinchProcessor.hpp:485
Opm::UgGridHelpers::globalCell
const int * globalCell(const UnstructuredGrid &grid)
Get the local to global index mapping.
Opm::PinchProcessor::PinchProcessor
PinchProcessor(const double minpvValue, const double thickness, const PinchMode::ModeEnum transMode, const PinchMode::ModeEnum multzMode)
Create a Pinch processor.
Definition: PinchProcessor.hpp:139
Opm::PinchProcessor
Definition: PinchProcessor.hpp:39
Opm::UgGridHelpers::cartDims
const int * cartDims(const UnstructuredGrid &grid)
Get the cartesion dimension of the underlying structured grid.
GridHelpers.hpp
grid
const UnstructuredGrid & grid
Definition: ColumnExtract.hpp:31
Opm::UgGridHelpers::cell2Faces
Cell2FacesTraits< UnstructuredGrid >::Type cell2Faces(const UnstructuredGrid &grid)
Get the cell to faces mapping of a grid.
TOP
Definition: preprocess.h:70
Units.hpp