opm-common
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 <algorithm>
29#include <array>
30#include <memory>
31#include <optional>
32#include <stdexcept>
33#include <string>
34#include <unordered_set>
35#include <vector>
36#include <map>
37
38namespace Opm {
39
40 class Deck;
41 namespace EclIO {
42 class EclFile;
43 class EclOutput;
44 }
45 struct NNCdata;
46 class UnitSystem;
47 class ZcornMapper;
48 class EclipseGridLGR;
49 class LgrCollection;
60
61 class EclipseGrid : public GridDims {
62 public:
63 EclipseGrid();
64 explicit EclipseGrid(const std::string& filename);
65
66 /*
67 These constructors will make a copy of the src grid, with
68 zcorn and or actnum have been adjustments.
69 */
70 EclipseGrid(const EclipseGrid& src, const std::vector<int>& actnum);
71 EclipseGrid(const EclipseGrid& src, const double* zcorn, const std::vector<int>& actnum);
72
73 EclipseGrid(size_t nx, size_t ny, size_t nz,
74 double dx = 1.0, double dy = 1.0, double dz = 1.0,
75 double top = 0.0);
76 explicit EclipseGrid(const GridDims& gd);
77
78 EclipseGrid(const std::array<int, 3>& dims ,
79 const std::vector<double>& coord ,
80 const std::vector<double>& zcorn ,
81 const int * actnum = nullptr);
82
83 virtual ~EclipseGrid() = default;
84
87 explicit EclipseGrid(const Deck& deck, const int * actnum = nullptr);
88
89 static bool hasGDFILE(const Deck& deck);
90 static bool hasRadialKeywords(const Deck& deck);
91 static bool hasSpiderKeywords(const Deck& deck);
92 static bool hasCylindricalKeywords(const Deck& deck);
93 static bool hasCornerPointKeywords(const Deck&);
94 static bool hasCartesianKeywords(const Deck&);
95 size_t getNumActive( ) const;
96 bool allActive() const;
97
98 size_t activeIndex(size_t i, size_t j, size_t k) const;
99 size_t activeIndex(size_t globalIndex) const;
100
101 size_t getTotalActiveLGR() const;
102 size_t getActiveIndexLGR(const std::string& label, size_t i, size_t j, size_t k) const;
103 size_t getActiveIndexLGR(const std::string& label, size_t localIndex) const;
104
105 size_t activeIndexLGR(const std::string& label, size_t i, size_t j, size_t k) const;
106 size_t activeIndexLGR(const std::string& label, size_t localIndex) const;
107
108 size_t getActiveIndex(size_t i, size_t j, size_t k) const {
109 return activeIndex(i, j, k);
110 }
111 const std::vector<std::size_t>& get_print_order_lgr () const {
112 return m_print_order_lgr_cells;
113 }
114
115 size_t get_lgr_cell_index(const std::string& lgr_tag) const {
116 const auto& labels = get_all_lgr_labels();
117
118 if (labels.empty()) {
119 throw std::runtime_error("No LGR labels available.");
120 }
121
122 auto it = std::find(labels.begin(), labels.end(), lgr_tag);
123 if (it == labels.end()) {
124 throw std::runtime_error("LGR tag not found: " + lgr_tag);
125 }
126
127 return static_cast<size_t>(std::distance(labels.begin(), it));
128 }
129
130 size_t getActiveIndex(size_t globalIndex) const {
131 return activeIndex(globalIndex);
132 }
133
134 std::vector<std::string> get_all_lgr_labels() const
135 {
136 if (this->all_lgr_labels.empty()) {
137 return {};
138 } else {
139 return {this->all_lgr_labels.begin() + 1, this->all_lgr_labels.end()};
140 }
141 }
142
143 const std::string& get_lgr_labels_by_number(std::size_t num) const
144 {
145 return all_lgr_labels[num];
146 }
147
148 const std::vector<std::string>& get_all_labels() const
149 {
150 return this->all_lgr_labels;
151 }
152
153 const std::string& get_lgr_tag() const {
154 return this->lgr_label;
155 }
156
157 std::vector<GridDims> get_lgr_children_gridim() const;
158
159
160 void assertIndexLGR(size_t localIndex) const;
161
162 void assertLabelLGR(const std::string& label) const;
163
164 void save(const std::string& filename, bool formatted, const std::vector<Opm::NNCdata>& nnc, const Opm::UnitSystem& units) const;
165 /*
166 Observe that the there is a getGlobalIndex(i,j,k)
167 implementation in the base class. This method - translating
168 from an active index to a global index must be implemented
169 in the current class.
170 */
171 void init_children_host_cells(bool logical = true);
172 void init_children_host_cells_logical(void);
173 void init_children_host_cells_geometrical(void);
174 std::array<int,3> getCellSubdivisionRatioLGR(const std::string& lgr_tag,
175 std::array<int,3> acum = {1,1,1}) const;
176 using GridDims::getGlobalIndex;
177 size_t getGlobalIndex(size_t active_index) const;
178
179 /*
180 For RADIAL grids you can *optionally* use the keyword
181 'CIRCLE' to denote that period boundary conditions should be
182 applied in the 'THETA' direction; this will only apply if
183 the theta keywords entered sum up to exactly 360 degrees!
184 */
185 bool circle() const;
186 bool isPinchActive() const;
187 double getPinchThresholdThickness() const;
188 PinchMode getPinchOption() const;
189 PinchMode getMultzOption() const;
190 PinchMode getPinchGapMode() const;
191 double getPinchMaxEmptyGap() const;
192 MinpvMode getMinpvMode() const;
193 const std::vector<double>& getMinpvVector( ) const;
194
195 /*
196 Will return a vector of nactive elements. The method will
197 behave differently depending on the lenght of the
198 input_vector:
199
200 nx*ny*nz: only the values corresponding to active cells
201 are copied out.
202
203 nactive: The input vector is copied straight out again.
204
205 ??? : Exception.
206 */
207
208 template<typename T>
209 std::vector<T> compressedVector(const std::vector<T>& input_vector) const {
210 if( input_vector.size() == this->getNumActive() ) {
211 return input_vector;
212 }
213
214 if (input_vector.size() != getCartesianSize())
215 throw std::invalid_argument("Input vector must have full size");
216
217 {
218 std::vector<T> compressed_vector( this->getNumActive() );
219 const auto& active_map = this->getActiveMap( );
220
221 for (size_t i = 0; i < this->getNumActive(); ++i)
222 compressed_vector[i] = input_vector[ active_map[i] ];
223
224 return compressed_vector;
225 }
226 }
227
228
231 const std::vector<int>& getActiveMap() const;
232
233 void init_lgr_cells(const LgrCollection& lgr_input);
234 void create_lgr_cells_tree(const LgrCollection& );
236 std::tuple<std::array<double, 3>,std::array<double, 3>,std::array<double, 3>>
237 getCellAndBottomCenterNormal(size_t globalIndex) const;
238 std::array<double, 3> getCellCenter(size_t i,size_t j, size_t k) const;
239 std::array<double, 3> getCellCenter(size_t globalIndex) const;
240 std::array<double, 3> getCornerPos(size_t i,size_t j, size_t k, size_t corner_index) const;
241 const std::vector<double>& activeVolume() const;
242 double getCellVolume(size_t globalIndex) const;
243 double getCellVolume(size_t i , size_t j , size_t k) const;
244 double getCellThickness(size_t globalIndex) const;
245 double getCellThickness(size_t i , size_t j , size_t k) const;
246 std::array<double, 3> getCellDims(size_t i,size_t j, size_t k) const;
247 std::array<double, 3> getCellDims(size_t globalIndex) const;
248 bool cellActive( size_t globalIndex ) const;
249 bool cellActive( size_t i , size_t j, size_t k ) const;
250 bool cellActiveAfterMINPV( size_t i , size_t j , size_t k, double cell_porv ) const;
251 bool is_lgr() const {return lgr_grid;};
252 std::array<double, 3> getCellDimensions(size_t i, size_t j, size_t k) const {
253 return getCellDims(i, j, k);
254 }
255 std::array<double,3> getCellDimensionsLGR(const std::size_t i,
256 const std::size_t j,
257 const std::size_t k,
258 const std::string& lgr_tag) const;
259 double getCellDepthLGR(size_t i, size_t j, size_t k, const std::string& lgr_tag) const;
260
261
262 bool isCellActive(size_t i, size_t j, size_t k) const {
263 return cellActive(i, j, k);
264 }
265
271 bool isValidCellGeomtry(const std::size_t globalIndex,
272 const UnitSystem& usys) const;
273
274 double getCellDepth(size_t i,size_t j, size_t k) const;
275 double getCellDepth(size_t globalIndex) const;
276 ZcornMapper zcornMapper() const;
277
278 const std::vector<double>& getCOORD() const;
279 const std::vector<double>& getZCORN() const;
280 const std::vector<int>& getACTNUM( ) const;
281
282 const std::optional<MapAxes>& getMapAxes() const;
283
284 const std::map<size_t, std::array<int,2>>& getAquiferCellTabnums() const;
285
286 /*
287 The fixupZCORN method is run as part of constructiong the grid. This will adjust the
288 z-coordinates to ensure that cells do not overlap. The return value is the number of
289 points which have been adjusted. The number of zcorn nodes that has ben fixed is
290 stored in private member zcorn_fixed.
291 */
292
293 size_t fixupZCORN();
294 size_t getZcornFixed() { return zcorn_fixed; };
295
296 // resetACTNUM with no arguments will make all cells in the grid active.
297
298 void resetACTNUM();
299 void resetACTNUM( const std::vector<int>& actnum);
300
302 void setMINPVV(const std::vector<double>& minpvv);
303
304 bool equal(const EclipseGrid& other) const;
305 static bool hasDVDEPTHZKeywords(const Deck&);
306
307 /*
308 For ALugrid we can *only* use the keyword <DXV, DXYV, DZV, DEPTHZ> so to
309 initialize a Regular Cartesian Grid; further we need equidistant mesh
310 spacing in each direction to initialize ALuGrid (mandatory for
311 mesh refinement!).
312 */
313
314 static bool hasEqualDVDEPTHZ(const Deck&);
315 static bool allEqual(const std::vector<double> &v);
316 EclipseGridLGR& getLGRCell(std::size_t index);
317 const EclipseGridLGR& getLGRCell(std::size_t index) const;
318 const EclipseGridLGR& getLGRCell(const std::string& lgr_tag) const;
319 int getLGR_global_father(std::size_t global_index, const std::string& lgr_tag) const;
320 int getLGR_father(std::size_t i, std::size_t j, std::size_t k, const std::string& lgr_tag) const;
321 int getLGR_father(std::size_t global_index, const std::string& lgr_tag) const;
322 std::array<int,3> getLGR_fatherIJK(std::size_t i, std::size_t j, std::size_t k, const std::string& lgr_tag) const;
323 std::vector<EclipseGridLGR> lgr_children_cells;
331 virtual void set_lgr_refinement(const std::string& lgr_tag, const std::vector<double>& coords, const std::vector<double>& zcorn);
332
333 protected:
334 std::size_t lgr_global_counter = 0;
335 std::string lgr_label = "GLOBAL";
336 int lgr_level = 0;
337 int lgr_level_father = 0;
338 std::vector<std::string> lgr_children_labels;
339 std::vector<std::size_t> lgr_active_index;
340 std::vector<std::size_t> lgr_level_active_map;
341 std::vector<std::string> all_lgr_labels;
342 std::map<std::vector<std::size_t>, std::size_t> num_lgr_children_cells;
343 std::vector<double> m_zcorn;
344 std::vector<double> m_coord;
345 std::vector<int> m_actnum;
346 std::vector<std::size_t> m_print_order_lgr_cells;
347 // Input grid data.
348 mutable std::optional<std::vector<double>> m_input_zcorn;
349 mutable std::optional<std::vector<double>> m_input_coord;
350
351 private:
352 std::vector<double> m_minpvVector;
353 MinpvMode m_minpvMode;
354 std::optional<double> m_pinch;
355
356 // Option 4 of PINCH (TOPBOT/ALL), how to calculate TRANS
357 PinchMode m_pinchoutMode;
358 // Option 5 of PINCH (TOP/ALL), how to apply MULTZ
359 PinchMode m_multzMode;
360 // Option 2 of PINCH (GAP/NOGAP)
361 PinchMode m_pinchGapMode;
362 double m_pinchMaxEmptyGap;
363 bool lgr_grid = false;
364 mutable std::optional<std::vector<double>> active_volume;
365
366 bool m_circle = false;
367 size_t zcorn_fixed = 0;
368 bool m_useActnumFromGdfile = false;
369
370 std::optional<MapAxes> m_mapaxes;
371
372 // Mapping to/from active cells.
373 int m_nactive {};
374 std::vector<int> m_active_to_global;
375 std::vector<int> m_global_to_active;
376 // Numerical aquifer cells, needs to be active
377 std::unordered_set<size_t> m_aquifer_cells;
378 // Keep track of aquifer cell depths and (pvtnum,satnum)
379 std::map<size_t, double> m_aquifer_cell_depths;
380 std::map<size_t, std::array<int,2>> m_aquifer_cell_tabnums;
381
382 // Radial grids need this for volume calculations.
383 std::optional<std::vector<double>> m_thetav;
384 std::optional<std::vector<double>> m_rv;
385 void parseGlobalReferenceToChildren(void);
386 int initializeLGRObjectIndices(int);
387 void initializeLGRTreeIndices(void);
388 void propagateParentIndicesToLGRChildren(int);
389 void updateNumericalAquiferCells(const Deck&);
390 double computeCellGeometricDepth(size_t globalIndex) const;
391
392 void initGridFromEGridFile(Opm::EclIO::EclFile& egridfile,
393 const std::string& fileName);
394 void resetACTNUM( const int* actnum);
395
396 void initBinaryGrid(const Deck& deck);
397
398 void initCornerPointGrid(const std::vector<double>& coord ,
399 const std::vector<double>& zcorn ,
400 const int * actnum);
401
402 bool keywInputBeforeGdfile(const Deck& deck, const std::string& keyword) const;
403
404 void initCylindricalGrid(const Deck&);
405 void initSpiderwebGrid(const Deck&);
406 void initSpiderwebOrCylindricalGrid(const Deck&, const bool);
407 void initCartesianGrid(const Deck&);
408 void initDTOPSGrid(const Deck&);
409 void initDVDEPTHZGrid(const Deck&);
410 void initGrid(const Deck&, const int* actnum);
411 void initCornerPointGrid(const Deck&);
412 void assertCornerPointKeywords(const Deck&);
413 void save_all_lgr_labels(const LgrCollection& );
414 static bool hasDTOPSKeywords(const Deck&);
415 static void assertVectorSize(const std::vector<double>& vector, size_t expectedSize, const std::string& msg);
416
417 static std::vector<double> createTOPSVector(const std::array<int, 3>& dims, const std::vector<double>& DZ, const Deck&);
418 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&);
419 static void scatterDim(const std::array<int, 3>& dims , size_t dim , const std::vector<double>& DV , std::vector<double>& D);
420
421
422 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;
423 std::vector<double> makeZcornDzTops(const std::vector<double>& dz, const std::vector<double>& tops) const ;
424 std::vector<double> makeZcornDzvDepthz(const std::vector<double>& dzv, const std::vector<double>& depthz) const;
425 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;
426
427 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;
428 void getCellCorners(const std::size_t globalIndex,
429 std::array<double,8>& X,
430 std::array<double,8>& Y,
431 std::array<double,8>& Z) const;
432
433 };
434
436 class EclipseGridLGR : public EclipseGrid
437 {
438 public:
439 using vec_size_t = std::vector<std::size_t>;
440 EclipseGridLGR() = default;
441 EclipseGridLGR(const std::string& self_label,
442 const std::string& father_label_,
443 size_t nx,
444 size_t ny,
445 size_t nz,
446 const vec_size_t& father_lgr_index,
447 const std::array<int, 3>& low_fatherIJK_,
448 const std::array<int, 3>& up_fatherIJK_);
449 const vec_size_t& getFatherGlobalID() const;
450 void save(Opm::EclIO::EclOutput&, const std::vector<Opm::NNCdata>&, const Opm::UnitSystem&) const;
451 void save_nnc(Opm::EclIO::EclOutput&) const;
452 void set_lgr_global_counter(std::size_t counter)
453 {
454 lgr_global_counter = counter;
455 }
456 const vec_size_t& get_father_global() const
457 {
458 return father_global;
459 }
460 std::optional<std::reference_wrapper<EclipseGridLGR>>
461 get_child_LGR_cell(const std::string& lgr_tag) const;
462 std::vector<int> save_hostnum(void) const;
463 int get_hostnum(std::size_t global_index) const {return(m_hostnum[global_index]);};
464
465 //parsing the father grid allows the global_father references to be given in terms of father_grid
466 std::vector<int> getLGRCell_global_father(const EclipseGrid& father_grid) const;
467 std::vector<double> getLGRCell_all_depth (const EclipseGrid& father_grid) const;
468
469 void get_label_child_to_top_father(std::vector<std::reference_wrapper<const std::string>>& list) const;
470 void get_global_index_child_to_top_father(std::vector<std::size_t> & list,
471 std::size_t i, std::size_t j, std::size_t k) const;
472 void get_global_index_child_to_top_father(std::vector<std::size_t> & list, std::size_t global_ind) const;
473
474 void set_hostnum(const std::vector<int>&);
475 const std::array<int,3>& get_low_fatherIJK() const{
476 return low_fatherIJK;
477 }
478 const std::string& get_father_label() const{
479 return father_label;
480 }
481
482 const std::array<int,3>& get_up_fatherIJK() const{
483 return up_fatherIJK;
484 }
485
486
494 void set_lgr_refinement(const std::string& lgr_tag,
495 const std::vector<double>& coord,
496 const std::vector<double>& zcorn) override;
497
498 void set_lgr_refinement(const std::vector<double>&, const std::vector<double>&);
499
500 private:
501 void init_father_global();
502 std::string father_label;
503 // references global on the father label
504 vec_size_t father_global;
505 std::array<int, 3> low_fatherIJK {};
506 std::array<int, 3> up_fatherIJK {};
507 std::vector<int> m_hostnum;
508 };
509
510
511 class CoordMapper {
512 public:
513 CoordMapper(size_t nx, size_t ny);
514 size_t size() const;
515
516
517 /*
518 dim = 0,1,2 for x, y and z coordinate respectively.
519 layer = 0,1 for k=0 and k=nz layers respectively.
520 */
521
522 size_t index(size_t i, size_t j, size_t dim, size_t layer) const;
523 private:
524 size_t nx;
525 size_t ny;
526 };
527
528
529
530 class ZcornMapper {
531 public:
532 ZcornMapper(size_t nx, size_t ny, size_t nz);
533 size_t index(size_t i, size_t j, size_t k, int c) const;
534 size_t index(size_t g, int c) const;
535 size_t size() const;
536
537 /*
538 The fixupZCORN method will take a zcorn vector as input and
539 run through it. If the following situation is detected:
540
541 /| /|
542 / | / |
543 / | / |
544 / | / |
545 / | ==> / |
546 / | / |
547 ---/------x /---------x
548 | /
549 |/
550
551 */
552 size_t fixupZCORN( std::vector<double>& zcorn);
553 bool validZCORN( const std::vector<double>& zcorn) const;
554 private:
555 std::array<size_t,3> dims;
556 std::array<size_t,3> stride;
557 std::array<size_t,8> cell_shift;
558 };
559}
560#endif // OPM_PARSER_ECLIPSE_GRID_HPP
Opm::Deck
Definition Deck.hpp:46
Opm::EclIO::EclFile
Definition EclFile.hpp:36
Opm::EclIO::EclOutput
Definition EclOutput.hpp:38
Opm::EclipseGridLGR
Specialized Class to describe LGR refined cells.
Definition EclipseGrid.hpp:437
Opm::EclipseGridLGR::set_lgr_refinement
void set_lgr_refinement(const std::string &lgr_tag, const std::vector< double > &coord, const std::vector< double > &zcorn) override
Sets Local Grid Refinement for the EclipseGridLGR.
Definition EclipseGrid.cpp:2768
Opm::EclipseGrid::setMINPVV
void setMINPVV(const std::vector< double > &minpvv)
Sets MINPVV if MINPV and MINPORV are not used.
Definition EclipseGrid.cpp:2487
Opm::EclipseGrid::getCellDimensionsLGR
std::array< double, 3 > getCellDimensionsLGR(const std::size_t i, const std::size_t j, const std::size_t k, const std::string &lgr_tag) const
Computes the dimensions of a local grid refinement (LGR) cell.
Definition EclipseGrid.cpp:2059
Opm::EclipseGrid::getCellAndBottomCenterNormal
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:1673
Opm::EclipseGrid::getActiveMap
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:2226
Opm::EclipseGrid::getGlobalIndex
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:588
Opm::EclipseGrid::set_lgr_refinement
virtual void set_lgr_refinement(const std::string &lgr_tag, const std::vector< double > &coords, const std::vector< double > &zcorn)
Sets Local Grid Refinement for the EclipseGrid.
Definition EclipseGrid.cpp:2149
Opm::EclipseGrid::isValidCellGeomtry
bool isValidCellGeomtry(const std::size_t globalIndex, const UnitSystem &usys) const
Whether or not given cell has a valid cell geometry.
Definition EclipseGrid.cpp:1782
Opm::LgrCollection
Definition LgrCollection.hpp:33
Opm::UnitSystem
Definition UnitSystem.hpp:34
Opm::ZcornMapper
Definition EclipseGrid.hpp:530
Opm
This class implements a small container which holds the transmissibility mulitpliers for all the face...
Definition Exceptions.hpp:30
Opm::NNCdata
Definition NNC.hpp:35