From 34410cbe417786f30256a475c5c4c33a0b5fef09 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Torbj=C3=B8rn=20Skille?= Date: Sat, 24 Aug 2019 14:59:52 +0200 Subject: [PATCH] ERT/libecl routines used in the EclipseGrid class has been replace by new member functions. Grid properties are calculated and stored as private data members in EclipseGrid. The API for the class is unchanged except for some minor changes for exportACTNUM, exportZCORN and exportCOORD. These changes have triggered some very few modifications in the opm-grid and opm-simulators repo. --- CMakeLists_files.cmake | 2 +- .../utility/numeric/calculateCellVol.hpp | 3 +- opm/output/eclipse/Summary.hpp | 1 + .../eclipse/EclipseState/Grid/EclipseGrid.hpp | 112 +- .../utility/numeric/calculateCellVol.cpp | 8 +- src/opm/io/eclipse/EclFile.cpp | 22 +- src/opm/output/eclipse/EclipseIO.cpp | 3 +- .../eclipse/EclipseState/EclipseState.cpp | 2 +- .../eclipse/EclipseState/Grid/EclipseGrid.cpp | 1367 +++++++++++--- tests/parser/ConnectionTests.cpp | 2 +- tests/parser/EclipseGridTests.cpp | 1638 +++++++++++++++-- tests/parser/ScheduleTests.cpp | 1 - .../integration/EclipseGridCreateFromDeck.cpp | 11 +- tests/test_AggregateConnectionData.cpp | 2 +- tests/test_EclIO.cpp | 4 +- tests/test_EclipseIO.cpp | 6 +- tests/test_calculateCellVol.cpp | 24 +- tests/test_restartwellinfo.cpp | 234 +++ tests/test_writenumwells.cpp | 173 -- tests/testblackoilstate3.DATA | 27 +- 20 files changed, 2953 insertions(+), 689 deletions(-) create mode 100644 tests/test_restartwellinfo.cpp delete mode 100644 tests/test_writenumwells.cpp diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 5f85c0aae..fb10c1716 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -316,7 +316,7 @@ if(ENABLE_ECL_OUTPUT) tests/test_Tables.cpp tests/test_Wells.cpp tests/test_WindowedArray.cpp - tests/test_writenumwells.cpp + tests/test_restartwellinfo.cpp ) endif() diff --git a/opm/common/utility/numeric/calculateCellVol.hpp b/opm/common/utility/numeric/calculateCellVol.hpp index aeb801a78..bc2533216 100755 --- a/opm/common/utility/numeric/calculateCellVol.hpp +++ b/opm/common/utility/numeric/calculateCellVol.hpp @@ -16,9 +16,10 @@ */ #include +#include #include -double calculateCellVol(const std::vector& X, const std::vector& Y, const std::vector& Z); +double calculateCellVol(const std::array& X, const std::array& Y, const std::array& Z); diff --git a/opm/output/eclipse/Summary.hpp b/opm/output/eclipse/Summary.hpp index 77ea9b127..6db1462b1 100644 --- a/opm/output/eclipse/Summary.hpp +++ b/opm/output/eclipse/Summary.hpp @@ -25,6 +25,7 @@ #include #include +#include #include #include diff --git a/opm/parser/eclipse/EclipseState/Grid/EclipseGrid.hpp b/opm/parser/eclipse/EclipseState/Grid/EclipseGrid.hpp index 5b481726e..3939a3cfe 100644 --- a/opm/parser/eclipse/EclipseState/Grid/EclipseGrid.hpp +++ b/opm/parser/eclipse/EclipseState/Grid/EclipseGrid.hpp @@ -26,9 +26,9 @@ #include #include #include +#include -#include -#include +#include #include #include @@ -59,10 +59,10 @@ namespace Opm { These constructors will make a copy of the src grid, with zcorn and or actnum have been adjustments. */ - EclipseGrid(const EclipseGrid& src, const double* zcorn , const std::vector& actnum); - EclipseGrid(const EclipseGrid& src, const std::vector& zcorn , const std::vector& actnum); + EclipseGrid(const EclipseGrid& src) = default; EclipseGrid(const EclipseGrid& src, const std::vector& actnum); - + EclipseGrid(const EclipseGrid& src, const double* zcorn, const std::vector& actnum); + EclipseGrid(size_t nx, size_t ny, size_t nz, double dx = 1.0, double dy = 1.0, double dz = 1.0); @@ -77,7 +77,6 @@ namespace Opm { /// explicitly. If a null pointer is passed, every cell is active. EclipseGrid(const Deck& deck, const int * actnum = nullptr); - static bool hasGDFILE(const Deck& deck); static bool hasCylindricalKeywords(const Deck& deck); static bool hasCornerPointKeywords(const Deck&); @@ -88,8 +87,7 @@ namespace Opm { size_t activeIndex(size_t i, size_t j, size_t k) const; size_t activeIndex(size_t globalIndex) const; - void save(const std::string& filename, UnitSystem::UnitType output_units) const; - void addNNC(const NNC& nnc); + void save(const std::string& filename, bool formatted, const Opm::NNC& nnc, const Opm::UnitSystem& units) const; /* Observe that the there is a getGlobalIndex(i,j,k) implementation in the base class. This method - translating @@ -105,6 +103,7 @@ namespace Opm { applied in the 'THETA' direction; this will only apply if the theta keywords entered sum up to exactly 360 degrees! */ + bool circle( ) const; bool isPinchActive( ) const; double getPinchThresholdThickness( ) const; @@ -114,7 +113,6 @@ namespace Opm { MinpvMode::ModeEnum getMinpvMode() const; const std::vector& getMinpvVector( ) const; - /* Will return a vector of nactive elements. The method will behave differently depending on the lenght of the @@ -127,6 +125,7 @@ namespace Opm { ??? : Exception. */ + template std::vector compressedVector(const std::vector& input_vector) const { if( input_vector.size() == this->getNumActive() ) { @@ -151,13 +150,13 @@ namespace Opm { /// Will return a vector a length num_active; where the value /// of each element is the corresponding global index. const std::vector& getActiveMap() const; - std::array getCellCenter(size_t i,size_t j, size_t k) const; - std::array getCellCenter(size_t globalIndex) const; + const std::array& getCellCenter(size_t i,size_t j, size_t k) const; + const std::array& getCellCenter(size_t globalIndex) const; std::array getCornerPos(size_t i,size_t j, size_t k, size_t corner_index) const; double getCellVolume(size_t globalIndex) const; double getCellVolume(size_t i , size_t j , size_t k) const; - double getCellThicknes(size_t globalIndex) const; - double getCellThicknes(size_t i , size_t j , size_t k) const; + double getCellThickness(size_t globalIndex) const; + double getCellThickness(size_t i , size_t j , size_t k) const; std::array getCellDims(size_t i,size_t j, size_t k) const; std::array getCellDims(size_t globalIndex) const; bool cellActive( size_t globalIndex ) const; @@ -166,19 +165,28 @@ namespace Opm { double getCellDepth(size_t globalIndex) const; ZcornMapper zcornMapper() const; + const std::vector& getCOORD() const; + const std::vector& getZCORN() const; + const std::vector& getACTNUM( ) const; + const std::vector& getMAPAXES() const; + const std::string& getMAPUNITS() const { return m_mapunits; } ; + /* - The exportZCORN method will adjust the z coordinates to ensure that cells do not - overlap. The return value is the number of points which have been adjusted. + The fixupZCORN method is run as part of constructiong the grid. This will adjust the + z-coordinates to ensure that cells do not overlap. The return value is the number of + points which have been adjusted. The number of zcorn nodes that has ben fixed is + stored in private member zcorn_fixed. */ - size_t exportZCORN( std::vector& zcorn) const; + + size_t fixupZCORN(); + size_t getZcornFixed() { return zcorn_fixed; }; + + // resetACTNUM with no arguments will make all cells in the grid active. + + void resetACTNUM(); + void resetACTNUM( const std::vector& actnum); - - void exportMAPAXES( std::vector& mapaxes) const; - void exportCOORD( std::vector& coord) const; - void exportACTNUM( std::vector& actnum) const; - void resetACTNUM( const int * actnum); bool equal(const EclipseGrid& other) const; - const ecl_grid_type * c_ptr() const; private: std::vector m_minpvVector; @@ -186,25 +194,35 @@ namespace Opm { Value m_pinch; PinchMode::ModeEnum m_pinchoutMode; PinchMode::ModeEnum m_multzMode; - mutable std::vector volume_cache; + mutable std::vector< int > activeMap; bool m_circle = false; - /* - The internal class grid_ptr is a a std::unique_ptr with - special copy semantics. The purpose of implementing this is - that the EclipseGrid class can now use the default - implementation for the copy and move constructors. - */ - using ert_ptr = ERT::ert_unique_ptr; - class grid_ptr : public ert_ptr { - public: - using ert_ptr::unique_ptr; - grid_ptr() = default; - grid_ptr(grid_ptr&&) = default; - grid_ptr(const grid_ptr& src) : - ert_ptr( ecl_grid_alloc_copy( src.get() ) ) {} - }; - grid_ptr m_grid; + + size_t zcorn_fixed = 0; + bool m_useActnumFromGdfile = false; + + // Input grid data. + std::vector m_zcorn; + std::vector m_coord; + std::vector m_mapaxes; + std::vector m_actnum; + std::string m_mapunits; + + // Mapping to/from active cells. + int m_nactive; + std::vector m_active_to_global; + std::vector m_global_to_active; + + // Geometry data. + std::vector m_volume; + std::vector m_depth; + std::vector m_dx; + std::vector m_dy; + std::vector m_dz; + std::vector> m_cellCenter; + + void initGridFromEGridFile(Opm::EclIO::EclFile& egridfile, std::string fileName); + void initBinaryGrid(const Deck& deck); void initCornerPointGrid(const std::array& dims , @@ -213,6 +231,8 @@ namespace Opm { const int * actnum, const double * mapaxes); + bool keywInputBeforeGdfile(const Deck& deck, const std::string keyword) const; + void initCylindricalGrid( const std::array&, const Deck&); void initCartesianGrid( const std::array&, const Deck&); void initCornerPointGrid( const std::array&, const Deck&); @@ -229,6 +249,20 @@ namespace Opm { static std::vector createDVector(const std::array& dims, size_t dim, const std::string& DKey, const std::string& DVKey, const Deck&); static void scatterDim(const std::array& dims , size_t dim , const std::vector& DV , std::vector& D); + + + std::vector makeCoordDxDyDzTops(const std::array& dims, const std::vector& dx, const std::vector& dy, const std::vector& dz, const std::vector& tops) const; + std::vector makeZcornDzTops(const std::array& dims, const std::vector& dz, const std::vector& tops) const ; + std::vector makeZcornDzvDepthz(const std::array& dims, const std::vector& dzv, const std::vector& depthz) const; + std::vector makeCoordDxvDyvDzvDepthz(const std::array& dims, const std::vector& dxv, const std::vector& dyv, const std::vector& dzv, const std::vector& depthz) const; + + double sumIdir(int j, int k, int i1, const std::array& dims, const std::vector& dx) const; + double sumJdir(int i, int k, int j1, const std::array& dims, const std::vector& dy) const; + double sumKdir(int i, int j, const std::array& dims, const std::vector& dz) const; + + void calculateGeometryData(); + + void getCellCorners(const std::array& ijk, const std::array& dims, std::array& X, std::array& Y, std::array& Z) const; }; class CoordMapper { diff --git a/src/opm/common/utility/numeric/calculateCellVol.cpp b/src/opm/common/utility/numeric/calculateCellVol.cpp index 87aceb473..748171a00 100755 --- a/src/opm/common/utility/numeric/calculateCellVol.cpp +++ b/src/opm/common/utility/numeric/calculateCellVol.cpp @@ -21,7 +21,6 @@ #include #include - /* Cell volume calculation based on following publication: @@ -89,7 +88,7 @@ double perm123sign(int i1, int i2, int i3){ if ((i1 ==1 ) && (i2==2) && (i3 == 3)){ temp = 1.0; } - + if ((i1 == 1) && (i2==3) && (i3 == 2)){ temp = -1.0; } @@ -113,7 +112,7 @@ double perm123sign(int i1, int i2, int i3){ return temp; } -double calculateCellVol(const std::vector& X, const std::vector& Y, const std::vector& Z){ +double calculateCellVol(const std::array& X, const std::array& Y, const std::array& Z){ double volume = 0.0; const double* vect[3]; @@ -134,7 +133,8 @@ double calculateCellVol(const std::vector& X, const std::vector& assert(false); vect[j] = 0; } - } + } + for (int pb = 0; pb < 2; ++pb){ for (int pg = 0; pg < 2; ++pg){ diff --git a/src/opm/io/eclipse/EclFile.cpp b/src/opm/io/eclipse/EclFile.cpp index 11c8f1714..928cab187 100644 --- a/src/opm/io/eclipse/EclFile.cpp +++ b/src/opm/io/eclipse/EclFile.cpp @@ -18,6 +18,7 @@ #include #include +#include #include #include @@ -28,24 +29,18 @@ #include #include #include - -#include - -#include - -//#include -//#include -//#include #include -//#include -//#include - // anonymous namespace for EclFile namespace { +bool fileExists(const std::string& filename){ + + std::ifstream fileH(filename.c_str()); + return fileH.good(); +} bool isFormatted(const std::string& filename) { @@ -471,6 +466,11 @@ namespace Opm { namespace EclIO { EclFile::EclFile(const std::string& filename) : inputFilename(filename) { + if (!fileExists(filename)){ + std::string message="Could not open EclFile: " + filename; + OPM_THROW(std::invalid_argument, message); + } + std::fstream fileH; formatted = isFormatted(filename); diff --git a/src/opm/output/eclipse/EclipseIO.cpp b/src/opm/output/eclipse/EclipseIO.cpp index 9ae36d0ec..8eac280f2 100644 --- a/src/opm/output/eclipse/EclipseIO.cpp +++ b/src/opm/output/eclipse/EclipseIO.cpp @@ -118,8 +118,7 @@ void EclipseIO::Impl::writeEGRIDFile( const NNC& nnc ) { ECL_EGRID_FILE, ioConfig.getFMTOUT() )); - this->grid.addNNC( nnc ); - this->grid.save( egridFile, this->es.getDeckUnitSystem().getType()); + this->grid.save( egridFile, ioConfig.getFMTOUT(), nnc, this->es.getDeckUnitSystem()); } /* diff --git a/src/opm/parser/eclipse/EclipseState/EclipseState.cpp b/src/opm/parser/eclipse/EclipseState/EclipseState.cpp index 43c88d34a..f6ae29e77 100644 --- a/src/opm/parser/eclipse/EclipseState/EclipseState.cpp +++ b/src/opm/parser/eclipse/EclipseState/EclipseState.cpp @@ -61,7 +61,7 @@ namespace Opm { m_simulationConfig( m_eclipseConfig.getInitConfig().restartRequested(), deck, m_eclipseProperties ), m_transMult( GridDims(deck), deck, m_eclipseProperties ) { - m_inputGrid.resetACTNUM(m_eclipseProperties.getIntGridProperty("ACTNUM").getData().data()); + m_inputGrid.resetACTNUM(m_eclipseProperties.getIntGridProperty("ACTNUM").getData()); if( this->runspec().phases().size() < 3 ) OpmLog::info("Only " + std::to_string( this->runspec().phases().size() ) diff --git a/src/opm/parser/eclipse/EclipseState/Grid/EclipseGrid.cpp b/src/opm/parser/eclipse/EclipseState/Grid/EclipseGrid.cpp index 5129a1ddf..1a8acadf9 100644 --- a/src/opm/parser/eclipse/EclipseState/Grid/EclipseGrid.cpp +++ b/src/opm/parser/eclipse/EclipseState/Grid/EclipseGrid.cpp @@ -19,6 +19,7 @@ #define _USE_MATH_DEFINES #include +#include #include #include @@ -27,6 +28,9 @@ #include #include +#include +#include + #include #include #include @@ -50,151 +54,254 @@ #include -#include namespace Opm { - EclipseGrid::EclipseGrid(std::array& dims , - const std::vector& coord , - const std::vector& zcorn , - const int * actnum, - const double * mapaxes) - : GridDims(dims), - m_minpvMode(MinpvMode::ModeEnum::Inactive), - m_pinch("PINCH"), - m_pinchoutMode(PinchMode::ModeEnum::TOPBOT), - m_multzMode(PinchMode::ModeEnum::TOP), - volume_cache(dims[0] * dims[1] * dims[2], -1.0) - { - initCornerPointGrid( dims, coord , zcorn , actnum , mapaxes ); +EclipseGrid::EclipseGrid(std::array& dims , + const std::vector& coord , + const std::vector& zcorn , + const int * actnum, + const double * mapaxes) + : GridDims(dims), + m_minpvMode(MinpvMode::ModeEnum::Inactive), + m_pinch("PINCH"), + m_pinchoutMode(PinchMode::ModeEnum::TOPBOT), + m_multzMode(PinchMode::ModeEnum::TOP) +{ + initCornerPointGrid( dims, coord , zcorn , actnum , mapaxes ); +} + +/** + Will create an EclipseGrid instance based on an existing + GRID/EGRID file. +*/ + + +EclipseGrid::EclipseGrid(const std::string& fileName ) + : GridDims(), + m_minpvMode(MinpvMode::ModeEnum::Inactive), + m_pinch("PINCH"), + m_pinchoutMode(PinchMode::ModeEnum::TOPBOT), + m_multzMode(PinchMode::ModeEnum::TOP) +{ + + Opm::EclIO::EclFile egridfile(fileName); + + m_useActnumFromGdfile=true; + + initGridFromEGridFile(egridfile, fileName); + + resetACTNUM(m_actnum); + + calculateGeometryData(); + +} + + +EclipseGrid::EclipseGrid(size_t nx, size_t ny , size_t nz, + double dx, double dy, double dz) + : GridDims(nx, ny, nz), + m_minpvMode(MinpvMode::ModeEnum::Inactive), + m_pinch("PINCH"), + m_pinchoutMode(PinchMode::ModeEnum::TOPBOT), + m_multzMode(PinchMode::ModeEnum::TOP) +{ + + m_coord.reserve((nx+1)*(ny+1)*6); + + for (size_t j = 0; j < (ny+1); j++) { + for (size_t i = 0; i < (nx+1); i++) { + m_coord.push_back(i*dx); + m_coord.push_back(j*dy); + m_coord.push_back(0.0); + m_coord.push_back(i*dx); + m_coord.push_back(j*dy); + m_coord.push_back(nz*dz); + } } + m_zcorn.assign(nx*ny*nz*8, 0); + for (size_t k = 0; k < nz ; k++) { + for (size_t j = 0; j < ny ; j++) { + for (size_t i = 0; i < nx ; i++) { - /** - Will create an EclipseGrid instance based on an existing - GRID/EGRID file. - */ - EclipseGrid::EclipseGrid(const std::string& filename ) - : GridDims(), - m_minpvMode(MinpvMode::ModeEnum::Inactive), - m_pinch("PINCH"), - m_pinchoutMode(PinchMode::ModeEnum::TOPBOT), - m_multzMode(PinchMode::ModeEnum::TOP) - { - ecl_grid_type * new_ptr = ecl_grid_load_case__( filename.c_str() , false ); - if (new_ptr) - m_grid.reset( new_ptr ); - else - throw std::invalid_argument("Could not load grid from binary file: " + filename); + // top face of cell + int zind = i*2 + j*nx*4 + k*nx*ny*8; - m_nx = ecl_grid_get_nx( c_ptr() ); - m_ny = ecl_grid_get_ny( c_ptr() ); - m_nz = ecl_grid_get_nz( c_ptr() ); + double zt = k*dz; + double zb = (k+1)*dz; - volume_cache.resize(m_nx * m_ny * m_nz, -1.0); - } + m_zcorn[zind] = zt; + m_zcorn[zind + 1] = zt; + zind = zind + nx*2; - EclipseGrid::EclipseGrid(size_t nx, size_t ny , size_t nz, - double dx, double dy, double dz) - : GridDims(nx, ny, nz), - m_minpvMode(MinpvMode::ModeEnum::Inactive), - m_pinch("PINCH"), - m_pinchoutMode(PinchMode::ModeEnum::TOPBOT), - m_multzMode(PinchMode::ModeEnum::TOP), - volume_cache(nx * ny * nz, -1.0), - m_grid( ecl_grid_alloc_rectangular(nx, ny, nz, dx, dy, dz, NULL) ) - { - } + m_zcorn[zind] = zt; + m_zcorn[zind + 1] = zt; - EclipseGrid::EclipseGrid(const EclipseGrid& src, const double* zcorn , const std::vector& actnum) - : GridDims(src.getNX(), src.getNY(), src.getNZ()), - m_minpvVector( src.m_minpvVector ), - m_minpvMode( src.m_minpvMode ), - m_pinch( src.m_pinch ), - m_pinchoutMode( src.m_pinchoutMode ), - m_multzMode( src.m_multzMode ), - volume_cache(src.volume_cache.size(), -1.0) - { - const int * actnum_data = (actnum.empty()) ? nullptr : actnum.data(); - m_grid.reset( ecl_grid_alloc_processed_copy( src.c_ptr(), zcorn , actnum_data )); - } + // bottom face of cell + zind = i*2 + j*nx*4 + k*nx*ny*8 + nx*ny*4; + m_zcorn[zind] = zb; + m_zcorn[zind + 1] = zb; - EclipseGrid::EclipseGrid(const EclipseGrid& src, const std::vector& zcorn , const std::vector& actnum) - : EclipseGrid( src , (zcorn.empty()) ? nullptr : zcorn.data(), actnum ) - { } + zind = zind + nx*2; - - EclipseGrid::EclipseGrid(const EclipseGrid& src, const std::vector& actnum) - : EclipseGrid( src , nullptr , actnum ) - { } - - - - /* - This is the main EclipseGrid constructor, it will inspect the - input Deck for grid keywords, either the corner point keywords - COORD and ZCORN, or the various rectangular keywords like DX,DY - and DZ. - - Actnum is treated specially: - - 1. If an actnum pointer is passed in that should be a pointer - to 0 and 1 values which will be used as actnum mask. - - 2. If the actnum pointer is not present the constructor will - look in the deck for an actnum keyword, and use that if it - is found. This is a best effort which will work in many - cases, but if the ACTNUM keyword is manipulated in the deck - those manipulations will be silently lost; if the ACTNUM - keyword has size different from nx*ny*nz it will also be - silently ignored. - - With a mutable EclipseGrid instance you can later call the - EclipseGrid::resetACTNUM() method when you have complete actnum - information. The EclipseState based construction of EclipseGrid - is a two-pass operation, which guarantees that actnum is handled - correctly. - */ - - - - - EclipseGrid::EclipseGrid(const Deck& deck, const int * actnum) - : GridDims(deck), - m_minpvMode(MinpvMode::ModeEnum::Inactive), - m_pinch("PINCH"), - m_pinchoutMode(PinchMode::ModeEnum::TOPBOT), - m_multzMode(PinchMode::ModeEnum::TOP), - volume_cache(m_nx * m_ny * m_nz, -1.0) - { - - const std::array dims = getNXYZ(); - initGrid(dims, deck); - - if (actnum != nullptr) - resetACTNUM(actnum); - else { - if (deck.hasKeyword()) { - const auto& actnumData = deck.getKeyword().getIntData(); - if (actnumData.size() == getCartesianSize()) - resetACTNUM( actnumData.data()); - else { - const std::string msg = "The ACTNUM keyword has " + std::to_string( actnumData.size() ) + " elements - expected : " + std::to_string( getCartesianSize()) + " - ignored."; - OpmLog::warning(msg); - } + m_zcorn[zind] = zb; + m_zcorn[zind + 1] = zb; } } } + resetACTNUM(); + calculateGeometryData(); +} + +EclipseGrid::EclipseGrid(const EclipseGrid& src, const double* zcorn, const std::vector& actnum) + : EclipseGrid(src) +{ + + if (zcorn != nullptr) { + std::array dims = getNXYZ(); + + size_t sizeZcorn = dims[0]*dims[1]*dims[2]*8; + + for (size_t n=0; n < sizeZcorn; n++) { + m_zcorn[n] = zcorn[n]; + } + + ZcornMapper mapper( getNX(), getNY(), getNZ()); + zcorn_fixed = mapper.fixupZCORN( m_zcorn ); + + calculateGeometryData(); + } + + resetACTNUM(actnum); +} + +EclipseGrid::EclipseGrid(const EclipseGrid& src, const std::vector& actnum) + : EclipseGrid( src , nullptr , actnum ) +{ } + +/* + This is the main EclipseGrid constructor, it will inspect the + input Deck for grid keywords, either the corner point keywords + COORD and ZCORN, or the various rectangular keywords like DX,DY + and DZ. + + Actnum is treated specially: + + 1. If an actnum pointer is passed in that should be a pointer + to 0 and 1 values which will be used as actnum mask. + + 2. If the actnum pointer is not present the constructor will + look in the deck for an actnum keyword, and use that if it + is found. This is a best effort which will work in many + cases, but if the ACTNUM keyword is manipulated in the deck + those manipulations will be silently lost; if the ACTNUM + keyword has size different from nx*ny*nz it will also be + silently ignored. + + With a mutable EclipseGrid instance you can later call the + EclipseGrid::resetACTNUM() method when you have complete actnum + information. The EclipseState based construction of EclipseGrid + is a two-pass operation, which guarantees that actnum is handled + correctly. +*/ + + + +EclipseGrid::EclipseGrid(const Deck& deck, const int * actnum) + : GridDims(deck), + m_minpvMode(MinpvMode::ModeEnum::Inactive), + m_pinch("PINCH"), + m_pinchoutMode(PinchMode::ModeEnum::TOPBOT), + m_multzMode(PinchMode::ModeEnum::TOP) +{ + + if (deck.hasKeyword("GDFILE")){ + + if (deck.hasKeyword("COORD")){ + throw std::invalid_argument("COORD can't be used together with GDFILE"); + } + + if (deck.hasKeyword("ZCORN")){ + throw std::invalid_argument("ZCORN can't be used together with GDFILE"); + } + + if (deck.hasKeyword("ACTNUM")){ + if (keywInputBeforeGdfile(deck, "ACTNUM")) { + m_useActnumFromGdfile = true; + } + } else { + m_useActnumFromGdfile = true; + } + + } + + + const std::array dims = getNXYZ(); + + initGrid(dims, deck); + + if (deck.hasKeyword("MAPUNITS")){ + if ((m_mapunits =="") || ( !keywInputBeforeGdfile(deck, "MAPUNITS"))) { + const auto& record = deck.getKeyword( ).getStringData(); + m_mapunits = record[0]; + } + } + + if (deck.hasKeyword("MAPAXES")){ + if ((m_mapaxes.size() == 0) || ( !keywInputBeforeGdfile(deck, "MAPAXES"))) { + const Opm::DeckKeyword& mapaxesKeyword = deck.getKeyword(); + m_mapaxes.resize(6); + for (size_t n=0; n < 6; n++){ + m_mapaxes[n] = mapaxesKeyword.getRecord(0).getItem(n).get(0); + } + } + } + + if (actnum != nullptr) { + + size_t nCells = dims[0]*dims[1]*dims[2]; + std::vector actnumVector(actnum, actnum + nCells); + + resetACTNUM( actnumVector); + } else { + + if (m_useActnumFromGdfile){ + // actnum already reset in initBinaryGrid + } else if (deck.hasKeyword()) { + + const auto& actnumData = deck.getKeyword().getIntData(); + + if (actnumData.size() == getCartesianSize()) { + resetACTNUM( actnumData); + } else { + const std::string msg = "The ACTNUM keyword has " + std::to_string( actnumData.size() ) + " elements - expected : " + std::to_string( getCartesianSize()) + " - ignored."; + OpmLog::warning(msg); + + resetACTNUM( ); + } + + } else { + + resetACTNUM(); + } + } + + calculateGeometryData(); +} + + bool EclipseGrid::circle( ) const{ return this->m_circle; } void EclipseGrid::initGrid( const std::array& dims, const Deck& deck) { + if (deck.hasKeyword()) { initCylindricalGrid( dims, deck ); } else { @@ -238,7 +345,6 @@ namespace Opm { m_minpvMode = MinpvMode::ModeEnum::EclSTD; } - if (deck.hasKeyword()) { const auto& record = deck.getKeyword( ).getRecord(0); const auto& item = record.getItem( ); @@ -246,29 +352,138 @@ namespace Opm { m_minpvMode = MinpvMode::ModeEnum::OpmFIL; } } + + void EclipseGrid::initGridFromEGridFile(Opm::EclIO::EclFile& egridfile, std::string fileName){ + if (!egridfile.hasKey("GRIDHEAD")) { + throw std::invalid_argument("file: " + fileName + " is not a valid egrid file, GRIDHEAD not found"); + } + if (!egridfile.hasKey("COORD")) { + throw std::invalid_argument("file: " + fileName + " is not a valid egrid file, COORD not found"); + } + + if (!egridfile.hasKey("ZCORN")) { + throw std::invalid_argument("file: " + fileName + " is not a valid egrid file, ZCORN not found"); + } + + if (!egridfile.hasKey("GRIDUNIT")) { + throw std::invalid_argument("file: " + fileName + " is not a valid egrid file, ZCORN not found"); + } + + const std::vector& gridunit = egridfile.get("GRIDUNIT"); + + const std::vector& gridhead = egridfile.get("GRIDHEAD"); + const std::array dims = {gridhead[1], gridhead[2], gridhead[3]}; + + m_nx=dims[0]; + m_ny=dims[1]; + m_nz=dims[2]; + + const std::vector& coord_f = egridfile.get("COORD"); + const std::vector& zcorn_f = egridfile.get("ZCORN"); + + m_coord.assign(coord_f.begin(), coord_f.end()); + m_zcorn.assign(zcorn_f.begin(), zcorn_f.end()); + + if (gridunit[0] != "METRES") { + + const auto length = ::Opm::UnitSystem::measure::length; + + if (gridunit[0] == "FEET"){ + Opm::UnitSystem units(Opm::UnitSystem::UnitType::UNIT_TYPE_FIELD ); + units.to_si(length, m_coord); + units.to_si(length, m_zcorn); + } else if (gridunit[0] == "CM"){ + Opm::UnitSystem units(Opm::UnitSystem::UnitType::UNIT_TYPE_LAB ); + units.to_si(length, m_coord); + units.to_si(length, m_zcorn); + } else { + std::string message = "gridunit '" + gridunit[0] + "' doesn't correspong to a valid unit system"; + throw std::invalid_argument(message); + } + } + + if ((egridfile.hasKey("ACTNUM")) && (m_useActnumFromGdfile)) { + const std::vector& actnum = egridfile.get("ACTNUM"); + resetACTNUM( actnum ); + } + + if (egridfile.hasKey("MAPAXES")) { + const std::vector& mapaxes_f = egridfile.get("MAPAXES"); + m_mapaxes.assign(mapaxes_f.begin(), mapaxes_f.end()); + } + + if (egridfile.hasKey("MAPUNITS")) { + const std::vector& mapunits = egridfile.get("MAPUNITS"); + m_mapunits=mapunits[0]; + } + + ZcornMapper mapper( getNX(), getNY(), getNZ()); + zcorn_fixed = mapper.fixupZCORN( m_zcorn ); + } + + bool EclipseGrid::keywInputBeforeGdfile(const Deck& deck, const std::string keyword) const { + + std::vector keywordList; + keywordList.reserve(deck.size()); + + for (size_t n=0;n( active_index ); + } + + return m_global_to_active[ globalIndex]; } /** Observe: the input argument is assumed to be in the space [0,num_active). */ + size_t EclipseGrid::getGlobalIndex(size_t active_index) const { - int global_index = ecl_grid_get_global_index1A( m_grid.get() , active_index ); - return static_cast(global_index); + + return m_active_to_global.at(active_index); } size_t EclipseGrid::getGlobalIndex(size_t i, size_t j, size_t k) const { + return GridDims::getGlobalIndex(i,j,k); } @@ -297,30 +512,30 @@ namespace Opm { return m_minpvVector; } - void EclipseGrid::initBinaryGrid(const Deck& deck) { + const DeckKeyword& gdfile_kw = deck.getKeyword("GDFILE"); const std::string& gdfile_arg = gdfile_kw.getRecord(0).getItem("filename").get(0); std::string filename = deck.makeDeckPath(gdfile_arg); - ecl_grid_type * grid_ptr_loc = ecl_grid_load_case__( filename.c_str(), false); - if (grid_ptr_loc) - this->m_grid.reset( grid_ptr_loc ); - else - throw std::invalid_argument("Failed to load grid from: " + filename); + Opm::EclIO::EclFile egridfile(filename); + + initGridFromEGridFile(egridfile, filename); } void EclipseGrid::initCartesianGrid(const std::array& dims , const Deck& deck) { - if (hasDVDEPTHZKeywords( deck )) + + if (hasDVDEPTHZKeywords( deck )) { initDVDEPTHZGrid( dims , deck ); - else if (hasDTOPSKeywords(deck)) + } else if (hasDTOPSKeywords(deck)) { initDTOPSGrid( dims ,deck ); - else + } else { throw std::invalid_argument("Tried to initialize cartesian grid without all required keywords"); + } } - void EclipseGrid::initDVDEPTHZGrid(const std::array& dims, const Deck& deck) { + const std::vector& DXV = deck.getKeyword().getSIDoubleData(); const std::vector& DYV = deck.getKeyword().getSIDoubleData(); const std::vector& DZV = deck.getKeyword().getSIDoubleData(); @@ -331,18 +546,435 @@ namespace Opm { assertVectorSize( DYV , static_cast( dims[1] ) , "DYV"); assertVectorSize( DZV , static_cast( dims[2] ) , "DZV"); - m_grid.reset( ecl_grid_alloc_dxv_dyv_dzv_depthz( dims[0] , dims[1] , dims[2] , DXV.data() , DYV.data() , DZV.data() , DEPTHZ.data() , nullptr ) ); + m_coord = makeCoordDxvDyvDzvDepthz(dims, DXV, DYV, DZV, DEPTHZ); + m_zcorn = makeZcornDzvDepthz(dims, DZV, DEPTHZ); + + ZcornMapper mapper( getNX(), getNY(), getNZ()); + zcorn_fixed = mapper.fixupZCORN( m_zcorn ); } - void EclipseGrid::initDTOPSGrid(const std::array& dims , const Deck& deck) { + std::vector DX = createDVector( dims , 0 , "DX" , "DXV" , deck); std::vector DY = createDVector( dims , 1 , "DY" , "DYV" , deck); std::vector DZ = createDVector( dims , 2 , "DZ" , "DZV" , deck); std::vector TOPS = createTOPSVector( dims , DZ , deck ); - m_grid.reset( ecl_grid_alloc_dx_dy_dz_tops( dims[0] , dims[1] , dims[2] , DX.data() , DY.data() , DZ.data() , TOPS.data() , nullptr ) ); + + m_coord = makeCoordDxDyDzTops(dims, DX, DY, DZ, TOPS); + m_zcorn = makeZcornDzTops(dims, DZ, TOPS); + + ZcornMapper mapper( getNX(), getNY(), getNZ()); + zcorn_fixed = mapper.fixupZCORN( m_zcorn ); } + void EclipseGrid::calculateGeometryData() { + + const std::array dims = getNXYZ(); + int nCells = getCartesianSize(); + + std::array X = {0.0}; + std::array Y = {0.0}; + std::array Z = {0.0}; + + m_volume.clear(); + m_cellCenter.clear(); + m_dx.clear(); + m_dy.clear(); + m_dz.clear(); + m_depth.clear(); + + m_volume.resize(nCells); + m_cellCenter.resize(nCells); + m_dx.resize(nCells); + m_dy.resize(nCells); + m_dz.resize(nCells); + m_depth.resize(nCells); + + std::array ijk; + size_t n = 0; + + for (size_t k = 0; k< getNZ(); k++) { + for (size_t j = 0; j< getNY(); j++) { + for (size_t i = 0; i< getNX(); i++) { + + ijk[0] = i; + ijk[1] = j; + ijk[2] = k; + + getCellCorners(ijk, dims, X, Y, Z ); + + m_volume[n] = calculateCellVol(X, Y, Z); + + std::array center; + + center[0]=std::accumulate(X.begin(), X.end(), 0.0) / 8.0; + center[1]=std::accumulate(Y.begin(), Y.end(), 0.0) / 8.0; + center[2]=std::accumulate(Z.begin(), Z.end(), 0.0) / 8.0; + + m_cellCenter[n] = center; + + // calculate dx + + double x2 = (X[1]+X[3]+X[5]+X[7])/4.0; + double y2 = (Y[1]+Y[3]+Y[5]+Y[7])/4.0; + + double x1 = (X[0]+X[2]+X[4]+X[6])/4.0; + double y1 = (Y[0]+Y[2]+Y[4]+Y[6])/4.0; + + double dx = sqrt(pow((x2-x1), 2.0) + pow((y2-y1), 2.0) ); + + // calculate dy + + x2 = (X[2]+X[3]+X[6]+X[7])/4.0; + y2 = (Y[2]+Y[3]+Y[6]+Y[7])/4.0; + + x1 = (X[0]+X[1]+X[4]+X[5])/4.0; + y1 = (Y[0]+Y[1]+Y[4]+Y[5])/4.0; + + double dy = sqrt(pow((x2-x1), 2.0) + pow((y2-y1), 2.0)); + + // calculate dz + + double z2 = (Z[4]+Z[5]+Z[6]+Z[7])/4.0; + double z1 = (Z[0]+Z[1]+Z[2]+Z[3])/4.0; + + double dz = z2-z1; + + m_dx[n] = dx; + m_dy[n] = dy; + m_dz[n] = dz; + + m_depth[n] = (z2 + z1) / 2.0; + + n++; + } + } + } + } + + void EclipseGrid::getCellCorners(const std::array& ijk, const std::array& dims, + std::array& X, + std::array& Y, + std::array& Z) const + { + + std::array zind; + std::array pind; + + // calculate indices for grid pillars in COORD arrray + const size_t p_offset = ijk[1]*(dims[0]+1)*6 + ijk[0]*6; + + pind[0] = p_offset; + pind[1] = p_offset + 6; + pind[2] = p_offset + (dims[0]+1)*6; + pind[3] = pind[2] + 6; + + // get depths from zcorn array in ZCORN array + const size_t z_offset = ijk[2]*dims[0]*dims[1]*8 + ijk[1]*dims[0]*4 + ijk[0]*2; + + zind[0] = z_offset; + zind[1] = z_offset + 1; + zind[2] = z_offset + dims[0]*2; + zind[3] = zind[2] + 1; + + for (int n = 0; n < 4; n++) + zind[n+4] = zind[n] + dims[0]*dims[1]*4; + + + for (int n = 0; n< 8; n++) + Z[n] = m_zcorn[zind[n]]; + + + for (int n=0; n<4; n++) { + double xt = m_coord[pind[n]]; + double yt = m_coord[pind[n] + 1]; + double zt = m_coord[pind[n] + 2]; + + double xb = m_coord[pind[n] + 3]; + double yb = m_coord[pind[n] + 4]; + double zb = m_coord[pind[n]+5]; + + X[n] = xt + (xb-xt) / (zt-zb) * (zt - Z[n]); + X[n+4] = xt + (xb-xt) / (zt-zb) * (zt-Z[n+4]); + + Y[n] = yt+(yb-yt)/(zt-zb)*(zt-Z[n]); + Y[n+4] = yt+(yb-yt)/(zt-zb)*(zt-Z[n+4]); + } + } + + std::vector EclipseGrid::makeCoordDxvDyvDzvDepthz(const std::array& dims, const std::vector& dxv, const std::vector& dyv, const std::vector& dzv, const std::vector& depthz) const { + + std::vector coord; + coord.reserve((dims[0]+1)*(dims[1]+1)*6); + + std::vector x(dims[0] + 1, 0.0); + std::partial_sum(dxv.begin(), dxv.end(), x.begin() + 1); + + std::vector y(dims[1] + 1, 0.0); + std::partial_sum(dyv.begin(), dyv.end(), y.begin() + 1); + + std::vector z(dims[2] + 1, 0.0); + std::partial_sum(dzv.begin(), dzv.end(), z.begin() + 1); + + + for (int j = 0; j < dims[1] + 1; j++) { + for (int i = 0; i EclipseGrid::makeZcornDzvDepthz(const std::array& dims, const std::vector& dzv, const std::vector& depthz) const { + + std::vector zcorn; + size_t sizeZcorn = dims[0]*dims[1]*dims[2]*8; + + zcorn.assign (sizeZcorn, 0.0); + + std::vector z(dims[2] + 1, 0.0); + std::partial_sum(dzv.begin(), dzv.end(), z.begin() + 1); + + + for (int k = 0; k < dims[2]; k++) { + for (int j = 0; j < dims[1]; j++) { + for (int i = 0; i < dims[0]; i++) { + + const double z0 = z[k]; + + // top face of cell + int zind = i*2 + j*dims[0]*4 + k*dims[0]*dims[1]*8; + + zcorn[zind] = depthz[i+j*(dims[0]+1)] + z0; + zcorn[zind + 1] = depthz[i+j*(dims[0]+1) +1 ] + z0; + + zind = zind + dims[0]*2; + + zcorn[zind] = depthz[i+(j+1)*(dims[0]+1)] + z0; + zcorn[zind + 1] = depthz[i+(j+1)*(dims[0]+1) + 1] + z0; + + // bottom face of cell + zind = i*2 + j*dims[0]*4 + k*dims[0]*dims[1]*8 + dims[0]*dims[1]*4; + + zcorn[zind] = depthz[i+j*(dims[0]+1)] + z0 + dzv[k]; + zcorn[zind + 1] = depthz[i+j*(dims[0]+1) +1 ] + z0 + dzv[k]; + + zind = zind + dims[0]*2; + + zcorn[zind] = depthz[i+(j+1)*(dims[0]+1)] + z0 + dzv[k]; + zcorn[zind + 1] = depthz[i+(j+1)*(dims[0]+1) + 1] + z0 + dzv[k]; + } + } + } + + return zcorn; + } + + std::vector EclipseGrid::makeCoordDxDyDzTops(const std::array& dims , const std::vector& dx, const std::vector& dy, const std::vector& dz, const std::vector& tops) const { + + std::vector coord; + coord.reserve((dims[0]+1)*(dims[1]+1)*6); + + for (int j = 0; j < dims[1]; j++) { + + double y0 = 0; + double zt = tops[0]; + double zb = zt + sumKdir(0, 0, dims, dz); + + if (j == 0) { + double x0 = 0.0; + + coord.push_back(x0); + coord.push_back(y0); + coord.push_back(zt); + coord.push_back(x0); + coord.push_back(y0); + coord.push_back(zb); + + for (int i = 0; i < dims[0]; i++) { + + size_t ind = i+j*dims[0]+1; + + if (i == (dims[0]-1)) { + ind=ind-1; + } + + zt = tops[ind]; + zb = zt + sumKdir(i, j, dims, dz); + + double xt = x0 + dx[i+j*dims[0]]; + double xb = sumIdir(i, j, dims[2]-1, dims, dx); + + coord.push_back(xt); + coord.push_back(y0); + coord.push_back(zt); + coord.push_back(xb); + coord.push_back(y0); + coord.push_back(zb); + + x0=xt; + } + } + + int ind = (j+1)*dims[0]; + + if (j == (dims[1]-1) ) { + ind = j*dims[0]; + } + + double x0 = 0.0; + + double yt = sumJdir(0, j, 0, dims, dy); + double yb = sumJdir(0, j, dims[2]-1, dims, dy); + + zt = tops[ind]; + zb = zt + sumKdir(0, j, dims, dz); + + coord.push_back(x0); + coord.push_back(yt); + coord.push_back(zt); + coord.push_back(x0); + coord.push_back(yb); + coord.push_back(zb); + + for (int i=0; i < dims[0]; i++) { + + ind = i+(j+1)*dims[0]+1; + + if (j == (dims[1]-1) ) { + ind = i+j*dims[0]+1; + } + + if (i == (dims[0] - 1) ) { + ind=ind-1; + } + + zt = tops[ind]; + zb = zt + sumKdir(i, j, dims, dz); + + double xt=-999; + double xb=-999; + + if (j == (dims[1]-1) ) { + xt = sumIdir(i, j, 0, dims, dx); + xb = sumIdir(i, j, dims[2]-1, dims, dx); + } else { + xt = sumIdir(i, j+1, 0, dims, dx); + xb = sumIdir(i, j+1, dims[2]-1, dims, dx); + } + + if (i == (dims[0] - 1) ) { + yt = sumJdir(i, j, 0, dims, dy); + yb = sumJdir(i, j, dims[2]-1, dims, dy); + } else { + yt = sumJdir(i+1, j, 0, dims, dy); + yb = sumJdir(i+1, j, dims[2]-1, dims, dy); + } + + coord.push_back(xt); + coord.push_back(yt); + coord.push_back(zt); + coord.push_back(xb); + coord.push_back(yb); + coord.push_back(zb); + + x0 = xt; + } + } + + return coord; + }; + + std::vector EclipseGrid::makeZcornDzTops(const std::array& dims, const std::vector& dz, const std::vector& tops) const { + + std::vector zcorn; + size_t sizeZcorn = dims[0]*dims[1]*dims[2]*8; + + zcorn.assign (sizeZcorn, 0.0); + + for (int j = 0; j < dims[1]; j++) { + for (int i = 0; i < dims[0]; i++) { + int ind = i + j*dims[0]; + double z = tops[ind]; + + for (int k = 0; k < dims[2]; k++) { + + // top face of cell + int zind = i*2 + j*dims[0]*4 + k*dims[0]*dims[1]*8; + + zcorn[zind] = z; + zcorn[zind + 1] = z; + + zind = zind + dims[0]*2; + + zcorn[zind] = z; + zcorn[zind + 1] = z; + + z = z + dz[i + j*dims[0] + k*dims[0]*dims[1]]; + + // bottom face of cell + zind = i*2 + j*dims[0]*4 + k*dims[0]*dims[1]*8 + dims[0]*dims[1]*4; + + zcorn[zind] = z; + zcorn[zind + 1] = z; + + zind = zind + dims[0]*2; + + zcorn[zind] = z; + zcorn[zind + 1] = z; + } + } + } + + return zcorn; + }; + + double EclipseGrid::sumIdir(int i1, int j, int k, const std::array& dims, const std::vector& dx) const { + + double sum = 0.0; + + for (int i = 0; i <= i1; i++) { + sum = sum + dx[i + j*dims[0] + k*dims[0]*dims[1]]; + } + + return sum; + } + + double EclipseGrid::sumJdir(int i, int j1, int k, const std::array& dims, const std::vector& dy) const { + + double sum = 0.0; + + for (int j = 0; j <= j1; j++) { + sum=sum + dy[i + j*dims[0] + k*dims[0]*dims[1]]; + } + + return sum; + } + + double EclipseGrid::sumKdir(int i, int j, const std::array& dims, const std::vector& dz) const { + + double sum = 0.0; + + for (int k = 0; k < dims[2]; k++) { + sum = sum + dz[i + j*dims[0] + k*dims[0]*dims[1]]; + } + + return sum; + } /* Limited implementaton - requires keywords: DRV, DTHETAV, DZV and TOPS. @@ -470,25 +1102,36 @@ namespace Opm { const std::vector& zcorn , const int * actnum, const double * mapaxes) + + { - float * mapaxes_float = nullptr; - if (mapaxes) { - mapaxes_float = new float[6]; - for (size_t i=0; i < 6; i++) - mapaxes_float[i] = mapaxes[i]; + size_t nCells=dims[0]*dims[1]*dims[2]; + + m_coord = coord; + m_zcorn = zcorn; + + ZcornMapper mapper( getNX(), getNY(), getNZ()); + zcorn_fixed = mapper.fixupZCORN( m_zcorn ); + + m_actnum.clear(); + + if (actnum != nullptr){ + m_actnum.resize(nCells); + + for (size_t n = 0; n < nCells; n++){ + m_actnum[n] = actnum[n]; + } + + } else { + m_actnum.assign(nCells, 1); } - m_grid.reset( ecl::ecl_grid_alloc_GRDECL_data(dims[0] , - dims[1] , - dims[2] , - zcorn.data() , - coord.data() , - actnum , - false, // We do not apply the MAPAXES transformations - mapaxes_float) ); - - if (mapaxes) - delete[] mapaxes_float; + if (mapaxes != nullptr){ + m_mapaxes.resize(6); + for (size_t n = 0; n < 6; n++){ + m_mapaxes[n] = mapaxes[n]; + } + } } void EclipseGrid::initCornerPointGrid(const std::array& dims, const Deck& deck) { @@ -496,26 +1139,23 @@ namespace Opm { { const auto& ZCORNKeyWord = deck.getKeyword(); const auto& COORDKeyWord = deck.getKeyword(); + const std::vector& zcorn = ZCORNKeyWord.getSIDoubleData(); const std::vector& coord = COORDKeyWord.getSIDoubleData(); - double * mapaxes = nullptr; - - if (deck.hasKeyword()) { - const auto& mapaxesKeyword = deck.getKeyword(); - const auto& record = mapaxesKeyword.getRecord(0); - mapaxes = new double[6]; - for (size_t i = 0; i < 6; i++) { - mapaxes[i] = record.getItem( i ).getSIDouble( 0 ); - } - } - initCornerPointGrid( dims, coord , zcorn , nullptr , mapaxes ); - if (mapaxes) - delete[] mapaxes; + + int * actnum = nullptr; + + if (deck.hasKeyword()) { + const auto& actnumKeyword = deck.getKeyword(); + std::vector actnumVector = actnumKeyword.getIntData(); + + actnum=actnumVector.data(); + } + + initCornerPointGrid( dims, coord , zcorn, actnum, nullptr ); } } - - bool EclipseGrid::hasCornerPointKeywords(const Deck& deck) { if (deck.hasKeyword() && deck.hasKeyword()) return true; @@ -599,8 +1239,6 @@ namespace Opm { return false; } - - void EclipseGrid::assertVectorSize(const std::vector& vector , size_t expectedSize , const std::string& vectorName) { if (vector.size() != expectedSize) throw std::invalid_argument("Wrong size for keyword: " + vectorName + ". Expected: " + std::to_string(expectedSize) + " got: " + std::to_string(vector.size())); @@ -722,22 +1360,45 @@ namespace Opm { } } - const ecl_grid_type * EclipseGrid::c_ptr() const { - return m_grid.get(); - } - - bool EclipseGrid::equal(const EclipseGrid& other) const { - bool status = (m_pinch.equal( other.m_pinch ) && (ecl_grid_compare( c_ptr() , other.c_ptr() , true , false , false )) && (m_minpvMode == other.getMinpvMode())); - if(m_minpvMode!=MinpvMode::ModeEnum::Inactive){ + + //double reltol = 1.0e-6; + + if (m_coord.size() != other.m_coord.size()) + return false; + + if (m_zcorn.size() != other.m_zcorn.size()) + return false; + + if (m_actnum.size() != other.m_actnum.size()) + return false; + + if (m_mapaxes.size() != other.m_mapaxes.size()) + return false; + + if (m_actnum != other.m_actnum) + return false; + + if (m_coord != other.m_coord) + return false; + + if (m_zcorn != other.m_zcorn) + return false; + + if (m_mapaxes != other.m_mapaxes) + return false; + + bool status = (m_pinch.equal( other.m_pinch ) && (m_minpvMode == other.getMinpvMode())); + + if(m_minpvMode!=MinpvMode::ModeEnum::Inactive) { status = status && (m_minpvVector == other.getMinpvVector()); } + return status; } - size_t EclipseGrid::getNumActive( ) const { - return static_cast(ecl_grid_get_nactive( c_ptr() )); + return m_nactive; } bool EclipseGrid::allActive( ) const { @@ -746,205 +1407,317 @@ namespace Opm { bool EclipseGrid::cellActive( size_t globalIndex ) const { assertGlobalIndex( globalIndex ); - return ecl_grid_cell_active1( c_ptr() , static_cast(globalIndex)); + + return m_actnum[globalIndex]>0; } bool EclipseGrid::cellActive( size_t i , size_t j , size_t k ) const { assertIJK(i,j,k); - return ecl_grid_cell_active3( c_ptr() , static_cast(i),static_cast(j),static_cast(k)); + + size_t globalIndex = getGlobalIndex(i,j,k); + + return m_actnum[globalIndex]>0; } double EclipseGrid::getCellVolume(size_t globalIndex) const { - assertGlobalIndex( globalIndex ); - if (volume_cache[globalIndex] < 0.0) { - // Calculate cell volume and put it in the cache. - std::vector x(8,0); - std::vector y(8,0); - std::vector z(8,0); - for (int i=0; i < 8; i++) { - ecl_grid_get_cell_corner_xyz1(c_ptr(), static_cast(globalIndex), i, &x[i], &y[i], &z[i]); - } - volume_cache[globalIndex] = calculateCellVol(x,y,z); - } - return volume_cache[globalIndex]; - } + assertGlobalIndex( globalIndex ); + + return m_volume[globalIndex]; + } double EclipseGrid::getCellVolume(size_t i , size_t j , size_t k) const { + assertIJK(i,j,k); - return getCellVolume(getGlobalIndex(i, j, k)); + + size_t globalIndex = getGlobalIndex( i,j,k ); + + return getCellVolume(globalIndex); } - double EclipseGrid::getCellThicknes(size_t i , size_t j , size_t k) const { + double EclipseGrid::getCellThickness(size_t i , size_t j , size_t k) const { assertIJK(i,j,k); - return ecl_grid_get_cell_thickness3( c_ptr() , static_cast(i),static_cast(j),static_cast(k)); + + const std::array dims = getNXYZ(); + size_t globalIndex = i + j*dims[0] + k*dims[0]*dims[1]; + + return getCellThickness(globalIndex); } - double EclipseGrid::getCellThicknes(size_t globalIndex) const { + double EclipseGrid::getCellThickness(size_t globalIndex) const { assertGlobalIndex( globalIndex ); - return ecl_grid_get_cell_thickness1( c_ptr() , static_cast(globalIndex)); - } + return m_dz[globalIndex]; + } std::array EclipseGrid::getCellDims(size_t globalIndex) const { assertGlobalIndex( globalIndex ); - { - double dx = ecl_grid_get_cell_dx1( c_ptr() , globalIndex); - double dy = ecl_grid_get_cell_dy1( c_ptr() , globalIndex); - double dz = ecl_grid_get_cell_thickness1( c_ptr() , globalIndex); - return std::array{ {dx , dy , dz }}; - } + return std::array {{m_dx[globalIndex] , m_dy[globalIndex] , m_dz[globalIndex] }}; } std::array EclipseGrid::getCellDims(size_t i , size_t j , size_t k) const { assertIJK(i,j,k); { size_t globalIndex = getGlobalIndex( i,j,k ); - double dx = ecl_grid_get_cell_dx1( c_ptr() , globalIndex); - double dy = ecl_grid_get_cell_dy1( c_ptr() , globalIndex); - double dz = ecl_grid_get_cell_thickness1( c_ptr() , globalIndex); - return std::array{ {dx , dy , dz }}; + return getCellDims(globalIndex); } } - std::array EclipseGrid::getCellCenter(size_t globalIndex) const { + const std::array& EclipseGrid::getCellCenter(size_t globalIndex) const { assertGlobalIndex( globalIndex ); - { - double x,y,z; - ecl_grid_get_xyz1( c_ptr() , static_cast(globalIndex) , &x , &y , &z); - return std::array{{x,y,z}}; - } + + return m_cellCenter[globalIndex]; + } + + const std::array& EclipseGrid::getCellCenter(size_t i,size_t j, size_t k) const { + assertIJK(i,j,k); + + const std::array dims = getNXYZ(); + + size_t globalIndex = i + j*dims[0] + k*dims[0]*dims[1]; + + return getCellCenter(globalIndex); } /* This is the numbering of the corners in the cell. - j + bottom j 6---7 /|\ | | | 4---5 | | - o----------> i + top o----------> i 2---3 | | 0---1 */ - std::array EclipseGrid::getCornerPos(size_t i,size_t j, size_t k, size_t corner_index) const { assertIJK(i,j,k); if (corner_index >= 8) throw std::invalid_argument("Invalid corner position"); { - double x,y,z; - ecl_grid_get_cell_corner_xyz3( c_ptr() , - static_cast(i), - static_cast(j), - static_cast(k), - static_cast(corner_index), - &x , &y , &z); - return std::array{{x,y,z}}; - } - } + const std::array dims = getNXYZ(); + std::array X = {0.0}; + std::array Y = {0.0}; + std::array Z = {0.0}; - std::array EclipseGrid::getCellCenter(size_t i,size_t j, size_t k) const { - assertIJK(i,j,k); - { - double x,y,z; - ecl_grid_get_xyz3( c_ptr() , static_cast(i),static_cast(j),static_cast(k), &x , &y , &z); - return std::array{{x,y,z}}; + std::array ijk; + + ijk[0] = i; + ijk[1] = j; + ijk[2] = k; + + getCellCorners(ijk, dims, X, Y, Z ); + + return std::array {{X[corner_index], Y[corner_index], Z[corner_index]}}; } } double EclipseGrid::getCellDepth(size_t globalIndex) const { assertGlobalIndex( globalIndex ); - return ecl_grid_get_cdepth1( c_ptr() , static_cast(globalIndex)); + + return m_depth[globalIndex]; } - - double EclipseGrid::getCellDepth(size_t i,size_t j, size_t k) const { + double EclipseGrid::getCellDepth(size_t i, size_t j, size_t k) const { assertIJK(i,j,k); - return ecl_grid_get_cdepth3( c_ptr() , static_cast(i),static_cast(j),static_cast(k)); + + size_t globalIndex = getGlobalIndex( i,j,k ); + + return getCellDepth(globalIndex); } + const std::vector& EclipseGrid::getACTNUM( ) const { - - void EclipseGrid::exportACTNUM( std::vector& actnum) const { - size_t volume = getNX() * getNY() * getNZ(); - if (getNumActive() == volume) - actnum.resize(0); - else { - actnum.resize( volume ); - ecl_grid_init_actnum_data( c_ptr() , actnum.data() ); - } + return m_actnum; } - void EclipseGrid::exportMAPAXES( std::vector& mapaxes) const { - if (ecl_grid_use_mapaxes( c_ptr())) { - mapaxes.resize(6); - ecl_grid_init_mapaxes_data_double( c_ptr() , mapaxes.data() ); - } else { - mapaxes.resize(0); - } + const std::vector& EclipseGrid::getMAPAXES( ) const { + + return m_mapaxes; } - void EclipseGrid::exportCOORD( std::vector& coord) const { - coord.resize( ecl_grid_get_coord_size( c_ptr() )); - ecl_grid_init_coord_data_double( c_ptr() , coord.data() ); + const std::vector& EclipseGrid::getCOORD() const { + + return m_coord; } - size_t EclipseGrid::exportZCORN( std::vector& zcorn) const { + size_t EclipseGrid::fixupZCORN() { + ZcornMapper mapper( getNX(), getNY(), getNZ()); - zcorn.resize( ecl_grid_get_zcorn_size( c_ptr() )); - ecl_grid_init_zcorn_data_double( c_ptr() , zcorn.data() ); - - return mapper.fixupZCORN( zcorn ); + return mapper.fixupZCORN( m_zcorn ); } + const std::vector& EclipseGrid::getZCORN( ) const { - void EclipseGrid::addNNC(const NNC& nnc) { - int idx = 0; - auto* ecl_grid = const_cast< ecl_grid_type* >( this->c_ptr() ); - for (const NNCdata& n : nnc.nncdata()) - ecl_grid_add_self_nnc( ecl_grid, n.cell1, n.cell2, idx++); + return m_zcorn; } + void EclipseGrid::save(const std::string& filename, bool formatted, const Opm::NNC& nnc, const Opm::UnitSystem& units) const { - void EclipseGrid::save(const std::string& filename, UnitSystem::UnitType output_units) const { - auto* ecl_grid = const_cast< ecl_grid_type* >( this->c_ptr() ); - ecl_grid_fwrite_EGRID2(ecl_grid, filename.c_str(), UnitSystem::ecl_units(output_units)); - } + Opm::UnitSystem::UnitType unitSystemType = units.getType(); + const auto length = ::Opm::UnitSystem::measure::length; + const std::array dims = getNXYZ(); - - const std::vector& EclipseGrid::getActiveMap() const { - if( !this->activeMap.empty() ) return this->activeMap; - - this->activeMap.resize( this->getNumActive() ); - const auto size = int(this->getCartesianSize()); - - for( int global_index = 0; global_index < size; global_index++) { - // Using the low level C function to get the active index, because the C++ - // version will throw for inactive cells. - int active_index = ecl_grid_get_active_index1( m_grid.get() , global_index ); - if (active_index >= 0) - this->activeMap[ active_index ] = global_index; + // Preparing vectors to be saved + + // create coord vector of floats with input units, converted from SI + std::vector coord_f; + coord_f.resize(m_coord.size()); + + for (size_t n=0; n< m_coord.size(); n++){ + coord_f[n] = static_cast(units.from_si(length, m_coord[n])); + } + + // create zcorn vector of floats with input units, converted from SI + std::vector zcorn_f; + zcorn_f.resize(m_zcorn.size()); + + for (size_t n=0; n< m_zcorn.size(); n++){ + zcorn_f[n] = static_cast(units.from_si(length, m_zcorn[n])); } - return this->activeMap; + std::vector mapaxes_f; + + std::vector filehead(100,0); + filehead[0] = 3; // version number + filehead[1] = 2007; // release year + filehead[6] = 1; // corner point grid + + std::vector gridhead(100,0); + gridhead[0] = 1; // corner point grid + gridhead[1] = dims[0]; // nI + gridhead[2] = dims[1]; // nJ + gridhead[3] = dims[2]; // nK + gridhead[24] = 1; // corner point grid + + std::vector nnchead(10, 0); + std::vector nnc1; + std::vector nnc2; + + for (const NNCdata& n : nnc.nncdata() ) { + nnc1.push_back(n.cell1 + 1); + nnc2.push_back(n.cell2 + 1); + } + + nnchead[0] = nnc1.size(); + + std::vector gridunits; + + switch (unitSystemType) { + case Opm::UnitSystem::UnitType::UNIT_TYPE_METRIC: + gridunits.push_back("METRES"); + break; + case Opm::UnitSystem::UnitType::UNIT_TYPE_FIELD: + gridunits.push_back("FEET"); + break; + case Opm::UnitSystem::UnitType::UNIT_TYPE_LAB: + gridunits.push_back("CM"); + break; + default: + OPM_THROW(std::runtime_error, "Unit system not supported when writing to EGRID file"); + break; + } + + gridunits.push_back(""); + + // map units is not dependent on deck units. A user may specify FIELD units for the model + // and metric units for the MAPAXES keyword (MAPUNITS) + + std::vector mapunits; + + if ((m_mapunits.size() > 0) && (m_mapaxes.size() > 0)) { + mapunits.push_back(m_mapunits); + } + + if (m_mapaxes.size() > 0){ + for (double dv : m_mapaxes){ + mapaxes_f.push_back(static_cast(dv)); + } + } + + std::vector endgrid = {}; + + // Writing vectors to egrid file + + Opm::EclIO::EclOutput egridfile(filename, formatted); + egridfile.write("FILEHEAD", filehead); + + if (mapunits.size() > 0) { + egridfile.write("MAPUNITS", mapunits); + } + + if (mapaxes_f.size() > 0){ + egridfile.write("MAPAXES", mapaxes_f); + } + + egridfile.write("GRIDUNIT", gridunits); + egridfile.write("GRIDHEAD", gridhead); + + egridfile.write("COORD", coord_f); + egridfile.write("ZCORN", zcorn_f); + + egridfile.write("ACTNUM", m_actnum); + egridfile.write("ENDGRID", endgrid); + + if (nnc1.size() > 0){ + egridfile.write("NNCHEAD", nnchead); + egridfile.write("NNC1", nnc1); + egridfile.write("NNC2", nnc2); + } } - void EclipseGrid::resetACTNUM( const int * actnum) { - ecl_grid_reset_actnum( m_grid.get() , actnum ); - /* re-build the active map cache */ - this->activeMap.clear(); - this->getActiveMap(); + const std::vector& EclipseGrid::getActiveMap() const { + + return m_active_to_global; + } + + void EclipseGrid::resetACTNUM() { + + const std::array dims = getNXYZ(); + m_nactive = dims[0]*dims[1]*dims[2]; + + std::vector all_active(m_nactive, 1); + + m_actnum.clear(); + m_actnum = all_active; + + m_global_to_active = all_active; + std::iota(m_global_to_active.begin(), m_global_to_active.end(), 0); + + m_active_to_global = m_global_to_active; + } + + void EclipseGrid::resetACTNUM( const std::vector& actnum) { + + if (actnum.size() != getCartesianSize()) { + throw std::runtime_error("resetACTNUM(): actnum vector size differs from logical cartesian size of grid."); + } + + m_actnum = actnum; + m_global_to_active.clear(); + m_active_to_global.clear(); + m_global_to_active.reserve(actnum.size()); + m_nactive = 0; + + for (size_t n = 0; n < actnum.size() ; n++) { + if (actnum[n] > 0) { + m_global_to_active.push_back(m_nactive); + m_active_to_global.push_back(n); + ++m_nactive; + } else { + m_global_to_active.push_back(-1); + } + } } ZcornMapper EclipseGrid::zcornMapper() const { @@ -1015,9 +1788,7 @@ namespace Opm { return true; } - - - + size_t ZcornMapper::fixupZCORN( std::vector& zcorn) { int sign = zcorn[ this->index(0,0,0,0) ] <= zcorn[this->index(0,0, this->dims[2] - 1,4)] ? 1 : -1; size_t cells_adjusted = 0; @@ -1037,7 +1808,6 @@ namespace Opm { } } - /* Cell internal */ { size_t index1 = this->index(i,j,k,c); @@ -1052,14 +1822,12 @@ namespace Opm { return cells_adjusted; } - CoordMapper::CoordMapper(size_t nx_, size_t ny_) : nx(nx_), ny(ny_) { } - size_t CoordMapper::size() const { return (this->nx + 1) * (this->ny + 1) * 6; } @@ -1083,3 +1851,4 @@ namespace Opm { } + diff --git a/tests/parser/ConnectionTests.cpp b/tests/parser/ConnectionTests.cpp index 545f54e59..f93fe0e91 100644 --- a/tests/parser/ConnectionTests.cpp +++ b/tests/parser/ConnectionTests.cpp @@ -134,7 +134,7 @@ BOOST_AUTO_TEST_CASE(ActiveCompletions) { std::vector actnum(grid.getCartesianSize(), 1); actnum[0] = 0; - grid.resetACTNUM( actnum.data() ); + grid.resetACTNUM( actnum); Opm::WellConnections active_completions(completions, grid); BOOST_CHECK_EQUAL( active_completions.size() , 2); diff --git a/tests/parser/EclipseGridTests.cpp b/tests/parser/EclipseGridTests.cpp index cfba38ae2..07490101c 100644 --- a/tests/parser/EclipseGridTests.cpp +++ b/tests/parser/EclipseGridTests.cpp @@ -22,13 +22,15 @@ #include #include #include - +#include +#include +#include #define BOOST_TEST_MODULE EclipseGridTests #include #include -#include +//#include #include #include @@ -40,9 +42,11 @@ #include #include #include +#include #include +#include BOOST_AUTO_TEST_CASE(CreateMissingDIMENS_throws) { Opm::Deck deck; @@ -115,14 +119,12 @@ BOOST_AUTO_TEST_CASE(MissingDimsThrows) { BOOST_CHECK_THROW( Opm::EclipseGrid{ deck }, std::invalid_argument); } - BOOST_AUTO_TEST_CASE(HasGridKeywords) { Opm::Deck deck = createDeckHeaders(); BOOST_CHECK( !Opm::EclipseGrid::hasCornerPointKeywords( deck )); BOOST_CHECK( !Opm::EclipseGrid::hasCartesianKeywords( deck )); } - BOOST_AUTO_TEST_CASE(CreateGridNoCells) { Opm::Deck deck = createDeckHeaders(); BOOST_CHECK_THROW( Opm::EclipseGrid{ deck }, std::invalid_argument); @@ -187,7 +189,6 @@ static Opm::Deck createCPDeck() { return parser.parseString( deckData) ; } - static Opm::Deck createPinchedCPDeck() { const char* deckData = "RUNSPEC\n" @@ -380,66 +381,59 @@ BOOST_AUTO_TEST_CASE(DEPTHZ_EQUAL_TOPS) { for (size_t k= 0; k < 10; k++) for (size_t j= 0; j < 10; j++) for (size_t i= 0; i < 10; i++) - BOOST_CHECK_CLOSE( grid1.getCellVolume(i,j,k) , 0.25*0.25*0.25 , 0.001 ); + BOOST_CHECK_CLOSE( grid1.getCellVolume(i, j, k) , 0.25*0.25*0.25 , 0.001 ); + } + { BOOST_CHECK_THROW( grid1.getCellCenter(1000) , std::invalid_argument); BOOST_CHECK_THROW( grid1.getCellCenter(10,0,0) , std::invalid_argument); BOOST_CHECK_THROW( grid1.getCellCenter(0,10,0) , std::invalid_argument); BOOST_CHECK_THROW( grid1.getCellCenter(0,0,10) , std::invalid_argument); - for (size_t k= 0; k < 10; k++) - for (size_t j= 0; j < 10; j++) + for (size_t k= 0; k < 10; k++) { + for (size_t j= 0; j < 10; j++) { for (size_t i= 0; i < 10; i++) { - auto pos = grid1.getCellCenter(i,j,k); + auto pos = grid1.getCellCenter(i, j, k); BOOST_CHECK_CLOSE( std::get<0>(pos) , i*0.25 + 0.125, 0.001); BOOST_CHECK_CLOSE( std::get<1>(pos) , j*0.25 + 0.125, 0.001); BOOST_CHECK_CLOSE( std::get<2>(pos) , k*0.25 + 0.125 + 0.25, 0.001); - } + } + } } } - - BOOST_AUTO_TEST_CASE(HasCPKeywords) { Opm::Deck deck = createCPDeck(); BOOST_CHECK( Opm::EclipseGrid::hasCornerPointKeywords( deck )); BOOST_CHECK( !Opm::EclipseGrid::hasCartesianKeywords( deck )); } - BOOST_AUTO_TEST_CASE(HasCartKeywords) { Opm::Deck deck = createCARTDeck(); BOOST_CHECK( !Opm::EclipseGrid::hasCornerPointKeywords( deck )); BOOST_CHECK( Opm::EclipseGrid::hasCartesianKeywords( deck )); } - BOOST_AUTO_TEST_CASE(HasCartKeywordsDEPTHZ) { Opm::Deck deck = createCARTDeckDEPTHZ(); BOOST_CHECK( !Opm::EclipseGrid::hasCornerPointKeywords( deck )); BOOST_CHECK( Opm::EclipseGrid::hasCartesianKeywords( deck )); } - BOOST_AUTO_TEST_CASE(HasINVALIDCartKeywords) { Opm::Deck deck = createCARTInvalidDeck(); BOOST_CHECK( !Opm::EclipseGrid::hasCornerPointKeywords( deck )); BOOST_CHECK( !Opm::EclipseGrid::hasCartesianKeywords( deck )); } - - - - BOOST_AUTO_TEST_CASE(CreateMissingGRID_throws) { auto deck= createDeckHeaders(); BOOST_CHECK_THROW( Opm::EclipseGrid{ deck }, std::invalid_argument); } - static Opm::Deck createInvalidDXYZCARTDeck() { const char* deckData = "RUNSPEC\n" @@ -462,14 +456,11 @@ static Opm::Deck createInvalidDXYZCARTDeck() { return parser.parseString( deckData) ; } - - BOOST_AUTO_TEST_CASE(CreateCartesianGRID) { auto deck = createInvalidDXYZCARTDeck(); BOOST_CHECK_THROW( Opm::EclipseGrid{ deck }, std::invalid_argument); } - static Opm::Deck createInvalidDXYZCARTDeckDEPTHZ() { const char* deckData = "RUNSPEC\n" @@ -492,14 +483,11 @@ static Opm::Deck createInvalidDXYZCARTDeckDEPTHZ() { return parser.parseString( deckData) ; } - - BOOST_AUTO_TEST_CASE(CreateCartesianGRIDDEPTHZ) { auto deck = createInvalidDXYZCARTDeckDEPTHZ(); BOOST_CHECK_THROW( Opm::EclipseGrid{ deck }, std::invalid_argument); } - static Opm::Deck createOnlyTopDZCartGrid() { const char* deckData = "RUNSPEC\n" @@ -545,13 +533,11 @@ static Opm::Deck createInvalidDEPTHZDeck1 () { return parser.parseString( deckData) ; } - BOOST_AUTO_TEST_CASE(CreateCartesianGRIDInvalidDEPTHZ1) { auto deck = createInvalidDEPTHZDeck1(); BOOST_CHECK_THROW( Opm::EclipseGrid{ deck }, std::invalid_argument); } - static Opm::Deck createInvalidDEPTHZDeck2 () { const char* deckData = "RUNSPEC\n" @@ -579,31 +565,24 @@ BOOST_AUTO_TEST_CASE(CreateCartesianGRIDInvalidDEPTHZ2) { BOOST_CHECK_THROW( Opm::EclipseGrid{ deck }, std::invalid_argument); } - - BOOST_AUTO_TEST_CASE(CreateCartesianGRIDOnlyTopLayerDZ) { Opm::Deck deck = createOnlyTopDZCartGrid(); Opm::EclipseGrid grid( deck ); - BOOST_CHECK_EQUAL( 10 , grid.getNX( )); BOOST_CHECK_EQUAL( 5 , grid.getNY( )); BOOST_CHECK_EQUAL( 20 , grid.getNZ( )); BOOST_CHECK_EQUAL( 1000 , grid.getNumActive()); } - - BOOST_AUTO_TEST_CASE(AllActiveExportActnum) { Opm::Deck deck = createOnlyTopDZCartGrid(); Opm::EclipseGrid grid( deck ); - std::vector actnum( 1, 100 ); - - grid.exportACTNUM( actnum ); - BOOST_CHECK_EQUAL( 0U , actnum.size()); + std::vector actnum = grid.getACTNUM(); + + BOOST_CHECK_EQUAL( 1000 , actnum.size()); } - BOOST_AUTO_TEST_CASE(CornerPointSizeMismatchCOORD) { const char* deckData = "RUNSPEC\n" @@ -650,7 +629,6 @@ BOOST_AUTO_TEST_CASE(CornerPointSizeMismatchZCORN) { BOOST_CHECK_THROW(Opm::EclipseGrid{ deck }, std::invalid_argument); } - BOOST_AUTO_TEST_CASE(ResetACTNUM) { const char* deckData = "RUNSPEC\n" @@ -675,19 +653,23 @@ BOOST_AUTO_TEST_CASE(ResetACTNUM) { actnum[2] = 1; actnum[4] = 1; actnum[6] = 1; - grid.resetACTNUM( actnum.data() ); + + grid.resetACTNUM( actnum ); + BOOST_CHECK_EQUAL( 4U , grid.getNumActive() ); { std::vector full(grid.getCartesianSize()); std::iota(full.begin(), full.end(), 0); auto compressed = grid.compressedVector( full ); + BOOST_CHECK_EQUAL( compressed.size() , 4U ); BOOST_CHECK_EQUAL( compressed[0] , 0 ); BOOST_CHECK_EQUAL( compressed[1] , 2 ); BOOST_CHECK_EQUAL( compressed[2] , 4 ); BOOST_CHECK_EQUAL( compressed[3] , 6 ); } + { const auto& activeMap = grid.getActiveMap( ); BOOST_CHECK_EQUAL( 4U , activeMap.size() ); @@ -697,9 +679,9 @@ BOOST_AUTO_TEST_CASE(ResetACTNUM) { BOOST_CHECK_EQUAL( 6 , activeMap[3] ); } - grid.resetACTNUM( NULL ); - BOOST_CHECK_EQUAL( 1000U , grid.getNumActive() ); + grid.resetACTNUM(); + BOOST_CHECK_EQUAL( 1000U , grid.getNumActive() ); { const auto& activeMap = grid.getActiveMap( ); BOOST_CHECK_EQUAL( 1000U , activeMap.size() ); @@ -708,6 +690,76 @@ BOOST_AUTO_TEST_CASE(ResetACTNUM) { BOOST_CHECK_EQUAL( 2 , activeMap[2] ); BOOST_CHECK_EQUAL( 999 , activeMap[999] ); } + + actnum.assign(1000, 1); + + actnum[0] = 0; + actnum[1] = 0; + actnum[2] = 0; + actnum[11] = 0; + actnum[21] = 0; + actnum[430] = 0; + actnum[431] = 0; + + grid.resetACTNUM( actnum ); + + std::vector actMap = grid.getActiveMap(); + + BOOST_CHECK_EQUAL(actMap.size(), 993); + BOOST_CHECK_THROW(grid.getGlobalIndex(993), std::out_of_range); + BOOST_CHECK_EQUAL(grid.getGlobalIndex(0), 3); + BOOST_CHECK_EQUAL(grid.getGlobalIndex(33), 38); + BOOST_CHECK_EQUAL(grid.getGlobalIndex(450), 457); + BOOST_CHECK_EQUAL(grid.getGlobalIndex(1,2,3), 321); +} + +BOOST_AUTO_TEST_CASE(TestCP_example) { + const char* deckData = + + "RUNSPEC\n" + "\n" + "DIMENS\n" + " 3 2 1 /\n" + "GRID\n" + "COORD\n" + " 2000.0000 2000.0000 2000.0000 1999.9476 2000.0000 2002.9995\n" + " 2049.9924 2000.0000 2000.8726 2049.9400 2000.0000 2003.8722 \n" + " 2099.9848 2000.0000 2001.7452 2099.9324 2000.0000 2004.7448 \n" + " 2149.9772 2000.0000 2002.6179 2149.9248 2000.0000 2005.6174 \n" + " 2000.0000 2050.0000 2000.0000 1999.9476 2050.0000 2002.9995 \n" + " 2049.9924 2050.0000 2000.8726 2049.9400 2050.0000 2003.8722 \n" + " 2099.9848 2050.0000 2001.7452 2099.9324 2050.0000 2004.7448 \n" + " 2149.9772 2050.0000 2002.6179 2149.9248 2050.0000 2005.6174 \n" + " 2000.0000 2100.0000 2000.0000 1999.9476 2100.0000 2002.9995 \n" + " 2049.9924 2100.0000 2000.8726 2049.9400 2100.0000 2003.8722 \n" + " 2099.9848 2100.0000 2001.7452 2099.9324 2100.0000 2004.7448 \n" + " 2149.9772 2100.0000 2002.6179 2149.9248 2100.0000 2005.6174 / \n" + "ZCORN\n" + " 2000.0000 2000.8726 2000.8726 2001.7452 2001.7452 2002.6179 \n" + " 2000.0000 2000.8726 2000.8726 2001.7452 2001.7452 2002.6179 \n" + " 2000.0000 2000.8726 2000.8726 2001.7452 2001.7452 2002.6179 \n" + " 2000.0000 2000.8726 2000.8726 2001.7452 2001.7452 2002.6179 \n" + " 2002.9995 2003.8722 2003.8722 2004.7448 2004.7448 2005.6174 \n" + " 2002.9995 2003.8722 2003.8722 2004.7448 2004.7448 2005.6174 \n" + " 2002.9995 2003.8722 2003.8722 2004.7448 2004.7448 2005.6174 \n" + " 2002.9995 2003.8722 2003.8722 2004.7448 2004.7448 2005.6174 / \n" + "EDIT\n" + "\n"; + + Opm::Parser parser; + auto deck = parser.parseString( deckData) ; + + Opm::EclipseGrid grid( deck); + BOOST_CHECK_EQUAL( 6U , grid.getNumActive()); + + std::vector actnum(6, 0); + actnum[0] = 1; + actnum[2] = 1; + actnum[4] = 1; + + grid.resetACTNUM( actnum ); + + BOOST_CHECK_EQUAL( 3U , grid.getNumActive() ); } @@ -756,15 +808,6 @@ BOOST_AUTO_TEST_CASE(ACTNUM_BEST_EFFORT) { -BOOST_AUTO_TEST_CASE(LoadFromBinary) { - BOOST_CHECK_THROW(Opm::EclipseGrid( "No/does/not/exist" ) , std::invalid_argument); -} - - - - - - BOOST_AUTO_TEST_CASE(ConstructorNORUNSPEC) { const char* deckData = "GRID\n" @@ -789,30 +832,6 @@ BOOST_AUTO_TEST_CASE(ConstructorNORUNSPEC) { BOOST_CHECK(grid1.equal( grid2 )); } -BOOST_AUTO_TEST_CASE(GDFILE) { - const char* gdfile_deck = - "GRID\n" - "GDFILE\n" - "'grid/CASE.EGRID' /\n" - "\n"; - - Opm::EclipseGrid grid1(createCPDeck()); - test_work_area_type * work_area = test_work_area_alloc("GDFILE"); - util_mkdir_p("ecl/grid"); - grid1.save("ecl/grid/CASE.EGRID", Opm::UnitSystem::UnitType::UNIT_TYPE_METRIC); - { - FILE * stream = fopen("ecl/DECK.DATA", "w"); - fputs(gdfile_deck, stream); - fclose(stream); - } - { - Opm::Parser parser; - Opm::EclipseGrid grid2(parser.parseFile("ecl/DECK.DATA" )); - BOOST_CHECK(grid1.equal(grid2)); - } - test_work_area_free( work_area ); -} - BOOST_AUTO_TEST_CASE(ConstructorNoSections) { const char* deckData = "DIMENS \n" @@ -835,8 +854,6 @@ BOOST_AUTO_TEST_CASE(ConstructorNoSections) { BOOST_CHECK(grid1.equal( grid2 )); } - - BOOST_AUTO_TEST_CASE(ConstructorNORUNSPEC_PINCH) { auto deck1 = createCPDeck(); auto deck2 = createPinchedCPDeck(); @@ -845,15 +862,13 @@ BOOST_AUTO_TEST_CASE(ConstructorNORUNSPEC_PINCH) { Opm::EclipseGrid grid2(deck2); BOOST_CHECK(!grid1.equal( grid2 )); + BOOST_CHECK(!grid1.isPinchActive()); BOOST_CHECK_THROW(grid1.getPinchThresholdThickness(), std::logic_error); BOOST_CHECK(grid2.isPinchActive()); BOOST_CHECK_EQUAL(grid2.getPinchThresholdThickness(), 0.2); } - - - BOOST_AUTO_TEST_CASE(ConstructorMINPV) { auto deck1 = createCPDeck(); auto deck2 = createMinpvDefaultCPDeck(); @@ -873,7 +888,6 @@ BOOST_AUTO_TEST_CASE(ConstructorMINPV) { BOOST_CHECK_EQUAL(grid4.getMinpvVector()[0], 20.0); } - static Opm::Deck createActnumDeck() { const char* deckData = "RUNSPEC\n" "\n" @@ -991,6 +1005,7 @@ BOOST_AUTO_TEST_CASE(GridBoxActnum) { } } + BOOST_AUTO_TEST_CASE(GridActnumVia3D) { auto deck = createActnumDeck(); @@ -999,14 +1014,20 @@ BOOST_AUTO_TEST_CASE(GridActnumVia3D) { const auto& grid = es.getInputGrid(); Opm::EclipseGrid grid2( grid ); + std::vector actnum = {1, 1, 0, 1, 1, 0, 1, 1}; + Opm::EclipseGrid grid3( grid , actnum); + BOOST_CHECK_NO_THROW(ep.getIntGridProperty("ACTNUM")); BOOST_CHECK_NO_THROW(grid.getNumActive()); BOOST_CHECK_EQUAL(grid.getNumActive(), 2 * 2 * 2 - 1); BOOST_CHECK_NO_THROW(grid2.getNumActive()); BOOST_CHECK_EQUAL(grid2.getNumActive(), 2 * 2 * 2 - 1); + + BOOST_CHECK_EQUAL(grid3.getNumActive(), 6); } + BOOST_AUTO_TEST_CASE(GridActnumViaState) { auto deck = createActnumDeck(); @@ -1024,6 +1045,7 @@ BOOST_AUTO_TEST_CASE(GridDimsSPECGRID) { BOOST_CHECK_EQUAL(gd.getNZ(), 19); } + BOOST_AUTO_TEST_CASE(GridDimsDIMENS) { auto deck = createDeckDIMENS(); auto gd = Opm::GridDims( deck ); @@ -1038,17 +1060,19 @@ BOOST_AUTO_TEST_CASE(ProcessedCopy) { std::vector zcorn; std::vector actnum; - gd.exportZCORN( zcorn ); - gd.exportACTNUM( actnum ); + zcorn = gd.getZCORN(); + actnum = gd.getACTNUM(); + Opm::EclipseGrid gd1(gd , actnum ); + BOOST_CHECK( gd.equal( gd1 )); { - Opm::EclipseGrid gd2(gd , zcorn , actnum ); + Opm::EclipseGrid gd2(gd , zcorn.data() , actnum ); BOOST_CHECK( gd.equal( gd2 )); } - zcorn[0] -= 1; + zcorn[0] -= 1.0; { - Opm::EclipseGrid gd2(gd , zcorn , actnum ); + Opm::EclipseGrid gd2(gd , zcorn.data() , actnum ); BOOST_CHECK( !gd.equal( gd2 )); } @@ -1066,50 +1090,101 @@ BOOST_AUTO_TEST_CASE(ProcessedCopy) { } } +BOOST_AUTO_TEST_CASE(regularCartGrid) { -BOOST_AUTO_TEST_CASE(ZcornMapper) { int nx = 3; int ny = 4; int nz = 5; - Opm::EclipseGrid grid(nx,ny,nz); + + double dx = 25; + double dy = 35; + double dz = 2; + + double ref_volume = dx* dy* dz; + + Opm::EclipseGrid grid(nx, ny, nz, dx, dy, dz); + + std::array dims = grid.getNXYZ(); + + int nCells = dims[0]*dims[1]*dims[2]; + + for (int n=0; n cc = grid.getCellCenter(i, j, k); + BOOST_CHECK_CLOSE(cc[0], ref_x, 1e-12); + BOOST_CHECK_CLOSE(cc[1], ref_y, 1e-12); + BOOST_CHECK_CLOSE(cc[2], ref_z, 1e-12); + } + } + } +} + +BOOST_AUTO_TEST_CASE(ZcornMapper) { + + int nx = 3; + int ny = 4; + int nz = 5; + + Opm::EclipseGrid grid(nx, ny, nz); Opm::ZcornMapper zmp = grid.zcornMapper( ); - const ecl_grid_type * ert_grid = grid.c_ptr(); - - + BOOST_CHECK_THROW(zmp.index(nx,1,1,0) , std::invalid_argument); BOOST_CHECK_THROW(zmp.index(0,ny,1,0) , std::invalid_argument); BOOST_CHECK_THROW(zmp.index(0,1,nz,0) , std::invalid_argument); BOOST_CHECK_THROW(zmp.index(0,1,2,8) , std::invalid_argument); - for (int k=0; k < nz; k++) - for (int j=0; j < ny; j++) - for (int i=0; i < nx; i++) - for (int c=0; c < 8; c++) { - size_t g = i + j*nx + k*nx*ny; - BOOST_CHECK_EQUAL( zmp.index(g , c) , zmp.index( i,j,k,c)); - BOOST_CHECK_EQUAL( zmp.index(i,j,k,c) , ecl_grid_zcorn_index( ert_grid, i , j , k, c)); - } + auto points_adjusted = grid.fixupZCORN(); - std::vector zcorn; - auto points_adjusted = grid.exportZCORN( zcorn ); + std::vector actnum = grid.getACTNUM(); + std::vector zcorn = grid.getZCORN(); + + zcorn[42] = zcorn[42] + 2.0; + zcorn[96] = zcorn[96] + 2.0; + + Opm::EclipseGrid grid2(grid , zcorn.data() , actnum ); + points_adjusted = grid2.getZcornFixed(); + BOOST_CHECK_EQUAL( points_adjusted , 4 ); + + points_adjusted = grid2.fixupZCORN(); BOOST_CHECK_EQUAL( points_adjusted , 0 ); + + zcorn = grid.getZCORN(); + BOOST_CHECK( zmp.validZCORN( zcorn )); - /* Manually destroy it - cell internal */ + // Manually destroy it - cell internal zcorn[ zmp.index(0,0,0,4) ] = zcorn[ zmp.index(0,0,0,0) ] - 0.1; BOOST_CHECK( !zmp.validZCORN( zcorn )); points_adjusted = zmp.fixupZCORN( zcorn ); BOOST_CHECK_EQUAL( points_adjusted , 1 ); BOOST_CHECK( zmp.validZCORN( zcorn )); - /* Manually destroy it - cell 2 cell */ + // Manually destroy it - cell 2 cell zcorn[ zmp.index(0,0,0,4) ] = zcorn[ zmp.index(0,0,1,0) ] + 0.1; BOOST_CHECK( !zmp.validZCORN( zcorn )); points_adjusted = zmp.fixupZCORN( zcorn ); BOOST_CHECK_EQUAL( points_adjusted , 1 ); BOOST_CHECK( zmp.validZCORN( zcorn )); - /* Manually destroy it - cell 2 cell and cell internal*/ + // Manually destroy it - cell 2 cell and cell internal zcorn[ zmp.index(0,0,0,4) ] = zcorn[ zmp.index(0,0,1,0) ] + 0.1; zcorn[ zmp.index(0,0,0,0) ] = zcorn[ zmp.index(0,0,0,4) ] + 0.1; BOOST_CHECK( !zmp.validZCORN( zcorn )); @@ -1118,8 +1193,6 @@ BOOST_AUTO_TEST_CASE(ZcornMapper) { BOOST_CHECK( zmp.validZCORN( zcorn )); } - - BOOST_AUTO_TEST_CASE(MoveTest) { int nx = 3; int ny = 4; @@ -1127,13 +1200,9 @@ BOOST_AUTO_TEST_CASE(MoveTest) { Opm::EclipseGrid grid1(nx,ny,nz); Opm::EclipseGrid grid2( std::move( grid1 )); // grid2 should be move constructed from grid1 - BOOST_CHECK( !grid1.c_ptr() ); // We peek at some internal details ... BOOST_CHECK( !grid1.circle( )); } - - - static Opm::Deck radial_missing_INRAD() { const char* deckData = "RUNSPEC\n" @@ -1205,7 +1274,6 @@ BOOST_AUTO_TEST_CASE(RadialTest) { BOOST_CHECK_THROW( Opm::EclipseGrid{ deck }, std::invalid_argument); } - BOOST_AUTO_TEST_CASE(RadialKeywordsOK) { Opm::Deck deck = radial_keywords_OK(); Opm::EclipseGrid grid( deck ); @@ -1218,8 +1286,6 @@ BOOST_AUTO_TEST_CASE(RadialKeywordsOK_CIRCLE) { BOOST_CHECK(grid.circle()); } - - static Opm::Deck radial_keywords_DRV_size_mismatch() { const char* deckData = "RUNSPEC\n" @@ -1293,10 +1359,8 @@ static Opm::Deck radial_keywords_DTHETAV_size_mismatch() { return parser.parseString( deckData); } -/* - This is stricter than the ECLIPSE implementation; we assume that - *only* the top layer is explicitly given. -*/ +// This is stricter than the ECLIPSE implementation; we assume that +// *only* the top layer is explicitly given. static Opm::Deck radial_keywords_TOPS_size_mismatch() { const char* deckData = @@ -1322,7 +1386,6 @@ static Opm::Deck radial_keywords_TOPS_size_mismatch() { return parser.parseString( deckData); } - static Opm::Deck radial_keywords_ANGLE_OVERFLOW() { const char* deckData = "RUNSPEC\n" @@ -1348,8 +1411,6 @@ static Opm::Deck radial_keywords_ANGLE_OVERFLOW() { } - - BOOST_AUTO_TEST_CASE(RadialKeywords_SIZE_ERROR) { BOOST_CHECK_THROW( Opm::EclipseGrid{ radial_keywords_DRV_size_mismatch() } , std::invalid_argument); BOOST_CHECK_THROW( Opm::EclipseGrid{ radial_keywords_DZV_size_mismatch() } , std::invalid_argument); @@ -1358,8 +1419,6 @@ BOOST_AUTO_TEST_CASE(RadialKeywords_SIZE_ERROR) { BOOST_CHECK_THROW( Opm::EclipseGrid{ radial_keywords_ANGLE_OVERFLOW() } , std::invalid_argument); } - - static Opm::Deck radial_details() { const char* deckData = "RUNSPEC\n" @@ -1384,8 +1443,6 @@ static Opm::Deck radial_details() { return parser.parseString( deckData); } - - BOOST_AUTO_TEST_CASE(RadialDetails) { Opm::Deck deck = radial_details(); Opm::EclipseGrid grid( deck ); @@ -1395,7 +1452,6 @@ BOOST_AUTO_TEST_CASE(RadialDetails) { auto pos0 = grid.getCellCenter(0,0,0); auto pos2 = grid.getCellCenter(0,2,0); - BOOST_CHECK_CLOSE( std::get<0>(pos0) , 0.75 , 0.0001); BOOST_CHECK_CLOSE( std::get<1>(pos0) , 0.75 , 0.0001); BOOST_CHECK_CLOSE( std::get<2>(pos0) , 1.50 , 0.0001); @@ -1407,7 +1463,6 @@ BOOST_AUTO_TEST_CASE(RadialDetails) { { const auto& p0 = grid.getCornerPos( 0,0,0 , 0 ); const auto& p6 = grid.getCornerPos( 0,0,0 , 6 ); - BOOST_CHECK_CLOSE( p0[0]*p0[0] + p0[1]*p0[1] , 1.0, 0.0001); BOOST_CHECK_CLOSE( p6[0]*p6[0] + p6[1]*p6[1] , 1.0, 0.0001); @@ -1415,7 +1470,6 @@ BOOST_AUTO_TEST_CASE(RadialDetails) { } } - BOOST_AUTO_TEST_CASE(CoordMapper) { size_t nx = 10; size_t ny = 7; @@ -1427,3 +1481,1333 @@ BOOST_AUTO_TEST_CASE(CoordMapper) { BOOST_CHECK_EQUAL( cmp.index(10,7,2,1) + 1 , cmp.size( )); } + +static Opm::Deck createCARTDeckTest3x4x2() { + const char* deckData = + "RUNSPEC\n" + "\n" + "DIMENS\n" + " 3 4 2 /\n" + "GRID\n" + "DX\n" + "100 120 110 100 120 110 100 120 110 100 120 110 /\n" + "DY\n" + "70 80 85 80 70 80 85 80 70 80 85 80 /\n" + "DZ\n" + "12*25 12*35 /\n" + "TOPS\n" + "2500 2510 2520 2520 2530 2540 2540 2550 2560 2560 2570 2580 /\n" + "EDIT\n" + "\n"; + + Opm::Parser parser; + return parser.parseString( deckData) ; +} + +BOOST_AUTO_TEST_CASE(CART_Deck_3x4x2) { + + Opm::Deck deck1 = createCARTDeckTest3x4x2(); + Opm::EclipseGrid grid1( deck1 ); + + const std::vector t_coord = grid1.getCOORD(); + const std::vector t_zcorn = grid1.getZCORN(); + + std::vector ref_coord = {0, 0, 2500, 0, 0, 2560, 100, 0, 2510, 100, 0, 2570, 220, 0, 2520, 220, 0, + 2580, 330, 0, 2520, 330, 0, 2580, 0, 70, 2520, 0, 70, 2580, 100, 80, 2530, 100, 80, 2590, 220, 85, + 2540, 220, 85, 2600, 330, 85, 2540, 330, 85, 2600, 0, 150, 2540, 0, 150, 2600, 100, 150, 2550, 100, + 150, 2610, 220, 165, 2560, 220, 165, 2620, 330, 165, 2560, 330, 165, 2620, 0, 235, 2560, 0, 235, + 2620, 100, 230, 2570, 100, 230, 2630, 220, 235, 2580, 220, 235, 2640, 330, 235, 2580, 330, 235, + 2640, 0, 315, 2560, 0, 315, 2620, 100, 315, 2570, 100, 315, 2630, 220, 315, 2580, 220, 315, 2640, + 330, 315, 2580, 330, 315, 2640}; + + std::vector ref_zcorn = {2500, 2500, 2510, 2510, 2520, 2520, 2500, 2500, 2510, 2510, 2520, 2520, + 2520, 2520, 2530, 2530, 2540, 2540, 2520, 2520, 2530, 2530, 2540, 2540, 2540, 2540, 2550, 2550, 2560, + 2560, 2540, 2540, 2550, 2550, 2560, 2560, 2560, 2560, 2570, 2570, 2580, 2580, 2560, 2560, 2570, 2570, + 2580, 2580, 2525, 2525, 2535, 2535, 2545, 2545, 2525, 2525, 2535, 2535, 2545, 2545, 2545, 2545, 2555, + 2555, 2565, 2565, 2545, 2545, 2555, 2555, 2565, 2565, 2565, 2565, 2575, 2575, 2585, 2585, 2565, 2565, + 2575, 2575, 2585, 2585, 2585, 2585, 2595, 2595, 2605, 2605, 2585, 2585, 2595, 2595, 2605, 2605, 2525, + 2525, 2535, 2535, 2545, 2545, 2525, 2525, 2535, 2535, 2545, 2545, 2545, 2545, 2555, 2555, 2565, 2565, + 2545, 2545, 2555, 2555, 2565, 2565, 2565, 2565, 2575, 2575, 2585, 2585, 2565, 2565, 2575, 2575, 2585, + 2585, 2585, 2585, 2595, 2595, 2605, 2605, 2585, 2585, 2595, 2595, 2605, 2605, 2560, 2560, 2570, 2570, + 2580, 2580, 2560, 2560, 2570, 2570, 2580, 2580, 2580, 2580, 2590, 2590, 2600, 2600, 2580, 2580, 2590, + 2590, 2600, 2600, 2600, 2600, 2610, 2610, 2620, 2620, 2600, 2600, 2610, 2610, 2620, 2620, 2620, 2620, + 2630, 2630, 2640, 2640, 2620, 2620, 2630, 2630, 2640, 2640}; + + BOOST_CHECK_EQUAL( t_coord.size() , ref_coord.size()); + + for (size_t i=0; i< t_coord.size(); i++) { + BOOST_CHECK_CLOSE( t_coord[i] , ref_coord[i], 1.0e-5); + } + + BOOST_CHECK_EQUAL( t_zcorn.size() , ref_zcorn.size()); + + for (size_t i=0; i< t_zcorn.size(); i++) { + BOOST_CHECK_CLOSE( t_zcorn[i] , ref_zcorn[i], 1.0e-5); + } +} + +static Opm::Deck createCARTDeckDEPTHZ_2x3x2() { + const char* deckData = + "RUNSPEC\n" + "\n" + "DIMENS\n" + " 2 3 2 /\n" + "GRID\n" + "DXV\n" + "100 120 /\n" + "DYV\n" + "70 80 85 /\n" + "DZV\n" + "25 35 /\n" + "DEPTHZ\n" + "2500 2510 2520 2502 2512 2522 2504 2514 2524 2505 2515 2525 /\n" + "EDIT\n" + "\n"; + + Opm::Parser parser; + return parser.parseString( deckData) ; +} + +BOOST_AUTO_TEST_CASE(CART_Deck_DEPTHZ_2x3x2) { + + Opm::Deck deck1 = createCARTDeckDEPTHZ_2x3x2(); + Opm::EclipseGrid grid1( deck1 ); + + std::vector ref_coord = { 0, 0, 2500, 0, 0, 2560, 100, 0, 2510, 100, 0, 2570, 220, 0, 2520, 220, 0, + 2580, 0, 70, 2502, 0, 70, 2562, 100, 70, 2512, 100, 70, 2572, 220, 70, 2522, 220, 70, 2582, 0, 150, + 2504, 0, 150, 2564, 100, 150, 2514, 100, 150, 2574, 220, 150, 2524, 220, 150, 2584, 0, 235, 2505, 0, + 235, 2565, 100, 235, 2515, 100, 235, 2575, 220, 235, 2525, 220, 235, 2585 }; + + std::vector ref_zcorn = { 2500, 2510, 2510, 2520, 2502, 2512, 2512, 2522, 2502, 2512, 2512, 2522, + 2504, 2514, 2514, 2524, 2504, 2514, 2514, 2524, 2505, 2515, 2515, 2525, 2525, 2535, 2535, 2545, 2527, + 2537, 2537, 2547, 2527, 2537, 2537, 2547, 2529, 2539, 2539, 2549, 2529, 2539, 2539, 2549, 2530, 2540, + 2540, 2550, 2525, 2535, 2535, 2545, 2527, 2537, 2537, 2547, 2527, 2537, 2537, 2547, 2529, 2539, 2539, + 2549, 2529, 2539, 2539, 2549, 2530, 2540, 2540, 2550, 2560, 2570, 2570, 2580, 2562, 2572, 2572, 2582, + 2562, 2572, 2572, 2582, 2564, 2574, 2574, 2584, 2564, 2574, 2574, 2584, 2565, 2575, 2575, 2585 }; + + const std::vector t_coord = grid1.getCOORD(); + const std::vector t_zcorn = grid1.getZCORN(); + + BOOST_CHECK_EQUAL( t_coord.size() , ref_coord.size()); + + for (size_t i=0; i< t_coord.size(); i++) { + BOOST_CHECK_CLOSE( t_coord[i] , ref_coord[i], 1.0e-5); + } + + BOOST_CHECK_EQUAL( t_zcorn.size() , ref_zcorn.size()); + + for (size_t i=0; i< t_zcorn.size(); i++) { + BOOST_CHECK_CLOSE( t_zcorn[i] , ref_zcorn[i], 1.0e-5); + } +} + +static Opm::Deck BAD_CP_GRID() { + const char* deckData = + "RUNSPEC\n" + "\n" + "DIMENS\n" + "2 2 2 /\n" + "GRID\n" + "SPECGRID\n" + " 2 2 2 1 F /\n" + "COORD\n" + " 2002.0000 2002.0000 100.0000 1999.8255 1999.9127 108.4935\n" + " 2011.9939 2000.0000 100.3490 2009.8194 1999.9127 108.8425\n" + " 2015.9878 2000.0000 100.6980 2019.8133 1999.9127 109.1915\n" + " 2000.0000 2009.9985 100.1745 1999.8255 2009.9112 108.6681 \n" + " 2010.9939 2011.9985 100.5235 2009.8194 2009.9112 109.0170\n" + " 2019.9878 2009.9985 100.8725 2019.8133 2009.9112 109.3660\n" + " 2005.0000 2019.9970 100.3490 1999.8255 2019.9097 108.8426\n" + " 2009.9939 2019.9970 100.6980 2009.8194 2019.9097 109.1916\n" + " 2016.9878 2019.9970 101.0470 2019.8133 2019.9097 109.5406 /\n" + "ZCORN\n" + " 98.0000 100.3490 97.3490 100.6980 100.1745 100.5235\n" + " 100.5235 100.8725 100.1745 100.5235 100.5235 100.8725\n" + " 100.3490 101.6980 101.6980 102.5470 102.4973 102.1463\n" + " 103.2463 104.1953 103.6719 104.0209 104.0209 104.3698\n" + " 103.6719 104.0209 104.0209 104.3698 103.8464 104.1954\n" + " 104.1954 104.5444 103.4973 103.8463 103.8463 104.1953\n" + " 103.6719 104.0209 104.0209 104.3698 103.6719 104.0209\n" + " 104.0209 104.3698 103.8464 104.1954 104.1954 104.5444\n" + " 108.4935 108.8425 108.8425 109.1915 108.6681 109.0170\n" + " 109.0170 109.3660 108.6681 109.0170 109.0170 109.3660\n" + " 108.8426 109.1916 109.1916 109.5406 /\n" + "\n" + "EDIT\n"; + + Opm::Parser parser; + return parser.parseString( deckData); +} + +static Opm::Deck BAD_CP_GRID_MAPAXES() { + const char* deckData = + "RUNSPEC\n" + "\n" + "DIMENS\n" + "2 2 2 /\n" + "GRID\n" + "MAPAXES\n" + " 0. 100. 0. 0. 100. 0. /\n" + "\n" + "SPECGRID\n" + " 2 2 2 1 F /\n" + "COORD\n" + " 2002.0000 2002.0000 100.0000 1999.8255 1999.9127 108.4935\n" + " 2011.9939 2000.0000 100.3490 2009.8194 1999.9127 108.8425\n" + " 2015.9878 2000.0000 100.6980 2019.8133 1999.9127 109.1915\n" + " 2000.0000 2009.9985 100.1745 1999.8255 2009.9112 108.6681 \n" + " 2010.9939 2011.9985 100.5235 2009.8194 2009.9112 109.0170\n" + " 2019.9878 2009.9985 100.8725 2019.8133 2009.9112 109.3660\n" + " 2005.0000 2019.9970 100.3490 1999.8255 2019.9097 108.8426\n" + " 2009.9939 2019.9970 100.6980 2009.8194 2019.9097 109.1916\n" + " 2016.9878 2019.9970 101.0470 2019.8133 2019.9097 109.5406 /\n" + "ZCORN\n" + " 98.0000 100.3490 97.3490 100.6980 100.1745 100.5235\n" + " 100.5235 100.8725 100.1745 100.5235 100.5235 100.8725\n" + " 100.3490 101.6980 101.6980 102.5470 102.4973 102.1463\n" + " 103.2463 104.1953 103.6719 104.0209 104.0209 104.3698\n" + " 103.6719 104.0209 104.0209 104.3698 103.8464 104.1954\n" + " 104.1954 104.5444 103.4973 103.8463 103.8463 104.1953\n" + " 103.6719 104.0209 104.0209 104.3698 103.6719 104.0209\n" + " 104.0209 104.3698 103.8464 104.1954 104.1954 104.5444\n" + " 108.4935 108.8425 108.8425 109.1915 108.6681 109.0170\n" + " 109.0170 109.3660 108.6681 109.0170 109.0170 109.3660\n" + " 108.8426 109.1916 109.1916 109.5406 /\n" + "\n" + "EDIT\n"; + + Opm::Parser parser; + return parser.parseString( deckData); +} + +BOOST_AUTO_TEST_CASE(SAVE_FIELD_UNITS) { + + const char* deckData = + + "RUNSPEC\n" + "\n" + "DIMENS\n" + " 4 4 3 /\n" + "FIELD\n" + "GRID\n" + "DX\n" + " 48*300 /\n" + "DY\n" + " 48*300 /\n" + "DZ\n" + " 16*20 16*30 16*50 / \n" + "TOPS\n" + " 16*8325 / \n" + "EDIT\n" + "\n"; + + const char* deckData2 = + + "RUNSPEC\n" + "\n" + "DIMENS\n" + " 4 4 3 /\n" + "FIELD\n" + "GRID\n" + "MAPUNITS\n" + " METRES /\n" + "MAPAXES\n" + " 0.0 101.1 0.0 0.0 102.2 0.0 /\n" + "DX\n" + " 48*300 /\n" + "DY\n" + " 48*300 /\n" + "DZ\n" + " 16*20 16*30 16*50 / \n" + "TOPS\n" + " 16*8325 / \n" + "EDIT\n" + "\n"; + + const char* deckData3 = + + "RUNSPEC\n" + "\n" + "DIMENS\n" + " 4 4 3 /\n" + "FIELD\n" + "GRID\n" + "MAPUNITS\n" + " FEET /\n" + "MAPAXES\n" + " 0.0 102.2 0.0 0.0 103.3 0.0 /\n" + "DX\n" + " 48*300 /\n" + "DY\n" + " 48*300 /\n" + "DZ\n" + " 16*20 16*30 16*50 / \n" + "TOPS\n" + " 16*8325 / \n" + "EDIT\n" + "\n"; + + std::vector ref2_mapaxes = {0.0, 101.1, 0.0, 0.0, 102.2, 0.0 }; + std::vector ref3_mapaxes = {0.0, 102.2, 0.0, 0.0, 103.3, 0.0 }; + + Opm::Parser parser; + auto deck = parser.parseString( deckData) ; + + Opm::EclipseState es(deck); + Opm::UnitSystem units = es.getDeckUnitSystem(); + const auto length = ::Opm::UnitSystem::measure::length; + + const auto& grid1 = es.getInputGrid(); + + Opm::NNC nnc( deck ); + bool formatted = false; + + time_t timer; + time(&timer); + + std::string cwd = boost::filesystem::current_path().c_str(); + std::string testDir = cwd + "/tmp_dir_" + std::to_string(timer); + + if ( boost::filesystem::exists( testDir )) { + boost::filesystem::remove_all(testDir); + } + + boost::filesystem::create_directory(testDir); + + std::string fileName = testDir + "/" + "TMP.EGRID"; + grid1.save(fileName, formatted, nnc, units); + + Opm::EclIO::EclFile file1(fileName); + + // Values getZCORNed from the grid needs to be converted from SI to Field units + // and then converted from double to single precissions before comparing with values saved to + // the EGRID file + + // check coord + const std::vector coord_egrid = file1.get("COORD"); + std::vector coord_input_si = grid1.getCOORD(); + + BOOST_CHECK( coord_egrid.size() == coord_input_si.size()); + + std::vector coord_input_f; + coord_input_f.reserve(coord_input_si.size()); + + for (size_t n =0; n< coord_egrid.size(); n++) { + coord_input_f.push_back( static_cast(units.from_si(length, coord_input_si[n]))); + BOOST_CHECK_CLOSE( coord_input_f[n] , coord_egrid[n], 1e-6 ); + } + + // check zcorn + const std::vector zcorn_egrid = file1.get("ZCORN"); + std::vector zcorn_input_si = grid1.getZCORN(); + + BOOST_CHECK( zcorn_egrid.size() == zcorn_input_si.size()); + + std::vector zcorn_input_f; + zcorn_input_f.reserve(zcorn_input_si.size()); + + for (size_t n =0; n< zcorn_egrid.size(); n++) { + zcorn_input_f.push_back( static_cast(units.from_si(length, zcorn_input_si[n]))); + BOOST_CHECK_CLOSE( zcorn_input_f[n] , zcorn_egrid[n], 1e-6 ); + } + + BOOST_CHECK( file1.hasKey("GRIDUNIT")); + const std::vector gridunits = file1.get("GRIDUNIT"); + + BOOST_CHECK( gridunits[0]=="FEET"); + + // input deck do not hold MAPAXES or MAPUNITS entries. Below keywords should not be written to EGRID file + BOOST_CHECK( !file1.hasKey("MAPAXES")); + BOOST_CHECK( !file1.hasKey("MAPUNITS")); + + // this deck do not have any nnc. Below keywords should not be written to EGRID file + BOOST_CHECK( !file1.hasKey("NNCHEAD")); + BOOST_CHECK( !file1.hasKey("NNC1")); + BOOST_CHECK( !file1.hasKey("NNC2")); + + // testing deck in field units and MAPUNITS in METRES + auto deck2 = parser.parseString( deckData2) ; + + Opm::EclipseState es2(deck2); + Opm::UnitSystem units2 = es.getDeckUnitSystem(); + Opm::NNC nnc2( deck2 ); + + const auto& grid2 = es2.getInputGrid(); + + std::string fileName2 = testDir + "/" + "TMP2.FEGRID"; + + grid2.save(fileName2, true, nnc2, units); + + Opm::EclIO::EclFile file2(fileName2); + + const std::vector& test_mapunits2 = file2.get("MAPUNITS"); + BOOST_CHECK( test_mapunits2[0] == "METRES"); + + const std::vector& test_mapaxes2 = file2.get("MAPAXES"); + + BOOST_CHECK( test_mapaxes2.size() == ref2_mapaxes.size()); + + for (size_t n =0; n< ref2_mapaxes.size(); n++) { + BOOST_CHECK( ref2_mapaxes[n] == test_mapaxes2[n]); + } + + // testing deck in field units and MAPUNITS in FEET + auto deck3 = parser.parseString( deckData3) ; + + Opm::EclipseState es3(deck3); + Opm::UnitSystem units3 = es.getDeckUnitSystem(); + Opm::NNC nnc3( deck3 ); + + const auto& grid3 = es3.getInputGrid(); + + std::string fileName3 = testDir + "/" + "TMP3.FEGRID"; + + grid3.save(fileName3, true, nnc3, units3); + + Opm::EclIO::EclFile file3(fileName3); + + const std::vector& test_mapunits3 = file3.get("MAPUNITS"); + BOOST_CHECK( test_mapunits3[0] == "FEET"); + + const std::vector& test_mapaxes3 = file3.get("MAPAXES"); + + BOOST_CHECK( test_mapaxes3.size() == ref3_mapaxes.size()); + + for (size_t n =0; n< ref3_mapaxes.size(); n++) { + BOOST_CHECK( ref3_mapaxes[n] == test_mapaxes3[n]); + } + + boost::filesystem::remove_all(testDir); +} + +BOOST_AUTO_TEST_CASE(SAVE_METRIC_UNITS) { + + const char* deckData1 = + + "RUNSPEC\n" + "\n" + "DIMENS\n" + " 4 4 3 /\n" + "GRID\n" + "MAPAXES\n" + " 0.0 45000.0 0.0 0.0 720000.0 0.0 / \n" + "MAPUNITS\n" + " METRES / \n" + "DX\n" + " 48*300 /\n" + "DY\n" + " 48*300 /\n" + "DZ\n" + " 16*20 16*30 16*50 / \n" + "TOPS\n" + " 16*8325 / \n" + "NNC\n" + " 2 2 1 2 3 2 0.95 / \n" + " 3 2 1 3 3 2 1.05 / \n" + " 4 2 1 4 3 2 1.15 / \n" + "/ \n" + "EDIT\n" + "\n"; + + const char* deckData2 = + + "RUNSPEC\n" + "\n" + "DIMENS\n" + " 4 4 3 /\n" + "GRID\n" + "MAPAXES\n" + " 0.0 450.0 0.0 0.0 7200.0 0.0 / \n" + "MAPUNITS\n" + " FEET / \n" + "DX\n" + " 48*300 /\n" + "DY\n" + " 48*300 /\n" + "DZ\n" + " 16*20 16*30 16*50 / \n" + "TOPS\n" + " 16*8325 / \n" + "NNC\n" + " 2 2 1 2 3 2 0.95 / \n" + " 3 2 1 3 3 2 1.05 / \n" + " 4 2 1 4 3 2 1.15 / \n" + "/ \n" + "EDIT\n" + "\n"; + + std::vector ref_mapaxes1 = { 0.0, 45000.0, 0.0, 0.0, 720000.0, 0.0 }; + std::vector ref_mapaxes2 = { 0.0, 450.0, 0.0, 0.0, 7200.0, 0.0 }; + + Opm::Parser parser; + auto deck1 = parser.parseString( deckData1) ; + + Opm::EclipseState es1(deck1); + Opm::UnitSystem units1 = es1.getDeckUnitSystem(); + const auto length = ::Opm::UnitSystem::measure::length; + + const auto& grid1 = es1.getInputGrid(); + Opm::NNC nnc( deck1 ); + + bool formatted = true; + + time_t timer; + time(&timer); + + std::string cwd = boost::filesystem::current_path().c_str(); + std::string testDir = cwd + "/tmp_dir_" + std::to_string(timer); + + if ( boost::filesystem::exists( testDir )) { + boost::filesystem::remove_all(testDir); + } + + boost::filesystem::create_directory(testDir); + + std::string fileName = testDir + "/" + "TMP.FEGRID"; + grid1.save(fileName, formatted, nnc, units1); + + Opm::EclIO::EclFile file1(fileName); + + // Values getZCORNed from the grid have same units as input deck (metric), however these needs to be + // converted from double to single precissions before comparing with values saved to the EGRID file + + // check coord + const std::vector coord_egrid = file1.get("COORD"); + std::vector coord_input_si = grid1.getCOORD(); + + BOOST_CHECK( coord_egrid.size() == coord_input_si.size()); + + std::vector coord_input_f; + coord_input_f.reserve(coord_input_si.size()); + + for (size_t n =0; n< coord_egrid.size(); n++) { + coord_input_f.push_back( static_cast(units1.from_si(length, coord_input_si[n]))); + BOOST_CHECK_CLOSE( coord_input_f[n] , coord_egrid[n], 1e-6 ); + } + + // check zcorn + const std::vector zcorn_egrid = file1.get("ZCORN"); + std::vector zcorn_input_si = grid1.getZCORN(); + + BOOST_CHECK( zcorn_egrid.size() == zcorn_input_si.size()); + + std::vector zcorn_input_f; + zcorn_input_f.reserve(zcorn_input_si.size()); + + for (size_t n =0; n< zcorn_egrid.size(); n++) { + zcorn_input_f.push_back( static_cast(units1.from_si(length, zcorn_input_si[n]))); + BOOST_CHECK_CLOSE( zcorn_input_f[n] , zcorn_egrid[n], 1e-6 ); + } + + BOOST_CHECK( file1.hasKey("GRIDUNIT")); + const std::vector gridunits = file1.get("GRIDUNIT"); + + BOOST_CHECK( gridunits[0]=="METRES"); + + BOOST_CHECK( file1.hasKey("MAPAXES")); + std::vector mapaxes = file1.get("MAPAXES"); + + for (size_t n = 0; n < 6; n++) { + BOOST_CHECK_CLOSE( mapaxes[n] , ref_mapaxes1[n], 1e-6 ); + } + + BOOST_CHECK( file1.hasKey("MAPUNITS")); + const std::vector mapunits = file1.get("MAPUNITS"); + BOOST_CHECK( gridunits[0]=="METRES"); + + BOOST_CHECK( file1.hasKey("NNCHEAD")); + const std::vector nnchead = file1.get("NNCHEAD"); + + BOOST_CHECK( nnchead[0] == static_cast(nnc.numNNC()) ); + + std::vector ref_nnc1 = { 6, 7, 8 }; + std::vector ref_nnc2 = { 26, 27, 28 }; + + BOOST_CHECK( file1.hasKey("NNC1")); + BOOST_CHECK( file1.hasKey("NNC2")); + + const std::vector nnc1 = file1.get("NNC1"); + const std::vector nnc2 = file1.get("NNC2"); + + BOOST_CHECK( nnc1.size() == nnc2.size() ); + + for (size_t n =0; n< nnc1.size(); n++) { + BOOST_CHECK( nnc1[n] == ref_nnc1[n] ); + } + + for (size_t n =0; n< nnc2.size(); n++) { + BOOST_CHECK( nnc2[n] == ref_nnc2[n] ); + } + + // testing deck in metric units with mapaxes in field units + auto deck2 = parser.parseString( deckData2) ; + + Opm::EclipseState es2(deck2); + Opm::UnitSystem units2 = es2.getDeckUnitSystem(); + + const auto& grid2 = es2.getInputGrid(); + //Opm::NNC nnc( deck2 ); + + std::string fileName2 = testDir + "/" + "TMP2.FEGRID"; + + grid2.save(fileName2, true, nnc, units2); + + Opm::EclIO::EclFile file2(fileName2); + + const std::vector& test_mapunits2 = file2.get("MAPUNITS"); + BOOST_CHECK( test_mapunits2[0] == "FEET"); + + const std::vector& test_mapaxes2 = file2.get("MAPAXES"); + + BOOST_CHECK( test_mapaxes2.size() == ref_mapaxes2.size()); + + for (size_t n =0; n< ref_mapaxes2.size(); n++) { + BOOST_CHECK( ref_mapaxes2[n] == test_mapaxes2[n]); + } + + + boost::filesystem::remove_all(testDir); +} + +BOOST_AUTO_TEST_CASE(CalcCellDims) { + + Opm::Deck deck = BAD_CP_GRID(); + Opm::EclipseGrid grid( deck ); + + std::array dims = grid.getNXYZ(); + + size_t nCells = dims[0]*dims[1]*dims[2]; + + std::vector dz_ref = { 0.33223500E+01, 0.40973248E+01, 0.32474000E+01, 0.28723750E+01, 0.49961748E+01, 0.49961748E+01, + 0.49961748E+01, 0.49961748E+01 + }; + + std::vector dx_ref = { 0.10309320E+02, 0.70301223E+01, 0.84377403E+01, 0.85725355E+01, 0.10140956E+02, 0.89693098E+01, + 0.94102650E+01, 0.94102678E+01 + }; + + std::vector dy_ref = { 0.99226236E+01, 0.10826077E+02, 0.93370037E+01, 0.93144703E+01, 0.10008223E+02, 0.10302064E+02, + 0.97221985E+01, 0.97221985E+01 + }; + + std::vector depth_ref = { 0.10142293E+03, 0.10190942E+03, 0.10230995E+03, 0.10284644E+03, 0.10625719E+03, + 0.10660616E+03, 0.10643174E+03, 0.10678072E+03 + }; + + for (size_t n=0; n cellDims = grid.getCellDims(n); + + BOOST_CHECK_CLOSE( cellDims[0] , dx_ref[n], 1e-5 ); + BOOST_CHECK_CLOSE( cellDims[1] , dy_ref[n], 1e-5 ); + BOOST_CHECK_CLOSE( cellDims[2] , dz_ref[n], 1e-5 ); + + BOOST_CHECK_CLOSE( grid.getCellDepth(n) , depth_ref[n], 1e-5 ); + } + + for (int k = 0; k < dims[2]; k++) { + for (int j = 0; j < dims[1]; j++) { + for (int i = 0; i < dims[0]; i++) { + size_t globInd = i + j*dims[0] + k*dims[0]*dims[1]; + BOOST_CHECK_CLOSE( grid.getCellThickness(i, j, k) , dz_ref[globInd], 1e-5 ); + + std::array cellDims = grid.getCellDims(i, j, k); + + BOOST_CHECK_CLOSE( cellDims[0] , dx_ref[globInd], 1e-5 ); + BOOST_CHECK_CLOSE( cellDims[1] , dy_ref[globInd], 1e-5 ); + BOOST_CHECK_CLOSE( cellDims[2] , dz_ref[globInd], 1e-5 ); + + BOOST_CHECK_CLOSE( grid.getCellDepth(i, j, k) , depth_ref[globInd], 1e-5 ); + } + } + } +} + +BOOST_AUTO_TEST_CASE(ExportMAPAXES_TEST) { + + Opm::Deck deck1 = BAD_CP_GRID_MAPAXES(); + Opm::EclipseGrid grid1( deck1 ); + + std::vector ref_mapaxes = { 0.0, 100.0, 0.0, 0.0, 100.0, 0.0 }; + + std::vector mapaxes = grid1.getMAPAXES(); + + for (size_t n=0; n< mapaxes.size(); n++ ) { + BOOST_CHECK_EQUAL( ref_mapaxes[n] , mapaxes[n]); + } + + Opm::Deck deck2 = BAD_CP_GRID(); + Opm::EclipseGrid grid2( deck2 ); + + BOOST_CHECK( !grid1.equal( grid2 )); + + std::vector coord = grid1.getCOORD(); + std::vector zcorn = grid1.getZCORN(); + std::vector actnum = grid1.getACTNUM(); + + std::array dims = grid1.getNXYZ(); + + Opm::EclipseGrid grid3(dims, coord, zcorn, actnum.data(), mapaxes.data()); + + BOOST_CHECK( grid3.equal( grid1 )); + + mapaxes[1] = 101; + Opm::EclipseGrid grid4(dims, coord, zcorn, actnum.data(), mapaxes.data()); + + BOOST_CHECK( !grid4.equal( grid1 )); +} + +BOOST_AUTO_TEST_CASE(TESTCP_ACTNUM_UPDATE) { + const char* deckData = + + "RUNSPEC\n" + "\n" + "DIMENS\n" + " 3 2 1 /\n" + "GRID\n" + "COORD\n" + " 2000.0000 2000.0000 2000.0000 1999.9476 2000.0000 2002.9995\n" + " 2049.9924 2000.0000 2000.8726 2049.9400 2000.0000 2003.8722 \n" + " 2099.9848 2000.0000 2001.7452 2099.9324 2000.0000 2004.7448 \n" + " 2149.9772 2000.0000 2002.6179 2149.9248 2000.0000 2005.6174 \n" + " 2000.0000 2050.0000 2000.0000 1999.9476 2050.0000 2002.9995 \n" + " 2049.9924 2050.0000 2000.8726 2049.9400 2050.0000 2003.8722 \n" + " 2099.9848 2050.0000 2001.7452 2099.9324 2050.0000 2004.7448 \n" + " 2149.9772 2050.0000 2002.6179 2149.9248 2050.0000 2005.6174 \n" + " 2000.0000 2100.0000 2000.0000 1999.9476 2100.0000 2002.9995 \n" + " 2049.9924 2100.0000 2000.8726 2049.9400 2100.0000 2003.8722 \n" + " 2099.9848 2100.0000 2001.7452 2099.9324 2100.0000 2004.7448 \n" + " 2149.9772 2100.0000 2002.6179 2149.9248 2100.0000 2005.6174 / \n" + "ZCORN\n" + " 2000.0000 2000.8726 2000.8726 2001.7452 2001.7452 2002.6179 \n" + " 2000.0000 2000.8726 2000.8726 2001.7452 2001.7452 2002.6179 \n" + " 2000.0000 2000.8726 2000.8726 2001.7452 2001.7452 2002.6179 \n" + " 2000.0000 2000.8726 2000.8726 2001.7452 2001.7452 2002.6179 \n" + " 2002.9995 2003.8722 2003.8722 2004.7448 2004.7448 2005.6174 \n" + " 2002.9995 2003.8722 2003.8722 2004.7448 2004.7448 2005.6174 \n" + " 2002.9995 2003.8722 2003.8722 2004.7448 2004.7448 2005.6174 \n" + " 2002.9995 2003.8722 2003.8722 2004.7448 2004.7448 2005.6174 / \n" + "ACTNUM\n" + " 0 1 1 1 0 1 / \n" + "EDIT\n" + "\n"; + + Opm::Parser parser; + auto deck = parser.parseString( deckData) ; + + std::vector actInDeck = {0, 1, 1, 1, 0, 1}; + std::vector newAct = {1, 0, 0, 0, 1}; + + Opm::EclipseGrid grid1( deck); + Opm::EclipseGrid grid2( deck, newAct.data()); + + std::vector actGrid1 = grid1.getACTNUM(); + std::vector actGrid2 = grid2.getACTNUM(); + + BOOST_CHECK( actGrid1.size() == actGrid2.size()); + + for (size_t n=0; n< actGrid1.size(); n++) { + BOOST_CHECK( actGrid1[n] == actInDeck[n]); + BOOST_CHECK( actGrid2[n] == newAct[n]); + } +} + + +BOOST_AUTO_TEST_CASE(TEST_altGridConstructors) { + + const char* deckData = + + "RUNSPEC\n" + "\n" + "DIMENS\n" + " 2 2 3 /\n" + "GRID\n" + "SPECGRID\n" + " 2 2 3 1 F / \n" + "\n" + "COORD\n" + " 2000.0000 2000.0000 1497.0000 1999.9127 2000.0000 1511.9977\n" + " 2049.9924 2000.0000 1500.8726 2049.9051 2000.0000 1515.8703\n" + " 2099.9848 2000.0000 1501.7452 2099.8975 2000.0000 1516.7430 \n" + " 2000.0000 2050.0000 1497.0000 1999.9127 2050.0000 1511.9977\n" + " 2049.9924 2050.0000 1500.8726 2049.9051 2050.0000 1515.8703\n" + " 2099.9848 2050.0000 1501.7452 2099.8975 2050.0000 1516.7430\n" + " 2000.0000 2100.0000 1497.0000 1999.9127 2100.0000 1511.9977\n" + " 2049.9924 2100.0000 1500.8726 2049.9051 2100.0000 1515.8703\n" + " 2099.9848 2100.0000 1501.7452 2099.8975 2100.0000 1516.7430 /\n" + "\n" + "ZCORN\n" + " 1497.0000 1497.8726 1500.8726 1501.7452 1497.0000 1497.8726\n" + " 1500.8726 1501.7452 1497.0000 1497.8726 1500.8726 1501.7452\n" + " 1497.0000 1497.8726 1500.8726 1501.7452 1501.9992 1502.8719\n" + " 1505.8719 1506.7445 1501.9992 1502.8719 1505.8719 1506.7445\n" + " 1501.9992 1502.8719 1505.8719 1506.7445 1501.9992 1502.8719\n" + " 1505.8719 1506.7445 1501.9992 1502.8719 1505.8719 1506.7445\n" + " 1501.9992 1502.8719 1505.8719 1506.7445 1501.9992 1502.8719\n" + " 1505.8719 1506.7445 1501.9992 1502.8719 1505.8719 1506.7445\n" + " 1506.9985 1507.8711 1510.8711 1511.7437 1506.9985 1507.8711\n" + " 1510.8711 1511.7437 1506.9985 1507.8711 1510.8711 1511.7437\n" + " 1506.9985 1507.8711 1510.8711 1511.7437 1506.9985 1507.8711\n" + " 1510.8711 1511.7437 1506.9985 1507.8711 1510.8711 1511.7437\n" + " 1506.9985 1507.8711 1510.8711 1511.7437 1506.9985 1507.8711\n" + " 1510.8711 1511.7437 1511.9977 1512.8703 1515.8703 1516.7430\n" + " 1511.9977 1512.8703 1515.8703 1516.7430 1511.9977 1512.8703\n" + " 1515.8703 1516.7430 1511.9977 1512.8703 1515.8703 1516.7430 /\n" + "\n" + "ACTNUM\n" + " 1 1 1 1 1 0 1 1 1 0 1 1 /\n" + + "EDIT\n" + "\n"; + + Opm::Parser parser; + auto deck = parser.parseString( deckData) ; + + Opm::EclipseGrid grid1( deck); + + std::vector actnum = grid1.getACTNUM(); + std::vector coord = grid1.getCOORD(); + std::vector zcorn = grid1.getZCORN(); + + Opm::EclipseGrid grid2( grid1 , zcorn.data(), actnum); + //Opm::EclipseGrid grid2( grid1 , zcorn, actnum); + + BOOST_CHECK( grid1.equal( grid2) ); + + std::vector emptyZcorn; + + Opm::EclipseGrid grid3( grid1 , emptyZcorn.data(), actnum); + BOOST_CHECK( grid1.equal( grid3) ); +} + +static Opm::Deck BAD_CP_GRID_ACTNUM() { + const char* deckData = + "RUNSPEC\n" + "\n" + "DIMENS\n" + "2 2 2 /\n" + "GRID\n" + "SPECGRID\n" + " 2 2 2 1 F /\n" + "COORD\n" + " 2002.0000 2002.0000 100.0000 1999.8255 1999.9127 108.4935\n" + " 2011.9939 2000.0000 100.3490 2009.8194 1999.9127 108.8425\n" + " 2015.9878 2000.0000 100.6980 2019.8133 1999.9127 109.1915\n" + " 2000.0000 2009.9985 100.1745 1999.8255 2009.9112 108.6681 \n" + " 2010.9939 2011.9985 100.5235 2009.8194 2009.9112 109.0170\n" + " 2019.9878 2009.9985 100.8725 2019.8133 2009.9112 109.3660\n" + " 2005.0000 2019.9970 100.3490 1999.8255 2019.9097 108.8426\n" + " 2009.9939 2019.9970 100.6980 2009.8194 2019.9097 109.1916\n" + " 2016.9878 2019.9970 101.0470 2019.8133 2019.9097 109.5406 /\n" + "ZCORN\n" + " 98.0000 100.3490 97.3490 100.6980 100.1745 100.5235\n" + " 100.5235 100.8725 100.1745 100.5235 100.5235 100.8725\n" + " 100.3490 101.6980 101.6980 102.5470 102.4973 102.1463\n" + " 103.2463 104.1953 103.6719 104.0209 104.0209 104.3698\n" + " 103.6719 104.0209 104.0209 104.3698 103.8464 104.1954\n" + " 104.1954 104.5444 103.4973 103.8463 103.8463 104.1953\n" + " 103.6719 104.0209 104.0209 104.3698 103.6719 104.0209\n" + " 104.0209 104.3698 103.8464 104.1954 104.1954 104.5444\n" + " 108.4935 108.8425 108.8425 109.1915 108.6681 109.0170\n" + " 109.0170 109.3660 108.6681 109.0170 109.0170 109.3660\n" + " 108.8426 109.1916 109.1916 109.5406 /\n" + "\n" + "ACTNUM\n" + " 1 1 1 1 0 1 0 1 /\n" + "EDIT\n"; + + Opm::Parser parser; + return parser.parseString( deckData ); +} + +BOOST_AUTO_TEST_CASE(TEST_getCellCenters) { + + Opm::Deck deck1 = BAD_CP_GRID_ACTNUM(); + Opm::EclipseGrid grid1( deck1 ); + + std::vector> ref_centers_prev = { + { 2.006104082026e+03, 2.005869733884e+03, 1.014229250000e+02 }, + { 2.014871600974e+03, 2.005382958337e+03, 1.019094125000e+02 }, + { 2.006149627699e+03, 2.015375548018e+03, 1.023099500000e+02 }, + { 2.014617670881e+03, 2.015373620907e+03, 1.028464375000e+02 }, + { 2.014794142285e+03, 2.005084683036e+03, 1.066061625000e+02 }, + { 2.014720614821e+03, 2.015083182885e+03, 1.067807125000e+02 } + }; + + std::vector> ref_dims_prev = { + { 1.030932076324e+01, 9.922624077310e+00, 3.322350000000e+00 }, + { 7.030122275667e+00, 1.082607766242e+01, 4.097325000000e+00 }, + { 8.437740561450e+00, 9.337003646693e+00, 3.247400000000e+00 }, + { 8.572535163904e+00, 9.314470420138e+00, 2.872375000000e+00 }, + { 8.969310273018e+00, 1.030206365401e+01, 4.996175000000e+00 }, + { 9.410267880374e+00, 9.722198202980e+00, 4.996175000000e+00 } + }; + + std::vector actMap = grid1.getActiveMap(); + + int n = 0; + for (auto ind : actMap) { + std::array cellC = grid1.getCellCenter(ind); + std::array cellD = grid1.getCellDims(ind); + + for (size_t i = 0; i < 3; i++) { + BOOST_CHECK_CLOSE( ref_centers_prev[n][i], cellC[i], 1e-10); + BOOST_CHECK_CLOSE( ref_dims_prev[n][i], cellD[i], 1e-10); + } + + n++; + } +} + +BOOST_AUTO_TEST_CASE(LoadFromBinary) { + BOOST_CHECK_THROW(Opm::EclipseGrid( "No/does/not/exist" ) , std::invalid_argument); +} + + +BOOST_AUTO_TEST_CASE(TEST_constructFromEgrid) { + + const char* deckData = + + "RUNSPEC\n" + "\n" + "DIMENS\n" + " 10 10 3 /\n" + "FIELD\n" + "GRID\n" + "DX\n" + " 300*1000 /\n" + "DY\n" + " 300*1000 /\n" + "DZ\n" + " 100*20 100*30 100*50 / \n" + "TOPS\n" + " 100*8325 / \n" + "ACTNUM\n" + " 44*1 3*0 7*1 3*0 243*1/\n" + "EDIT\n" + "\n"; + + Opm::Parser parser; + auto deck = parser.parseString( deckData) ; + + Opm::EclipseGrid grid1( deck); + Opm::EclipseGrid grid2( "SPE1CASE1.EGRID"); + + // compare actnum + std::vector actGrid1 = grid1.getACTNUM(); + std::vector actGrid2 = grid2.getACTNUM(); + + BOOST_CHECK( actGrid1.size() == actGrid2.size() ); + + for (size_t n= 0; n< actGrid1.size(); n++) { + BOOST_CHECK( actGrid1[n] == actGrid2[n] ); + } + + // compare coord + std::vector coordGrid1 = grid1.getCOORD(); + std::vector coordGrid2 = grid2.getCOORD(); + + BOOST_CHECK( coordGrid1.size() == coordGrid2.size() ); + + std::vector zcornGrid1 = grid1.getZCORN(); + std::vector zcornGrid2 = grid2.getZCORN(); + + BOOST_CHECK( zcornGrid1.size() == zcornGrid2.size() ); + + for (size_t n = 0; n < zcornGrid1.size(); n++){ + BOOST_CHECK_CLOSE( zcornGrid1[n], zcornGrid2[n], 1.0e-6 ); + } + + BOOST_CHECK( grid1.getCartesianSize() == grid2.getCartesianSize() ); + + for (size_t n=0; n < grid1.getCartesianSize(); n++){ + BOOST_CHECK_CLOSE( grid1.getCellVolume(n), grid2.getCellVolume(n), 1e-6); + + } +} + + + +BOOST_AUTO_TEST_CASE(TEST_GDFILE_1) { + + const char* deckData1 = + "RUNSPEC\n" + "DIMENS\n" + "1 1 2 /\n" + "GRID\n" + "COORD\n" + "10.0000 10.0000 2000.0000 9.8255 10.0000 2014.9977\n" + "109.9848 10.0000 2001.7452 109.8102 10.0000 2016.7430\n" + "10.0000 110.0000 2000.0000 9.8255 110.0000 2014.9977\n" + "109.9848 110.0000 2001.7452 109.8102 110.0000 2016.7430 /\n" + "ZCORN\n" + "2000.0000 2001.7452 2000.0000 2001.7452 2004.9992 2006.7445\n" + "2004.9992 2006.7445 2004.9992 2006.7445 2004.9992 2006.7445\n" + "2014.9977 2016.7430 2014.9977 2016.7430 /\n"; + + const char* deckData2 = + "RUNSPEC\n" + "DIMENS\n" + "1 1 2 /\n" + "GRID\n" + "GDFILE\n" + " 'BAD_CP_M.EGRID' /\n" + "COORD\n" + "10.0000 10.0000 2000.0000 9.8255 10.0000 2014.9977\n" + "109.9848 10.0000 2001.7452 109.8102 10.0000 2016.7430\n" + "10.0000 110.0000 2000.0000 9.8255 110.0000 2014.9977\n" + "109.9848 110.0000 2001.7452 109.8102 110.0000 2016.7430 /\n"; + + const char* deckData3 = + "RUNSPEC\n" + "DIMENS\n" + "1 1 2 /\n" + "GRID\n" + "GDFILE\n" + " 'BAD_CP_M.EGRID' /\n" + "ZCORN\n" + "2000.0000 2001.7452 2000.0000 2001.7452 2004.9992 2006.7445\n" + "2004.9992 2006.7445 2004.9992 2006.7445 2004.9992 2006.7445\n" + "2014.9977 2016.7430 2014.9977 2016.7430 /\n"; + + Opm::Parser parser; + auto deck1 = parser.parseString( deckData1) ; + auto deck2 = parser.parseString( deckData2) ; + auto deck3 = parser.parseString( deckData3) ; + + BOOST_CHECK_NO_THROW( Opm::EclipseGrid grid1(deck1) ); + BOOST_CHECK_THROW(Opm::EclipseGrid grid2(deck2), std::invalid_argument); + BOOST_CHECK_THROW(Opm::EclipseGrid grid3(deck3), std::invalid_argument); +} + + +BOOST_AUTO_TEST_CASE(TEST_GDFILE_2) { + + const char* deckData1 = + "RUNSPEC\n" + "\n" + "DIMENS\n" + "2 2 2 /\n" + "GRID\n" + "SPECGRID\n" + " 2 2 2 1 F /\n" + "COORD\n" + " 2002.0000 2002.0000 100.0000 1999.8255 1999.9127 108.4935\n" + " 2011.9939 2000.0000 100.3490 2009.8194 1999.9127 108.8425\n" + " 2015.9878 2000.0000 100.6980 2019.8133 1999.9127 109.1915\n" + " 2000.0000 2009.9985 100.1745 1999.8255 2009.9112 108.6681 \n" + " 2010.9939 2011.9985 100.5235 2009.8194 2009.9112 109.0170\n" + " 2019.9878 2009.9985 100.8725 2019.8133 2009.9112 109.3660\n" + " 2005.0000 2019.9970 100.3490 1999.8255 2019.9097 108.8426\n" + " 2009.9939 2019.9970 100.6980 2009.8194 2019.9097 109.1916\n" + " 2016.9878 2019.9970 101.0470 2019.8133 2019.9097 109.5406 /\n" + "ZCORN\n" + " 98.0000 100.3490 97.3490 100.6980 100.1745 100.5235\n" + " 100.5235 100.8725 100.1745 100.5235 100.5235 100.8725\n" + " 100.3490 101.6980 101.6980 102.5470 102.4973 102.1463\n" + " 103.2463 104.1953 103.6719 104.0209 104.0209 104.3698\n" + " 103.6719 104.0209 104.0209 104.3698 103.8464 104.1954\n" + " 104.1954 104.5444 103.4973 103.8463 103.8463 104.1953\n" + " 103.6719 104.0209 104.0209 104.3698 103.6719 104.0209\n" + " 104.0209 104.3698 103.8464 104.1954 104.1954 104.5444\n" + " 108.4935 108.8425 108.8425 109.1915 108.6681 109.0170\n" + " 109.0170 109.3660 108.6681 109.0170 109.0170 109.3660\n" + " 108.8426 109.1916 109.1916 109.5406 /\n" + "\n" + "EDIT\n"; + + const char* deckData1a = + "RUNSPEC\n" + "\n" + "DIMENS\n" + "2 2 2 /\n" + "GRID\n" + "SPECGRID\n" + " 2 2 2 1 F /\n" + "COORD\n" + " 2002.0000 2002.0000 100.0000 1999.8255 1999.9127 108.4935\n" + " 2011.9939 2000.0000 100.3490 2009.8194 1999.9127 108.8425\n" + " 2015.9878 2000.0000 100.6980 2019.8133 1999.9127 109.1915\n" + " 2000.0000 2009.9985 100.1745 1999.8255 2009.9112 108.6681 \n" + " 2010.9939 2011.9985 100.5235 2009.8194 2009.9112 109.0170\n" + " 2019.9878 2009.9985 100.8725 2019.8133 2009.9112 109.3660\n" + " 2005.0000 2019.9970 100.3490 1999.8255 2019.9097 108.8426\n" + " 2009.9939 2019.9970 100.6980 2009.8194 2019.9097 109.1916\n" + " 2016.9878 2019.9970 101.0470 2019.8133 2019.9097 109.5406 /\n" + "ZCORN\n" + " 98.0000 100.3490 97.3490 100.6980 100.1745 100.5235\n" + " 100.5235 100.8725 100.1745 100.5235 100.5235 100.8725\n" + " 100.3490 101.6980 101.6980 102.5470 102.4973 102.1463\n" + " 103.2463 104.1953 103.6719 104.0209 104.0209 104.3698\n" + " 103.6719 104.0209 104.0209 104.3698 103.8464 104.1954\n" + " 104.1954 104.5444 103.4973 103.8463 103.8463 104.1953\n" + " 103.6719 104.0209 104.0209 104.3698 103.6719 104.0209\n" + " 104.0209 104.3698 103.8464 104.1954 104.1954 104.5444\n" + " 108.4935 108.8425 108.8425 109.1915 108.6681 109.0170\n" + " 109.0170 109.3660 108.6681 109.0170 109.0170 109.3660\n" + " 108.8426 109.1916 109.1916 109.5406 /\n" + "\n" + "ACTNUM\n" + " 1 1 1 1 0 1 0 1 /\n" + "EDIT\n"; + + const char* deckData1b = + "RUNSPEC\n" + "\n" + "DIMENS\n" + "2 2 2 /\n" + "FIELD\n" + "GRID\n" + "MAPUNITS\n" + " METRES /\n" + "MAPAXES\n" + " 0. 100. 0. 0. 100. 0. /\n" + "SPECGRID\n" + " 2 2 2 1 F /\n" + "COORD\n" + " 2002.0000 2002.0000 100.0000 1999.8255 1999.9127 108.4935\n" + " 2011.9939 2000.0000 100.3490 2009.8194 1999.9127 108.8425\n" + " 2015.9878 2000.0000 100.6980 2019.8133 1999.9127 109.1915\n" + " 2000.0000 2009.9985 100.1745 1999.8255 2009.9112 108.6681 \n" + " 2010.9939 2011.9985 100.5235 2009.8194 2009.9112 109.0170\n" + " 2019.9878 2009.9985 100.8725 2019.8133 2009.9112 109.3660\n" + " 2005.0000 2019.9970 100.3490 1999.8255 2019.9097 108.8426\n" + " 2009.9939 2019.9970 100.6980 2009.8194 2019.9097 109.1916\n" + " 2016.9878 2019.9970 101.0470 2019.8133 2019.9097 109.5406 /\n" + "ZCORN\n" + " 98.0000 100.3490 97.3490 100.6980 100.1745 100.5235\n" + " 100.5235 100.8725 100.1745 100.5235 100.5235 100.8725\n" + " 100.3490 101.6980 101.6980 102.5470 102.4973 102.1463\n" + " 103.2463 104.1953 103.6719 104.0209 104.0209 104.3698\n" + " 103.6719 104.0209 104.0209 104.3698 103.8464 104.1954\n" + " 104.1954 104.5444 103.4973 103.8463 103.8463 104.1953\n" + " 103.6719 104.0209 104.0209 104.3698 103.6719 104.0209\n" + " 104.0209 104.3698 103.8464 104.1954 104.1954 104.5444\n" + " 108.4935 108.8425 108.8425 109.1915 108.6681 109.0170\n" + " 109.0170 109.3660 108.6681 109.0170 109.0170 109.3660\n" + " 108.8426 109.1916 109.1916 109.5406 /\n" + "\n" + "ACTNUM\n" + " 1 1 1 1 0 1 0 1 /\n" + "EDIT\n"; + + const char* deckData2 = + "RUNSPEC\n" + "\n" + "DIMENS\n" + "2 2 2 /\n" + "GRID\n" + "GDFILE\n" + " 'BAD_CP_M.EGRID' /\n" + "EDIT\n"; + + const char* deckData3a = + "RUNSPEC\n" + "\n" + "DIMENS\n" + "2 2 2 /\n" + "GRID\n" + "ACTNUM\n" + " 1 0 1 0 1 1 1 1 /\n" + "MAPUNITS\n" + " FEET /\n" + "MAPAXES\n" + " 0. 200. 0. 0. 200. 0. /\n" + "GDFILE\n" + " 'BAD_CP_M.EGRID' /\n" + "EDIT\n"; + + const char* deckData3b = + "RUNSPEC\n" + "\n" + "DIMENS\n" + "2 2 2 /\n" + "GRID\n" + "MAPUNITS\n" + " FEET /\n" + "MAPAXES\n" + " 0. 200. 0. 0. 200. 0. /\n" + "GDFILE\n" + " 'BAD_CP_F.EGRID' /\n" + "ACTNUM\n" + " 1 0 1 0 1 1 1 1 /\n" + "EDIT\n"; + + const char* deckData3c = + "RUNSPEC\n" + "\n" + "DIMENS\n" + "2 2 2 /\n" + "GRID\n" + "GDFILE\n" + " 'BAD_CP_F.EGRID' /\n" + "MAPUNITS\n" + " FEET /\n" + "MAPAXES\n" + " 0. 200. 0. 0. 200. 0. /\n" + "ACTNUM\n" + " 1 0 1 0 1 1 1 1 /\n" + "EDIT\n"; + + Opm::Parser parser; + + std::vector ref_act_egrid = {1, 1, 1, 1, 0, 1, 0, 1}; + std::vector ref_act_deck3 = {1, 0, 1, 0, 1, 1, 1, 1}; + + std::vector ref_mapaxes_egrid = { 0.0, 100.0, 0.0, 0.0, 100.0, 0.0 }; + std::vector ref_mapaxes_deck = { 0.0, 200.0, 0.0, 0.0, 200.0, 0.0 }; + + + // egrid file in si units, no conversion requied by grid constructor + std::vector refDepthGrid3a = {101.42292, 101.90941, 102.30995, 102.84644, 106.25719, 106.60616, 106.43174, 106.78071 }; + + // egrid file in field units, depths converted to SI units when loaded by grid constructor + std::vector refDepthGrid3b = {30.913707, 31.061988, 31.184072, 31.347594, 32.38719, 32.493558, 32.440393 }; + + auto deck1a = parser.parseString( deckData1a) ; + + Opm::EclipseState es1a( deck1a ); + Opm::UnitSystem units1a = es1a.getDeckUnitSystem(); + + const auto& grid1a = es1a.getInputGrid(); + Opm::NNC nnc( deck1a ); + + grid1a.save("BAD_CP_M.EGRID", false, nnc, units1a); + + auto deck1b = parser.parseString( deckData1b) ; + Opm::EclipseState es1b( deck1b ); + Opm::UnitSystem units1b = es1b.getDeckUnitSystem(); + const auto& grid1b = es1b.getInputGrid(); + + grid1b.save("BAD_CP_F.EGRID", false, nnc, units1b); + + auto deck1 = parser.parseString( deckData1) ; + Opm::EclipseGrid grid1( deck1); + + Opm::EclIO::EclFile file1("BAD_CP_M.EGRID"); + + // actnum not defined in deck. keyword GDFILE not present in the DECK + // check that coord and zcorn from deck-grid identical to coord and zcorn + // from egrid - grid + + std::vector coordGrid1 = grid1.getCOORD(); + std::vector zcornGrid1 = grid1.getZCORN(); + std::vector coordGrid1_f(coordGrid1.begin(), coordGrid1.end() ); + std::vector zcornGrid1_f(zcornGrid1.begin(), zcornGrid1.end() ); + + const std::vector coord_egrid_f = file1.get("COORD"); + const std::vector zcorn_egrid_f = file1.get("ZCORN"); + + BOOST_CHECK( coordGrid1.size() == coord_egrid_f.size() ); + BOOST_CHECK( zcornGrid1.size() == zcorn_egrid_f.size() ); + + for (size_t n = 0; n < coordGrid1.size(); n++){ + BOOST_CHECK( coordGrid1_f[n] == coord_egrid_f[n] ); + } + + for (size_t n = 0; n < zcornGrid1.size(); n++){ + BOOST_CHECK( zcornGrid1_f[n] == zcorn_egrid_f[n] ); + } + + // all cells are active, since ACTNUM not present + std::vector actGrid1 = grid1.getACTNUM(); + for (size_t n = 0; n < actGrid1.size(); n++){ + BOOST_CHECK( actGrid1[n] == 1 ); + } + + BOOST_CHECK( grid1.getMAPUNITS() == "" ); + + std::vector mapaxes = grid1.getMAPAXES(); + BOOST_CHECK( mapaxes.size() == 0 ); + + + auto deck2 = parser.parseString( deckData2) ; + Opm::EclipseGrid grid2( deck2); + + std::vector actGrid2 = grid2.getACTNUM(); + + // check that actnum is reset from gdfile + + for (size_t n = 0; n < actGrid2.size(); n++){ + BOOST_CHECK( actGrid2[n] == ref_act_egrid[n] ); + } + + BOOST_CHECK( grid2.getMAPUNITS() == "" ); + + mapaxes = grid2.getMAPAXES(); + BOOST_CHECK( mapaxes.size() == 0 ); + + + auto deck3a = parser.parseString( deckData3a) ; + Opm::EclipseGrid grid3a( deck3a); + + // mapunits and mapaxes define in deck (only) + + BOOST_CHECK( grid3a.getMAPUNITS() == "FEET" ); + + mapaxes = grid3a.getMAPAXES(); + BOOST_CHECK( mapaxes.size() == 6 ); + + for (size_t n = 0; n < mapaxes.size(); n++){ + BOOST_CHECK( mapaxes[n] == ref_mapaxes_deck[n] ); + } + + std::vector actGrid3 = grid3a.getACTNUM(); + + // check that actnum is reset from gdfile, ACTNUM input in deck + // but before keyword GDFILE + + for (size_t n = 0; n < actGrid3.size(); n++){ + BOOST_CHECK( actGrid3[n] == ref_act_egrid[n] ); + } + + // check that depth values are in SI units + for (size_t n = 0; n < refDepthGrid3a.size(); n++){ + BOOST_CHECK_CLOSE( grid3a.getCellDepth(n), refDepthGrid3a[n], 1e-3 ); + } + + auto deck3b = parser.parseString( deckData3b) ; + Opm::EclipseGrid grid3b( deck3b); + + // mapunits and mapaxes both in egrid and deck. Uses properties + // from the egrid keyword gdfile input after MAPUNITS and MAPAXES + + BOOST_CHECK( grid3b.getMAPUNITS() == "METRES" ); + + mapaxes = grid3b.getMAPAXES(); + BOOST_CHECK( mapaxes.size() == 6 ); + + for (size_t n = 0; n < mapaxes.size(); n++){ + BOOST_CHECK( mapaxes[n] == ref_mapaxes_egrid[n] ); + } + + actGrid3 = grid3b.getACTNUM(); + + // check that actnum is reset from deck since input after keyword GDFILE + for (size_t n = 0; n < actGrid3.size(); n++){ + BOOST_CHECK( actGrid3[n] == ref_act_deck3[n] ); + } + + // check that depth values are converted from Field to SI units + for (size_t n = 0; n < refDepthGrid3b.size(); n++){ + BOOST_CHECK_CLOSE( grid3b.getCellDepth(n), refDepthGrid3b[n], 1e-3 ); + } + + // mapunits and mapaxes both in egrid and deck. Uses properties + // from the deck sinze these are input after GDfile + + auto deck3c = parser.parseString( deckData3c) ; + Opm::EclipseGrid grid3c( deck3c); + + BOOST_CHECK( grid3c.getMAPUNITS() == "FEET" ); + + mapaxes = grid3c.getMAPAXES(); + BOOST_CHECK( mapaxes.size() == 6 ); + + for (size_t n = 0; n < mapaxes.size(); n++){ + BOOST_CHECK( mapaxes[n] == ref_mapaxes_deck[n] ); + } +} + diff --git a/tests/parser/ScheduleTests.cpp b/tests/parser/ScheduleTests.cpp index 19d2f6927..22c20ad53 100644 --- a/tests/parser/ScheduleTests.cpp +++ b/tests/parser/ScheduleTests.cpp @@ -2960,7 +2960,6 @@ BOOST_AUTO_TEST_CASE(FilterCompletions2) { - BOOST_AUTO_TEST_CASE(VFPINJ_TEST) { const char *deckData = "\ START\n \ diff --git a/tests/parser/integration/EclipseGridCreateFromDeck.cpp b/tests/parser/integration/EclipseGridCreateFromDeck.cpp index f3366772f..ba330ce1e 100644 --- a/tests/parser/integration/EclipseGridCreateFromDeck.cpp +++ b/tests/parser/integration/EclipseGridCreateFromDeck.cpp @@ -75,9 +75,8 @@ BOOST_AUTO_TEST_CASE(ExportFromCPGridAllActive) { std::vector actnum; - actnum.push_back(100); - grid.exportACTNUM( actnum ); - BOOST_CHECK_EQUAL( actnum.size() , 0U ); + actnum = grid.getACTNUM(); + BOOST_CHECK_EQUAL( actnum.size() , 500U ); } @@ -95,13 +94,13 @@ BOOST_AUTO_TEST_CASE(ExportFromCPGridACTNUM) { std::vector actnum; size_t volume = grid.getNX()*grid.getNY()*grid.getNZ(); - grid.exportCOORD( coord ); + coord = grid.getCOORD(); BOOST_CHECK_EQUAL( coord.size() , (grid.getNX() + 1) * (grid.getNY() + 1) * 6); - grid.exportZCORN( zcorn ); + zcorn = grid.getZCORN(); BOOST_CHECK_EQUAL( zcorn.size() , volume * 8); - grid.exportACTNUM( actnum ); + actnum = grid.getACTNUM(); BOOST_CHECK_EQUAL( actnum.size() , volume ); { diff --git a/tests/test_AggregateConnectionData.cpp b/tests/test_AggregateConnectionData.cpp index f96be2e39..4b84020cd 100644 --- a/tests/test_AggregateConnectionData.cpp +++ b/tests/test_AggregateConnectionData.cpp @@ -723,7 +723,7 @@ BOOST_AUTO_TEST_CASE(InactiveCell) { std::vector actnum(500, 1); // Here we deactive the cell holding connection number 2. actnum[simCase.grid.getGlobalIndex(2,4,1)] = 0; - simCase.grid.resetACTNUM(actnum.data()); + simCase.grid.resetACTNUM(actnum); auto conn1 = Opm::RestartIO::Helpers::AggregateConnectionData{ih.value}; conn1.captureDeclaredConnData(simCase.sched, diff --git a/tests/test_EclIO.cpp b/tests/test_EclIO.cpp index 39cad9576..07860c8c8 100644 --- a/tests/test_EclIO.cpp +++ b/tests/test_EclIO.cpp @@ -73,7 +73,7 @@ BOOST_AUTO_TEST_CASE(TestEclFile_BINARY) { // check that exception is thrown when input file doesn't exist - BOOST_CHECK_THROW(EclFile file1("DUMMY.DAT") , std::runtime_error ); + BOOST_CHECK_THROW(EclFile file1("DUMMY.DAT") , std::invalid_argument ); EclFile file1(testFile); @@ -331,4 +331,4 @@ BOOST_AUTO_TEST_CASE(TestEcl_Write_CHAR) { }; -} \ No newline at end of file +} diff --git a/tests/test_EclipseIO.cpp b/tests/test_EclipseIO.cpp index 8fb402829..25ec18052 100644 --- a/tests/test_EclipseIO.cpp +++ b/tests/test_EclipseIO.cpp @@ -120,19 +120,19 @@ void checkEgridFile( const EclipseGrid& eclGrid ) { std::string keywordName(ecl_kw_get_header(eclKeyword)); if (keywordName == "COORD") { std::vector< double > sourceData; - eclGrid.exportCOORD( sourceData ); + sourceData = eclGrid.getCOORD(); auto resultData = getErtData< float >( eclKeyword ); compareErtData(sourceData, resultData, 1e-6); } else if (keywordName == "ZCORN") { std::vector< double > sourceData; - eclGrid.exportZCORN(sourceData); + sourceData = eclGrid.getZCORN(); auto resultData = getErtData< float >( eclKeyword ); compareErtData(sourceData, resultData, /*percentTolerance=*/1e-6); } else if (keywordName == "ACTNUM") { std::vector< int > sourceData( numCells ); - eclGrid.exportACTNUM(sourceData); + sourceData = eclGrid.getACTNUM(); auto resultData = getErtData< int >( eclKeyword ); if( sourceData.empty() ) diff --git a/tests/test_calculateCellVol.cpp b/tests/test_calculateCellVol.cpp index 364f00367..101bb99e3 100755 --- a/tests/test_calculateCellVol.cpp +++ b/tests/test_calculateCellVol.cpp @@ -38,21 +38,21 @@ BOOST_AUTO_TEST_CASE (calc_cellvol) /* This is four example cells with different geometries. */ - std::vector x1 {488100.140035, 488196.664549, 488085.584866, 488182.365605, 488099.065709, 488195.880889, 488084.559409, 488181.633495}; - std::vector y1 {6692539.945578, 6692550.834909, 6692638.574346, 6692650.086244, 6692538.810649, 6692550.080826, 6692637.628127, 6692649.429649}; - std::vector z1 {2841.856000, 2840.138000, 2842.042000, 2839.816000, 2846.142000, 2844.252000, 2846.244000, 2843.868000}; + std::array x1 {488100.140035, 488196.664549, 488085.584866, 488182.365605, 488099.065709, 488195.880889, 488084.559409, 488181.633495}; + std::array y1 {6692539.945578, 6692550.834909, 6692638.574346, 6692650.086244, 6692538.810649, 6692550.080826, 6692637.628127, 6692649.429649}; + std::array z1 {2841.856000, 2840.138000, 2842.042000, 2839.816000, 2846.142000, 2844.252000, 2846.244000, 2843.868000}; - std::vector x2 {489539.050892, 489638.562319, 489527.025803, 489627.914272, 489537.173839, 489638.562319, 489524.119410, 489627.914272}; - std::vector y2 {6695907.426615, 6695923.399538, 6696003.688945, 6696020.073168, 6695908.526577, 6695923.399538, 6696005.533276, 6696020.073168}; - std::vector z2 {2652.859000, 2651.765000, 2652.381000, 2652.608000, 2655.276000, 2651.765000, 2656.381000, 2652.608000}; + std::array x2 {489539.050892, 489638.562319, 489527.025803, 489627.914272, 489537.173839, 489638.562319, 489524.119410, 489627.914272}; + std::array y2 {6695907.426615, 6695923.399538, 6696003.688945, 6696020.073168, 6695908.526577, 6695923.399538, 6696005.533276, 6696020.073168}; + std::array z2 {2652.859000, 2651.765000, 2652.381000, 2652.608000, 2655.276000, 2651.765000, 2656.381000, 2652.608000}; - std::vector x3 {486819.009056, 486904.877179, 486812.975440, 486900.527416, 486818.450563, 486903.978383, 486812.975440, 486900.527416}; - std::vector y3 {6694966.253697, 6694998.360834, 6695061.104896, 6695086.717018, 6694967.135567, 6695000.197284, 6695061.104896, 6695086.717018}; - std::vector z3 {2795.193000, 2799.997000, 2792.240000, 2793.972000, 2795.671000, 2800.927000, 2792.240000, 2793.972000}; + std::array x3 {486819.009056, 486904.877179, 486812.975440, 486900.527416, 486818.450563, 486903.978383, 486812.975440, 486900.527416}; + std::array y3 {6694966.253697, 6694998.360834, 6695061.104896, 6695086.717018, 6694967.135567, 6695000.197284, 6695061.104896, 6695086.717018}; + std::array z3 {2795.193000, 2799.997000, 2792.240000, 2793.972000, 2795.671000, 2800.927000, 2792.240000, 2793.972000}; - std::vector x4 {489194.711544, 489259.232449, 489206.220219, 489294.099099, 489194.711544, 489254.725623, 489201.723535, 489289.423088}; - std::vector y4 {6694510.504054, 6694525.862296, 6694608.410134, 6694621.029382, 6694510.504054, 6694526.272988, 6694608.819829, 6694621.467114}; - std::vector z4 {2740.7340, 2754.7810, 2730.0150, 2726.1080, 2740.7340, 2758.3530, 2733.9490, 2730.1080}; + std::array x4 {489194.711544, 489259.232449, 489206.220219, 489294.099099, 489194.711544, 489254.725623, 489201.723535, 489289.423088}; + std::array y4 {6694510.504054, 6694525.862296, 6694608.410134, 6694621.029382, 6694510.504054, 6694526.272988, 6694608.819829, 6694621.467114}; + std::array z4 {2740.7340, 2754.7810, 2730.0150, 2726.1080, 2740.7340, 2758.3530, 2733.9490, 2730.1080}; BOOST_REQUIRE_CLOSE (calculateCellVol(x1,y1,z1), 40368.7852157, 1e-9); BOOST_REQUIRE_CLOSE (calculateCellVol(x2,y2,z2), 15766.9187847524, 1e-9); diff --git a/tests/test_restartwellinfo.cpp b/tests/test_restartwellinfo.cpp new file mode 100644 index 000000000..c0f40539c --- /dev/null +++ b/tests/test_restartwellinfo.cpp @@ -0,0 +1,234 @@ + +/* + Copyright 2014 Statoil IT + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#include "config.h" + +#define BOOST_TEST_MODULE EclipseWriter +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include + +#include + + +void verifyWellState(const std::string& rst_filename, const Opm::Schedule& schedule) { + + /* + reference variables ref_wellList, ref_wellHead and ref_wellConn are + based on WELSPECS and COMPDAT in data deck -> testblackoilstate3.DATA. + + */ + + int step = std::stoi( rst_filename.substr(rst_filename.size()-1)) -1; + Opm::EclIO::EclFile rstFile(rst_filename); + + std::vector> ref_wellList = { + {}, + {"OP_1", "OP_2", "OP_3"}, + {"OP_1", "OP_2", "OP_3", "OP_4", "OP_5", "OP_6"}, + {"OP_1", "OP_2", "OP_3", "OP_4", "OP_5", "OP_6", "WI_1", "WI_2", "GI_1"} + }; + + std::vector>> ref_wellHead = { + {}, + {{9,9}, {8,8}, {7,7}}, + {{9,9}, {8,8}, {7,7}, {2,2}, {5,4}, {8,2}}, + {{9,9}, {8,8}, {7,7}, {2,2}, {5,4}, {8,2}, {3,3}, {3,9}, {3,6}} + }; + + std::vector>>> ref_wellConn = { + {{{}}}, // <- timestep 0 + { {{9,9,1},{9,9,2},{9,9,3},{9,9,4},{9,9,5},{9,9,6},{9,9,7},{9,9,8},{9,9,9},{9,9,10}}, // OP_1 + {{8,8,1},{8,8,2},{8,8,3},{8,7,3},{8,7,4},{8,7,5},{8,7,6}}, // OP_2 + {{7,7,1},{7,7,2},{7,7,3},{7,7,4},{7,7,5}} // OP_3 + }, // <- timestep 1 + { {{9,9,1},{9,9,2},{9,9,3},{9,9,4},{9,9,5},{9,9,6},{9,9,7},{9,9,8},{9,9,9},{9,9,10}}, + {{8,8,1},{8,8,2},{8,8,3},{8,7,3},{8,7,4},{8,7,5},{8,7,6}}, // OP_2 + {{7,7,1},{7,7,2},{7,7,3},{7,7,4},{7,7,5}}, // OP_3 + {{2,2,1},{2,2,2},{2,2,3},{2,2,4},{2,2,5},{2,2,6},{2,2,7},{2,2,8},{2,2,9},{2,2,10}}, // OP_4 + {{5,4,1},{5,4,2},{5,4,3}}, // OP_5 + {{8,2,1},{8,2,2},{8,2,3},{8,3,3},{8,4,3},{8,5,3},{8,5,4},{8,5,5},{8,5,6}} // OP_6 + }, // <- timestep 2 + { {{9,9,1},{9,9,2},{9,9,3},{9,9,4},{9,9,5},{9,9,6},{9,9,7},{9,9,8},{9,9,9},{9,9,10}}, + {{8,8,1},{8,8,2},{8,8,3},{8,7,3},{8,7,4},{8,7,5},{8,7,6}}, // OP_2 + {{7,7,1},{7,7,2},{7,7,3},{7,7,4},{7,7,5}}, // OP_3 + {{2,2,1},{2,2,2},{2,2,3},{2,2,4},{2,2,5},{2,2,6},{2,2,7},{2,2,8},{2,2,9},{2,2,10}}, // OP_4 + {{5,4,1},{5,4,2},{5,4,3}}, // OP_5 + {{8,2,1},{8,2,2},{8,2,3},{8,3,3},{8,4,3},{8,5,3},{8,5,4},{8,5,5},{8,5,6}}, // OP_6 + {{3,3,1},{3,3,2},{3,3,3},{3,3,4},{3,3,5},{3,3,6},{3,3,7},{3,3,8},{3,3,9},{3,3,10}}, // WI_1 + {{3,9,1},{3,9,2},{3,9,3},{3,9,4},{3,9,5},{3,9,6},{3,9,7}}, // WI_2 + {{3,6,1},{3,6,2},{3,6,3}} // WI_3 + } // <- timestep 3 + + }; + + std::vector intehead = rstFile.get("INTEHEAD"); + std::vector zwel; + std::vector iwel; + std::vector icon; + + int ncwmax = intehead[17]; + int niwelz = intehead[24]; + int niconz = intehead[32]; + + if (rstFile.hasKey("ZWEL")) { + zwel = rstFile.get("ZWEL"); + } + + if (rstFile.hasKey("IWEL")) { + iwel = rstFile.get("IWEL"); + } + + if (rstFile.hasKey("ICON")) { + icon = rstFile.get("ICON"); + } + + std::vector wellList = schedule.wellNames(step); + + //Verify number of active wells + BOOST_CHECK_EQUAL( wellList.size(), intehead[16]); + + for (size_t i=0; i< wellList.size(); i++) { + + // Verify wellname + BOOST_CHECK_EQUAL(zwel[i*3], ref_wellList[step][i]); + BOOST_CHECK_EQUAL(zwel[i*3], wellList[i]); + + // Verify well I, J head + + BOOST_CHECK_EQUAL(iwel[i*niwelz], std::get<0>(ref_wellHead[step][i])); + BOOST_CHECK_EQUAL(iwel[i*niwelz + 1], std::get<1>(ref_wellHead[step][i])); + + Opm::Well2 sched_well2 = schedule.getWell2(wellList[i], step); + + BOOST_CHECK_EQUAL(iwel[i*niwelz], sched_well2.getHeadI() +1 ); + BOOST_CHECK_EQUAL(iwel[i*niwelz + 1], sched_well2.getHeadJ() +1 ); + + int sched_wtype = -99; + + if (sched_well2.isProducer()) { + sched_wtype = 1; + } else { + switch( sched_well2.getInjectionProperties( ).injectorType ) { + case Opm::Well2::InjectorType::WATER: + sched_wtype = 3; + break; + case Opm::Well2::InjectorType::GAS: + sched_wtype = 4; + break; + case Opm::Well2::InjectorType::OIL: + sched_wtype = 2; + break; + default: + break; + } + } + + // Verify well type + // 1 = producer, 2 = oil injection, 3 = water injector, 4 = gas injector + + BOOST_CHECK_EQUAL(iwel[i*niwelz + 6], sched_wtype ); + + const auto& connections_set = sched_well2.getConnections(); + + // Verify number of connections + + BOOST_CHECK_EQUAL(iwel[i*niwelz + 4], connections_set.size() ); + BOOST_CHECK_EQUAL(ref_wellConn[step][i].size(), connections_set.size() ); + + + for (size_t n=0; n< connections_set.size(); n++) { + const auto& completion = connections_set.get(n); + + // Verify I, J and K indices for each connection + + BOOST_CHECK_EQUAL(completion.getI()+1 , std::get<0>(ref_wellConn[step][i][n])); + BOOST_CHECK_EQUAL(completion.getJ()+1 , std::get<1>(ref_wellConn[step][i][n])); + BOOST_CHECK_EQUAL(completion.getK()+1 , std::get<2>(ref_wellConn[step][i][n])); + + size_t ind = i*ncwmax*niconz+n*niconz; + + BOOST_CHECK_EQUAL(completion.getI()+1 , icon[ind+1]); + BOOST_CHECK_EQUAL(completion.getJ()+1 , icon[ind+2]); + BOOST_CHECK_EQUAL(completion.getK()+1 , icon[ind+3]); + } + } +} + + +BOOST_AUTO_TEST_CASE(EclipseWriteRestartWellInfo) { + + std::string eclipse_data_filename = "testblackoilstate3.DATA"; + + Opm::Parser parser; + Opm::Deck deck( parser.parseFile( eclipse_data_filename )); + Opm::EclipseState es(deck); + const Opm::EclipseGrid& grid = es.getInputGrid(); + Opm::Schedule schedule( deck, es); + Opm::SummaryConfig summary_config( deck, schedule, es.getTableManager( )); + const auto num_cells = grid.getCartesianSize(); + Opm::EclipseIO eclipseWriter( es, grid , schedule, summary_config); + int countTimeStep = schedule.getTimeMap().numTimesteps(); + Opm::SummaryState st(std::chrono::system_clock::from_time_t(schedule.getStartTime())); + + Opm::data::Solution solution; + solution.insert( "PRESSURE", Opm::UnitSystem::measure::pressure , std::vector< double >( num_cells, 1 ) , Opm::data::TargetType::RESTART_SOLUTION); + solution.insert( "SWAT" , Opm::UnitSystem::measure::identity , std::vector< double >( num_cells, 1 ) , Opm::data::TargetType::RESTART_SOLUTION); + solution.insert( "SGAS" , Opm::UnitSystem::measure::identity , std::vector< double >( num_cells, 1 ) , Opm::data::TargetType::RESTART_SOLUTION); + Opm::data::Wells wells; + + for(int timestep = 0; timestep <= countTimeStep; ++timestep) { + + eclipseWriter.writeTimeStep( st, + timestep, + false, + schedule.seconds(timestep), + Opm::RestartValue(solution, wells)); + } + + for (int i=1; i <=4; i++) { + verifyWellState("TESTBLACKOILSTATE3.X000" + std::to_string(i), schedule); + } + + // cleaning up after test + + for (int i=1; i <=4; i++) { + std::string fileName = "TESTBLACKOILSTATE3.X000" + std::to_string(i); + remove(fileName.c_str()); + + fileName = "testblackoilstate3.s000" + std::to_string(i); + remove(fileName.c_str()); + } + + remove("testblackoilstate3.smspec"); +} diff --git a/tests/test_writenumwells.cpp b/tests/test_writenumwells.cpp deleted file mode 100644 index f95c50282..000000000 --- a/tests/test_writenumwells.cpp +++ /dev/null @@ -1,173 +0,0 @@ - -/* - Copyright 2014 Statoil IT - This file is part of the Open Porous Media project (OPM). - - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ -#include "config.h" - -#define BOOST_TEST_MODULE EclipseWriter -#include -#include -#include - -#include -#include -#include -#include -#include -#include -#include -#include - -#include -#include - -// ERT stuff -#include -#include -#include -#include -#include -#include - -using namespace Opm; - -void verifyWellState(const std::string& rst_filename, - const EclipseGrid& ecl_grid, - const Schedule& schedule) { - - well_info_type* well_info = well_info_alloc(ecl_grid.c_ptr()); - well_info_load_rstfile(well_info, rst_filename.c_str(), false); - - //Verify numwells - int numwells = well_info_get_num_wells(well_info); - BOOST_CHECK_EQUAL( numwells, schedule.numWells() ); - - auto wells = schedule.getWells2atEnd(); - - for (int i = 0; i < numwells; ++i) { - - //Verify wellnames - const char * wellname = well_info_iget_well_name(well_info, i); - well_ts_type* well_ts = well_info_get_ts(well_info , wellname); - well_state_type* well_state = well_ts_iget_state(well_ts, 0); - { - const auto& sched_well = wells.at(i); - BOOST_CHECK_EQUAL( wellname, sched_well.name() ); - - // Verify well-head position data - const well_conn_type* well_head = well_state_get_wellhead(well_state, ECL_GRID_GLOBAL_GRID); - BOOST_CHECK_EQUAL(well_conn_get_i(well_head), sched_well.getHeadI()); - BOOST_CHECK_EQUAL(well_conn_get_j(well_head), sched_well.getHeadJ()); - } - - for (int j = 0; j < well_ts_get_size(well_ts); ++j) { - const auto& well_at_end = schedule.getWell2atEnd(wellname); - if (!well_at_end.hasBeenDefined(j)) - continue; - - const auto& well = schedule.getWell2(wellname, j); - well_state = well_ts_iget_state(well_ts, j); - - //Verify welltype - int well_type = ERT_UNDOCUMENTED_ZERO; - if( well.isProducer( ) ) { - well_type = ERT_PRODUCER; - } - else { - switch( well.getInjectionProperties( ).injectorType ) { - case Well2::InjectorType::WATER: - well_type = ERT_WATER_INJECTOR; - break; - case Well2::InjectorType::GAS: - well_type = ERT_GAS_INJECTOR; - break; - case Well2::InjectorType::OIL: - well_type = ERT_OIL_INJECTOR; - break; - default: - break; - } - } - - int ert_well_type = well_state_get_type( well_state ); - BOOST_CHECK_EQUAL( ert_well_type, well_type ); - - //Verify wellstatus - int ert_well_status = well_state_is_open(well_state) ? 1 : 0; - int wellstatus = well.getStatus( ) == Well2::Status::OPEN ? 1 : 0; - - BOOST_CHECK_EQUAL(ert_well_status, wellstatus); - - //Verify number of completion connections - const well_conn_collection_type * well_connections = well_state_get_global_connections( well_state ); - size_t num_wellconnections = well_conn_collection_get_size(well_connections); - - const auto& connections_set = well.getConnections(); - - BOOST_CHECK_EQUAL(num_wellconnections, connections_set.size()); - - //Verify coordinates for each completion connection - for (size_t k = 0; k < num_wellconnections; ++k) { - const well_conn_type * well_connection = well_conn_collection_iget_const(well_connections , k); - - const auto& completion = connections_set.get(k); - - BOOST_CHECK_EQUAL(well_conn_get_i(well_connection), completion.getI()); - BOOST_CHECK_EQUAL(well_conn_get_j(well_connection), completion.getJ()); - BOOST_CHECK_EQUAL(well_conn_get_k(well_connection), completion.getK()); - } - } - } - - well_info_free(well_info); -} - -BOOST_AUTO_TEST_CASE(EclipseWriteRestartWellInfo) { - std::string eclipse_data_filename = "testblackoilstate3.DATA"; - std::string eclipse_restart_filename = "TESTBLACKOILSTATE3.X0004"; - test_work_area_type * test_area = test_work_area_alloc("TEST_EclipseWriteNumWells"); - test_work_area_copy_file(test_area, eclipse_data_filename.c_str()); - - Parser parser; - Deck deck( parser.parseFile( eclipse_data_filename )); - EclipseState es(deck); - const EclipseGrid& grid = es.getInputGrid(); - Schedule schedule( deck, es); - SummaryConfig summary_config( deck, schedule, es.getTableManager( )); - const auto num_cells = grid.getCartesianSize(); - EclipseIO eclipseWriter( es, grid , schedule, summary_config); - int countTimeStep = schedule.getTimeMap().numTimesteps(); - SummaryState st(std::chrono::system_clock::now()); - - data::Solution solution; - solution.insert( "PRESSURE",UnitSystem::measure::pressure , std::vector< double >( num_cells, 1 ) , data::TargetType::RESTART_SOLUTION); - solution.insert( "SWAT" ,UnitSystem::measure::identity , std::vector< double >( num_cells, 1 ) , data::TargetType::RESTART_SOLUTION); - solution.insert( "SGAS" ,UnitSystem::measure::identity , std::vector< double >( num_cells, 1 ) , data::TargetType::RESTART_SOLUTION); - data::Wells wells; - - for(int timestep = 0; timestep <= countTimeStep; ++timestep){ - eclipseWriter.writeTimeStep( st, - timestep, - false, - timestep, - RestartValue(solution, wells)); - } - - verifyWellState(eclipse_restart_filename, grid, schedule); - - test_work_area_free(test_area); -} diff --git a/tests/testblackoilstate3.DATA b/tests/testblackoilstate3.DATA index 84a53fa85..6a408aff2 100644 --- a/tests/testblackoilstate3.DATA +++ b/tests/testblackoilstate3.DATA @@ -79,6 +79,13 @@ COMPDAT 'OP_3' 7 7 5 5 'OPEN' 1* 4.728 0.311 420.067 1* 1* 'Y' 18.162 / / +WCONPROD + 'OP_1' OPEN ORAT 1000.0 4* 90.0 / + 'OP_2' OPEN ORAT 1000.0 4* 90.0 / + 'OP_3' SHUT ORAT 1000.0 4* 90.0 / +/ + + DATES 1 JUN 1980/ / @@ -100,24 +107,34 @@ COMPDAT 'OP_6' 8 5 5 6 'OPEN' 1* 4.728 0.311 420.067 1* 1* 'Y' 18.162 / / +WCONPROD + 'OP_4' OPEN ORAT 1000.0 4* 90.0 / + 'OP_5' OPEN ORAT 1000.0 4* 90.0 / + 'OP_6' OPEN ORAT 1000.0 4* 90.0 / +/ + DATES 1 NOV 1980/ / + WELSPECS 'WI_1' 'WI' 3 3 1* 'OIL' 1* 1* 1* 1* 1* 1* 1* / - 'WI_2' 'WI' 3 9 1* 'OIL' 1* 1* 1* 1* 1* 1* 1* / - 'WI_3' 'WI' 3 6 1* 'OIL' 1* 1* 1* 1* 1* 1* 1* / + 'WI_2' 'WI' 3 9 1* 'OIL' 1* 1* 1* 1* 1* 1* 1* / + 'GI_1' 'WI' 3 6 1* 'GAS' 1* 1* 1* 1* 1* 1* 1* / / COMPDAT 'WI_1' 3 3 1 10 'OPEN' 1* 6.642 0.311 615.914 1* 1* 'X' 22.372 / 'WI_2' 3 9 1 7 'OPEN' 1* 1.168 0.311 107.872 1* 1* 'Y' 21.925 / - 'WI_3' 3 6 1 3 'OPEN' 1* 27.412 0.311 2445.337 1* 1* 'Y' 18.521 / + 'GI_1' 3 6 1 3 'OPEN' 1* 27.412 0.311 2445.337 1* 1* 'Y' 18.521 / / - - +WCONINJE + 'WI_1' WATER OPEN RATE 1000.0 / + 'WI_2' WATER OPEN RATE 1000.0 / + 'GI_1' GAS OPEN RATE 1.0E6 / +/ DATES 1 NOV 1982/