My Project
EclipseGrid.hpp
1 /*
2  Copyright 2014 Statoil ASA.
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_PARSER_ECLIPSE_GRID_HPP
21 #define OPM_PARSER_ECLIPSE_GRID_HPP
22 
23 #include <opm/input/eclipse/EclipseState/Grid/GridDims.hpp>
24 #include <opm/input/eclipse/EclipseState/Grid/MapAxes.hpp>
25 #include <opm/input/eclipse/EclipseState/Grid/MinpvMode.hpp>
26 #include <opm/input/eclipse/EclipseState/Grid/PinchMode.hpp>
27 
28 #include <array>
29 #include <memory>
30 #include <optional>
31 #include <stdexcept>
32 #include <unordered_set>
33 #include <vector>
34 
35 namespace Opm {
36 
37  class Deck;
38  namespace EclIO { class EclFile; }
39  struct NNCdata;
40  class UnitSystem;
41  class ZcornMapper;
42 
54  class EclipseGrid : public GridDims {
55  public:
56  EclipseGrid() = default;
57  explicit EclipseGrid(const std::string& filename);
58 
59  /*
60  These constructors will make a copy of the src grid, with
61  zcorn and or actnum have been adjustments.
62  */
63  EclipseGrid(const EclipseGrid& src, const std::vector<int>& actnum);
64  EclipseGrid(const EclipseGrid& src, const double* zcorn, const std::vector<int>& actnum);
65 
66  EclipseGrid(size_t nx, size_t ny, size_t nz,
67  double dx = 1.0, double dy = 1.0, double dz = 1.0);
68  explicit EclipseGrid(const GridDims& gd);
69 
70  EclipseGrid(const std::array<int, 3>& dims ,
71  const std::vector<double>& coord ,
72  const std::vector<double>& zcorn ,
73  const int * actnum = nullptr);
74 
75 
78  EclipseGrid(const Deck& deck, const int * actnum = nullptr);
79 
80  static bool hasGDFILE(const Deck& deck);
81  static bool hasCylindricalKeywords(const Deck& deck);
82  static bool hasCornerPointKeywords(const Deck&);
83  static bool hasCartesianKeywords(const Deck&);
84  size_t getNumActive( ) const;
85  bool allActive() const;
86 
87  size_t activeIndex(size_t i, size_t j, size_t k) const;
88  size_t activeIndex(size_t globalIndex) const;
89 
90  size_t getActiveIndex(size_t i, size_t j, size_t k) const {
91  return activeIndex(i, j, k);
92  }
93 
94  void save(const std::string& filename, bool formatted, const std::vector<Opm::NNCdata>& nnc, const Opm::UnitSystem& units) const;
95  /*
96  Observe that the there is a getGlobalIndex(i,j,k)
97  implementation in the base class. This method - translating
98  from an active index to a global index must be implemented
99  in the current class.
100  */
101  size_t getGlobalIndex(size_t active_index) const;
102  size_t getGlobalIndex(size_t i, size_t j, size_t k) const;
103 
104  /*
105  For RADIAL grids you can *optionally* use the keyword
106  'CIRCLE' to denote that period boundary conditions should be
107  applied in the 'THETA' direction; this will only apply if
108  the theta keywords entered sum up to exactly 360 degrees!
109  */
110 
111  bool circle( ) const;
112  bool isPinchActive( ) const;
113  double getPinchThresholdThickness( ) const;
114  PinchMode::ModeEnum getPinchOption( ) const;
115  PinchMode::ModeEnum getMultzOption( ) const;
116  PinchMode::ModeEnum getPinchGapMode( ) const;
117 
118  MinpvMode::ModeEnum getMinpvMode() const;
119  const std::vector<double>& getMinpvVector( ) const;
120 
121  /*
122  Will return a vector of nactive elements. The method will
123  behave differently depending on the lenght of the
124  input_vector:
125 
126  nx*ny*nz: only the values corresponding to active cells
127  are copied out.
128 
129  nactive: The input vector is copied straight out again.
130 
131  ??? : Exception.
132  */
133 
134  template<typename T>
135  std::vector<T> compressedVector(const std::vector<T>& input_vector) const {
136  if( input_vector.size() == this->getNumActive() ) {
137  return input_vector;
138  }
139 
140  if (input_vector.size() != getCartesianSize())
141  throw std::invalid_argument("Input vector must have full size");
142 
143  {
144  std::vector<T> compressed_vector( this->getNumActive() );
145  const auto& active_map = this->getActiveMap( );
146 
147  for (size_t i = 0; i < this->getNumActive(); ++i)
148  compressed_vector[i] = input_vector[ active_map[i] ];
149 
150  return compressed_vector;
151  }
152  }
153 
154 
157  const std::vector<int>& getActiveMap() const;
158  std::array<double, 3> getCellCenter(size_t i,size_t j, size_t k) const;
159  std::array<double, 3> getCellCenter(size_t globalIndex) const;
160  std::array<double, 3> getCornerPos(size_t i,size_t j, size_t k, size_t corner_index) const;
161  const std::vector<double>& activeVolume() const;
162  double getCellVolume(size_t globalIndex) const;
163  double getCellVolume(size_t i , size_t j , size_t k) const;
164  double getCellThickness(size_t globalIndex) const;
165  double getCellThickness(size_t i , size_t j , size_t k) const;
166  std::array<double, 3> getCellDims(size_t i,size_t j, size_t k) const;
167  std::array<double, 3> getCellDims(size_t globalIndex) const;
168  bool cellActive( size_t globalIndex ) const;
169  bool cellActive( size_t i , size_t j, size_t k ) const;
170 
171  std::array<double, 3> getCellDimensions(size_t i, size_t j, size_t k) const {
172  return getCellDims(i, j, k);
173  }
174 
175  bool isCellActive(size_t i, size_t j, size_t k) const {
176  return cellActive(i, j, k);
177  }
178 
184  bool isValidCellGeomtry(const std::size_t globalIndex,
185  const UnitSystem& usys) const;
186 
187  double getCellDepth(size_t i,size_t j, size_t k) const;
188  double getCellDepth(size_t globalIndex) const;
189  ZcornMapper zcornMapper() const;
190 
191  const std::vector<double>& getCOORD() const;
192  const std::vector<double>& getZCORN() const;
193  const std::vector<int>& getACTNUM( ) const;
194 
195  /*
196  The fixupZCORN method is run as part of constructiong the grid. This will adjust the
197  z-coordinates to ensure that cells do not overlap. The return value is the number of
198  points which have been adjusted. The number of zcorn nodes that has ben fixed is
199  stored in private member zcorn_fixed.
200  */
201 
202  size_t fixupZCORN();
203  size_t getZcornFixed() { return zcorn_fixed; };
204 
205  // resetACTNUM with no arguments will make all cells in the grid active.
206 
207  void resetACTNUM();
208  void resetACTNUM( const std::vector<int>& actnum);
209 
210  bool equal(const EclipseGrid& other) const;
211  static bool hasDVDEPTHZKeywords(const Deck&);
212 
213  /*
214  For ALugrid we can *only* use the keyword <DXV, DXYV, DZV, DEPTHZ> so to
215  initialize a Regular Cartesian Grid; further we need equidistant mesh
216  spacing in each direction to initialize ALuGrid (mandatory for
217  mesh refinement!).
218  */
219 
220  static bool hasEqualDVDEPTHZ(const Deck&);
221  static bool allEqual(const std::vector<double> &v);
222 
223  private:
224  std::vector<double> m_minpvVector;
225  MinpvMode::ModeEnum m_minpvMode;
226  std::optional<double> m_pinch;
227  PinchMode::ModeEnum m_pinchoutMode;
228  PinchMode::ModeEnum m_multzMode;
229  PinchMode::ModeEnum m_pinchGapMode;
230 
231  mutable std::optional<std::vector<double>> active_volume;
232 
233  bool m_circle = false;
234 
235  size_t zcorn_fixed = 0;
236  bool m_useActnumFromGdfile = false;
237 
238  // Input grid data.
239  std::vector<double> m_zcorn;
240  std::vector<double> m_coord;
241  std::vector<int> m_actnum;
242  std::optional<MapAxes> m_mapaxes;
243 
244  // Mapping to/from active cells.
245  int m_nactive;
246  std::vector<int> m_active_to_global;
247  std::vector<int> m_global_to_active;
248  // Numerical aquifer cells, needs to be active
249  std::unordered_set<size_t> m_aquifer_cells;
250 
251  // Radial grids need this for volume calculations.
252  std::optional<std::vector<double>> m_thetav;
253  std::optional<std::vector<double>> m_rv;
254 
255  void updateNumericalAquiferCells(const Deck&);
256 
257  void initGridFromEGridFile(Opm::EclIO::EclFile& egridfile, std::string fileName);
258  void resetACTNUM( const int* actnum);
259 
260  void initBinaryGrid(const Deck& deck);
261 
262  void initCornerPointGrid(const std::vector<double>& coord ,
263  const std::vector<double>& zcorn ,
264  const int * actnum);
265 
266  bool keywInputBeforeGdfile(const Deck& deck, const std::string& keyword) const;
267 
268  void initCylindricalGrid(const Deck&);
269  void initSpiderwebGrid(const Deck&);
270  void initSpiderwebOrCylindricalGrid(const Deck&, const bool);
271  void initCartesianGrid(const Deck&);
272  void initDTOPSGrid(const Deck&);
273  void initDVDEPTHZGrid(const Deck&);
274  void initGrid(const Deck&, const int* actnum);
275  void initCornerPointGrid(const Deck&);
276  void assertCornerPointKeywords(const Deck&);
277 
278  static bool hasDTOPSKeywords(const Deck&);
279  static void assertVectorSize(const std::vector<double>& vector, size_t expectedSize, const std::string& msg);
280 
281  static std::vector<double> createTOPSVector(const std::array<int, 3>& dims, const std::vector<double>& DZ, const Deck&);
282  static std::vector<double> createDVector(const std::array<int, 3>& dims, std::size_t dim, const std::string& DKey, const std::string& DVKey, const Deck&);
283  static void scatterDim(const std::array<int, 3>& dims , size_t dim , const std::vector<double>& DV , std::vector<double>& D);
284 
285 
286  std::vector<double> makeCoordDxDyDzTops(const std::vector<double>& dx, const std::vector<double>& dy, const std::vector<double>& dz, const std::vector<double>& tops) const;
287  std::vector<double> makeZcornDzTops(const std::vector<double>& dz, const std::vector<double>& tops) const ;
288  std::vector<double> makeZcornDzvDepthz(const std::vector<double>& dzv, const std::vector<double>& depthz) const;
289  std::vector<double> makeCoordDxvDyvDzvDepthz(const std::vector<double>& dxv, const std::vector<double>& dyv, const std::vector<double>& dzv, const std::vector<double>& depthz) const;
290 
291  void getCellCorners(const std::array<int, 3>& ijk, const std::array<int, 3>& dims, std::array<double,8>& X, std::array<double,8>& Y, std::array<double,8>& Z) const;
292  void getCellCorners(const std::size_t globalIndex,
293  std::array<double,8>& X,
294  std::array<double,8>& Y,
295  std::array<double,8>& Z) const;
296 
297  };
298 
299  class CoordMapper {
300  public:
301  CoordMapper(size_t nx, size_t ny);
302  size_t size() const;
303 
304 
305  /*
306  dim = 0,1,2 for x, y and z coordinate respectively.
307  layer = 0,1 for k=0 and k=nz layers respectively.
308  */
309 
310  size_t index(size_t i, size_t j, size_t dim, size_t layer) const;
311  private:
312  size_t nx;
313  size_t ny;
314  };
315 
316 
317 
318  class ZcornMapper {
319  public:
320  ZcornMapper(size_t nx, size_t ny, size_t nz);
321  size_t index(size_t i, size_t j, size_t k, int c) const;
322  size_t index(size_t g, int c) const;
323  size_t size() const;
324 
325  /*
326  The fixupZCORN method will take a zcorn vector as input and
327  run through it. If the following situation is detected:
328 
329  /| /|
330  / | / |
331  / | / |
332  / | / |
333  / | ==> / |
334  / | / |
335  ---/------x /---------x
336  | /
337  |/
338 
339  */
340  size_t fixupZCORN( std::vector<double>& zcorn);
341  bool validZCORN( const std::vector<double>& zcorn) const;
342  private:
343  std::array<size_t,3> dims;
344  std::array<size_t,3> stride;
345  std::array<size_t,8> cell_shift;
346  };
347 }
348 
349 #endif // OPM_PARSER_ECLIPSE_GRID_HPP
Definition: EclipseGrid.hpp:299
Definition: Deck.hpp:63
Definition: EclFile.hpp:35
About cell information and dimension: The actual grid information is held in a pointer to an ERT ecl_...
Definition: EclipseGrid.hpp:54
EclipseGrid(const Deck &deck, const int *actnum=nullptr)
EclipseGrid ignores ACTNUM in Deck, and therefore needs ACTNUM explicitly.
bool isValidCellGeomtry(const std::size_t globalIndex, const UnitSystem &usys) const
Whether or not given cell has a valid cell geometry.
const std::vector< int > & getActiveMap() const
Will return a vector a length num_active; where the value of each element is the corresponding global...
Definition: GridDims.hpp:31
Definition: UnitSystem.hpp:33
Definition: EclipseGrid.hpp:318
This class implements a small container which holds the transmissibility mulitpliers for all the face...
Definition: Exceptions.hpp:29