20 #ifndef OPM_MINPVPROCESSOR_HEADER_INCLUDED
21 #define OPM_MINPVPROCESSOR_HEADER_INCLUDED
24 #include <opm/grid/utility/ErrorMacros.hpp>
25 #include <opm/grid/utility/OpmParserIncludes.hpp>
38 std::vector<std::size_t> removed_cells;
39 std::map<int,int> nnc;
42 void add_nnc(
int cell1,
int cell2) {
43 auto key = std::min(cell1, cell2);
44 auto value = std::max(cell1,cell2);
46 this->nnc.insert({key, value});
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,
76 bool pinchNOGAP =
false)
const;
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_;
86 dims_( {{nx,ny,nz}} ),
87 delta_( {{1 , 2*nx , 4*nx*ny}} )
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,
99 bool pinchNOGAP)
const
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.");
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.");
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])) {
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];
147 setCellZcorn(ii, jj, kk, cz, zcorn);
150 int kk_iter = kk + 1;
151 if (kk_iter == dims_[2]) {
152 result.removed_cells.push_back(c);
156 int c_below = ii + dims_[0] * (jj + dims_[1] * (kk_iter));
158 while ( ((actnum.empty() || !actnum[c_below]) && (thickness[c_below] <= z_tolerance)) ){
161 setCellZcorn(ii, jj, kk_iter, cz, zcorn);
163 if (kk_iter == dims_[2])
166 c_below = ii + dims_[0] * (jj + dims_[1] * (kk_iter));
169 if (kk_iter == dims_[2]) {
170 result.removed_cells.push_back(c);
175 if (!mergeMinPVCells) {
179 result.removed_cells.push_back(c);
183 int c_above = ii + dims_[0] * (jj + dims_[1] * (kk - 1));
186 auto above_active = actnum.empty() || actnum[c_above];
187 auto above_inactive = 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;
199 if ( (above_active && (above_significant_pv || (pinchNOGAP && above_broad) ) ) || (above_inactive && above_broad)) {
206 auto below_active = actnum.empty() || actnum[c_below];
207 auto below_inactive = 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];
216 auto below_significant_pv = pv[c_below] > minpvv[c_below];
217 auto below_broad = thickness[c_above] > z_tolerance;
219 if ( (below_active && (below_significant_pv || (pinchNOGAP && below_broad) ) ) || (below_inactive && below_broad)) {
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);
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];
238 setCellZcorn(ii, jj, kk_iter, cz_below, zcorn);
241 result.removed_cells.push_back(c);
252 inline std::array<int,8> MinpvProcessor::cornerIndices(
const int i,
const int j,
const int k)
const
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] }};
269 inline std::array<double, 8> MinpvProcessor::getCellZcorn(
const int i,
const int j,
const int k,
const double* z)
const
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]];
281 inline void MinpvProcessor::setCellZcorn(
const int i,
const int j,
const int k,
const std::array<double, 8>& cellz,
double* z)
const
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];
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