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 <array>
26 
27 
28 namespace Opm
29 {
30 
33  {
34  public:
39  MinpvProcessor(const int nx, const int ny, const int nz);
51  int process(const std::vector<double>& pv,
52  const double minpv,
53  const std::vector<int>& actnum,
54  const bool mergeMinPVCells,
55  double* zcorn) const;
56  private:
57  std::array<int,8> cornerIndices(const int i, const int j, const int k) const;
58  std::array<double, 8> getCellZcorn(const int i, const int j, const int k, const double* z) const;
59  void setCellZcorn(const int i, const int j, const int k, const std::array<double, 8>& cellz, double* z) const;
60  std::array<int, 3> dims_;
61  std::array<int, 3> delta_;
62  };
63 
64  inline MinpvProcessor::MinpvProcessor(const int nx, const int ny, const int nz) :
65  dims_( {{nx,ny,nz}} ),
66  delta_( {{1 , 2*nx , 4*nx*ny}} )
67  { }
68 
69 
70 
71  inline int MinpvProcessor::process(const std::vector<double>& pv,
72  const double minpv,
73  const std::vector<int>& actnum,
74  const bool mergeMinPVCells,
75  double* zcorn) const
76  {
77  // Algorithm:
78  // 1. Process each column of cells (with same i and j
79  // coordinates) from top (low k) to bottom (high k).
80  // 2. For each cell 'c' visited, check if its pore volume
81  // pv[c] is less than minpv.
82  // 3. If below the minpv threshold, move the lower four
83  // zcorn associated with the cell c to coincide with
84  // the upper four (so it becomes degenerate). If mergeMinPVcells
85  // is true, the higher four zcorn associated with the cell below
86  // is moved to these values (so it gains the deleted volume).
87 
88  int cells_modified = 0;
89 
90  // Check for sane input sizes.
91  const size_t log_size = dims_[0] * dims_[1] * dims_[2];
92  if (pv.size() != log_size) {
93  OPM_THROW(std::runtime_error, "Wrong size of PORV input, must have one element per logical cartesian cell.");
94  }
95  if (!actnum.empty() && actnum.size() != log_size) {
96  OPM_THROW(std::runtime_error, "Wrong size of ACTNUM input, must have one element per logical cartesian cell.");
97  }
98 
99  // Main loop.
100  for (int kk = 0; kk < dims_[2]; ++kk) {
101  for (int jj = 0; jj < dims_[1]; ++jj) {
102  for (int ii = 0; ii < dims_[0]; ++ii) {
103  const int c = ii + dims_[0] * (jj + dims_[1] * kk);
104  if (pv[c] < minpv && (actnum.empty() || actnum[c])) {
105  // Move deeper (higher k) coordinates to lower k coordinates.
106  std::array<double, 8> cz = getCellZcorn(ii, jj, kk, zcorn);
107  for (int count = 0; count < 4; ++count) {
108  cz[count + 4] = cz[count];
109  }
110  setCellZcorn(ii, jj, kk, cz, zcorn);
111 
112 
113  // optionally add removed volume to the cell below.
114  if (mergeMinPVCells) {
115  // Check if there is a cell below.
116  if (pv[c] > 0.0 && kk < dims_[2] - 1) {
117  // Set lower k coordinates of cell below to upper cells's coordinates.
118  std::array<double, 8> cz_below = getCellZcorn(ii, jj, kk + 1, zcorn);
119  for (int count = 0; count < 4; ++count) {
120  cz_below[count] = cz[count];
121  }
122  setCellZcorn(ii, jj, kk + 1, cz_below, zcorn);
123  }
124  }
125  ++cells_modified;
126  }
127  }
128  }
129  }
130  return cells_modified;
131  }
132 
133 
134 
135  inline std::array<int,8> MinpvProcessor::cornerIndices(const int i, const int j, const int k) const
136  {
137  const int ix = 2*(i*delta_[0] + j*delta_[1] + k*delta_[2]);
138  std::array<int, 8> ixs = {{ ix, ix + delta_[0],
139  ix + delta_[1], ix + delta_[1] + delta_[0],
140  ix + delta_[2], ix + delta_[2] + delta_[0],
141  ix + delta_[2] + delta_[1], ix + delta_[2] + delta_[1] + delta_[0] }};
142 
143  return ixs;
144  }
145 
146 
147 
148  // Returns the eight z-values associated with a given cell.
149  // The ordering is such that i runs fastest. That is, with
150  // L = low and H = high:
151  // {LLL, HLL, LHL, HHL, LLH, HLH, LHH, HHH }.
152  inline std::array<double, 8> MinpvProcessor::getCellZcorn(const int i, const int j, const int k, const double* z) const
153  {
154  const std::array<int, 8> ixs = cornerIndices(i, j, k);
155  std::array<double, 8> cellz;
156  for (int count = 0; count < 8; ++count) {
157  cellz[count] = z[ixs[count]];
158  }
159  return cellz;
160  }
161 
162 
163 
164  inline void MinpvProcessor::setCellZcorn(const int i, const int j, const int k, const std::array<double, 8>& cellz, double* z) const
165  {
166  const std::array<int, 8> ixs = cornerIndices(i, j, k);
167  for (int count = 0; count < 8; ++count) {
168  z[ixs[count]] = cellz[count];
169  }
170  }
171 
172 
173 
174 } // namespace Opm
175 
176 #endif // OPM_MINPVPROCESSOR_HEADER_INCLUDED
Definition: CellQuadrature.hpp:28
MinpvProcessor(const int nx, const int ny, const int nz)
Create a processor.
Definition: MinpvProcessor.hpp:64
int process(const std::vector< double > &pv, const double minpv, const std::vector< int > &actnum, const bool mergeMinPVCells, double *zcorn) const
Change zcorn so that it respects the minpv property.
Definition: MinpvProcessor.hpp:71
Transform a corner-point grid ZCORN field to account for MINPV processing.
Definition: MinpvProcessor.hpp:32