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