OilfieldToolsNet/opm-common
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Support cylindrical grids as a spiderweb with corrected volumes

Browse Source
This commit is contained in:
Peter Verveer
2020-12-09 15:51:01 +01:00
parent 7541b5e503
commit 81c7b53449
6 changed files with 179 additions and 13 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
opm/common/utility/numeric/calculateCellVol.hpp
View File

@@ -21,5 +21,6 @@
double calculateCellVol(const std::array<double,8>& X, const std::array<double,8>& Y, const std::array<double,8>& Z);
double calculateCylindricalCellVol(const double R1, const double R2, const double dTheta, const double dZ);

6
opm/parser/eclipse/EclipseState/Grid/EclipseGrid.hpp
View File

@@ -81,7 +81,6 @@ namespace Opm {
static bool hasGDFILE(const Deck& deck);
static bool hasCylindricalKeywords(const Deck& deck);
static bool hasSpiderwebKeywords(const Deck& deck);
static bool hasCornerPointKeywords(const Deck&);
static bool hasCartesianKeywords(const Deck&);
size_t getNumActive( ) const;
@@ -220,6 +219,10 @@ namespace Opm {
// Numerical aquifer cells, needs to be active
std::unordered_set<size_t> m_aquifer_cells;
// Radial grids need this for volume calculations.
std::optional<std::vector<double>> m_thetav;
std::optional<std::vector<double>> m_rv;
void updateNumericalAquiferCells(const Deck&);
void initGridFromEGridFile(Opm::EclIO::EclFile& egridfile, std::string fileName);
@@ -236,6 +239,7 @@ namespace Opm {
void initCylindricalGrid(const Deck&);
void initSpiderwebGrid(const Deck&);
void initSpiderwebOrCylindricalGrid(const Deck&, const bool);
void initCartesianGrid(const Deck&);
void initDTOPSGrid(const Deck&);
void initDVDEPTHZGrid(const Deck&);

9
src/opm/common/utility/numeric/calculateCellVol.cpp
View File

@@ -131,4 +131,11 @@ double calculateCellVol(const std::array<double,8>& X, const std::array<double,8
}
/*
Cell volume calculation for a cell from a cylindrical grid, given by the
inner and outer radius of the cell, and its spans in the angle and Z.
*/
double calculateCylindricalCellVol(const double r_inner, const double r_outer, const double delta_theta, const double delta_z)
{
return M_PI * std::abs((std::pow(r_outer,2) - std::pow(r_inner,2)) * delta_theta * delta_z) / 360.0;
}

50
src/opm/parser/eclipse/EclipseState/Grid/EclipseGrid.cpp
View File

@@ -923,24 +923,31 @@ EclipseGrid::EclipseGrid(const Deck& deck, const int * actnum)
return zcorn;
}
void EclipseGrid::initCylindricalGrid([[maybe_unused]] const Deck& deck)
void EclipseGrid::initCylindricalGrid(const Deck& deck)
{
throw std::invalid_argument("Cylindrical grid not implemented yet, use SPIDER web grid keyword instead for radial flow modeling");
initSpiderwebOrCylindricalGrid(deck, true);
}
/*
Limited implementaton - requires keywords: DRV, DTHETAV, DZV and TOPS.
*/
void EclipseGrid::initSpiderwebGrid(const Deck& deck)
{
// The hasSpiderKeywords( ) checks according to the
initSpiderwebOrCylindricalGrid(deck, false);
}
/*
Limited implementaton - requires keywords: DRV, DTHETAV, DZV or DZ, and TOPS.
*/
void EclipseGrid::initSpiderwebOrCylindricalGrid(const Deck& deck, const bool is_cylindrical)
{
const std::string kind = is_cylindrical ? "cylindrical" : "spiderweb";
// The hasCylindricalKeywords( ) checks according to the
// eclipse specification for RADIAL grid. 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 spiderweb grids present");
throw std::invalid_argument("Not all keywords required for " + kind + " grids present");
if (!deck.hasKeyword<ParserKeywords::DTHETAV>())
throw std::logic_error("The current implementation *must* have theta values specified using the DTHETAV keyword");
@@ -954,7 +961,7 @@ EclipseGrid::EclipseGrid(const Deck& deck, const int * actnum)
const std::vector<double>& drv = deck.getKeyword<ParserKeywords::DRV>().getSIDoubleData();
const std::vector<double>& dthetav = deck.getKeyword<ParserKeywords::DTHETAV>().getSIDoubleData();
const std::vector<double>& tops = deck.getKeyword<ParserKeywords::TOPS>().getSIDoubleData();
OpmLog::info(fmt::format("\nCreating spiderweb grid from keywords DRV, DTHETAV, DZV and TOPS"));
OpmLog::info(fmt::format("\nCreating {} grid from keywords DRV, DTHETAV, DZV and TOPS", kind));
if (drv.size() != this->getNX())
throw std::invalid_argument("DRV keyword should have exactly " + std::to_string( this->getNX() ) + " elements");
@@ -1055,6 +1062,13 @@ EclipseGrid::EclipseGrid(const Deck& deck, const int * actnum)
coord[ cm.index(i,j,2,1) ] = z2;
}
}
// Save angles, used by the cylindrical grid to calculate volumes.
if (is_cylindrical) {
m_rv = ri;
m_thetav = dthetav;
}
}
initCornerPointGrid( coord, zcorn, nullptr, nullptr);
}
@@ -1381,7 +1395,14 @@ std::vector<double> EclipseGrid::createDVector(const std::array<int,3>& dims, st
std::array<double,8> Z;
auto global_index = this->m_active_to_global[active_index];
this->getCellCorners(global_index, X, Y, Z );
active_volume[active_index] = calculateCellVol(X, Y, Z);
if (m_rv && m_thetav) {
const auto[i,j,k] = this->getIJK(global_index);
auto& r = *m_rv;
auto& t = *m_thetav;
active_volume[active_index] = calculateCylindricalCellVol(r[i], r[i+1], t[j], Z[4] - Z[4]);
} else {
active_volume[active_index] = calculateCellVol(X, Y, Z);
}
}
return active_volume;
@@ -1394,7 +1415,14 @@ std::vector<double> EclipseGrid::createDVector(const std::array<int,3>& dims, st
std::array<double,8> Y;
std::array<double,8> Z;
this->getCellCorners(globalIndex, X, Y, Z );
return calculateCellVol(X, Y, Z);
if (m_rv && m_thetav) {
const auto[i,j,k] = this->getIJK(globalIndex);
auto& r = *m_rv;
auto& t = *m_thetav;
return calculateCylindricalCellVol(r[i], r[i+1], t[j], Z[4] - Z[0]);
} else {
return calculateCellVol(X, Y, Z);
}
}
double EclipseGrid::getCellVolume(size_t i , size_t j , size_t k) const {

106
tests/parser/EclipseGridTests.cpp
View File

@@ -1465,6 +1465,112 @@ BOOST_AUTO_TEST_CASE(SpiderDetailsDZ) {
}
}
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"
"PORO \n"
" 10*0.15 /"
"\n";
Opm::Parser parser;
return parser.parseString( deckData);
}
BOOST_AUTO_TEST_CASE(RadialDetails) {
Opm::Deck deck = radial_details();
Opm::EclipseGrid grid( deck );
BOOST_CHECK_CLOSE( grid.getCellVolume( 0 , 0 , 0 ) , 0.75 * M_PI, 0.0001);
BOOST_CHECK_CLOSE( grid.getCellVolume( 0 , 3 , 0 ) , 0.5 * M_PI , 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);
{
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);
BOOST_CHECK_THROW( grid.getCornerPos( 0,0,0 , 8 ) , std::invalid_argument);
}
}
static Opm::Deck radial_details_dz() {
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"
"DZ\n"
"10*1 /\n"
"TOPS\n"
"5*1.0 /\n"
"PORO \n"
" 10*0.15 /"
"\n";
Opm::Parser parser;
return parser.parseString( deckData);
}
BOOST_AUTO_TEST_CASE(RadialDetailsDZ) {
Opm::Deck deck = radial_details_dz();
Opm::EclipseGrid grid( deck );
BOOST_CHECK_CLOSE( grid.getCellVolume( 0 , 0 , 0 ) , 0.75 * M_PI, 0.0001);
BOOST_CHECK_CLOSE( grid.getCellVolume( 0 , 3 , 0 ) , 0.5 * M_PI , 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);
{
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);
BOOST_CHECK_THROW( grid.getCornerPos( 0,0,0 , 8 ) , std::invalid_argument);
}
}
BOOST_AUTO_TEST_CASE(CoordMapper) {
size_t nx = 10;
size_t ny = 7;

20
tests/test_calculateCellVol.cpp
View File

@@ -59,4 +59,24 @@ BOOST_AUTO_TEST_CASE (calc_cellvol)
BOOST_REQUIRE_CLOSE (calculateCellVol(x4,y4,z4), 23391.4917234564, 1e-9);
}
BOOST_AUTO_TEST_CASE (calc_cellvol_cylindric)
{
{
double r_inner = 0.0, r_outer = 1.0, delta_theta = 360.0, delta_z = 1.0;
BOOST_REQUIRE_CLOSE(calculateCylindricalCellVol(r_inner, r_outer, delta_theta, delta_z), M_PI, 1e-9);
}
{
double r_inner = 1.0, r_outer = 2.0, delta_theta = 360.0, delta_z = 1.0;
BOOST_REQUIRE_CLOSE(calculateCylindricalCellVol(r_inner, r_outer, delta_theta, delta_z), 3.0 * M_PI, 1e-9);
}
{
double r_inner = 1.0, r_outer = 2.0, delta_theta = 120.0, delta_z = 1.0;
BOOST_REQUIRE_CLOSE(calculateCylindricalCellVol(r_inner, r_outer, delta_theta, delta_z), M_PI, 1e-9);
}
{
double r_inner = 1.0, r_outer = 2.0, delta_theta = 120.0, delta_z = 1.5;
BOOST_REQUIRE_CLOSE(calculateCylindricalCellVol(r_inner, r_outer, delta_theta, delta_z), 1.5 * M_PI, 1e-9);
}
}
BOOST_AUTO_TEST_SUITE_END()