20 #ifndef OPM_CORNERPOINTCHOPPER_HEADER_INCLUDED 21 #define OPM_CORNERPOINTCHOPPER_HEADER_INCLUDED 23 #include <opm/parser/eclipse/Parser/Parser.hpp> 24 #include <opm/parser/eclipse/Parser/ParseContext.hpp> 25 #include <opm/parser/eclipse/Units/UnitSystem.hpp> 26 #include <opm/parser/eclipse/Deck/Deck.hpp> 27 #include <opm/parser/eclipse/Deck/DeckKeyword.hpp> 28 #include <opm/parser/eclipse/Deck/DeckRecord.hpp> 29 #include <opm/parser/eclipse/Deck/DeckItem.hpp> 45 Opm::ParseContext parseContext;
47 deck_ = parser.parseFile(file , parseContext);
49 metricUnits_ = Opm::UnitSystem::newMETRIC();
51 const auto& specgridRecord = deck_.getKeyword(
"SPECGRID").getRecord(0);
52 dims_[0] = specgridRecord.getItem(
"NX").get<
int >(0);
53 dims_[1] = specgridRecord.getItem(
"NY").get<
int >(0);
54 dims_[2] = specgridRecord.getItem(
"NZ").get<
int >(0);
56 int layersz = 8*dims_[0]*dims_[1];
57 const std::vector<double>& ZCORN = deck_.getKeyword(
"ZCORN").getRawDoubleData();
58 botmax_ = *std::max_element(ZCORN.begin(), ZCORN.begin() + layersz/2);
59 topmin_ = *std::min_element(ZCORN.begin() + dims_[2]*layersz - layersz/2,
60 ZCORN.begin() + dims_[2]*layersz);
62 abszmax_ = *std::max_element(ZCORN.begin(), ZCORN.end());
63 abszmin_ = *std::min_element(ZCORN.begin(), ZCORN.end());
65 std::cout <<
"Parsed grdecl file with dimensions (" 66 << dims_[0] <<
", " << dims_[1] <<
", " << dims_[2] <<
")" << std::endl;
72 const int* dimensions()
const 80 const int* newDimensions()
const 88 const std::pair<double, double> zLimits()
const 90 return std::make_pair(botmax_, topmin_);
93 const std::pair<double, double> abszLimits()
const 95 return std::make_pair(abszmin_, abszmax_);
99 void verifyInscribedShoebox(
int imin,
int ilen,
int imax,
100 int jmin,
int jlen,
int jmax,
101 double zmin,
double zlen,
double zmax)
104 std::cerr <<
"Error! imin < 0 (imin = " << imin <<
")\n";
105 throw std::runtime_error(
"Inconsistent user input.");
107 if (ilen > dims_[0]) {
108 std::cerr <<
"Error! ilen larger than grid (ilen = " << ilen <<
")\n";
109 throw std::runtime_error(
"Inconsistent user input.");
111 if (imax > dims_[0]) {
112 std::cerr <<
"Error! imax larger than input grid (imax = " << imax <<
")\n";
113 throw std::runtime_error(
"Inconsistent user input.");
116 std::cerr <<
"Error! jmin < 0 (jmin = " << jmin <<
")\n";
117 throw std::runtime_error(
"Inconsistent user input.");
119 if (jlen > dims_[1]) {
120 std::cerr <<
"Error! jlen larger than grid (jlen = " << jlen <<
")\n";
121 throw std::runtime_error(
"Inconsistent user input.");
123 if (jmax > dims_[1]) {
124 std::cerr <<
"Error! jmax larger than input grid (jmax = " << jmax <<
")\n";
125 throw std::runtime_error(
"Inconsistent user input.");
127 if (zmin < abszmin_) {
128 std::cerr <<
"Error! zmin ("<< zmin <<
") less than minimum ZCORN value ("<< abszmin_ <<
")\n";
129 throw std::runtime_error(
"Inconsistent user input.");
131 if (zmax > abszmax_) {
132 std::cerr <<
"Error! zmax ("<< zmax <<
") larger than maximal ZCORN value ("<< abszmax_ <<
")\n";
133 throw std::runtime_error(
"Inconsistent user input.");
135 if (zlen > (abszmax_ - abszmin_)) {
136 std::cerr <<
"Error! zlen ("<< zlen <<
") larger than maximal ZCORN (" << abszmax_ <<
") minus minimal ZCORN ("<< abszmin_ <<
")\n";
137 throw std::runtime_error(
"Inconsistent user input.");
141 void chop(
int imin,
int imax,
int jmin,
int jmax,
double zmin,
double zmax,
bool resettoorigin=
true)
143 new_dims_[0] = imax - imin;
144 new_dims_[1] = jmax - jmin;
147 const std::vector<double>& COORD = deck_.getKeyword(
"COORD").getRawDoubleData();
148 int num_coord = COORD.size();
149 if (num_coord != 6*(dims_[0] + 1)*(dims_[1] + 1)) {
150 std::cerr <<
"Error! COORD size (" << COORD.size() <<
") not consistent with SPECGRID\n";
151 throw std::runtime_error(
"Inconsistent COORD and SPECGRID.");
153 int num_new_coord = 6*(new_dims_[0] + 1)*(new_dims_[1] + 1);
154 double x_correction = COORD[6*((dims_[0] + 1)*jmin + imin)];
155 double y_correction = COORD[6*((dims_[0] + 1)*jmin + imin) + 1];
156 new_COORD_.resize(num_new_coord, 1e100);
157 for (
int j = jmin; j < jmax + 1; ++j) {
158 for (
int i = imin; i < imax + 1; ++i) {
159 int pos = (dims_[0] + 1)*j + i;
160 int new_pos = (new_dims_[0] + 1)*(j-jmin) + (i-imin);
162 std::copy(COORD.begin() + 6*pos, COORD.begin() + 6*(pos + 1), new_COORD_.begin() + 6*new_pos);
165 new_COORD_[6*new_pos] -= x_correction;
166 new_COORD_[6*new_pos + 1] -= y_correction;
167 new_COORD_[6*new_pos + 2] = 0;
168 new_COORD_[6*new_pos + 3] -= x_correction;
169 new_COORD_[6*new_pos + 4] -= y_correction;
170 new_COORD_[6*new_pos + 5] = zmax-zmin;
179 int layersz = 8*dims_[0]*dims_[1];
180 const std::vector<double>& ZCORN = deck_.getKeyword(
"ZCORN").getRawDoubleData();
181 int num_zcorn = ZCORN.size();
182 if (num_zcorn != layersz*dims_[2]) {
183 std::cerr <<
"Error! ZCORN size (" << ZCORN.size() <<
") not consistent with SPECGRID\n";
184 throw std::runtime_error(
"Inconsistent ZCORN and SPECGRID.");
187 zmin = std::max(zmin, botmax_);
188 zmax = std::min(zmax, topmin_);
190 std::cerr <<
"Error: zmin >= zmax (zmin = " << zmin <<
", zmax = " << zmax <<
")\n";
191 throw std::runtime_error(
"zmin >= zmax");
193 std::cout <<
"Chopping subsample, i: (" << imin <<
"--" << imax <<
") j: (" << jmin <<
"--" << jmax <<
") z: (" << zmin <<
"--" << zmax <<
")" << std::endl;
198 for (
int k = 0; k < dims_[2]; ++k) {
199 double layer_max = *std::max_element(ZCORN.begin() + k*layersz, ZCORN.begin() + (k + 1)*layersz);
200 if (layer_max > zmin) {
207 for (
int k = dims_[2]; k > 0; --k) {
208 double layer_min = *std::min_element(ZCORN.begin() + (k - 1)*layersz, ZCORN.begin() + k*layersz);
209 if (layer_min < zmax) {
214 new_dims_[2] = kmax - kmin;
217 double z_origin_correction = 0.0;
219 z_origin_correction = zmin;
221 new_ZCORN_.resize(8*new_dims_[0]*new_dims_[1]*new_dims_[2], 1e100);
222 new_to_old_cell_.resize(new_dims_[0]*new_dims_[1]*new_dims_[2], -1);
224 int delta[3] = { 1, 2*dims_[0], 4*dims_[0]*dims_[1] };
225 int new_delta[3] = { 1, 2*new_dims_[0], 4*new_dims_[0]*new_dims_[1] };
226 for (
int k = kmin; k < kmax; ++k) {
227 for (
int j = jmin; j < jmax; ++j) {
228 for (
int i = imin; i < imax; ++i) {
229 new_to_old_cell_[cellcount++] = dims_[0]*dims_[1]*k + dims_[0]*j + i;
230 int old_ix = 2*(i*delta[0] + j*delta[1] + k*delta[2]);
231 int new_ix = 2*((i-imin)*new_delta[0] + (j-jmin)*new_delta[1] + (k-kmin)*new_delta[2]);
232 int old_indices[8] = { old_ix, old_ix + delta[0],
233 old_ix + delta[1], old_ix + delta[1] + delta[0],
234 old_ix + delta[2], old_ix + delta[2] + delta[0],
235 old_ix + delta[2] + delta[1], old_ix + delta[2] + delta[1] + delta[0] };
236 int new_indices[8] = { new_ix, new_ix + new_delta[0],
237 new_ix + new_delta[1], new_ix + new_delta[1] + new_delta[0],
238 new_ix + new_delta[2], new_ix + new_delta[2] + new_delta[0],
239 new_ix + new_delta[2] + new_delta[1], new_ix + new_delta[2] + new_delta[1] + new_delta[0] };
240 for (
int cc = 0; cc < 8; ++cc) {
241 new_ZCORN_[new_indices[cc]] = std::min(zmax, std::max(zmin, ZCORN[old_indices[cc]])) - z_origin_correction;
247 filterIntegerField(
"ACTNUM", new_ACTNUM_);
248 filterDoubleField(
"PORO", new_PORO_);
249 filterDoubleField(
"NTG", new_NTG_);
250 filterDoubleField(
"SWCR", new_SWCR_);
251 filterDoubleField(
"SOWCR", new_SOWCR_);
252 filterDoubleField(
"PERMX", new_PERMX_);
253 filterDoubleField(
"PERMY", new_PERMY_);
254 filterDoubleField(
"PERMZ", new_PERMZ_);
255 filterIntegerField(
"SATNUM", new_SATNUM_);
263 Opm::DeckKeyword specGridKw(
"SPECGRID");
264 Opm::DeckRecord specGridRecord;
266 auto nxItem = Opm::DeckItem(
"NX",
int());
267 auto nyItem = Opm::DeckItem(
"NY",
int());
268 auto nzItem = Opm::DeckItem(
"NZ",
int());
269 auto numresItem = Opm::DeckItem(
"NUMRES",
int());
270 auto coordTypeItem = Opm::DeckItem(
"COORD_TYPE", std::string());
272 nxItem.push_back(new_dims_[0]);
273 nyItem.push_back(new_dims_[1]);
274 nzItem.push_back(new_dims_[2]);
275 numresItem.push_back(1);
276 coordTypeItem.push_back(
"F");
278 specGridRecord.addItem(std::move(nxItem));
279 specGridRecord.addItem(std::move(nyItem));
280 specGridRecord.addItem(std::move(nzItem));
281 specGridRecord.addItem(std::move(numresItem));
282 specGridRecord.addItem(std::move(coordTypeItem));
284 specGridKw.addRecord(std::move(specGridRecord));
286 subDeck.addKeyword(std::move(specGridKw));
287 addDoubleKeyword_(
subDeck,
"COORD",
"Length", new_COORD_);
288 addDoubleKeyword_(
subDeck,
"ZCORN",
"Length", new_ZCORN_);
289 addIntKeyword_(
subDeck,
"ACTNUM", new_ACTNUM_);
290 addDoubleKeyword_(
subDeck,
"PORO",
"1", new_PORO_);
291 addDoubleKeyword_(
subDeck,
"NTG",
"1", new_NTG_);
292 addDoubleKeyword_(
subDeck,
"SWCR",
"1", new_SWCR_);
293 addDoubleKeyword_(
subDeck,
"SOWCR",
"1", new_SOWCR_);
294 addDoubleKeyword_(
subDeck,
"PERMX",
"Permeability", new_PERMX_);
295 addDoubleKeyword_(
subDeck,
"PERMY",
"Permeability", new_PERMY_);
296 addDoubleKeyword_(
subDeck,
"PERMZ",
"Permeability", new_PERMZ_);
297 addIntKeyword_(
subDeck,
"SATNUM", new_SATNUM_);
300 void writeGrdecl(
const std::string& filename)
303 std::ofstream out(filename.c_str());
305 std::cerr <<
"Could not open file " << filename <<
"\n";
306 throw std::runtime_error(
"Could not open output file.");
308 out <<
"SPECGRID\n" << new_dims_[0] <<
' ' << new_dims_[1] <<
' ' << new_dims_[2]
312 out.setf(std::ios::scientific);
314 outputField(out, new_COORD_,
"COORD", 3);
315 outputField(out, new_ZCORN_,
"ZCORN", 4);
316 outputField(out, new_ACTNUM_,
"ACTNUM");
317 outputField(out, new_PORO_,
"PORO", 4);
318 if (hasNTG()) {outputField(out, new_NTG_,
"NTG", 4);}
319 if (hasSWCR()) {outputField(out, new_SWCR_,
"SWCR", 4);}
320 if (hasSOWCR()) {outputField(out, new_SOWCR_,
"SOWCR", 4);}
321 outputField(out, new_PERMX_,
"PERMX", 4);
322 outputField(out, new_PERMY_,
"PERMY", 4);
323 outputField(out, new_PERMZ_,
"PERMZ", 4);
324 outputField(out, new_SATNUM_,
"SATNUM");
326 bool hasNTG()
const {
return !new_NTG_.empty(); }
327 bool hasSWCR()
const {
return !new_SWCR_.empty(); }
328 bool hasSOWCR()
const {
return !new_SOWCR_.empty(); }
332 Opm::UnitSystem metricUnits_;
338 std::vector<double> new_COORD_;
339 std::vector<double> new_ZCORN_;
340 std::vector<int> new_ACTNUM_;
341 std::vector<double> new_PORO_;
342 std::vector<double> new_NTG_;
343 std::vector<double> new_SWCR_;
344 std::vector<double> new_SOWCR_;
345 std::vector<double> new_PERMX_;
346 std::vector<double> new_PERMY_;
347 std::vector<double> new_PERMZ_;
348 std::vector<int> new_SATNUM_;
351 std::vector<int> new_to_old_cell_;
353 void addDoubleKeyword_(Deck&
subDeck,
354 const std::string& keywordName,
355 const std::string& dimensionString,
356 const std::vector<double>& data)
361 Opm::DeckKeyword dataKw(keywordName);
362 Opm::DeckRecord dataRecord;
363 auto dataItem = Opm::DeckItem(
"DATA",
double());
365 for (
size_t i = 0; i < data.size(); ++i) {
366 dataItem.push_back(data[i]);
369 auto dimension = metricUnits_.parse(dimensionString);
370 dataItem.push_backDimension(dimension, dimension);
372 dataRecord.addItem(std::move(dataItem));
373 dataKw.addRecord(std::move(dataRecord));
374 subDeck.addKeyword(std::move(dataKw));
377 void addIntKeyword_(Deck&
subDeck,
378 const std::string& keywordName,
379 const std::vector<int>& data)
384 Opm::DeckKeyword dataKw(keywordName);
385 Opm::DeckRecord dataRecord;
386 auto dataItem = Opm::DeckItem(
"DATA",
int());
388 for (
size_t i = 0; i < data.size(); ++i) {
389 dataItem.push_back(data[i]);
392 dataRecord.addItem(std::move(dataItem));
393 dataKw.addRecord(std::move(dataRecord));
394 subDeck.addKeyword(std::move(dataKw));
397 template <
typename T>
398 void outputField(std::ostream& os,
399 const std::vector<T>& field,
400 const std::string& keyword,
401 const typename std::vector<T>::size_type nl = 20)
403 if (field.empty())
return;
405 os << keyword <<
'\n';
407 typedef typename std::vector<T>::size_type sz_t;
409 const sz_t n = field.size();
410 for (sz_t i = 0; i < n; ++i) {
412 << (((i + 1) % nl == 0) ?
'\n' :
' ');
422 template <
typename T>
423 void filterField(
const std::vector<T>& field,
424 std::vector<T>& output_field)
426 int sz = new_to_old_cell_.size();
427 output_field.resize(sz);
428 for (
int i = 0; i < sz; ++i) {
429 output_field[i] = field[new_to_old_cell_[i]];
433 void filterDoubleField(
const std::string& keyword, std::vector<double>& output_field)
435 if (deck_.hasKeyword(keyword)) {
436 const std::vector<double>& field = deck_.getKeyword(keyword).getRawDoubleData();
437 filterField(field, output_field);
441 void filterIntegerField(
const std::string& keyword, std::vector<int>& output_field)
443 if (deck_.hasKeyword(keyword)) {
444 const std::vector<int>& field = deck_.getKeyword(keyword).getIntData();
445 filterField(field, output_field);
456 #endif // OPM_CORNERPOINTCHOPPER_HEADER_INCLUDED Class for immiscible dead oil and dry gas.
Definition: applier.hpp:18
Definition: CornerpointChopper.hpp:40
Opm::Deck subDeck()
Return a sub-deck with fields corresponding to the selected subset.
Definition: CornerpointChopper.hpp:259