20 #ifndef OPM_PINCHPROCESSOR_HEADER_INCLUDED 21 #define OPM_PINCHPROCESSOR_HEADER_INCLUDED 23 #include <opm/grid/utility/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/parser/eclipse/Units/Units.hpp> 32 #include <unordered_map> 48 const double thickness,
49 const PinchMode::ModeEnum transMode,
50 const PinchMode::ModeEnum multzMode);
63 const std::vector<double>& htrans,
64 const std::vector<int>& actnum,
65 const std::vector<double>& multz,
66 const std::vector<double>& pv,
72 PinchMode::ModeEnum transMode_;
73 PinchMode::ModeEnum multzMode_;
76 std::vector<int> getMinpvCells_(
const std::vector<int>& actnum,
77 const std::vector<double>& pv);
80 int interface_(
const Grid& grid,
85 int interface_(
const Grid& grid,
87 const Opm::FaceDir::DirEnum& faceDir);
90 std::vector<std::vector<int> >
91 getPinchoutsColumn_(
const Grid& grid,
92 const std::vector<int>& actnum,
93 const std::vector<double>& pv);
96 int getGlobalIndex_(
const int i,
const int j,
const int k,
const int* dims);
99 std::array<int, 3> getCartIndex_(
const int idx,
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);
109 std::vector<int> getHfIdxMap_(
const Grid& grid);
112 int getActiveCellIdx_(
const Grid& grid,
113 const int globalIdx);
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,
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);
131 void applyMultz_(std::vector<double>& trans,
132 const std::unordered_multimap<int, double>& multzmap);
138 template <
class Gr
id>
140 const double thickness,
141 const PinchMode::ModeEnum transMode,
142 const PinchMode::ModeEnum multzMode)
145 thickness_ = thickness;
146 transMode_ = transMode;
147 multzMode_ = multzMode;
152 template <
class Gr
id>
155 return i + dims[0] * (j + dims[1] * k);
160 template <
class Gr
id>
161 inline std::array<int, 3> PinchProcessor<Grid>::getCartIndex_(
const int idx,
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]);
175 inline int PinchProcessor<Grid>::interface_(
const Grid& grid,
179 const auto cell_faces = Opm::UgGridHelpers::cell2Faces(grid);
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) {
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);
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]<<
")");
210 inline int PinchProcessor<Grid>::interface_(
const Grid& grid,
212 const Opm::FaceDir::DirEnum& faceDir)
214 const auto actCellIdx = getActiveCellIdx_(grid, cellIdx);
215 const auto cell_faces = Opm::UgGridHelpers::cell2Faces(grid);
216 const auto cellFacesRange = cell_faces[actCellIdx];
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;
226 OPM_THROW(std::logic_error,
"Couldn't find the face for cell ." << cellIdx);
235 inline std::vector<int> PinchProcessor<Grid>::getMinpvCells_(
const std::vector<int>& actnum,
236 const std::vector<double>& pv)
238 std::vector<int> minpvCells(pv.size(), 0);
239 for (
int idx = 0; idx < static_cast<int>(pv.size()); ++idx) {
241 if (pv[idx] < minpvValue_) {
252 inline std::vector<int> PinchProcessor<Grid>::getHfIdxMap_(
const Grid& grid)
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);
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++;
270 inline int PinchProcessor<Grid>::getActiveCellIdx_(
const Grid& grid,
273 const int nc = Opm::UgGridHelpers::numCells(grid);
274 const int* global_cell = Opm::UgGridHelpers::globalCell(grid);
276 for (
int i = 0; i < nc; ++i) {
277 if (global_cell[i] == globalIdx) {
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)
293 const int nc = Opm::UgGridHelpers::numCells(grid);
294 const int nf = Opm::UgGridHelpers::numFaces(grid);
295 std::vector<double> trans(nf, 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];
308 const int idx1 = std::distance(std::begin(pinFaces), pos);
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];
323 for (
auto f = 0; f < nf; ++f) {
324 trans[f] = 1. / trans[f];
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)
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;
348 for (
int zz = z+1; zz < dims[2]; ++zz) {
349 const int cc = getGlobalIndex_(x, y, zz, dims);
350 if (minpvCells[cc]) {
357 segment.push_back(seg);
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,
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);
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 for (
int topk = ijk1[2]-2; topk > 0; --topk) {
393 topCell = getGlobalIndex_(ijk1[0], ijk1[1], topk, dims);
394 if (actnum[topCell]) {
399 pinFaces.push_back(interface_(grid, topCell, Opm::FaceDir::ZPlus));
400 pinCells.push_back(topCell);
402 tmp.insert(tmp.begin(), topCell);
403 newSeg.push_back(tmp);
404 if (!actnum[botCell]) {
405 for (
int botk = ijk2[2]+2; botk < dims[2]; ++botk) {
406 botCell = getGlobalIndex_(ijk2[0], ijk2[1], botk, dims);
407 if (actnum[botCell]) {
413 pinFaces.push_back(interface_(grid, botCell, Opm::FaceDir::ZMinus));
414 pinCells.push_back(botCell);
419 auto faceTrans = transCompute_(grid, htrans, pinCells, pinFaces);
420 auto multzmap = multzOptions_(grid, pinCells, pinFaces, multz, newSeg);
421 applyMultz_(faceTrans, multzmap);
422 for (
int i = 0; i < static_cast<int>(pinCells.size())/2; ++i) {
423 nnc.addNNC(static_cast<int>(pinCells[2*i]), static_cast<int>(pinCells[2*i+1]), faceTrans[pinFaces[2*i]]);
431 inline std::unordered_multimap<int, double> PinchProcessor<Grid>::multzOptions_(
const Grid& grid,
432 const std::vector<int>& pinCells,
433 const std::vector<int>& pinFaces,
434 const std::vector<double>& multz,
435 const std::vector<std::vector<int> >& segs)
437 std::unordered_multimap<int, double> multzmap;
439 for (
int i = 0; i < static_cast<int>(pinFaces.size())/2; ++i) {
440 multzmap.insert(std::make_pair(pinFaces[2*i], multz[getActiveCellIdx_(grid, pinCells[2*i])]));
441 multzmap.insert(std::make_pair(pinFaces[2*i+1],multz[getActiveCellIdx_(grid, pinCells[2*i])]));
443 }
else if (multzMode_ == PinchMode::ModeEnum::ALL) {
444 for (
auto& seg : segs) {
446 auto multzValue = std::numeric_limits<double>::max();
447 for (
auto& cellIdx : seg) {
448 auto activeIdx = getActiveCellIdx_(grid, cellIdx);
449 if (activeIdx != -1) {
450 multzValue = std::min(multzValue, multz[activeIdx]);
454 auto index = std::distance(std::begin(pinCells), std::find(pinCells.begin(), pinCells.end(), seg.front()));
455 multzmap.insert(std::make_pair(pinFaces[index], multzValue));
456 multzmap.insert(std::make_pair(pinFaces[index+1], multzValue));
466 inline void PinchProcessor<Grid>::applyMultz_(std::vector<double>& trans,
467 const std::unordered_multimap<int, double>& multzmap)
469 for (
auto& x : multzmap) {
470 trans[x.first] *= x.second;
478 const std::vector<double>& htrans,
479 const std::vector<int>& actnum,
480 const std::vector<double>& multz,
481 const std::vector<double>& pv,
484 transTopbot_(grid, htrans, actnum, multz, pv, nnc);
488 #endif // OPM_PINCHPROCESSOR_HEADER_INCLUDED Definition: PinchProcessor.hpp:39
Definition: CellQuadrature.hpp:28
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)
Generate NNCs for cells which pv is less than MINPV.
Definition: PinchProcessor.hpp:477
Connection topologically parallel to I-J plane.
Definition: preprocess.h:70
PinchProcessor(const double minpvValue, const double thickness, const PinchMode::ModeEnum transMode, const PinchMode::ModeEnum multzMode)
Create a Pinch processor.
Definition: PinchProcessor.hpp:139