From 1c2b04f1501719deeae0d54385c03c6beba752a0 Mon Sep 17 00:00:00 2001 From: Liu Ming Date: Thu, 15 Oct 2015 12:59:10 +0800 Subject: [PATCH] add comments and some minor changes. --- opm/core/grid/PinchProcessor.hpp | 74 +++++++++++++++++++------------- 1 file changed, 45 insertions(+), 29 deletions(-) diff --git a/opm/core/grid/PinchProcessor.hpp b/opm/core/grid/PinchProcessor.hpp index dd746115..2d01f4fc 100755 --- a/opm/core/grid/PinchProcessor.hpp +++ b/opm/core/grid/PinchProcessor.hpp @@ -51,9 +51,9 @@ namespace Opm /// Generate NNCs for cells which pv is less than MINPV. /// \param[in] Grid cpgrid or unstructured grid /// \param[in] htrans half cell transmissibility - /// \param[in] multz Z+ transmissibility multiplier - /// \param[in] pv pore volume for all the cells - /// \param[in] dz dz for all the cells + /// \param[in] multz Z+ transmissibility multiplier for all active cells + /// \param[in] pv pore volume for all the cartesian cells + /// \param[in] dz dz for all the cartesian cells /// \param[in] nnc non-neighbor connection class /// Algorithm: /// 1. Mark all the cells which pv and dz less than minpvValue and thickness. @@ -94,7 +94,8 @@ namespace Opm std::vector > getPinchoutsColumn_(const Grid& grid, const std::vector& actnum, - std::vector& minpvCells); + const std::vector& pv, + const std::vector& dz); /// Get global cell index. int getGlobalIndex_(const int i, const int j, const int k, const int* dims); @@ -115,7 +116,7 @@ namespace Opm /// Get active cell index. int getActiveCellIdx_(const Grid& grid, - const int cellIdx); + const int globalIdx); /// Item 4 in PINCH keyword. void transTopbot_(const Grid& grid, @@ -140,6 +141,7 @@ namespace Opm }; + template inline PinchProcessor::PinchProcessor(const double minpv, const double thickness, @@ -153,12 +155,15 @@ namespace Opm } + template inline int PinchProcessor::getGlobalIndex_(const int i, const int j, const int k, const int* dims) { return i + dims[0] * (j + dims[1] * k); } + + template inline std::array PinchProcessor::getCartIndex_(const int idx, const int* dims) @@ -171,6 +176,8 @@ namespace Opm return ijk; } + + template inline int PinchProcessor::interface_(const Grid& grid, const int cellIdx1, @@ -191,17 +198,21 @@ namespace Opm } } - const auto dims = Opm::UgGridHelpers::cartDims(grid); - auto ijk1 = getCartIndex_(cellIdx1, dims); - auto ijk2 = getCartIndex_(cellIdx2, dims); - if (commonFace == -1) { - OPM_THROW(std::logic_error, "Couldn't find the common face for cell " << cellIdx1<< "("< inline int PinchProcessor::interface_(const Grid& grid, const int cellIdx, @@ -226,17 +237,15 @@ namespace Opm } + template inline std::vector PinchProcessor::getMinpvCells_(const Grid& grid, const std::vector& actnum, const std::vector& pv, const std::vector& dz) { - const int nc = Opm::UgGridHelpers::numCells(grid); std::vector minpvCells(pv.size(), 0); - const int* global_cell = Opm::UgGridHelpers::globalCell(grid); - for (int cellIdx = 0; cellIdx < nc; ++cellIdx) { - const int idx = global_cell[cellIdx]; + for (int idx = 0; idx < pv.size(); ++idx) { if (actnum[idx]) { if (pv[idx] < minpvValue_) { minpvCells[idx] = 1; @@ -266,13 +275,13 @@ namespace Opm template inline int PinchProcessor::getActiveCellIdx_(const Grid& grid, - const int cellIdx) + const int globalIdx) { const int nc = Opm::UgGridHelpers::numCells(grid); const int* global_cell = Opm::UgGridHelpers::globalCell(grid); int idx = -1; for (int i = 0; i < nc; ++i) { - if (global_cell[i] == cellIdx) { + if (global_cell[i] == globalIdx) { idx = i; break; } @@ -298,11 +307,12 @@ namespace Opm for (int cellIdx = 0; cellIdx < nc; ++cellIdx) { auto cellFacesRange = cell_faces[cellIdx]; for (auto cellFaceIter = cellFacesRange.begin(); cellFaceIter != cellFacesRange.end(); ++cellFaceIter, ++cellFaceIdx) { - const int faceIdx = *cellFaceIter; - if (std::find(pinFaces.begin(), pinFaces.end(), faceIdx) == pinFaces.end()) { + const int faceIdx = *cellFaceIter; + const auto pos = std::find(pinFaces.begin(), pinFaces.end(), faceIdx); + if (pos == pinFaces.end()) { trans[faceIdx] += 1. / htrans[cellFaceIdx]; } else { - const int idx1 = std::distance(std::begin(pinFaces), std::find(pinFaces.begin(), pinFaces.end(), faceIdx)); + const int idx1 = std::distance(std::begin(pinFaces), pos); int idx2; if (idx1 % 2 == 0) { idx2 = idx1 + 1; @@ -326,10 +336,12 @@ namespace Opm template inline std::vector> PinchProcessor::getPinchoutsColumn_(const Grid& grid, - const std::vector& actnum, - std::vector& minpvCells) + const std::vector& actnum, + const std::vector& pv, + const std::vector& dz) { const int* dims = Opm::UgGridHelpers::cartDims(grid); + std::vector minpvCells = getMinpvCells_(grid, actnum, pv, dz); std::vector> segment; for (int z = 0; z < dims[2]; ++z) { for (int y = 0; y < dims[1]; ++y) { @@ -371,8 +383,7 @@ namespace Opm std::vector pinFaces; std::vector pinCells; std::vector > newSeg; - std::vector minpvCells = getMinpvCells_(grid, actnum, pv, dz); - auto minpvSeg = getPinchoutsColumn_(grid, actnum, minpvCells); + auto minpvSeg = getPinchoutsColumn_(grid, actnum, pv, dz); for (auto& seg : minpvSeg) { std::array ijk1 = getCartIndex_(seg.front(), dims); std::array ijk2 = getCartIndex_(seg.back(), dims); @@ -380,8 +391,12 @@ namespace Opm if ((ijk1[2]-1) >= 0 && (ijk2[2]+1) < dims[2]) { int topCell = getGlobalIndex_(ijk1[0], ijk1[1], ijk1[2]-1, dims); int botCell = getGlobalIndex_(ijk2[0], ijk2[1], ijk2[2]+1, dims); + /// for any segments, we need to find the active top and bottom cells. + /// if the original segment's top and bottom is inactive, we need to lookup + /// the column until they're found otherwise just ignore this segment. if (!actnum[topCell]) { - for (auto topk = ijk1[2]-2; topk > 0; --topk) { + seg.insert(seg.begin(), topCell); + for (int topk = ijk1[2]-2; topk > 0; --topk) { topCell = getGlobalIndex_(ijk1[0], ijk1[1], topk, dims); if (actnum[topCell]) { break; @@ -398,7 +413,8 @@ namespace Opm newSeg.push_back(tmp); pinCells.push_back(topCell); if (!actnum[botCell]) { - for (auto botk = ijk2[2]+2; botk < dims[2]; ++botk) { + seg.push_back(botCell); + for (int botk = ijk2[2]+2; botk < dims[2]; ++botk) { botCell = getGlobalIndex_(ijk2[0], ijk2[1], botk, dims); if (actnum[botCell]) { break; @@ -417,7 +433,7 @@ namespace Opm auto faceTrans = transCompute_(grid, htrans, pinCells, pinFaces, multz); auto multzmap = multzOptions_(grid, pinCells, pinFaces, multz, newSeg); applyMultz_(faceTrans, multzmap); - for (int i = 0; i < pinCells.size()/2; ++i) { + for (int i = 0; i < static_cast(pinCells.size())/2; ++i) { nnc.addNNC(static_cast(pinCells[2*i]), static_cast(pinCells[2*i+1]), faceTrans[pinFaces[2*i]]); } } @@ -434,14 +450,12 @@ namespace Opm { std::unordered_multimap multzmap; if (multzMode_ == PinchMode::ModeEnum::TOP) { - for (int i = 0; i < pinFaces.size()/2; ++i) { + for (int i = 0; i < static_cast(pinFaces.size())/2; ++i) { multzmap.insert(std::make_pair(pinFaces[2*i], multz[getActiveCellIdx_(grid, pinCells[2*i])])); multzmap.insert(std::make_pair(pinFaces[2*i+1],multz[getActiveCellIdx_(grid, pinCells[2*i])])); } } else if (multzMode_ == PinchMode::ModeEnum::ALL) { for (auto& seg : segs) { - //find the right face. - auto index = std::distance(std::begin(pinCells), std::find(pinCells.begin(), pinCells.end(), seg.front())); //find the min multz in seg cells. auto multzValue = std::numeric_limits::max(); for (auto& cellIdx : seg) { @@ -450,6 +464,8 @@ namespace Opm multzValue = std::min(multzValue, multz[activeIdx]); } } + //find the right face. + auto index = std::distance(std::begin(pinCells), std::find(pinCells.begin(), pinCells.end(), seg.front())); multzmap.insert(std::make_pair(pinFaces[index], multzValue)); multzmap.insert(std::make_pair(pinFaces[index+1], multzValue)); }