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