diff --git a/opm/parser/eclipse/EclipseState/Grid/EclipseGrid.hpp b/opm/parser/eclipse/EclipseState/Grid/EclipseGrid.hpp index 472ef0771..3f59f4897 100644 --- a/opm/parser/eclipse/EclipseState/Grid/EclipseGrid.hpp +++ b/opm/parser/eclipse/EclipseState/Grid/EclipseGrid.hpp @@ -78,6 +78,7 @@ namespace Opm { /// explicitly. If a null pointer is passed, every cell is active. EclipseGrid(const Deck& deck, const int * actnum = nullptr); + static bool hasCylindricalKeywords(const Deck& deck); static bool hasCornerPointKeywords(const Deck&); static bool hasCartesianKeywords(const Deck&); size_t getNumActive( ) const; @@ -210,6 +211,7 @@ namespace Opm { const int * actnum, const double * mapaxes); + void initCylindricalGrid( const std::array&, const Deck&); void initCartesianGrid( const std::array&, const Deck&); void initCornerPointGrid( const std::array&, const Deck&); void initDTOPSGrid( const std::array&, const Deck&); @@ -227,11 +229,31 @@ namespace Opm { static void scatterDim(const std::array& dims , size_t dim , const std::vector& DV , std::vector& D); }; + class CoordMapper { + public: + CoordMapper(size_t nx, size_t ny); + size_t size() const; + + + /* + dim = 0,1,2 for x, y and z coordinate respectively. + layer = 0,1 for k=0 and k=nz layers respectively. + */ + + size_t index(size_t i, size_t j, size_t dim, size_t layer) const; + private: + size_t nx; + size_t ny; + }; + + + class ZcornMapper { public: ZcornMapper(size_t nx, size_t ny, size_t nz); size_t index(size_t i, size_t j, size_t k, int c) const; size_t index(size_t g, int c) const; + size_t size() const; /* The fixupZCORN method will take a zcorn vector as input and diff --git a/src/opm/parser/eclipse/EclipseState/Grid/EclipseGrid.cpp b/src/opm/parser/eclipse/EclipseState/Grid/EclipseGrid.cpp index 787b57163..22282ffc3 100644 --- a/src/opm/parser/eclipse/EclipseState/Grid/EclipseGrid.cpp +++ b/src/opm/parser/eclipse/EclipseState/Grid/EclipseGrid.cpp @@ -17,12 +17,12 @@ along with OPM. If not, see . */ +#define _USE_MATH_DEFINES +#include #include #include -#include - -#include +#include #include #include @@ -35,8 +35,10 @@ #include #include #include +#include #include #include +#include #include #include #include @@ -186,12 +188,16 @@ namespace Opm { } void EclipseGrid::initGrid( const std::array& dims, const Deck& deck) { - if (hasCornerPointKeywords(deck)) { - initCornerPointGrid(dims , deck); - } else if (hasCartesianKeywords(deck)) { - initCartesianGrid(dims , deck); + if (deck.hasKeyword()) { + initCylindricalGrid( dims, deck ); } else { - throw std::invalid_argument("EclipseGrid needs cornerpoint or cartesian keywords."); + if (hasCornerPointKeywords(deck)) { + initCornerPointGrid(dims , deck); + } else if (hasCartesianKeywords(deck)) { + initCartesianGrid(dims , deck); + } else { + throw std::invalid_argument("EclipseGrid needs cornerpoint or cartesian keywords."); + } } if (deck.hasKeyword()) { @@ -311,6 +317,126 @@ namespace Opm { } + /* + Limited implementaton - requires keywords: DRV, DTHETAV, DZV and TOPS. + */ + + void EclipseGrid::initCylindricalGrid(const std::array& dims , const Deck& deck) + { + // The hasCyindricalKeywords( ) checks according to the + // eclipse specification. We currently do not support all + // aspects of cylindrical grids, we therefor have an + // additional test here, which checks if we have the keywords + // required by the current implementation. + if (!hasCylindricalKeywords(deck)) + throw std::invalid_argument("Not all keywords required for cylindrical grids present"); + + if (!deck.hasKeyword()) + throw std::logic_error("The current implementation *must* have theta values specified using the DTHETAV keyword"); + + if (!deck.hasKeyword()) + throw std::logic_error("The current implementation *must* have radial values specified using the DRV keyword"); + + if (!deck.hasKeyword() || !deck.hasKeyword()) + throw std::logic_error("The current implementation *must* have vertical cell size specified using the DZV and TOPS keywords"); + + const std::vector& drv = deck.getKeyword().getSIDoubleData(); + const std::vector& dthetav = deck.getKeyword().getSIDoubleData(); + const std::vector& dzv = deck.getKeyword().getSIDoubleData(); + const std::vector& tops = deck.getKeyword().getSIDoubleData(); + + if (drv.size() != dims[0]) + throw std::invalid_argument("DRV keyword should have exactly " + std::to_string( dims[0] ) + " elements"); + + if (dthetav.size() != dims[1]) + throw std::invalid_argument("DTHETAV keyword should have exactly " + std::to_string( dims[1] ) + " elements"); + + if (dzv.size() != dims[2]) + throw std::invalid_argument("DZV keyword should have exactly " + std::to_string( dims[2] ) + " elements"); + + if (tops.size() != dims[0] * dims[1]) + throw std::invalid_argument("TOPS keyword should have exactly " + std::to_string( dims[0] * dims[1] ) + " elements"); + + { + double total_angle = 0; + for (auto theta : dthetav) + total_angle += theta; + + if (std::abs( total_angle - 360 ) < 0.01) + m_circle = deck.hasKeyword(); + else { + if (total_angle > 360) + throw std::invalid_argument("More than 360 degrees rotation - cells will be double covered"); + } + } + + /* + Now the data has been validated, now we continue to create + ZCORN and COORD vectors, and we are done. + */ + { + ZcornMapper zm( dims[0], dims[1] , dims[2]); + CoordMapper cm(dims[0] , dims[1]); + std::vector zcorn( zm.size() ); + std::vector coord( cm.size() ); + { + std::vector zk(dims[2]); + zk[0] = 0; + for (size_t k = 1; k < dims[2]; k++) + zk[k] = zk[k - 1] + dzv[k - 1]; + + for (size_t k= 0; k < dims[2]; k++) { + for (size_t j=0; j < dims[1]; j++) { + for (size_t i = 0; i < dims[0]; i++) { + size_t tops_value = tops[ i + dims[0] * j]; + for (size_t c=0; c < 4; c++) { + zcorn[ zm.index(i,j,k,c) ] = zk[k] + tops_value; + zcorn[ zm.index(i,j,k,c + 4) ] = zk[k] + tops_value + dzv[k]; + } + } + } + } + } + { + std::vector ri(dims[0] + 1); + std::vector tj(dims[1] + 1); + double z1 = *std::min_element( zcorn.begin() , zcorn.end()); + double z2 = *std::max_element( zcorn.begin() , zcorn.end()); + ri[0] = deck.getKeyword().getRecord(0).getItem(0).getSIDouble( 0 ); + for (size_t i = 1; i <= dims[0]; i++) + ri[i] = ri[i - 1] + drv[i - 1]; + + tj[0] = 0; + for (size_t j = 1; j <= dims[1]; j++) + tj[j] = tj[j - 1] + dthetav[j - 1]; + + + for (size_t j = 0; j <= dims[1]; j++) { + /* + The theta value is supposed to go counterclockwise, starting at 'twelve o clock'. + */ + double t = M_PI * (90 - tj[j]) / 180; + double c = cos( t ); + double s = sin( t ); + for (size_t i=0; i <= dims[0]; i++) { + double r = ri[i]; + double x = r*c; + double y = r*s; + + coord[ cm.index(i,j,0,0) ] = x; + coord[ cm.index(i,j,1,0) ] = y; + coord[ cm.index(i,j,2,0) ] = z1; + + coord[ cm.index(i,j,0,1) ] = x; + coord[ cm.index(i,j,1,1) ] = y; + coord[ cm.index(i,j,2,1) ] = z2; + } + } + } + initCornerPointGrid( dims , coord, zcorn, nullptr, nullptr); + } + } + void EclipseGrid::initCornerPointGrid(const std::array& dims , const std::vector& coord , @@ -413,6 +539,17 @@ namespace Opm { return hasDTOPSKeywords(deck); } + bool EclipseGrid::hasCylindricalKeywords(const Deck& deck) { + if (deck.hasKeyword() && + deck.hasKeyword() && + (deck.hasKeyword() || deck.hasKeyword()) && + (deck.hasKeyword() || deck.hasKeyword()) && + (deck.hasKeyword() || deck.hasKeyword())) + return true; + else + return false; + } + bool EclipseGrid::hasDVDEPTHZKeywords(const Deck& deck) { if (deck.hasKeyword() && @@ -762,7 +899,9 @@ namespace Opm { return i*stride[0] + j*stride[1] + k*stride[2] + cell_shift[c]; } - + size_t ZcornMapper::size() const { + return dims[0] * dims[1] * dims[2] * 8; + } size_t ZcornMapper::index(size_t g, int c) const { int k = g / (dims[0] * dims[1]); @@ -838,6 +977,35 @@ 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; + } + + + size_t CoordMapper::index(size_t i, size_t j, size_t dim, size_t layer) const { + if (i > this->nx) + throw std::invalid_argument("Out of range"); + + if (j > this->ny) + throw std::invalid_argument("Out of range"); + + if (dim > 2) + throw std::invalid_argument("Out of range"); + + if (layer > 1) + throw std::invalid_argument("Out of range"); + + return 6*( i + j*(this->nx + 1) ) + layer * 3 + dim; + } } diff --git a/tests/parser/EclipseGridTests.cpp b/tests/parser/EclipseGridTests.cpp index 71beaea6d..1c768a141 100644 --- a/tests/parser/EclipseGridTests.cpp +++ b/tests/parser/EclipseGridTests.cpp @@ -1098,3 +1098,290 @@ BOOST_AUTO_TEST_CASE(MoveTest) { 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" + "\n" + "DIMENS\n" + " 10 10 10 /\n" + "RADIAL\n" + "\n"; + + Opm::Parser parser; + return parser.parseString( deckData, Opm::ParseContext()); +} + + + +static Opm::Deck radial_keywords_OK() { + const char* deckData = + "RUNSPEC\n" + "\n" + "DIMENS\n" + " 10 6 10 /\n" + "RADIAL\n" + "GRID\n" + "INRAD\n" + "1 /\n" + "DRV\n" + "10*1 /\n" + "DTHETAV\n" + "6*60 /\n" + "DZV\n" + "10*0.25 /\n" + "TOPS\n" + "60*0.0 /\n" + "\n"; + + Opm::Parser parser; + return parser.parseString( deckData, Opm::ParseContext()); +} + +static Opm::Deck radial_keywords_OK_CIRCLE() { + const char* deckData = + "RUNSPEC\n" + "\n" + "DIMENS\n" + " 10 6 10 /\n" + "RADIAL\n" + "GRID\n" + "CIRCLE\n" + "INRAD\n" + "1 /\n" + "DRV\n" + "10*1 /\n" + "DTHETAV\n" + "6*60 /\n" + "DZV\n" + "10*0.25 /\n" + "TOPS\n" + "60*0.0 /\n" + "\n"; + + Opm::Parser parser; + return parser.parseString( deckData, Opm::ParseContext()); +} + + + +BOOST_AUTO_TEST_CASE(RadialTest) { + Opm::Deck deck = radial_missing_INRAD(); + BOOST_CHECK_THROW( Opm::EclipseGrid{ deck }, std::invalid_argument); +} + + +BOOST_AUTO_TEST_CASE(RadialKeywordsOK) { + Opm::Deck deck = radial_keywords_OK(); + Opm::EclipseGrid grid( deck ); + BOOST_CHECK(!grid.circle()); +} + +BOOST_AUTO_TEST_CASE(RadialKeywordsOK_CIRCLE) { + Opm::Deck deck = radial_keywords_OK_CIRCLE(); + Opm::EclipseGrid grid( deck ); + BOOST_CHECK(grid.circle()); +} + + + +static Opm::Deck radial_keywords_DRV_size_mismatch() { + const char* deckData = + "RUNSPEC\n" + "\n" + "DIMENS\n" + "10 6 12 /\n" + "RADIAL\n" + "GRID\n" + "INRAD\n" + "1 /\n" + "DRV\n" + "9*1 /\n" + "DTHETAV\n" + "6*60 /\n" + "DZV\n" + "12*0.25 /\n" + "TOPS\n" + "60*0.0 /\n" + "\n"; + + Opm::Parser parser; + return parser.parseString( deckData, Opm::ParseContext()); +} + + +static Opm::Deck radial_keywords_DZV_size_mismatch() { + const char* deckData = + "RUNSPEC\n" + "\n" + "DIMENS\n" + "10 6 12 /\n" + "RADIAL\n" + "GRID\n" + "INRAD\n" + "1 /\n" + "DRV\n" + "10*1 /\n" + "DTHETAV\n" + "6*60 /\n" + "DZV\n" + "11*0.25 /\n" + "TOPS\n" + "60*0.0 /\n" + "\n"; + + Opm::Parser parser; + return parser.parseString( deckData, Opm::ParseContext()); +} + +static Opm::Deck radial_keywords_DTHETAV_size_mismatch() { + const char* deckData = + "RUNSPEC\n" + "\n" + "DIMENS\n" + "10 6 12 /\n" + "RADIAL\n" + "GRID\n" + "INRAD\n" + "1 /\n" + "DRV\n" + "10*1 /\n" + "DTHETAV\n" + "5*60 /\n" + "DZV\n" + "12*0.25 /\n" + "TOPS\n" + "60*0.0 /\n" + "\n"; + + Opm::Parser parser; + return parser.parseString( deckData, Opm::ParseContext()); +} + +/* + 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 = + "RUNSPEC\n" + "\n" + "DIMENS\n" + "10 6 12 /\n" + "RADIAL\n" + "GRID\n" + "INRAD\n" + "1 /\n" + "DRV\n" + "10*1 /\n" + "DTHETAV\n" + "6*60 /\n" + "DZV\n" + "12*0.25 /\n" + "TOPS\n" + "65*0.0 /\n" + "\n"; + + Opm::Parser parser; + return parser.parseString( deckData, Opm::ParseContext()); +} + + +static Opm::Deck radial_keywords_ANGLE_OVERFLOW() { + const char* deckData = + "RUNSPEC\n" + "\n" + "DIMENS\n" + "10 6 12 /\n" + "RADIAL\n" + "GRID\n" + "INRAD\n" + "1 /\n" + "DRV\n" + "10*1 /\n" + "DTHETAV\n" + "6*70 /\n" + "DZV\n" + "12*0.25 /\n" + "TOPS\n" + "60*0.0 /\n" + "\n"; + + Opm::Parser parser; + return parser.parseString( deckData, Opm::ParseContext()); +} + + + + +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); + BOOST_CHECK_THROW( Opm::EclipseGrid{ radial_keywords_TOPS_size_mismatch() } , std::invalid_argument); + BOOST_CHECK_THROW( Opm::EclipseGrid{ radial_keywords_DTHETAV_size_mismatch() } , std::invalid_argument); + BOOST_CHECK_THROW( Opm::EclipseGrid{ radial_keywords_ANGLE_OVERFLOW() } , std::invalid_argument); +} + + + +static Opm::Deck radial_details() { + const char* deckData = + "RUNSPEC\n" + "\n" + "DIMENS\n" + "1 5 2 /\n" + "RADIAL\n" + "GRID\n" + "INRAD\n" + "1 /\n" + "DRV\n" + "1 /\n" + "DTHETAV\n" + "3*90 60 30/\n" + "DZV\n" + "2*1 /\n" + "TOPS\n" + "5*1.0 /\n" + "\n"; + + Opm::Parser parser; + return parser.parseString( deckData, Opm::ParseContext()); +} + + + +BOOST_AUTO_TEST_CASE(RadialDetails) { + Opm::Deck deck = radial_details(); + Opm::EclipseGrid grid( deck ); + + BOOST_CHECK_CLOSE( grid.getCellVolume( 0 , 0 , 0 ) , 0.5*(2*2 - 1)*1, 0.0001); + BOOST_CHECK_CLOSE( grid.getCellVolume( 0 , 3 , 0 ) , sqrt(3.0)*0.25*( 4 - 1 ) , 0.0001); + 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); + + BOOST_CHECK_CLOSE( std::get<0>(pos2) , -0.75 , 0.0001); + BOOST_CHECK_CLOSE( std::get<1>(pos2) , -0.75 , 0.0001); + BOOST_CHECK_CLOSE( std::get<2>(pos2) , 1.50 , 0.0001); +} + + +BOOST_AUTO_TEST_CASE(CoordMapper) { + size_t nx = 10; + size_t ny = 7; + Opm::CoordMapper cmp = Opm::CoordMapper( nx , ny ); + BOOST_CHECK_THROW( cmp.index(12,6,0,0), std::invalid_argument ); + BOOST_CHECK_THROW( cmp.index(10,8,0,0), std::invalid_argument ); + BOOST_CHECK_THROW( cmp.index(10,7,5,0), std::invalid_argument ); + BOOST_CHECK_THROW( cmp.index(10,5,1,2), std::invalid_argument ); + + BOOST_CHECK_EQUAL( cmp.index(10,7,2,1) + 1 , cmp.size( )); +}