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