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