My Project
Loading...
Searching...
No Matches
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#include <map>
35
36namespace 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:
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, const std::vector<int>& actnum);
65 EclipseGrid(const EclipseGrid& src, const double* zcorn, const std::vector<int>& actnum);
66
67 EclipseGrid(size_t nx, size_t ny, size_t nz,
68 double dx = 1.0, double dy = 1.0, double dz = 1.0,
69 double top = 0.0);
70 explicit EclipseGrid(const GridDims& gd);
71
72 EclipseGrid(const std::array<int, 3>& dims ,
73 const std::vector<double>& coord ,
74 const std::vector<double>& zcorn ,
75 const int * actnum = nullptr);
76
77
80 explicit EclipseGrid(const Deck& deck, const int * actnum = nullptr);
81
82 static bool hasGDFILE(const Deck& deck);
83 static bool hasRadialKeywords(const Deck& deck);
84 static bool hasSpiderKeywords(const Deck& deck);
85 static bool hasCylindricalKeywords(const Deck& deck);
86 static bool hasCornerPointKeywords(const Deck&);
87 static bool hasCartesianKeywords(const Deck&);
88 size_t getNumActive( ) const;
89 bool allActive() const;
90
91 size_t activeIndex(size_t i, size_t j, size_t k) const;
92 size_t activeIndex(size_t globalIndex) const;
93
94 size_t getActiveIndex(size_t i, size_t j, size_t k) const {
95 return activeIndex(i, j, k);
96 }
97
98 void save(const std::string& filename, bool formatted, const std::vector<Opm::NNCdata>& nnc, const Opm::UnitSystem& units) const;
99 /*
100 Observe that the there is a getGlobalIndex(i,j,k)
101 implementation in the base class. This method - translating
102 from an active index to a global index must be implemented
103 in the current class.
104 */
105 using GridDims::getGlobalIndex;
106 size_t getGlobalIndex(size_t active_index) const;
107
108 /*
109 For RADIAL grids you can *optionally* use the keyword
110 'CIRCLE' to denote that period boundary conditions should be
111 applied in the 'THETA' direction; this will only apply if
112 the theta keywords entered sum up to exactly 360 degrees!
113 */
114
115 bool circle() const;
116 bool isPinchActive() const;
117 double getPinchThresholdThickness() const;
118 PinchMode getPinchOption() const;
119 PinchMode getMultzOption() const;
120 PinchMode getPinchGapMode() const;
121 double getPinchMaxEmptyGap() const;
122
123 MinpvMode getMinpvMode() const;
124 const std::vector<double>& getMinpvVector( ) const;
125
126 /*
127 Will return a vector of nactive elements. The method will
128 behave differently depending on the lenght of the
129 input_vector:
130
131 nx*ny*nz: only the values corresponding to active cells
132 are copied out.
133
134 nactive: The input vector is copied straight out again.
135
136 ??? : Exception.
137 */
138
139 template<typename T>
140 std::vector<T> compressedVector(const std::vector<T>& input_vector) const {
141 if( input_vector.size() == this->getNumActive() ) {
142 return input_vector;
143 }
144
145 if (input_vector.size() != getCartesianSize())
146 throw std::invalid_argument("Input vector must have full size");
147
148 {
149 std::vector<T> compressed_vector( this->getNumActive() );
150 const auto& active_map = this->getActiveMap( );
151
152 for (size_t i = 0; i < this->getNumActive(); ++i)
153 compressed_vector[i] = input_vector[ active_map[i] ];
154
155 return compressed_vector;
156 }
157 }
158
159
162 const std::vector<int>& getActiveMap() const;
164 std::tuple<std::array<double, 3>,std::array<double, 3>,std::array<double, 3>>
165 getCellAndBottomCenterNormal(size_t globalIndex) const;
166 std::array<double, 3> getCellCenter(size_t i,size_t j, size_t k) const;
167 std::array<double, 3> getCellCenter(size_t globalIndex) const;
168 std::array<double, 3> getCornerPos(size_t i,size_t j, size_t k, size_t corner_index) const;
169 const std::vector<double>& activeVolume() const;
170 double getCellVolume(size_t globalIndex) const;
171 double getCellVolume(size_t i , size_t j , size_t k) const;
172 double getCellThickness(size_t globalIndex) const;
173 double getCellThickness(size_t i , size_t j , size_t k) const;
174 std::array<double, 3> getCellDims(size_t i,size_t j, size_t k) const;
175 std::array<double, 3> getCellDims(size_t globalIndex) const;
176 bool cellActive( size_t globalIndex ) const;
177 bool cellActive( size_t i , size_t j, size_t k ) const;
178 bool cellActiveAfterMINPV( size_t i , size_t j , size_t k, double cell_porv ) const;
179
180 std::array<double, 3> getCellDimensions(size_t i, size_t j, size_t k) const {
181 return getCellDims(i, j, k);
182 }
183
184 bool isCellActive(size_t i, size_t j, size_t k) const {
185 return cellActive(i, j, k);
186 }
187
193 bool isValidCellGeomtry(const std::size_t globalIndex,
194 const UnitSystem& usys) const;
195
196 double getCellDepth(size_t i,size_t j, size_t k) const;
197 double getCellDepth(size_t globalIndex) const;
198 ZcornMapper zcornMapper() const;
199
200 const std::vector<double>& getCOORD() const;
201 const std::vector<double>& getZCORN() const;
202 const std::vector<int>& getACTNUM( ) const;
203
204
205 const std::map<size_t, std::array<int,2>>& getAquiferCellTabnums() const;
206
207 /*
208 The fixupZCORN method is run as part of constructiong the grid. This will adjust the
209 z-coordinates to ensure that cells do not overlap. The return value is the number of
210 points which have been adjusted. The number of zcorn nodes that has ben fixed is
211 stored in private member zcorn_fixed.
212 */
213
214 size_t fixupZCORN();
215 size_t getZcornFixed() { return zcorn_fixed; };
216
217 // resetACTNUM with no arguments will make all cells in the grid active.
218
219 void resetACTNUM();
220 void resetACTNUM( const std::vector<int>& actnum);
221
223 void setMINPVV(const std::vector<double>& minpvv);
224
225 bool equal(const EclipseGrid& other) const;
226 static bool hasDVDEPTHZKeywords(const Deck&);
227
228 /*
229 For ALugrid we can *only* use the keyword <DXV, DXYV, DZV, DEPTHZ> so to
230 initialize a Regular Cartesian Grid; further we need equidistant mesh
231 spacing in each direction to initialize ALuGrid (mandatory for
232 mesh refinement!).
233 */
234
235 static bool hasEqualDVDEPTHZ(const Deck&);
236 static bool allEqual(const std::vector<double> &v);
237
238 private:
239 std::vector<double> m_minpvVector;
240 MinpvMode m_minpvMode;
241 std::optional<double> m_pinch;
242 // Option 4 of PINCH (TOPBOT/ALL), how to calculate TRANS
243 PinchMode m_pinchoutMode;
244 // Option 5 of PINCH (TOP/ALL), how to apply MULTZ
245 PinchMode m_multzMode;
246 // Option 2 of PINCH (GAP/NOGAP)
247 PinchMode m_pinchGapMode;
248 double m_pinchMaxEmptyGap;
249
250 mutable std::optional<std::vector<double>> active_volume;
251
252 bool m_circle = false;
253
254 size_t zcorn_fixed = 0;
255 bool m_useActnumFromGdfile = false;
256
257 // Input grid data.
258 mutable std::optional<std::vector<double>> m_input_zcorn;
259 mutable std::optional<std::vector<double>> m_input_coord;
260
261 std::vector<double> m_zcorn;
262 std::vector<double> m_coord;
263
264
265 std::vector<int> m_actnum;
266 std::optional<MapAxes> m_mapaxes;
267
268 // Mapping to/from active cells.
269 int m_nactive {};
270 std::vector<int> m_active_to_global;
271 std::vector<int> m_global_to_active;
272 // Numerical aquifer cells, needs to be active
273 std::unordered_set<size_t> m_aquifer_cells;
274 // Keep track of aquifer cell depths and (pvtnum,satnum)
275 std::map<size_t, double> m_aquifer_cell_depths;
276 std::map<size_t, std::array<int,2>> m_aquifer_cell_tabnums;
277
278 // Radial grids need this for volume calculations.
279 std::optional<std::vector<double>> m_thetav;
280 std::optional<std::vector<double>> m_rv;
281
282 void updateNumericalAquiferCells(const Deck&);
283 double computeCellGeometricDepth(size_t globalIndex) const;
284
285 void initGridFromEGridFile(Opm::EclIO::EclFile& egridfile,
286 const std::string& fileName);
287 void resetACTNUM( const int* actnum);
288
289 void initBinaryGrid(const Deck& deck);
290
291 void initCornerPointGrid(const std::vector<double>& coord ,
292 const std::vector<double>& zcorn ,
293 const int * actnum);
294
295 bool keywInputBeforeGdfile(const Deck& deck, const std::string& keyword) const;
296
297 void initCylindricalGrid(const Deck&);
298 void initSpiderwebGrid(const Deck&);
299 void initSpiderwebOrCylindricalGrid(const Deck&, const bool);
300 void initCartesianGrid(const Deck&);
301 void initDTOPSGrid(const Deck&);
302 void initDVDEPTHZGrid(const Deck&);
303 void initGrid(const Deck&, const int* actnum);
304 void initCornerPointGrid(const Deck&);
305 void assertCornerPointKeywords(const Deck&);
306
307 static bool hasDTOPSKeywords(const Deck&);
308 static void assertVectorSize(const std::vector<double>& vector, size_t expectedSize, const std::string& msg);
309
310 static std::vector<double> createTOPSVector(const std::array<int, 3>& dims, const std::vector<double>& DZ, const Deck&);
311 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&);
312 static void scatterDim(const std::array<int, 3>& dims , size_t dim , const std::vector<double>& DV , std::vector<double>& D);
313
314
315 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;
316 std::vector<double> makeZcornDzTops(const std::vector<double>& dz, const std::vector<double>& tops) const ;
317 std::vector<double> makeZcornDzvDepthz(const std::vector<double>& dzv, const std::vector<double>& depthz) const;
318 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;
319
320 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;
321 void getCellCorners(const std::size_t globalIndex,
322 std::array<double,8>& X,
323 std::array<double,8>& Y,
324 std::array<double,8>& Z) const;
325
326 };
327
329 public:
330 CoordMapper(size_t nx, size_t ny);
331 size_t size() const;
332
333
334 /*
335 dim = 0,1,2 for x, y and z coordinate respectively.
336 layer = 0,1 for k=0 and k=nz layers respectively.
337 */
338
339 size_t index(size_t i, size_t j, size_t dim, size_t layer) const;
340 private:
341 size_t nx;
342 size_t ny;
343 };
344
345
346
348 public:
349 ZcornMapper(size_t nx, size_t ny, size_t nz);
350 size_t index(size_t i, size_t j, size_t k, int c) const;
351 size_t index(size_t g, int c) const;
352 size_t size() const;
353
354 /*
355 The fixupZCORN method will take a zcorn vector as input and
356 run through it. If the following situation is detected:
357
358 /| /|
359 / | / |
360 / | / |
361 / | / |
362 / | ==> / |
363 / | / |
364 ---/------x /---------x
365 | /
366 |/
367
368 */
369 size_t fixupZCORN( std::vector<double>& zcorn);
370 bool validZCORN( const std::vector<double>& zcorn) const;
371 private:
372 std::array<size_t,3> dims;
373 std::array<size_t,3> stride;
374 std::array<size_t,8> cell_shift;
375 };
376}
377
378#endif // OPM_PARSER_ECLIPSE_GRID_HPP
Definition EclipseGrid.hpp:328
Definition Deck.hpp:49
Definition EclFile.hpp:36
About cell information and dimension: The actual grid information is held in a pointer to an ERT ecl_...
Definition EclipseGrid.hpp:55
void setMINPVV(const std::vector< double > &minpvv)
Sets MINPVV if MINPV and MINPORV are not used.
Definition EclipseGrid.cpp:2006
size_t getGlobalIndex(size_t active_index) const
Observe: the input argument is assumed to be in the space [0,num_active).
Definition EclipseGrid.cpp:576
std::tuple< std::array< double, 3 >, std::array< double, 3 >, std::array< double, 3 > > getCellAndBottomCenterNormal(size_t globalIndex) const
get cell center, and center and normal of bottom face
Definition EclipseGrid.cpp:1661
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 EclipseGrid.cpp:1961
bool isValidCellGeomtry(const std::size_t globalIndex, const UnitSystem &usys) const
Whether or not given cell has a valid cell geometry.
Definition EclipseGrid.cpp:1770
Definition GridDims.hpp:31
Definition UnitSystem.hpp:34
Definition EclipseGrid.hpp:347
This class implements a small container which holds the transmissibility mulitpliers for all the face...
Definition Exceptions.hpp:30