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 
32 namespace 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