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