-
Notifications
You must be signed in to change notification settings - Fork 73
/
MinpvProcessor.hpp
272 lines (227 loc) · 13.8 KB
/
MinpvProcessor.hpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
/*
Copyright 2014 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_MINPVPROCESSOR_HEADER_INCLUDED
#define OPM_MINPVPROCESSOR_HEADER_INCLUDED
#include <opm/grid/utility/ErrorMacros.hpp>
#include <opm/grid/utility/OpmParserIncludes.hpp>
#include <array>
namespace Opm
{
/// \brief Transform a corner-point grid ZCORN field to account for MINPV processing.
class MinpvProcessor
{
public:
/// \brief Create a processor.
/// \param[in] nx logical cartesian number of cells in I-direction
/// \param[in] ny logical cartesian number of cells in J-direction
/// \param[in] nz logical cartesian number of cells in K-direction
MinpvProcessor(const int nx, const int ny, const int nz);
/// Change zcorn so that it respects the minpv property.
/// \param[in] thickness thickness of the cell
/// \param[in] z_tolerance cells with thickness below z_tolerance will be bypassed in the minpv process.
/// \param[in] pv pore volumes of all logical cartesian cells
/// \param[in] minpvv minimum pore volume to accept a cell
/// \param[in] actnum active cells, inactive cells are not considered
/// \param[in] mergeMinPVCells flag to determine whether cells below minpv
/// should be included in the cell below
/// \param[in, out] zcorn ZCORN array to be manipulated
/// After processing, all cells that have lower pore volume than minpv
/// will have the zcorn numbers changed so they are zero-thickness. Any
/// cell below will be changed to include the deleted volume if mergeMinPCCells is true
/// els the volume will be lost
std::map<int,int> 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;
private:
std::array<int,8> cornerIndices(const int i, const int j, const int k) const;
std::array<double, 8> getCellZcorn(const int i, const int j, const int k, const double* z) const;
void setCellZcorn(const int i, const int j, const int k, const std::array<double, 8>& cellz, double* z) const;
std::array<int, 3> dims_;
std::array<int, 3> delta_;
};
inline MinpvProcessor::MinpvProcessor(const int nx, const int ny, const int nz) :
dims_( {{nx,ny,nz}} ),
delta_( {{1 , 2*nx , 4*nx*ny}} )
{ }
inline std::map<int,int> MinpvProcessor::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) const
{
// Algorithm:
// 1. Process each column of cells (with same i and j
// coordinates) from top (low k) to bottom (high k).
// 2. For each cell 'c' visited, check if its pore volume
// pv[c] is less than minpvv[c] .
// 3. If below the minpv threshold, move the lower four
// zcorn associated with the cell c to coincide with
// the upper four (so it becomes degenerate).
// 4. Look for the next active cell by skipping
// inactive cells with thickness below the z_tolerance.
// 5. If mergeMinPVcells:
// is true, the higher four zcorn associated with the cell below
// is moved to these values (so it gains the deleted volume).
// is false, a nnc is created between the cell above the removed
// cell and the cell below it. Note that the connection is only
// created if the cell below and above are active
// Inactive cells with thickness below z_tolerance and cells with porv<minpv
// are bypassed.
// 6. If pinchNOGAP (only has an effect if mergeMinPVcells==false holds):
// is true active cells with porevolume less than minpvv will only be disregarded
// if their thickness is below z_tolerance and nncs will be created in this case.
// return a list of the non-neighbor connection.
std::map<int,int> nnc;
// Check for sane input sizes.
const size_t log_size = dims_[0] * dims_[1] * dims_[2];
if (pv.size() != log_size) {
OPM_THROW(std::runtime_error, "Wrong size of PORV input, must have one element per logical cartesian cell.");
}
if (!actnum.empty() && actnum.size() != log_size) {
OPM_THROW(std::runtime_error, "Wrong size of ACTNUM input, must have one element per logical cartesian cell.");
}
// Main loop.
for (int kk = 0; kk < dims_[2]; ++kk) {
for (int jj = 0; jj < dims_[1]; ++jj) {
for (int ii = 0; ii < dims_[0]; ++ii) {
const int c = ii + dims_[0] * (jj + dims_[1] * kk);
if (pv[c] < minpvv[c] && (actnum.empty() || actnum[c])) {
// Move deeper (higher k) coordinates to lower k coordinates.
// i.e remove the cell
std::array<double, 8> cz = getCellZcorn(ii, jj, kk, zcorn);
for (int count = 0; count < 4; ++count) {
cz[count + 4] = cz[count];
}
setCellZcorn(ii, jj, kk, cz, zcorn);
// Find the next cell
int kk_iter = kk + 1;
if (kk_iter == dims_[2]) // we are at the end of the pillar.
continue;
int c_below = ii + dims_[0] * (jj + dims_[1] * (kk_iter));
// bypass inactive cells with thickness less then the tolerance
while ( ((actnum.empty() || !actnum[c_below]) && (thickness[c_below] <= z_tolerance)) ){
// move these cell to the posistion of the first cell to make the
// coordinates strictly sorted
setCellZcorn(ii, jj, kk_iter, cz, zcorn);
kk_iter ++;
if (kk_iter == dims_[2])
break;
c_below = ii + dims_[0] * (jj + dims_[1] * (kk_iter));
}
if (kk_iter == dims_[2]) // we have come to the end of the pillar.
continue;
// create nnc if false or merge the cells if true
if (!mergeMinPVCells) {
// We are at the top, so no nnc is created.
if (kk == 0)
continue;
int c_above = ii + dims_[0] * (jj + dims_[1] * (kk - 1));
// Bypass inactive cells with thickness below tolerance and active cells with volume below minpv
auto above_active = actnum.empty() || actnum[c_above];
auto above_inactive = actnum.empty() || !actnum[c_above]; // \todo Kept original, but should be !actnum.empty() && !actnum[c_above]
auto above_thin = thickness[c_above] < z_tolerance;
auto above_small_pv = pv[c_above] < minpvv[c_above];
if ((above_inactive && above_thin) || (above_active && above_small_pv
&& (!pinchNOGAP || above_thin) ) ) {
for (int topk = kk - 2; topk > 0; --topk) {
c_above = ii + dims_[0] * (jj + dims_[1] * (topk));
above_active = actnum.empty() || actnum[c_above];
above_inactive = actnum.empty() || !actnum[c_above];
auto above_significant_pv = pv[c_above] > minpvv[c_above];
auto above_broad = thickness[c_above] > z_tolerance;
// \todo if condition seems wrong and should be the negation of above?
if ( (above_active && (above_significant_pv || (pinchNOGAP && above_broad) ) ) || (above_inactive && above_broad)) {
break;
}
}
}
// Bypass inactive cells with thickness below tolerance and active cells with volume below minpv
auto below_active = actnum.empty() || actnum[c_below];
auto below_inactive = actnum.empty() || !actnum[c_below]; // \todo Kept original, but should be !actnum.empty() && !actnum[c_below]
auto below_thin = thickness[c_below] < z_tolerance;
auto below_small_pv = pv[c_below] < minpvv[c];
if ((below_inactive && below_thin) || (below_active && below_small_pv
&& (!pinchNOGAP || below_thin ) ) ) {
for (int botk = kk_iter + 1; botk < dims_[2]; ++botk) {
c_below = ii + dims_[0] * (jj + dims_[1] * (botk));
below_active = actnum.empty() || actnum[c_below];
below_inactive = actnum.empty() || !actnum[c_below]; // \todo Kept original, but should be !actnum.empty() && !actnum[c_below]
auto below_significant_pv = pv[c_below] > minpvv[c_below];
auto below_broad = thickness[c_above] > z_tolerance;
// \todo if condition seems wrong and should be the negation of above?
if ( (below_active && (below_significant_pv || (pinchNOGAP && below_broad) ) ) || (below_inactive && below_broad)) {
break;
}
}
}
// Add a connection if the cell above and below is active and has porv > minpv
if ((actnum.empty() || (actnum[c_above] && actnum[c_below])) && pv[c_above] > minpvv[c_above] && pv[c_below] > minpvv[c_below]) {
nnc.insert(std::make_pair(c_above, c_below));
}
} else {
// Set lower k coordinates of cell below to upper cells's coordinates.
// i.e fill the void using the cell below
std::array<double, 8> cz_below = getCellZcorn(ii, jj, kk_iter, zcorn);
for (int count = 0; count < 4; ++count) {
cz_below[count] = cz[count];
}
setCellZcorn(ii, jj, kk_iter, cz_below, zcorn);
}
}
}
}
}
return nnc;
}
inline std::array<int,8> MinpvProcessor::cornerIndices(const int i, const int j, const int k) const
{
const int ix = 2*(i*delta_[0] + j*delta_[1] + k*delta_[2]);
std::array<int, 8> ixs = {{ ix, ix + delta_[0],
ix + delta_[1], ix + delta_[1] + delta_[0],
ix + delta_[2], ix + delta_[2] + delta_[0],
ix + delta_[2] + delta_[1], ix + delta_[2] + delta_[1] + delta_[0] }};
return ixs;
}
// Returns the eight z-values associated with a given cell.
// The ordering is such that i runs fastest. That is, with
// L = low and H = high:
// {LLL, HLL, LHL, HHL, LLH, HLH, LHH, HHH }.
inline std::array<double, 8> MinpvProcessor::getCellZcorn(const int i, const int j, const int k, const double* z) const
{
const std::array<int, 8> ixs = cornerIndices(i, j, k);
std::array<double, 8> cellz;
for (int count = 0; count < 8; ++count) {
cellz[count] = z[ixs[count]];
}
return cellz;
}
inline void MinpvProcessor::setCellZcorn(const int i, const int j, const int k, const std::array<double, 8>& cellz, double* z) const
{
const std::array<int, 8> ixs = cornerIndices(i, j, k);
for (int count = 0; count < 8; ++count) {
z[ixs[count]] = cellz[count];
}
}
} // namespace Opm
#endif // OPM_MINPVPROCESSOR_HEADER_INCLUDED