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