All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
CornerpointChopper.hpp
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/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>
30 #include <iostream>
31 #include <fstream>
32 #include <string>
33 #include <vector>
34 #include <stdexcept>
35 #include <memory>
36 
37 namespace Opm
38 {
39 
41  {
42  public:
43  CornerPointChopper(const std::string& file)
44  {
45  Opm::ParseContext parseContext;
46  Opm::Parser parser;
47  deck_ = parser.parseFile(file , parseContext);
48 
49  metricUnits_ = Opm::UnitSystem::newMETRIC();
50 
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);
55 
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);
61 
62  abszmax_ = *std::max_element(ZCORN.begin(), ZCORN.end());
63  abszmin_ = *std::min_element(ZCORN.begin(), ZCORN.end());
64 
65  std::cout << "Parsed grdecl file with dimensions ("
66  << dims_[0] << ", " << dims_[1] << ", " << dims_[2] << ")" << std::endl;
67  }
68 
69 
70 
71 
72  const int* dimensions() const
73  {
74  return dims_;
75  }
76 
77 
78 
79 
80  const int* newDimensions() const
81  {
82  return new_dims_;
83  }
84 
85 
86 
87 
88  const std::pair<double, double> zLimits() const
89  {
90  return std::make_pair(botmax_, topmin_);
91  }
92 
93  const std::pair<double, double> abszLimits() const
94  {
95  return std::make_pair(abszmin_, abszmax_);
96  }
97 
98 
99  void verifyInscribedShoebox(int imin, int ilen, int imax,
100  int jmin, int jlen, int jmax,
101  double zmin, double zlen, double zmax)
102  {
103  if (imin < 0) {
104  std::cerr << "Error! imin < 0 (imin = " << imin << ")\n";
105  throw std::runtime_error("Inconsistent user input.");
106  }
107  if (ilen > dims_[0]) {
108  std::cerr << "Error! ilen larger than grid (ilen = " << ilen <<")\n";
109  throw std::runtime_error("Inconsistent user input.");
110  }
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.");
114  }
115  if (jmin < 0) {
116  std::cerr << "Error! jmin < 0 (jmin = " << jmin << ")\n";
117  throw std::runtime_error("Inconsistent user input.");
118  }
119  if (jlen > dims_[1]) {
120  std::cerr << "Error! jlen larger than grid (jlen = " << jlen <<")\n";
121  throw std::runtime_error("Inconsistent user input.");
122  }
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.");
126  }
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.");
130  }
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.");
134  }
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.");
138  }
139  }
140 
141  void chop(int imin, int imax, int jmin, int jmax, double zmin, double zmax, bool resettoorigin=true)
142  {
143  new_dims_[0] = imax - imin;
144  new_dims_[1] = jmax - jmin;
145 
146  // Filter the coord field
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.");
152  }
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);
161  // Copy all 6 coordinates for a pillar.
162  std::copy(COORD.begin() + 6*pos, COORD.begin() + 6*(pos + 1), new_COORD_.begin() + 6*new_pos);
163  if (resettoorigin) {
164  // Substract lowest x value from all X-coords, similarly for y, and truncate in z-direction
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;
171  }
172  }
173  }
174 
175  // Get the z limits, check if they must be changed to make a shoe-box.
176  // This means that zmin must be greater than or equal to the highest
177  // coordinate of the bottom surface, while zmax must be less than or
178  // equal to the lowest coordinate of the top surface.
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.");
185  }
186 
187  zmin = std::max(zmin, botmax_);
188  zmax = std::min(zmax, topmin_);
189  if (zmin >= zmax) {
190  std::cerr << "Error: zmin >= zmax (zmin = " << zmin << ", zmax = " << zmax << ")\n";
191  throw std::runtime_error("zmin >= zmax");
192  }
193  std::cout << "Chopping subsample, i: (" << imin << "--" << imax << ") j: (" << jmin << "--" << jmax << ") z: (" << zmin << "--" << zmax << ")" << std::endl;
194 
195  // We must find the maximum and minimum k value for the given z limits.
196  // First, find the first layer with a z-coordinate strictly above zmin.
197  int kmin = -1;
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) {
201  kmin = k;
202  break;
203  }
204  }
205  // Then, find the last layer with a z-coordinate strictly below zmax.
206  int kmax = -1;
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) {
210  kmax = k;
211  break;
212  }
213  }
214  new_dims_[2] = kmax - kmin;
215 
216  // Filter the ZCORN field, build mapping from new to old cells.
217  double z_origin_correction = 0.0;
218  if (resettoorigin) {
219  z_origin_correction = zmin;
220  }
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);
223  int cellcount = 0;
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;
242  }
243  }
244  }
245  }
246 
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_);
256  }
257 
259  Opm::Deck subDeck()
260  {
261  Opm::Deck subDeck;
262 
263  Opm::DeckKeyword specGridKw("SPECGRID");
264  Opm::DeckRecord specGridRecord;
265 
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());
271 
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");
277 
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));
283 
284  specGridKw.addRecord(std::move(specGridRecord));
285 
286  subDeck.addKeyword(std::move(specGridKw));
287  addDoubleKeyword_(subDeck, "COORD", /*dimension=*/"Length", new_COORD_);
288  addDoubleKeyword_(subDeck, "ZCORN", /*dimension=*/"Length", new_ZCORN_);
289  addIntKeyword_(subDeck, "ACTNUM", new_ACTNUM_);
290  addDoubleKeyword_(subDeck, "PORO", /*dimension=*/"1", new_PORO_);
291  addDoubleKeyword_(subDeck, "NTG", /*dimension=*/"1", new_NTG_);
292  addDoubleKeyword_(subDeck, "SWCR", /*dimension=*/"1", new_SWCR_);
293  addDoubleKeyword_(subDeck, "SOWCR", /*dimension=*/"1", new_SOWCR_);
294  addDoubleKeyword_(subDeck, "PERMX", /*dimension=*/"Permeability", new_PERMX_);
295  addDoubleKeyword_(subDeck, "PERMY", /*dimension=*/"Permeability", new_PERMY_);
296  addDoubleKeyword_(subDeck, "PERMZ", /*dimension=*/"Permeability", new_PERMZ_);
297  addIntKeyword_(subDeck, "SATNUM", new_SATNUM_);
298  return subDeck;
299  }
300  void writeGrdecl(const std::string& filename)
301  {
302  // Output new versions of SPECGRID, COORD, ZCORN, ACTNUM, PERMX, PORO, SATNUM.
303  std::ofstream out(filename.c_str());
304  if (!out) {
305  std::cerr << "Could not open file " << filename << "\n";
306  throw std::runtime_error("Could not open output file.");
307  }
308  out << "SPECGRID\n" << new_dims_[0] << ' ' << new_dims_[1] << ' ' << new_dims_[2]
309  << " 1 F\n/\n\n";
310 
311  out.precision(15);
312  out.setf(std::ios::scientific);
313 
314  outputField(out, new_COORD_, "COORD", /* nl = */ 3);
315  outputField(out, new_ZCORN_, "ZCORN", /* nl = */ 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");
325  }
326  bool hasNTG() const {return !new_NTG_.empty(); }
327  bool hasSWCR() const {return !new_SWCR_.empty(); }
328  bool hasSOWCR() const {return !new_SOWCR_.empty(); }
329 
330  private:
331  Opm::Deck deck_;
332  Opm::UnitSystem metricUnits_;
333 
334  double botmax_;
335  double topmin_;
336  double abszmin_;
337  double abszmax_;
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_;
349  int dims_[3];
350  int new_dims_[3];
351  std::vector<int> new_to_old_cell_;
352 
353  void addDoubleKeyword_(Deck& subDeck,
354  const std::string& keywordName,
355  const std::string& dimensionString,
356  const std::vector<double>& data)
357  {
358  if (data.empty())
359  return;
360 
361  Opm::DeckKeyword dataKw(keywordName);
362  Opm::DeckRecord dataRecord;
363  auto dataItem = Opm::DeckItem("DATA", double());
364 
365  for (size_t i = 0; i < data.size(); ++i) {
366  dataItem.push_back(data[i]);
367  }
368 
369  auto dimension = metricUnits_.parse(dimensionString);
370  dataItem.push_backDimension(/*active=*/dimension, /*default=*/dimension);
371 
372  dataRecord.addItem(std::move(dataItem));
373  dataKw.addRecord(std::move(dataRecord));
374  subDeck.addKeyword(std::move(dataKw));
375  }
376 
377  void addIntKeyword_(Deck& subDeck,
378  const std::string& keywordName,
379  const std::vector<int>& data)
380  {
381  if (data.empty())
382  return;
383 
384  Opm::DeckKeyword dataKw(keywordName);
385  Opm::DeckRecord dataRecord;
386  auto dataItem = Opm::DeckItem("DATA", int());
387 
388  for (size_t i = 0; i < data.size(); ++i) {
389  dataItem.push_back(data[i]);
390  }
391 
392  dataRecord.addItem(std::move(dataItem));
393  dataKw.addRecord(std::move(dataRecord));
394  subDeck.addKeyword(std::move(dataKw));
395  }
396 
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)
402  {
403  if (field.empty()) return;
404 
405  os << keyword << '\n';
406 
407  typedef typename std::vector<T>::size_type sz_t;
408 
409  const sz_t n = field.size();
410  for (sz_t i = 0; i < n; ++i) {
411  os << field[i]
412  << (((i + 1) % nl == 0) ? '\n' : ' ');
413  }
414  if (n % nl != 0) {
415  os << '\n';
416  }
417  os << "/\n\n";
418  }
419 
420 
421 
422  template <typename T>
423  void filterField(const std::vector<T>& field,
424  std::vector<T>& output_field)
425  {
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]];
430  }
431  }
432 
433  void filterDoubleField(const std::string& keyword, std::vector<double>& output_field)
434  {
435  if (deck_.hasKeyword(keyword)) {
436  const std::vector<double>& field = deck_.getKeyword(keyword).getRawDoubleData();
437  filterField(field, output_field);
438  }
439  }
440 
441  void filterIntegerField(const std::string& keyword, std::vector<int>& output_field)
442  {
443  if (deck_.hasKeyword(keyword)) {
444  const std::vector<int>& field = deck_.getKeyword(keyword).getIntData();
445  filterField(field, output_field);
446  }
447  }
448 
449  };
450 
451 }
452 
453 
454 
455 
456 #endif // OPM_CORNERPOINTCHOPPER_HEADER_INCLUDED
Definition: CornerpointChopper.hpp:40
Opm::Deck subDeck()
Return a sub-deck with fields corresponding to the selected subset.
Definition: CornerpointChopper.hpp:259