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
35namespace 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
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
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.
const std::vector< int > & getActiveMap() const
Will return a vector a length num_active; where the value of each element is the corresponding global...
bool isValidCellGeomtry(const std::size_t globalIndex, const UnitSystem &usys) const
Whether or not given cell has a valid cell geometry.
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