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