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

Limited support for cylindrical grids.

Browse Source 
For cylindrical grids based on the keyword set DRV, DTHETAV, DZV and
TOPS the EclipseGrid implementation will create a cornerpoint
representation of the grid.
This commit is contained in:
Joakim Hove
2017-05-19 11:09:49 +02:00
parent fbb2ce73cd
commit 494cb74124
3 changed files with 486 additions and 9 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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

@@ -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<int, 3>&, const Deck&);
void initCartesianGrid( const std::array<int, 3>&, const Deck&);
void initCornerPointGrid( const std::array<int, 3>&, const Deck&);
void initDTOPSGrid( const std::array<int, 3>&, const Deck&);
@@ -227,11 +229,31 @@ namespace Opm {
static void scatterDim(const std::array<int, 3>& dims , size_t dim , const std::vector<double>& DV , std::vector<double>& 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

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

@@ -17,12 +17,12 @@
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#define _USE_MATH_DEFINES
#include <cmath>
#include <iostream>
#include <tuple>
#include <cmath>
#include <boost/lexical_cast.hpp>
#include <functional>
#include <opm/parser/eclipse/Deck/Section.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
@@ -35,8 +35,10 @@
#include <opm/parser/eclipse/Parser/ParserKeywords/A.hpp>
#include <opm/parser/eclipse/Parser/ParserKeywords/C.hpp>
#include <opm/parser/eclipse/Parser/ParserKeywords/D.hpp>
#include <opm/parser/eclipse/Parser/ParserKeywords/I.hpp>
#include <opm/parser/eclipse/Parser/ParserKeywords/M.hpp>
#include <opm/parser/eclipse/Parser/ParserKeywords/P.hpp>
#include <opm/parser/eclipse/Parser/ParserKeywords/R.hpp>
#include <opm/parser/eclipse/Parser/ParserKeywords/S.hpp>
#include <opm/parser/eclipse/Parser/ParserKeywords/T.hpp>
#include <opm/parser/eclipse/Parser/ParserKeywords/Z.hpp>
@@ -186,12 +188,16 @@ namespace Opm {
}
void EclipseGrid::initGrid( const std::array<int, 3>& dims, const Deck& deck) {
if (hasCornerPointKeywords(deck)) {
initCornerPointGrid(dims , deck);
} else if (hasCartesianKeywords(deck)) {
initCartesianGrid(dims , deck);
if (deck.hasKeyword<ParserKeywords::RADIAL>()) {
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<ParserKeywords::PINCH>()) {
@@ -311,6 +317,126 @@ namespace Opm {
}
/*
Limited implementaton - requires keywords: DRV, DTHETAV, DZV and TOPS.
*/
void EclipseGrid::initCylindricalGrid(const std::array<int, 3>& 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<ParserKeywords::DTHETAV>())
throw std::logic_error("The current implementation *must* have theta values specified using the DTHETAV keyword");
if (!deck.hasKeyword<ParserKeywords::DRV>())
throw std::logic_error("The current implementation *must* have radial values specified using the DRV keyword");
if (!deck.hasKeyword<ParserKeywords::DZV>() || !deck.hasKeyword<ParserKeywords::TOPS>())
throw std::logic_error("The current implementation *must* have vertical cell size specified using the DZV and TOPS keywords");
const std::vector<double>& drv = deck.getKeyword<ParserKeywords::DRV>().getSIDoubleData();
const std::vector<double>& dthetav = deck.getKeyword<ParserKeywords::DTHETAV>().getSIDoubleData();
const std::vector<double>& dzv = deck.getKeyword<ParserKeywords::DZV>().getSIDoubleData();
const std::vector<double>& tops = deck.getKeyword<ParserKeywords::TOPS>().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<ParserKeywords::CIRCLE>();
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<double> zcorn( zm.size() );
std::vector<double> coord( cm.size() );
{
std::vector<double> 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<double> ri(dims[0] + 1);
std::vector<double> 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<ParserKeywords::INRAD>().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<int,3>& dims ,
const std::vector<double>& coord ,
@@ -413,6 +539,17 @@ namespace Opm {
return hasDTOPSKeywords(deck);
}
bool EclipseGrid::hasCylindricalKeywords(const Deck& deck) {
if (deck.hasKeyword<ParserKeywords::INRAD>() &&
deck.hasKeyword<ParserKeywords::TOPS>() &&
(deck.hasKeyword<ParserKeywords::DZ>() || deck.hasKeyword<ParserKeywords::DZV>()) &&
(deck.hasKeyword<ParserKeywords::DRV>() || deck.hasKeyword<ParserKeywords::DR>()) &&
(deck.hasKeyword<ParserKeywords::DTHETA>() || deck.hasKeyword<ParserKeywords::DTHETAV>()))
return true;
else
return false;
}
bool EclipseGrid::hasDVDEPTHZKeywords(const Deck& deck) {
if (deck.hasKeyword<ParserKeywords::DXV>() &&
@@ -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;
}
}

287
tests/parser/EclipseGridTests.cpp
View File

@@ -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( ));
}