CornerpointChopper.hpp
Go to the documentation of this file.
1/*
2 Copyright 2010 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_CORNERPOINTCHOPPER_HEADER_INCLUDED
21#define OPM_CORNERPOINTCHOPPER_HEADER_INCLUDED
22
23#include <opm/input/eclipse/Parser/Parser.hpp>
24#include <opm/input/eclipse/Units/UnitSystem.hpp>
25#include <opm/input/eclipse/Deck/Deck.hpp>
26#include <opm/input/eclipse/Deck/DeckKeyword.hpp>
27#include <opm/input/eclipse/Deck/DeckRecord.hpp>
28#include <opm/input/eclipse/Deck/DeckItem.hpp>
29#include <iostream>
30#include <fstream>
31#include <string>
32#include <vector>
33#include <stdexcept>
34#include <memory>
35
36namespace Opm
37{
38
39 class CornerPointChopper
40 {
41 public:
42 CornerPointChopper(const std::string& file) :
43 parser_(Parser{}),
44 deck_(parser_.parseFile(file)),
45 metricUnits_(Opm::UnitSystem::newMETRIC())
46 {
47 const auto& specgridRecord = deck_["SPECGRID"].back().getRecord(0);
48 dims_[0] = specgridRecord.getItem("NX").get< int >(0);
49 dims_[1] = specgridRecord.getItem("NY").get< int >(0);
50 dims_[2] = specgridRecord.getItem("NZ").get< int >(0);
51
52 int layersz = 8*dims_[0]*dims_[1];
53 const std::vector<double>& ZCORN = deck_["ZCORN"].back().getRawDoubleData();
54 botmax_ = *std::max_element(ZCORN.begin(), ZCORN.begin() + layersz/2);
55 topmin_ = *std::min_element(ZCORN.begin() + dims_[2]*layersz - layersz/2,
56 ZCORN.begin() + dims_[2]*layersz);
57
58 abszmax_ = *std::max_element(ZCORN.begin(), ZCORN.end());
59 abszmin_ = *std::min_element(ZCORN.begin(), ZCORN.end());
60
61 std::cout << "Parsed grdecl file with dimensions ("
62 << dims_[0] << ", " << dims_[1] << ", " << dims_[2] << ")" << std::endl;
63 }
64
65
66
67
68 const int* dimensions() const
69 {
70 return dims_;
71 }
72
73
74
75
76 const int* newDimensions() const
77 {
78 return new_dims_;
79 }
80
81
82
83
84 const std::pair<double, double> zLimits() const
85 {
86 return std::make_pair(botmax_, topmin_);
87 }
88
89 const std::pair<double, double> abszLimits() const
90 {
91 return std::make_pair(abszmin_, abszmax_);
92 }
93
94
95 void verifyInscribedShoebox(int imin, int ilen, int imax,
96 int jmin, int jlen, int jmax,
97 double zmin, double zlen, double zmax)
98 {
99 if (imin < 0) {
100 std::cerr << "Error! imin < 0 (imin = " << imin << ")\n";
101 throw std::runtime_error("Inconsistent user input.");
102 }
103 if (ilen > dims_[0]) {
104 std::cerr << "Error! ilen larger than grid (ilen = " << ilen <<")\n";
105 throw std::runtime_error("Inconsistent user input.");
106 }
107 if (imax > dims_[0]) {
108 std::cerr << "Error! imax larger than input grid (imax = " << imax << ")\n";
109 throw std::runtime_error("Inconsistent user input.");
110 }
111 if (jmin < 0) {
112 std::cerr << "Error! jmin < 0 (jmin = " << jmin << ")\n";
113 throw std::runtime_error("Inconsistent user input.");
114 }
115 if (jlen > dims_[1]) {
116 std::cerr << "Error! jlen larger than grid (jlen = " << jlen <<")\n";
117 throw std::runtime_error("Inconsistent user input.");
118 }
119 if (jmax > dims_[1]) {
120 std::cerr << "Error! jmax larger than input grid (jmax = " << jmax << ")\n";
121 throw std::runtime_error("Inconsistent user input.");
122 }
123 if (zmin < abszmin_) {
124 std::cerr << "Error! zmin ("<< zmin << ") less than minimum ZCORN value ("<< abszmin_ << ")\n";
125 throw std::runtime_error("Inconsistent user input.");
126 }
127 if (zmax > abszmax_) {
128 std::cerr << "Error! zmax ("<< zmax << ") larger than maximal ZCORN value ("<< abszmax_ << ")\n";
129 throw std::runtime_error("Inconsistent user input.");
130 }
131 if (zlen > (abszmax_ - abszmin_)) {
132 std::cerr << "Error! zlen ("<< zlen <<") larger than maximal ZCORN (" << abszmax_ << ") minus minimal ZCORN ("<< abszmin_ <<")\n";
133 throw std::runtime_error("Inconsistent user input.");
134 }
135 }
136
137 void chop(int imin, int imax, int jmin, int jmax, double zmin, double zmax, bool resettoorigin=true)
138 {
139 new_dims_[0] = imax - imin;
140 new_dims_[1] = jmax - jmin;
141
142 // Filter the coord field
143 const std::vector<double>& COORD = deck_["COORD"].back().getRawDoubleData();
144 int num_coord = COORD.size();
145 if (num_coord != 6*(dims_[0] + 1)*(dims_[1] + 1)) {
146 std::cerr << "Error! COORD size (" << COORD.size() << ") not consistent with SPECGRID\n";
147 throw std::runtime_error("Inconsistent COORD and SPECGRID.");
148 }
149 int num_new_coord = 6*(new_dims_[0] + 1)*(new_dims_[1] + 1);
150 double x_correction = COORD[6*((dims_[0] + 1)*jmin + imin)];
151 double y_correction = COORD[6*((dims_[0] + 1)*jmin + imin) + 1];
152 new_COORD_.resize(num_new_coord, 1e100);
153 for (int j = jmin; j < jmax + 1; ++j) {
154 for (int i = imin; i < imax + 1; ++i) {
155 int pos = (dims_[0] + 1)*j + i;
156 int new_pos = (new_dims_[0] + 1)*(j-jmin) + (i-imin);
157 // Copy all 6 coordinates for a pillar.
158 std::copy(COORD.begin() + 6*pos, COORD.begin() + 6*(pos + 1), new_COORD_.begin() + 6*new_pos);
159 if (resettoorigin) {
160 // Substract lowest x value from all X-coords, similarly for y, and truncate in z-direction
161 new_COORD_[6*new_pos] -= x_correction;
162 new_COORD_[6*new_pos + 1] -= y_correction;
163 new_COORD_[6*new_pos + 2] = 0;
164 new_COORD_[6*new_pos + 3] -= x_correction;
165 new_COORD_[6*new_pos + 4] -= y_correction;
166 new_COORD_[6*new_pos + 5] = zmax-zmin;
167 }
168 }
169 }
170
171 // Get the z limits, check if they must be changed to make a shoe-box.
172 // This means that zmin must be greater than or equal to the highest
173 // coordinate of the bottom surface, while zmax must be less than or
174 // equal to the lowest coordinate of the top surface.
175 int layersz = 8*dims_[0]*dims_[1];
176 const std::vector<double>& ZCORN = deck_["ZCORN"].back().getRawDoubleData();
177 int num_zcorn = ZCORN.size();
178 if (num_zcorn != layersz*dims_[2]) {
179 std::cerr << "Error! ZCORN size (" << ZCORN.size() << ") not consistent with SPECGRID\n";
180 throw std::runtime_error("Inconsistent ZCORN and SPECGRID.");
181 }
182
183 zmin = std::max(zmin, botmax_);
184 zmax = std::min(zmax, topmin_);
185 if (zmin >= zmax) {
186 std::cerr << "Error: zmin >= zmax (zmin = " << zmin << ", zmax = " << zmax << ")\n";
187 throw std::runtime_error("zmin >= zmax");
188 }
189 std::cout << "Chopping subsample, i: (" << imin << "--" << imax << ") j: (" << jmin << "--" << jmax << ") z: (" << zmin << "--" << zmax << ")" << std::endl;
190
191 // We must find the maximum and minimum k value for the given z limits.
192 // First, find the first layer with a z-coordinate strictly above zmin.
193 int kmin = -1;
194 for (int k = 0; k < dims_[2]; ++k) {
195 double layer_max = *std::max_element(ZCORN.begin() + k*layersz, ZCORN.begin() + (k + 1)*layersz);
196 if (layer_max > zmin) {
197 kmin = k;
198 break;
199 }
200 }
201 // Then, find the last layer with a z-coordinate strictly below zmax.
202 int kmax = -1;
203 for (int k = dims_[2]; k > 0; --k) {
204 double layer_min = *std::min_element(ZCORN.begin() + (k - 1)*layersz, ZCORN.begin() + k*layersz);
205 if (layer_min < zmax) {
206 kmax = k;
207 break;
208 }
209 }
210 new_dims_[2] = kmax - kmin;
211
212 // Filter the ZCORN field, build mapping from new to old cells.
213 double z_origin_correction = 0.0;
214 if (resettoorigin) {
215 z_origin_correction = zmin;
216 }
217 new_ZCORN_.resize(8*new_dims_[0]*new_dims_[1]*new_dims_[2], 1e100);
218 new_to_old_cell_.resize(new_dims_[0]*new_dims_[1]*new_dims_[2], -1);
219 int cellcount = 0;
220 int delta[3] = { 1, 2*dims_[0], 4*dims_[0]*dims_[1] };
221 int new_delta[3] = { 1, 2*new_dims_[0], 4*new_dims_[0]*new_dims_[1] };
222 for (int k = kmin; k < kmax; ++k) {
223 for (int j = jmin; j < jmax; ++j) {
224 for (int i = imin; i < imax; ++i) {
225 new_to_old_cell_[cellcount++] = dims_[0]*dims_[1]*k + dims_[0]*j + i;
226 int old_ix = 2*(i*delta[0] + j*delta[1] + k*delta[2]);
227 int new_ix = 2*((i-imin)*new_delta[0] + (j-jmin)*new_delta[1] + (k-kmin)*new_delta[2]);
228 int old_indices[8] = { old_ix, old_ix + delta[0],
229 old_ix + delta[1], old_ix + delta[1] + delta[0],
230 old_ix + delta[2], old_ix + delta[2] + delta[0],
231 old_ix + delta[2] + delta[1], old_ix + delta[2] + delta[1] + delta[0] };
232 int new_indices[8] = { new_ix, new_ix + new_delta[0],
233 new_ix + new_delta[1], new_ix + new_delta[1] + new_delta[0],
234 new_ix + new_delta[2], new_ix + new_delta[2] + new_delta[0],
235 new_ix + new_delta[2] + new_delta[1], new_ix + new_delta[2] + new_delta[1] + new_delta[0] };
236 for (int cc = 0; cc < 8; ++cc) {
237 new_ZCORN_[new_indices[cc]] = std::min(zmax, std::max(zmin, ZCORN[old_indices[cc]])) - z_origin_correction;
238 }
239 }
240 }
241 }
242
243 filterIntegerField("ACTNUM", new_ACTNUM_);
244 filterDoubleField("PORO", new_PORO_);
245 filterDoubleField("NTG", new_NTG_);
246 filterDoubleField("SWCR", new_SWCR_);
247 filterDoubleField("SOWCR", new_SOWCR_);
248 filterDoubleField("PERMX", new_PERMX_);
249 filterDoubleField("PERMY", new_PERMY_);
250 filterDoubleField("PERMZ", new_PERMZ_);
251 filterIntegerField("SATNUM", new_SATNUM_);
252 }
253
255 Opm::Deck subDeck()
256 {
257 Opm::Deck subDeck;
258
259 Opm::DeckKeyword specGridKw(this->parser_.getKeyword("SPECGRID"));
260 Opm::DeckRecord specGridRecord;
261
262 auto nxItem = Opm::DeckItem("NX", int());
263 auto nyItem = Opm::DeckItem("NY", int());
264 auto nzItem = Opm::DeckItem("NZ", int());
265 auto numresItem = Opm::DeckItem("NUMRES", int());
266 auto coordTypeItem = Opm::DeckItem("COORD_TYPE", std::string());
267
268 nxItem.push_back(new_dims_[0]);
269 nyItem.push_back(new_dims_[1]);
270 nzItem.push_back(new_dims_[2]);
271 numresItem.push_back(1);
272 coordTypeItem.push_back("F");
273
274 specGridRecord.addItem(std::move(nxItem));
275 specGridRecord.addItem(std::move(nyItem));
276 specGridRecord.addItem(std::move(nzItem));
277 specGridRecord.addItem(std::move(numresItem));
278 specGridRecord.addItem(std::move(coordTypeItem));
279
280 specGridKw.addRecord(std::move(specGridRecord));
281
282 subDeck.addKeyword(std::move(specGridKw));
283 addDoubleKeyword_(subDeck, this->parser_.getKeyword("COORD"), /*dimension=*/"Length", new_COORD_);
284 addDoubleKeyword_(subDeck, this->parser_.getKeyword("ZCORN"), /*dimension=*/"Length", new_ZCORN_);
285 addIntKeyword_(subDeck, this->parser_.getKeyword("ACTNUM"), new_ACTNUM_);
286 addDoubleKeyword_(subDeck, this->parser_.getKeyword("PORO"), /*dimension=*/"1", new_PORO_);
287 addDoubleKeyword_(subDeck, this->parser_.getKeyword("NTG"), /*dimension=*/"1", new_NTG_);
288 addDoubleKeyword_(subDeck, this->parser_.getKeyword("SWCR"), /*dimension=*/"1", new_SWCR_);
289 addDoubleKeyword_(subDeck, this->parser_.getKeyword("SOWCR"), /*dimension=*/"1", new_SOWCR_);
290 addDoubleKeyword_(subDeck, this->parser_.getKeyword("PERMX"), /*dimension=*/"Permeability", new_PERMX_);
291 addDoubleKeyword_(subDeck, this->parser_.getKeyword("PERMY"), /*dimension=*/"Permeability", new_PERMY_);
292 addDoubleKeyword_(subDeck, this->parser_.getKeyword("PERMZ"), /*dimension=*/"Permeability", new_PERMZ_);
293 addIntKeyword_(subDeck, this->parser_.getKeyword("SATNUM"), new_SATNUM_);
294 return subDeck;
295 }
296 void writeGrdecl(const std::string& filename)
297 {
298 // Output new versions of SPECGRID, COORD, ZCORN, ACTNUM, PERMX, PORO, SATNUM.
299 std::ofstream out(filename.c_str());
300 if (!out) {
301 std::cerr << "Could not open file " << filename << "\n";
302 throw std::runtime_error("Could not open output file.");
303 }
304 out << "SPECGRID\n" << new_dims_[0] << ' ' << new_dims_[1] << ' ' << new_dims_[2]
305 << " 1 F\n/\n\n";
306
307 out.precision(15);
308 out.setf(std::ios::scientific);
309
310 outputField(out, new_COORD_, "COORD", /* nl = */ 3);
311 outputField(out, new_ZCORN_, "ZCORN", /* nl = */ 4);
312 outputField(out, new_ACTNUM_, "ACTNUM");
313 outputField(out, new_PORO_, "PORO", 4);
314 if (hasNTG()) {outputField(out, new_NTG_, "NTG", 4);}
315 if (hasSWCR()) {outputField(out, new_SWCR_, "SWCR", 4);}
316 if (hasSOWCR()) {outputField(out, new_SOWCR_, "SOWCR", 4);}
317 outputField(out, new_PERMX_, "PERMX", 4);
318 outputField(out, new_PERMY_, "PERMY", 4);
319 outputField(out, new_PERMZ_, "PERMZ", 4);
320 outputField(out, new_SATNUM_, "SATNUM");
321 }
322 bool hasNTG() const {return !new_NTG_.empty(); }
323 bool hasSWCR() const {return !new_SWCR_.empty(); }
324 bool hasSOWCR() const {return !new_SOWCR_.empty(); }
325
326 private:
327 Opm::Parser parser_;
328 Opm::Deck deck_;
329 Opm::UnitSystem metricUnits_;
330
331 double botmax_;
332 double topmin_;
333 double abszmin_;
334 double abszmax_;
335 std::vector<double> new_COORD_;
336 std::vector<double> new_ZCORN_;
337 std::vector<int> new_ACTNUM_;
338 std::vector<double> new_PORO_;
339 std::vector<double> new_NTG_;
340 std::vector<double> new_SWCR_;
341 std::vector<double> new_SOWCR_;
342 std::vector<double> new_PERMX_;
343 std::vector<double> new_PERMY_;
344 std::vector<double> new_PERMZ_;
345 std::vector<int> new_SATNUM_;
346 int dims_[3];
347 int new_dims_[3];
348 std::vector<int> new_to_old_cell_;
349
350
351 void addDoubleKeyword_(Deck& subDeck,
352 const ParserKeyword& keyword,
353 const std::string& dimensionString,
354 const std::vector<double>& data)
355 {
356 if (data.empty())
357 return;
358
359 Opm::DeckKeyword dataKw(keyword);
360 Opm::DeckRecord dataRecord;
361 auto dimension = metricUnits_.parse(dimensionString);
362 auto dataItem = Opm::DeckItem("DATA", double(), { dimension }, { dimension });
363
364 for (size_t i = 0; i < data.size(); ++i) {
365 dataItem.push_back(data[i]);
366 }
367
368 dataRecord.addItem(std::move(dataItem));
369 dataKw.addRecord(std::move(dataRecord));
370 subDeck.addKeyword(std::move(dataKw));
371 }
372
373 void addIntKeyword_(Deck& subDeck,
374 const ParserKeyword& keyword,
375 const std::vector<int>& data)
376 {
377 if (data.empty())
378 return;
379
380 Opm::DeckKeyword dataKw(keyword);
381 Opm::DeckRecord dataRecord;
382 auto dataItem = Opm::DeckItem("DATA", int());
383
384 for (size_t i = 0; i < data.size(); ++i) {
385 dataItem.push_back(data[i]);
386 }
387
388 dataRecord.addItem(std::move(dataItem));
389 dataKw.addRecord(std::move(dataRecord));
390 subDeck.addKeyword(std::move(dataKw));
391 }
392
393 template <typename T>
394 void outputField(std::ostream& os,
395 const std::vector<T>& field,
396 const std::string& keyword,
397 const typename std::vector<T>::size_type nl = 20)
398 {
399 if (field.empty()) return;
400
401 os << keyword << '\n';
402
403 typedef typename std::vector<T>::size_type sz_t;
404
405 const sz_t n = field.size();
406 for (sz_t i = 0; i < n; ++i) {
407 os << field[i]
408 << (((i + 1) % nl == 0) ? '\n' : ' ');
409 }
410 if (n % nl != 0) {
411 os << '\n';
412 }
413 os << "/\n\n";
414 }
415
416
417
418 template <typename T>
419 void filterField(const std::vector<T>& field,
420 std::vector<T>& output_field)
421 {
422 int sz = new_to_old_cell_.size();
423 output_field.resize(sz);
424 for (int i = 0; i < sz; ++i) {
425 output_field[i] = field[new_to_old_cell_[i]];
426 }
427 }
428
429 void filterDoubleField(const std::string& keyword, std::vector<double>& output_field)
430 {
431 if (deck_.hasKeyword(keyword)) {
432 const std::vector<double>& field = deck_[keyword].back().getRawDoubleData();
433 filterField(field, output_field);
434 }
435 }
436
437 void filterIntegerField(const std::string& keyword, std::vector<int>& output_field)
438 {
439 if (deck_.hasKeyword(keyword)) {
440 const std::vector<int>& field = deck_[keyword].back().getIntData();
441 filterField(field, output_field);
442 }
443 }
444
445 };
446
447}
448
449
450
451
452#endif // OPM_CORNERPOINTCHOPPER_HEADER_INCLUDED
Opm::CornerPointChopper
Definition: CornerpointChopper.hpp:40
Opm::CornerPointChopper::writeGrdecl
void writeGrdecl(const std::string &filename)
Definition: CornerpointChopper.hpp:296
Opm::CornerPointChopper::newDimensions
const int * newDimensions() const
Definition: CornerpointChopper.hpp:76
Opm::CornerPointChopper::abszLimits
const std::pair< double, double > abszLimits() const
Definition: CornerpointChopper.hpp:89
Opm::CornerPointChopper::hasNTG
bool hasNTG() const
Definition: CornerpointChopper.hpp:322
Opm::CornerPointChopper::zLimits
const std::pair< double, double > zLimits() const
Definition: CornerpointChopper.hpp:84
Opm::CornerPointChopper::hasSWCR
bool hasSWCR() const
Definition: CornerpointChopper.hpp:323
Opm::CornerPointChopper::subDeck
Opm::Deck subDeck()
Return a sub-deck with fields corresponding to the selected subset.
Definition: CornerpointChopper.hpp:255
Opm::CornerPointChopper::dimensions
const int * dimensions() const
Definition: CornerpointChopper.hpp:68
Opm::CornerPointChopper::verifyInscribedShoebox
void verifyInscribedShoebox(int imin, int ilen, int imax, int jmin, int jlen, int jmax, double zmin, double zlen, double zmax)
Definition: CornerpointChopper.hpp:95
Opm::CornerPointChopper::CornerPointChopper
CornerPointChopper(const std::string &file)
Definition: CornerpointChopper.hpp:42
Opm::CornerPointChopper::chop
void chop(int imin, int imax, int jmin, int jmax, double zmin, double zmax, bool resettoorigin=true)
Definition: CornerpointChopper.hpp:137
Opm::CornerPointChopper::hasSOWCR
bool hasSOWCR() const
Definition: CornerpointChopper.hpp:324
Opm::Elasticity::min
min[0]
Definition: elasticity_upscale_impl.hpp:146
Opm::Elasticity::j
int j
Definition: elasticity_upscale_impl.hpp:658
Opm
Definition: ImplicitAssembly.hpp:43