diff --git a/opm/parser/eclipse/EclipseState/Grid/EclipseGrid.hpp b/opm/parser/eclipse/EclipseState/Grid/EclipseGrid.hpp index a7a0266ab..0f96cf50d 100644 --- a/opm/parser/eclipse/EclipseState/Grid/EclipseGrid.hpp +++ b/opm/parser/eclipse/EclipseState/Grid/EclipseGrid.hpp @@ -33,6 +33,7 @@ #include #include #include +#include namespace Opm { @@ -216,6 +217,10 @@ namespace Opm { int m_nactive; std::vector m_active_to_global; std::vector m_global_to_active; + // Numerical aquifer cells, needs to be active + std::unordered_set m_aquifer_cells; + + void updateNumericalAquiferCells(const Deck&); void initGridFromEGridFile(Opm::EclIO::EclFile& egridfile, std::string fileName); void resetACTNUM( const int* actnum); diff --git a/src/opm/parser/eclipse/EclipseState/Grid/EclipseGrid.cpp b/src/opm/parser/eclipse/EclipseState/Grid/EclipseGrid.cpp index ce1c8ae54..c06993ac5 100644 --- a/src/opm/parser/eclipse/EclipseState/Grid/EclipseGrid.cpp +++ b/src/opm/parser/eclipse/EclipseState/Grid/EclipseGrid.cpp @@ -223,6 +223,8 @@ EclipseGrid::EclipseGrid(const Deck& deck, const int * actnum) } + updateNumericalAquiferCells(deck); + initGrid(deck); if (deck.hasKeyword("MAPUNITS")){ @@ -1092,7 +1094,7 @@ EclipseGrid::EclipseGrid(const Deck& deck, const int * actnum) actnum = actnumVector.data(); OpmLog::info(fmt::format("\nCreating cornerpoint grid from keywords ZCORN, COORD and ACTNUM")); - } else + } else OpmLog::info(fmt::format("\nCreating cornerpoint grid from keywords ZCORN and COORD")); initCornerPointGrid( coord , zcorn, actnum, nullptr ); @@ -1700,7 +1702,11 @@ std::vector EclipseGrid::createDVector(const std::array& dims, st for (size_t n = 0; n < global_size; n++) { this->m_actnum[n] = actnum[n]; - if (actnum[n] > 0) { + // numerical aquifer cells need to be active + if (this->m_aquifer_cells.count(n) > 0) { + this->m_actnum[n] = 1; + } + if (this->m_actnum[n] > 0) { this->m_global_to_active.push_back(this->m_nactive); this->m_active_to_global.push_back(n); this->m_nactive++; @@ -1723,6 +1729,23 @@ std::vector EclipseGrid::createDVector(const std::array& dims, st return ZcornMapper( getNX() , getNY(), getNZ() ); } + void EclipseGrid::updateNumericalAquiferCells(const Deck& deck) { + using AQUNUM =ParserKeywords::AQUNUM; + if ( !deck.hasKeyword() ) { + return; + } + const auto &aqunum_keywords = deck.getKeywordList(); + for (const auto &keyword : aqunum_keywords) { + for (const auto &record : *keyword) { + const size_t i = record.getItem().get(0) - 1; + const size_t j = record.getItem().get(0) - 1; + const size_t k = record.getItem().get(0) - 1; + const size_t global_index = this->getGlobalIndex(i, j, k); + this->m_aquifer_cells.insert(global_index); + } + } + } + ZcornMapper::ZcornMapper(size_t nx , size_t ny, size_t nz) : dims( {{nx,ny,nz}} ), stride( {{2 , 4*nx, 8*nx*ny}} ), diff --git a/tests/parser/EclipseGridTests.cpp b/tests/parser/EclipseGridTests.cpp index 070b93702..1c4a8cb59 100644 --- a/tests/parser/EclipseGridTests.cpp +++ b/tests/parser/EclipseGridTests.cpp @@ -2148,6 +2148,69 @@ BOOST_AUTO_TEST_CASE(TESTCP_ACTNUM_UPDATE) { } +BOOST_AUTO_TEST_CASE(TESTCP_ACTNUM_AQUNUM) { + const std::string deckData = R"( + RUNSPEC + DIMENS + 3 2 1 / + GRID + COORD + 2000.0000 2000.0000 2000.0000 1999.9476 2000.0000 2002.9995 + 2049.9924 2000.0000 2000.8726 2049.9400 2000.0000 2003.8722 + 2099.9848 2000.0000 2001.7452 2099.9324 2000.0000 2004.7448 + 2149.9772 2000.0000 2002.6179 2149.9248 2000.0000 2005.6174 + 2000.0000 2050.0000 2000.0000 1999.9476 2050.0000 2002.9995 + 2049.9924 2050.0000 2000.8726 2049.9400 2050.0000 2003.8722 + 2099.9848 2050.0000 2001.7452 2099.9324 2050.0000 2004.7448 + 2149.9772 2050.0000 2002.6179 2149.9248 2050.0000 2005.6174 + 2000.0000 2100.0000 2000.0000 1999.9476 2100.0000 2002.9995 + 2049.9924 2100.0000 2000.8726 2049.9400 2100.0000 2003.8722 + 2099.9848 2100.0000 2001.7452 2099.9324 2100.0000 2004.7448 + 2149.9772 2100.0000 2002.6179 2149.9248 2100.0000 2005.6174 / + ZCORN + 2000.0000 2000.8726 2000.8726 2001.7452 2001.7452 2002.6179 + 2000.0000 2000.8726 2000.8726 2001.7452 2001.7452 2002.6179 + 2000.0000 2000.8726 2000.8726 2001.7452 2001.7452 2002.6179 + 2000.0000 2000.8726 2000.8726 2001.7452 2001.7452 2002.6179 + 2002.9995 2003.8722 2003.8722 2004.7448 2004.7448 2005.6174 + 2002.9995 2003.8722 2003.8722 2004.7448 2004.7448 2005.6174 + 2002.9995 2003.8722 2003.8722 2004.7448 2004.7448 2005.6174 + 2002.9995 2003.8722 2003.8722 2004.7448 2004.7448 2005.6174 / + ACTNUM + 0 0 1 1 0 1 / + PORO + 6*0.15 / + AQUNUM + 1 1 1 1 1000000.0 10000 0.25 400 2585.00 285.00 1 1 / + 1 2 1 1 1000000.0 10000 0.25 400 2585.00 285.00 1 1 / + / + AQUCON + 1 3 3 2 2 1 1 'J+' 1.00 1 / + / + EDIT)"; + + Opm::Parser parser; + const auto deck = parser.parseString( deckData) ; + Opm::EclipseGrid grid( deck); + + const std::vector& grid_actnum = grid.getACTNUM(); + const std::vector desired_actnum = {1, 1, 1, 1, 0, 1}; + + BOOST_CHECK( grid_actnum.size() == desired_actnum.size() ); + + for (size_t n=0; n< grid_actnum.size(); n++) { + BOOST_CHECK_EQUAL( grid_actnum[n], desired_actnum[n] ); + } + + Opm::EclipseState es(deck); + const auto& grid_actnum2 = es.getInputGrid().getACTNUM(); + + BOOST_CHECK( desired_actnum.size() == grid_actnum2.size()); + for (size_t n=0; n< grid_actnum.size(); n++) { + BOOST_CHECK_EQUAL( grid_actnum2[n], desired_actnum[n] ); + } +} + BOOST_AUTO_TEST_CASE(TEST_altGridConstructors) { const char* deckData =