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/