My Project
MinpvProcessor.hpp
1/*
2 Copyright 2014 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_MINPVPROCESSOR_HEADER_INCLUDED
21#define OPM_MINPVPROCESSOR_HEADER_INCLUDED
22
23
24#include <opm/grid/utility/ErrorMacros.hpp>
25#include <opm/grid/utility/OpmParserIncludes.hpp>
26
27#include <array>
28#include <cstddef>
29#include <map>
30#include <vector>
31
32namespace Opm
33{
34
37 {
38 public:
39
40 struct Result {
41 std::vector<std::size_t> removed_cells;
42 std::map<int,int> nnc;
43
44
45 void add_nnc(int cell1, int cell2) {
46 auto key = std::min(cell1, cell2);
47 auto value = std::max(cell1,cell2);
48
49 this->nnc.insert({key, value});
50 }
51 };
52
53
58 MinpvProcessor(const int nx, const int ny, const int nz);
72 Result process(const std::vector<double>& thickness,
73 const double z_tolerance,
74 const std::vector<double>& pv,
75 const std::vector<double>& minpvv,
76 const std::vector<int>& actnum,
77 const bool mergeMinPVCells,
78 double* zcorn,
79 bool pinchNOGAP = false) const;
80 private:
81 std::array<int,8> cornerIndices(const int i, const int j, const int k) const;
82 std::array<double, 8> getCellZcorn(const int i, const int j, const int k, const double* z) const;
83 void setCellZcorn(const int i, const int j, const int k, const std::array<double, 8>& cellz, double* z) const;
84 std::array<int, 3> dims_;
85 std::array<int, 3> delta_;
86 };
87
88 inline MinpvProcessor::MinpvProcessor(const int nx, const int ny, const int nz) :
89 dims_( {{nx,ny,nz}} ),
90 delta_( {{1 , 2*nx , 4*nx*ny}} )
91 { }
92
93
94
95 inline MinpvProcessor::Result MinpvProcessor::process(const std::vector<double>& thickness,
96 const double z_tolerance,
97 const std::vector<double>& pv,
98 const std::vector<double>& minpvv,
99 const std::vector<int>& actnum,
100 const bool mergeMinPVCells,
101 double* zcorn,
102 bool pinchNOGAP) const
103 {
104 // Algorithm:
105 // 1. Process each column of cells (with same i and j
106 // coordinates) from top (low k) to bottom (high k).
107 // 2. For each cell 'c' visited, check if its pore volume
108 // pv[c] is less than minpvv[c] .
109 // 3. If below the minpv threshold, move the lower four
110 // zcorn associated with the cell c to coincide with
111 // the upper four (so it becomes degenerate).
112 // 4. Look for the next active cell by skipping
113 // inactive cells with thickness below the z_tolerance.
114 // 5. If mergeMinPVcells:
115 // is true, the higher four zcorn associated with the cell below
116 // is moved to these values (so it gains the deleted volume).
117 // is false, a nnc is created between the cell above the removed
118 // cell and the cell below it. Note that the connection is only
119 // created if the cell below and above are active
120 // Inactive cells with thickness below z_tolerance and cells with porv<minpv
121 // are bypassed.
122 // 6. If pinchNOGAP (only has an effect if mergeMinPVcells==false holds):
123 // is true active cells with porevolume less than minpvv will only be disregarded
124 // if their thickness is below z_tolerance and nncs will be created in this case.
125
126
127 Result result;
128
129 // Check for sane input sizes.
130 const size_t log_size = dims_[0] * dims_[1] * dims_[2];
131 if (pv.size() != log_size) {
132 OPM_THROW(std::runtime_error, "Wrong size of PORV input, must have one element per logical cartesian cell.");
133 }
134 if (!actnum.empty() && actnum.size() != log_size) {
135 OPM_THROW(std::runtime_error, "Wrong size of ACTNUM input, must have one element per logical cartesian cell.");
136 }
137
138 // Main loop.
139 for (int kk = 0; kk < dims_[2]; ++kk) {
140 for (int jj = 0; jj < dims_[1]; ++jj) {
141 for (int ii = 0; ii < dims_[0]; ++ii) {
142 const int c = ii + dims_[0] * (jj + dims_[1] * kk);
143 if (pv[c] < minpvv[c] && (actnum.empty() || actnum[c])) {
144 // Move deeper (higher k) coordinates to lower k coordinates.
145 // i.e remove the cell
146 std::array<double, 8> cz = getCellZcorn(ii, jj, kk, zcorn);
147 for (int count = 0; count < 4; ++count) {
148 cz[count + 4] = cz[count];
149 }
150 setCellZcorn(ii, jj, kk, cz, zcorn);
151
152 // Find the next cell
153 int kk_iter = kk + 1;
154 if (kk_iter == dims_[2]) { // we are at the end of the pillar.
155 result.removed_cells.push_back(c);
156 continue;
157 }
158
159 int c_below = ii + dims_[0] * (jj + dims_[1] * (kk_iter));
160 // bypass inactive cells with thickness less then the tolerance
161 while ( ((actnum.empty() || !actnum[c_below]) && (thickness[c_below] <= z_tolerance)) ){
162 // move these cell to the posistion of the first cell to make the
163 // coordinates strictly sorted
164 setCellZcorn(ii, jj, kk_iter, cz, zcorn);
165 kk_iter ++;
166 if (kk_iter == dims_[2])
167 break;
168
169 c_below = ii + dims_[0] * (jj + dims_[1] * (kk_iter));
170 }
171
172 if (kk_iter == dims_[2]) { // we have come to the end of the pillar.
173 result.removed_cells.push_back(c);
174 continue;
175 }
176
177 // create nnc if false or merge the cells if true
178 if (!mergeMinPVCells) {
179
180 // We are at the top, so no nnc is created.
181 if (kk == 0) {
182 result.removed_cells.push_back(c);
183 continue;
184 }
185
186 int c_above = ii + dims_[0] * (jj + dims_[1] * (kk - 1));
187
188 // Bypass inactive cells with thickness below tolerance and active cells with volume below minpv
189 auto above_active = actnum.empty() || actnum[c_above];
190 auto above_inactive = actnum.empty() || !actnum[c_above]; // \todo Kept original, but should be !actnum.empty() && !actnum[c_above]
191 auto above_thin = thickness[c_above] < z_tolerance;
192 auto above_small_pv = pv[c_above] < minpvv[c_above];
193 if ((above_inactive && above_thin) || (above_active && above_small_pv
194 && (!pinchNOGAP || above_thin) ) ) {
195 for (int topk = kk - 2; topk > 0; --topk) {
196 c_above = ii + dims_[0] * (jj + dims_[1] * (topk));
197 above_active = actnum.empty() || actnum[c_above];
198 above_inactive = actnum.empty() || !actnum[c_above];
199 auto above_significant_pv = pv[c_above] > minpvv[c_above];
200 auto above_broad = thickness[c_above] > z_tolerance;
201 // \todo if condition seems wrong and should be the negation of above?
202 if ( (above_active && (above_significant_pv || (pinchNOGAP && above_broad) ) ) || (above_inactive && above_broad)) {
203 break;
204 }
205 }
206 }
207
208 // Bypass inactive cells with thickness below tolerance and active cells with volume below minpv
209 auto below_active = actnum.empty() || actnum[c_below];
210 auto below_inactive = actnum.empty() || !actnum[c_below]; // \todo Kept original, but should be !actnum.empty() && !actnum[c_below]
211 auto below_thin = thickness[c_below] < z_tolerance;
212 auto below_small_pv = pv[c_below] < minpvv[c];
213 if ((below_inactive && below_thin) || (below_active && below_small_pv
214 && (!pinchNOGAP || below_thin ) ) ) {
215 for (int botk = kk_iter + 1; botk < dims_[2]; ++botk) {
216 c_below = ii + dims_[0] * (jj + dims_[1] * (botk));
217 below_active = actnum.empty() || actnum[c_below];
218 below_inactive = actnum.empty() || !actnum[c_below]; // \todo Kept original, but should be !actnum.empty() && !actnum[c_below]
219 auto below_significant_pv = pv[c_below] > minpvv[c_below];
220 auto below_broad = thickness[c_above] > z_tolerance;
221 // \todo if condition seems wrong and should be the negation of above?
222 if ( (below_active && (below_significant_pv || (pinchNOGAP && below_broad) ) ) || (below_inactive && below_broad)) {
223 break;
224 }
225 }
226 }
227
228 // Add a connection if the cell above and below is active and has porv > minpv
229 if ((actnum.empty() || (actnum[c_above] && actnum[c_below])) && pv[c_above] > minpvv[c_above] && pv[c_below] > minpvv[c_below]) {
230 result.add_nnc(c_above, c_below);
231 }
232 }
233 else {
234 // Set lower k coordinates of cell below to upper cells's coordinates.
235 // i.e fill the void using the cell below
236 std::array<double, 8> cz_below = getCellZcorn(ii, jj, kk_iter, zcorn);
237 for (int count = 0; count < 4; ++count) {
238 cz_below[count] = cz[count];
239 }
240
241 setCellZcorn(ii, jj, kk_iter, cz_below, zcorn);
242 }
243
244 result.removed_cells.push_back(c);
245 }
246 }
247 }
248 }
249
250 return result;
251 }
252
253
254
255 inline std::array<int,8> MinpvProcessor::cornerIndices(const int i, const int j, const int k) const
256 {
257 const int ix = 2*(i*delta_[0] + j*delta_[1] + k*delta_[2]);
258 std::array<int, 8> ixs = {{ ix, ix + delta_[0],
259 ix + delta_[1], ix + delta_[1] + delta_[0],
260 ix + delta_[2], ix + delta_[2] + delta_[0],
261 ix + delta_[2] + delta_[1], ix + delta_[2] + delta_[1] + delta_[0] }};
262
263 return ixs;
264 }
265
266
267
268 // Returns the eight z-values associated with a given cell.
269 // The ordering is such that i runs fastest. That is, with
270 // L = low and H = high:
271 // {LLL, HLL, LHL, HHL, LLH, HLH, LHH, HHH }.
272 inline std::array<double, 8> MinpvProcessor::getCellZcorn(const int i, const int j, const int k, const double* z) const
273 {
274 const std::array<int, 8> ixs = cornerIndices(i, j, k);
275 std::array<double, 8> cellz;
276 for (int count = 0; count < 8; ++count) {
277 cellz[count] = z[ixs[count]];
278 }
279 return cellz;
280 }
281
282
283
284 inline void MinpvProcessor::setCellZcorn(const int i, const int j, const int k, const std::array<double, 8>& cellz, double* z) const
285 {
286 const std::array<int, 8> ixs = cornerIndices(i, j, k);
287 for (int count = 0; count < 8; ++count) {
288 z[ixs[count]] = cellz[count];
289 }
290 }
291
292
293
294} // namespace Opm
295
296#endif // OPM_MINPVPROCESSOR_HEADER_INCLUDED
Transform a corner-point grid ZCORN field to account for MINPV processing.
Definition: MinpvProcessor.hpp:37
Result process(const std::vector< double > &thickness, const double z_tolerance, const std::vector< double > &pv, const std::vector< double > &minpvv, const std::vector< int > &actnum, const bool mergeMinPVCells, double *zcorn, bool pinchNOGAP=false) const
Change zcorn so that it respects the minpv property.
Definition: MinpvProcessor.hpp:95
MinpvProcessor(const int nx, const int ny, const int nz)
Create a processor.
Definition: MinpvProcessor.hpp:88
Holds the implementation of the CpGrid as a pimple.
Definition: CellQuadrature.hpp:29
Definition: MinpvProcessor.hpp:40