Updating class EGrid
- now supporting parsing/reading of lgrs - detect if radial grid is used - return xyz coordinates also for radial grids - support for getting nnc data
This commit is contained in:
Torbjørn Skille
parent
e1334adbf6
commit
a7314ab5a0
@ -479,6 +479,7 @@ if(ENABLE_ECL_INPUT)
|
||||
list (APPEND TEST_DATA_FILES
|
||||
tests/ECLFILE.INIT
|
||||
tests/ECLFILE.FINIT
|
||||
tests/LGR_TESTMOD.EGRID
|
||||
tests/LGR_TESTMOD.INIT
|
||||
tests/MODEL1_IX.INIT
|
||||
tests/MODEL1_IX.SMSPEC
|
||||
|
||||
@ -20,6 +20,7 @@
|
||||
#define OPM_IO_EGRID_HPP
|
||||
|
||||
#include <opm/io/eclipse/EclFile.hpp>
|
||||
#include <opm/common/utility/FileSystem.hpp>
|
||||
|
||||
#include <array>
|
||||
#include <iostream>
|
||||
@ -34,7 +35,7 @@ namespace Opm { namespace EclIO {
|
||||
class EGrid : public EclFile
|
||||
{
|
||||
public:
|
||||
explicit EGrid(const std::string& filename);
|
||||
explicit EGrid(const std::string& filename, std::string grid_name = "global");
|
||||
|
||||
int global_index(int i, int j, int k) const;
|
||||
int active_index(int i, int j, int k) const;
|
||||
@ -44,20 +45,54 @@ public:
|
||||
std::array<int, 3> ijk_from_active_index(int actInd) const;
|
||||
std::array<int, 3> ijk_from_global_index(int globInd) const;
|
||||
|
||||
void getCellCorners(int globindex, std::array<double,8>& X, std::array<double,8>& Y, std::array<double,8>& Z) const;
|
||||
void getCellCorners(const std::array<int, 3>& ijk, std::array<double,8>& X, std::array<double,8>& Y, std::array<double,8>& Z) const;
|
||||
void getCellCorners(int globindex, std::array<double,8>& X, std::array<double,8>& Y, std::array<double,8>& Z);
|
||||
void getCellCorners(const std::array<int, 3>& ijk, std::array<double,8>& X, std::array<double,8>& Y, std::array<double,8>& Z);
|
||||
|
||||
int activeCells() const { return nactive; }
|
||||
int totalNumberOfCells() const { return nijk[0] * nijk[1] * nijk[2]; }
|
||||
|
||||
void load_grid_data();
|
||||
void load_nnc_data();
|
||||
bool is_radial() const { return m_radial; }
|
||||
|
||||
const std::vector<int>& hostCellsGlobalIndex() const { return host_cells; }
|
||||
std::vector<std::array<int, 3>> hostCellsIJK();
|
||||
|
||||
// zero based: i1, j1,k1, i2,j2,k2, transmisibility
|
||||
using NNCentry = std::tuple<int, int, int, int, int, int, float>;
|
||||
std::vector<NNCentry> get_nnc_ijk();
|
||||
|
||||
const std::vector<std::string>& list_of_lgrs() const { return lgr_names; }
|
||||
|
||||
private:
|
||||
Opm::filesystem::path inputFileName, initFileName;
|
||||
std::string m_grid_name;
|
||||
bool m_radial;
|
||||
|
||||
std::array<int, 3> nijk;
|
||||
int nactive;
|
||||
std::array<int, 3> host_nijk;
|
||||
|
||||
int nactive, m_nncs;
|
||||
mutable bool m_nncs_loaded;
|
||||
|
||||
std::vector<int> act_index;
|
||||
std::vector<int> glob_index;
|
||||
|
||||
std::vector<float> coord_array;
|
||||
std::vector<float> zcorn_array;
|
||||
|
||||
std::vector<int> nnc1_array;
|
||||
std::vector<int> nnc2_array;
|
||||
std::vector<float> transnnc_array;
|
||||
std::vector<int> host_cells;
|
||||
|
||||
std::vector<std::string> lgr_names;
|
||||
|
||||
int zcorn_array_index;
|
||||
int coord_array_index;
|
||||
int actnum_array_index;
|
||||
int nnc1_array_index;
|
||||
int nnc2_array_index;
|
||||
};
|
||||
|
||||
}} // namespace Opm::EclIO
|
||||
|
||||
@ -17,6 +17,8 @@
|
||||
*/
|
||||
|
||||
#include <opm/io/eclipse/EGrid.hpp>
|
||||
#include <opm/io/eclipse/EInit.hpp>
|
||||
#include <opm/io/eclipse/EclUtil.hpp>
|
||||
|
||||
#include <opm/common/ErrorMacros.hpp>
|
||||
|
||||
@ -28,43 +30,202 @@
|
||||
#include <string>
|
||||
#include <sstream>
|
||||
|
||||
#include <math.h>
|
||||
|
||||
namespace Opm { namespace EclIO {
|
||||
|
||||
EGrid::EGrid(const std::string &filename) : EclFile(filename)
|
||||
using NNCentry = std::tuple<int, int, int, int, int,int, float>;
|
||||
|
||||
EGrid::EGrid(const std::string &filename, std::string grid_name) :
|
||||
EclFile(filename), inputFileName { filename }, m_grid_name {grid_name}
|
||||
{
|
||||
EclFile file(filename);
|
||||
initFileName = inputFileName.parent_path() / inputFileName.stem();
|
||||
|
||||
auto gridhead = get<int>("GRIDHEAD");
|
||||
if (this->formattedInput())
|
||||
initFileName += ".FINIT";
|
||||
else
|
||||
initFileName += ".INIT";
|
||||
|
||||
nijk[0] = gridhead[1];
|
||||
nijk[1] = gridhead[2];
|
||||
nijk[2] = gridhead[3];
|
||||
std::string lgrname = "global";
|
||||
|
||||
if (file.hasKey("ACTNUM")) {
|
||||
auto actnum = get<int>("ACTNUM");
|
||||
m_nncs_loaded = false;
|
||||
actnum_array_index = -1;
|
||||
nnc1_array_index = -1;
|
||||
nnc2_array_index = -1;
|
||||
|
||||
nactive = 0;
|
||||
for (unsigned int i = 0; i < actnum.size(); i++) {
|
||||
if (actnum[i] > 0) {
|
||||
act_index.push_back(nactive);
|
||||
glob_index.push_back(i);
|
||||
nactive++;
|
||||
} else {
|
||||
act_index.push_back(-1);
|
||||
}
|
||||
int hostnum_index = -1;
|
||||
|
||||
for (size_t n = 0; n < array_name.size(); n++) {
|
||||
|
||||
if (array_name[n] == "ENDLGR")
|
||||
lgrname = "global";
|
||||
|
||||
if (array_name[n] == "LGR") {
|
||||
auto lgr = this->get<std::string>(n);
|
||||
lgrname = lgr[0];
|
||||
lgr_names.push_back(lgr[0]);
|
||||
}
|
||||
} else {
|
||||
int nCells = nijk[0] * nijk[1] * nijk[2];
|
||||
act_index.resize(nCells);
|
||||
glob_index.resize(nCells);
|
||||
std::iota(act_index.begin(), act_index.end(), 0);
|
||||
std::iota(glob_index.begin(), glob_index.end(), 0);
|
||||
}
|
||||
|
||||
coord_array = get<float>("COORD");
|
||||
zcorn_array = get<float>("ZCORN");
|
||||
if (array_name[n] == "NNCHEAD"){
|
||||
auto nnchead = this->get<int>(n);
|
||||
|
||||
if (nnchead[1] == 0)
|
||||
lgrname = "global";
|
||||
else
|
||||
lgrname = lgr_names[nnchead[1] - 1];
|
||||
}
|
||||
|
||||
if (lgrname == grid_name) {
|
||||
if (array_name[n] == "GRIDHEAD") {
|
||||
auto gridhead = get<int>(n);
|
||||
nijk[0] = gridhead[1];
|
||||
nijk[1] = gridhead[2];
|
||||
nijk[2] = gridhead[3];
|
||||
|
||||
m_radial = gridhead[26] > 0 ? true: false;
|
||||
}
|
||||
|
||||
if (array_name[n] == "COORD")
|
||||
coord_array_index= n;
|
||||
else if (array_name[n] == "ZCORN")
|
||||
zcorn_array_index= n;
|
||||
else if (array_name[n] == "ACTNUM")
|
||||
actnum_array_index= n;
|
||||
else if (array_name[n] == "NNC1")
|
||||
nnc1_array_index= n;
|
||||
else if (array_name[n] == "NNC2")
|
||||
nnc2_array_index= n;
|
||||
else if (array_name[n] == "HOSTNUM")
|
||||
hostnum_index= n;
|
||||
}
|
||||
|
||||
if ((lgrname == "global") && (array_name[n] == "GRIDHEAD")) {
|
||||
auto gridhead = get<int>(n);
|
||||
host_nijk[0] = gridhead[1];
|
||||
host_nijk[1] = gridhead[2];
|
||||
host_nijk[2] = gridhead[3];
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
if (actnum_array_index != -1) {
|
||||
auto actnum = this->get<int>(actnum_array_index);
|
||||
nactive = 0;
|
||||
for (size_t i = 0; i < actnum.size(); i++) {
|
||||
if (actnum[i] > 0) {
|
||||
act_index.push_back(nactive);
|
||||
glob_index.push_back(i);
|
||||
nactive++;
|
||||
} else {
|
||||
act_index.push_back(-1);
|
||||
}
|
||||
}
|
||||
} else {
|
||||
int nCells = nijk[0] * nijk[1] * nijk[2];
|
||||
act_index.resize(nCells);
|
||||
glob_index.resize(nCells);
|
||||
std::iota(act_index.begin(), act_index.end(), 0);
|
||||
std::iota(glob_index.begin(), glob_index.end(), 0);
|
||||
}
|
||||
|
||||
if (hostnum_index > -1){
|
||||
auto hostnum = getImpl(hostnum_index, INTE, inte_array, "integer");
|
||||
host_cells.reserve(hostnum.size());
|
||||
|
||||
for (auto val : hostnum)
|
||||
host_cells.push_back(val -1);
|
||||
}
|
||||
}
|
||||
|
||||
std::vector<std::array<int, 3>> EGrid::hostCellsIJK()
|
||||
{
|
||||
std::vector<std::array<int, 3>> res_vect;
|
||||
res_vect.reserve(host_cells.size());
|
||||
|
||||
for (auto val : host_cells){
|
||||
std::array<int, 3> tmp;
|
||||
tmp[2] = val / (host_nijk[0] * host_nijk[1]);
|
||||
int rest = val % (host_nijk[0] * host_nijk[1]);
|
||||
|
||||
tmp[1] = rest / host_nijk[0];
|
||||
tmp[0] = rest % host_nijk[0];
|
||||
|
||||
res_vect.push_back(tmp);
|
||||
}
|
||||
|
||||
return res_vect;
|
||||
}
|
||||
|
||||
std::vector<NNCentry> EGrid::get_nnc_ijk()
|
||||
{
|
||||
if (!m_nncs_loaded)
|
||||
load_nnc_data();
|
||||
|
||||
std::vector<NNCentry> res_vect;
|
||||
res_vect.reserve(nnc1_array.size());
|
||||
|
||||
for (size_t n=0; n< nnc1_array.size(); n++){
|
||||
auto ijk1 = ijk_from_global_index(nnc1_array[n] - 1);
|
||||
auto ijk2 = ijk_from_global_index(nnc2_array[n] - 1);
|
||||
|
||||
if (transnnc_array.size() > 0)
|
||||
res_vect.push_back({ijk1[0], ijk1[1], ijk1[2], ijk2[0], ijk2[1], ijk2[2], transnnc_array[n]});
|
||||
else
|
||||
res_vect.push_back({ijk1[0], ijk1[1], ijk1[2], ijk2[0], ijk2[1], ijk2[2], -1.0 });
|
||||
}
|
||||
|
||||
return res_vect;
|
||||
}
|
||||
|
||||
|
||||
void EGrid::load_grid_data()
|
||||
{
|
||||
coord_array = getImpl(coord_array_index, REAL, real_array, "float");
|
||||
zcorn_array = getImpl(zcorn_array_index, REAL, real_array, "float");
|
||||
}
|
||||
|
||||
void EGrid::load_nnc_data()
|
||||
{
|
||||
if ((nnc1_array_index > -1) && (nnc2_array_index > -1)) {
|
||||
|
||||
nnc1_array = getImpl(nnc1_array_index, Opm::EclIO::INTE, inte_array, "inte");
|
||||
nnc2_array = getImpl(nnc2_array_index, Opm::EclIO::INTE, inte_array, "inte");
|
||||
|
||||
if ((Opm::filesystem::exists(initFileName)) && (nnc1_array.size() > 0)){
|
||||
Opm::EclIO::EInit init(initFileName);
|
||||
|
||||
auto init_dims = init.grid_dimension(m_grid_name);
|
||||
int init_nactive = init.activeCells(m_grid_name);
|
||||
|
||||
if (init_dims != nijk){
|
||||
std::string message = "Dimensions of Egrid differ from dimensions found in init file. ";
|
||||
std::string grid_str = std::to_string(nijk[0]) + "x" + std::to_string(nijk[1]) + "x" + std::to_string(nijk[2]);
|
||||
std::string init_str = std::to_string(init_dims[0]) + "x" + std::to_string(init_dims[1]) + "x" + std::to_string(init_dims[2]);
|
||||
message = message + "Egrid: " + grid_str + ". INIT file: " + init_str;
|
||||
OPM_THROW(std::invalid_argument, message);
|
||||
}
|
||||
|
||||
if (init_nactive != nactive){
|
||||
std::string message = "Number of active cells are different in Egrid and Init file.";
|
||||
message = message + " Egrid: " + std::to_string(nactive) + ". INIT file: " + std::to_string(init_nactive);
|
||||
OPM_THROW(std::invalid_argument, message);
|
||||
}
|
||||
|
||||
auto trans_data = init.getInitData<float>("TRANNNC", m_grid_name);
|
||||
|
||||
if (trans_data.size() != nnc1_array.size()){
|
||||
std::string message = "inconsistent size of array TRANNNC in init file. ";
|
||||
message = message + " Size of NNC1 and NNC2: " + std::to_string(nnc1_array.size());
|
||||
message = message + " Size of TRANNNC: " + std::to_string(trans_data.size());
|
||||
OPM_THROW(std::invalid_argument, message);
|
||||
}
|
||||
|
||||
transnnc_array = trans_data;
|
||||
}
|
||||
|
||||
m_nncs_loaded = true;
|
||||
}
|
||||
}
|
||||
|
||||
int EGrid::global_index(int i, int j, int k) const
|
||||
{
|
||||
@ -129,52 +290,65 @@ std::array<int, 3> EGrid::ijk_from_global_index(int globInd) const
|
||||
void EGrid::getCellCorners(const std::array<int, 3>& ijk,
|
||||
std::array<double,8>& X,
|
||||
std::array<double,8>& Y,
|
||||
std::array<double,8>& Z) const
|
||||
std::array<double,8>& Z)
|
||||
{
|
||||
if (coord_array.empty())
|
||||
load_grid_data();
|
||||
|
||||
std::vector<int> zind;
|
||||
std::vector<int> pind;
|
||||
|
||||
|
||||
// calculate indices for grid pillars in COORD arrray
|
||||
pind.push_back(ijk[1]*(nijk[0]+1)*6 + ijk[0]*6);
|
||||
pind.push_back(pind[0] + 6);
|
||||
pind.push_back(pind[0] + (nijk[0]+1)*6);
|
||||
pind.push_back(pind[2] + 6);
|
||||
pind.push_back(ijk[1]*(nijk[0]+1)*6 + ijk[0]*6);
|
||||
pind.push_back(pind[0] + 6);
|
||||
pind.push_back(pind[0] + (nijk[0]+1)*6);
|
||||
pind.push_back(pind[2] + 6);
|
||||
|
||||
// get depths from zcorn array in ZCORN array
|
||||
zind.push_back(ijk[2]*nijk[0]*nijk[1]*8 + ijk[1]*nijk[0]*4 + ijk[0]*2);
|
||||
zind.push_back(zind[0] + 1);
|
||||
zind.push_back(zind[0] + nijk[0]*2);
|
||||
zind.push_back(zind[2] + 1);
|
||||
// get depths from zcorn array in ZCORN array
|
||||
zind.push_back(ijk[2]*nijk[0]*nijk[1]*8 + ijk[1]*nijk[0]*4 + ijk[0]*2);
|
||||
zind.push_back(zind[0] + 1);
|
||||
zind.push_back(zind[0] + nijk[0]*2);
|
||||
zind.push_back(zind[2] + 1);
|
||||
|
||||
for (int n = 0; n < 4; n++) {
|
||||
zind.push_back(zind[n] + nijk[0]*nijk[1]*4);
|
||||
}
|
||||
for (int n = 0; n < 4; n++)
|
||||
zind.push_back(zind[n] + nijk[0]*nijk[1]*4);
|
||||
|
||||
for (int n = 0; n< 8; n++){
|
||||
Z[n] = zcorn_array[zind[n]];
|
||||
}
|
||||
for (int n = 0; n< 8; n++)
|
||||
Z[n] = zcorn_array[zind[n]];
|
||||
|
||||
for (int n = 0; n < 4; n++) {
|
||||
double xt = coord_array[pind[n]];
|
||||
double yt = coord_array[pind[n] + 1];
|
||||
double zt = coord_array[pind[n] + 2];
|
||||
for (int n = 0; n < 4; n++) {
|
||||
double xt;
|
||||
double yt;
|
||||
double xb;
|
||||
double yb;
|
||||
|
||||
double xb = coord_array[pind[n] + 3];
|
||||
double yb = coord_array[pind[n] + 4];
|
||||
double zb = coord_array[pind[n]+5];
|
||||
double zt = coord_array[pind[n] + 2];
|
||||
double zb = coord_array[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]);
|
||||
if (m_radial) {
|
||||
xt = coord_array[pind[n]] * cos(coord_array[pind[n] + 1] / 180.0 * M_PI);
|
||||
yt = coord_array[pind[n]] * sin(coord_array[pind[n] + 1] / 180.0 * M_PI);
|
||||
xb = coord_array[pind[n]+3] * cos(coord_array[pind[n] + 4] / 180.0 * M_PI);
|
||||
yb = coord_array[pind[n]+3] * sin(coord_array[pind[n] + 4] / 180.0 * M_PI);
|
||||
} else {
|
||||
xt = coord_array[pind[n]];
|
||||
yt = coord_array[pind[n] + 1];
|
||||
xb = coord_array[pind[n] + 3];
|
||||
yb = coord_array[pind[n] + 4];
|
||||
}
|
||||
|
||||
Y[n] = yt+(yb-yt)/(zt-zb)*(zt-Z[n]);
|
||||
Y[n+4] = yt+(yb-yt)/(zt-zb)*(zt-Z[n+4]);
|
||||
}
|
||||
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]);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
void EGrid::getCellCorners(int globindex, std::array<double,8>& X,
|
||||
std::array<double,8>& Y, std::array<double,8>& Z) const
|
||||
std::array<double,8>& Y, std::array<double,8>& Z)
|
||||
{
|
||||
return getCellCorners(ijk_from_global_index(globindex),X,Y,Z);
|
||||
}
|
||||
|
||||
BIN
tests/LGR_TESTMOD.EGRID
Normal file
BIN
tests/LGR_TESTMOD.EGRID
Normal file
Binary file not shown.
@ -20,6 +20,9 @@
|
||||
#include "config.h"
|
||||
|
||||
#include <opm/io/eclipse/EGrid.hpp>
|
||||
#include <opm/io/eclipse/EInit.hpp>
|
||||
#include <opm/common/utility/numeric/calculateCellVol.hpp>
|
||||
|
||||
|
||||
#define BOOST_TEST_MODULE Test EGrid
|
||||
#include <boost/test/unit_test.hpp>
|
||||
@ -71,22 +74,22 @@ bool operator==(const std::vector<T> & t1, const std::vector<T> & t2)
|
||||
}
|
||||
|
||||
// test is using SPE1CASE1, with minor modifications in order to test API for EGrid class
|
||||
// -> 6 cells made inactive, box: 5 7 5 6 1 1
|
||||
|
||||
// -> 6 cells made inactive, box: 5 7 5 6 1 1
|
||||
|
||||
// first pilar at x=0.0, y=0.0 and z=0.0
|
||||
// dx = 1000 ft, dy = 1000 ft. dz = 20 ft for layer 1, 30 ft for layer 2 and 50 ft for layer 3.
|
||||
// size of grid is 10x10x3
|
||||
// size of grid is 10x10x3
|
||||
|
||||
|
||||
|
||||
BOOST_AUTO_TEST_CASE(DimensAndIndices) {
|
||||
|
||||
|
||||
std::string testFile="SPE1CASE1.EGRID";
|
||||
|
||||
EGrid grid1(testFile);
|
||||
|
||||
|
||||
int nAct=grid1.activeCells();
|
||||
int nTot=grid1.totalNumberOfCells();
|
||||
|
||||
|
||||
BOOST_CHECK_EQUAL(nAct,294);
|
||||
BOOST_CHECK_EQUAL(nTot,300);
|
||||
|
||||
@ -96,7 +99,7 @@ BOOST_AUTO_TEST_CASE(DimensAndIndices) {
|
||||
BOOST_CHECK_EQUAL(dim[1],10);
|
||||
BOOST_CHECK_EQUAL(dim[2],3);
|
||||
|
||||
int globInd = grid1.global_index(3, 2, 1);
|
||||
int globInd = grid1.global_index(3, 2, 1);
|
||||
|
||||
BOOST_CHECK_EQUAL(globInd, 123); // 10*10*1 + 10*2 + 3 = 100+20+3 = 123
|
||||
|
||||
@ -104,7 +107,7 @@ BOOST_AUTO_TEST_CASE(DimensAndIndices) {
|
||||
BOOST_CHECK_EQUAL(grid1.global_index(dim[0] - 1, dim[1] - 1, dim[2] - 1), nTot - 1);
|
||||
|
||||
// check global_index valid range, should throw exception if outside
|
||||
|
||||
|
||||
BOOST_CHECK_THROW(grid1.global_index(-1, -1, -1), std::invalid_argument);
|
||||
BOOST_CHECK_THROW(grid1.global_index(0, -1, 0), std::invalid_argument);
|
||||
BOOST_CHECK_THROW(grid1.global_index(0, 0, -1), std::invalid_argument);
|
||||
@ -112,11 +115,11 @@ BOOST_AUTO_TEST_CASE(DimensAndIndices) {
|
||||
BOOST_CHECK_THROW(grid1.global_index(dim[0], 0, 0), std::invalid_argument);
|
||||
BOOST_CHECK_THROW(grid1.global_index(dim[0], dim[1], 0), std::invalid_argument);
|
||||
BOOST_CHECK_THROW(grid1.global_index(dim[0], dim[1], dim[2]), std::invalid_argument);
|
||||
|
||||
// 6 inactive cells in first layer, actInd not same as glogInd
|
||||
// actInd and globInd are zero based indices
|
||||
|
||||
int actInd = grid1.active_index(3,2,1);
|
||||
// 6 inactive cells in first layer, actInd not same as glogInd
|
||||
// actInd and globInd are zero based indices
|
||||
|
||||
int actInd = grid1.active_index(3,2,1);
|
||||
|
||||
BOOST_CHECK_EQUAL(actInd, 117); // global index 123, - 6 inactive
|
||||
BOOST_CHECK_EQUAL(grid1.active_index(0, 0, 0), 0);
|
||||
@ -131,18 +134,18 @@ BOOST_AUTO_TEST_CASE(DimensAndIndices) {
|
||||
BOOST_CHECK_THROW(grid1.active_index(dim[0], 0, 0), std::invalid_argument);
|
||||
BOOST_CHECK_THROW(grid1.active_index(dim[0], dim[1], 0), std::invalid_argument);
|
||||
BOOST_CHECK_THROW(grid1.active_index(dim[0], dim[1], dim[2]), std::invalid_argument);
|
||||
|
||||
|
||||
auto res = grid1.ijk_from_active_index(actInd);
|
||||
|
||||
|
||||
BOOST_CHECK_EQUAL(res[0], 3);
|
||||
BOOST_CHECK_EQUAL(res[1], 2);
|
||||
BOOST_CHECK_EQUAL(res[2], 1);
|
||||
|
||||
BOOST_CHECK_THROW(grid1.ijk_from_active_index(-1), std::invalid_argument);
|
||||
BOOST_CHECK_THROW(grid1.ijk_from_active_index(nAct), std::invalid_argument);
|
||||
|
||||
|
||||
res = grid1.ijk_from_global_index(globInd);
|
||||
|
||||
|
||||
BOOST_CHECK_EQUAL(res[0], 3);
|
||||
BOOST_CHECK_EQUAL(res[1], 2);
|
||||
BOOST_CHECK_EQUAL(res[2], 1);
|
||||
@ -153,7 +156,7 @@ BOOST_AUTO_TEST_CASE(DimensAndIndices) {
|
||||
|
||||
|
||||
BOOST_AUTO_TEST_CASE(getCellCorners) {
|
||||
|
||||
|
||||
std::string testFile="SPE1CASE1.EGRID";
|
||||
|
||||
std::array<double,8> ref_X = {3000,4000,3000,4000,3000,4000,3000,4000};
|
||||
@ -165,7 +168,7 @@ BOOST_AUTO_TEST_CASE(getCellCorners) {
|
||||
std::array<double,8> X = {0.0};
|
||||
std::array<double,8> Y = {0.0};
|
||||
std::array<double,8> Z = {0.0};
|
||||
|
||||
|
||||
// cell 4,3,2 => zero based 3,2,1
|
||||
grid1.getCellCorners({3, 2, 1}, X, Y, Z);
|
||||
|
||||
@ -177,7 +180,7 @@ BOOST_AUTO_TEST_CASE(getCellCorners) {
|
||||
Y = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
|
||||
Z = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
|
||||
|
||||
int globInd=grid1.global_index(3,2,1);
|
||||
int globInd=grid1.global_index(3,2,1);
|
||||
|
||||
grid1.getCellCorners(globInd,X,Y,Z);
|
||||
|
||||
@ -185,3 +188,159 @@ BOOST_AUTO_TEST_CASE(getCellCorners) {
|
||||
BOOST_CHECK_EQUAL(Y == ref_Y, true);
|
||||
BOOST_CHECK_EQUAL(Z == ref_Z, true);
|
||||
}
|
||||
|
||||
|
||||
BOOST_AUTO_TEST_CASE(lgr_1) {
|
||||
|
||||
std::string testEgridFile = "LGR_TESTMOD.EGRID";
|
||||
std::string testInitFile = "LGR_TESTMOD.INIT";
|
||||
|
||||
std::vector<EGrid::NNCentry> nnc_global_ref = {{0,0,1,0,1,0,8.52637}, {0,0,2,0,1,1,8.52637},
|
||||
{0,0,3,0,1,2,8.52637}, {0,0,4,0,1,3,8.52637}, {1,0,1,1,1,0,8.52637}, {1,0,2,1,1,1,8.52637},
|
||||
{1,0,3,1,1,2,8.52637}, {1,0,4,1,1,3,8.52637} };
|
||||
|
||||
std::vector<EGrid::NNCentry> nnc_lgr1_ref = {{0,3,1,0,4,0,7.67373}, {0,3,2,0,4,1,7.67373},
|
||||
{0,3,3,0,4,2,7.67373}, {1,3,1,1,4,0,7.67373}, {1,3,2,1,4,1,7.67373}, {1,3,3,1,4,2,7.67373},
|
||||
{2,3,1,2,4,0,7.67373}, {2,3,2,2,4,1,7.67373}, {2,3,3,2,4,2,7.67373}, {3,3,1,3,4,0,7.67373},
|
||||
{3,3,2,3,4,1,7.67373}, {3,3,3,3,4,2,7.67373}};
|
||||
|
||||
const std::vector<std::string> ref_lgr_list = {"LGR1", "LGR2"};
|
||||
|
||||
const std::array<int, 3> ref_dim_global = {2,3,5};
|
||||
const std::array<int, 3> ref_dim_lgr1 = {4,8,4};
|
||||
const std::array<int, 3> ref_dim_lgr2 = {6,8,4};
|
||||
|
||||
EGrid grid1(testEgridFile);
|
||||
|
||||
auto nijk_global = grid1.dimension();
|
||||
BOOST_CHECK_EQUAL(nijk_global == ref_dim_global, true);
|
||||
|
||||
auto lgr_list = grid1.list_of_lgrs();
|
||||
|
||||
BOOST_CHECK_EQUAL(lgr_list.size(), 2);
|
||||
BOOST_CHECK_EQUAL(lgr_list[0], "LGR1");
|
||||
BOOST_CHECK_EQUAL(lgr_list[1], "LGR2");
|
||||
|
||||
BOOST_CHECK_EQUAL(grid1.is_radial(), false);
|
||||
|
||||
BOOST_CHECK_EQUAL(grid1.activeCells(), 30);
|
||||
BOOST_CHECK_EQUAL(grid1.totalNumberOfCells(), 30);
|
||||
|
||||
std::array<double,8> global_X_1_1_0 = { 2099.985, 2199.969, 2099.985, 2199.969, 2099.974, 2199.958, 2099.974, 2199.958 };
|
||||
std::array<double,8> global_Y_1_1_0 = { 2100.000, 2100.000, 2200.000, 2200.000, 2100.000, 2100.000, 2200.000, 2200.000 };
|
||||
std::array<double,8> global_Z_1_1_0 = { 2002.745, 2004.490, 2002.745, 2004.490, 2005.245, 2006.990, 2005.245, 2006.990 };
|
||||
|
||||
std::array<double,8> X = {0.0};
|
||||
std::array<double,8> Y = {0.0};
|
||||
std::array<double,8> Z = {0.0};
|
||||
|
||||
// cell 2,2,1 => zero based 1,1,0
|
||||
grid1.getCellCorners({1, 1, 0}, X, Y, Z);
|
||||
|
||||
for (size_t n = 0; n < 8; n++){
|
||||
BOOST_REQUIRE_CLOSE(global_X_1_1_0[n], X[n], 1e-3);
|
||||
BOOST_REQUIRE_CLOSE(global_Y_1_1_0[n], Y[n], 1e-3);
|
||||
BOOST_REQUIRE_CLOSE(global_Z_1_1_0[n], Z[n], 1e-3);
|
||||
}
|
||||
|
||||
Opm::EclIO::EInit init1(testInitFile);
|
||||
|
||||
{
|
||||
auto porv = init1.getInitData<float>("PORV");
|
||||
auto poro = init1.getInitData<float>("PORO");
|
||||
|
||||
BOOST_REQUIRE_CLOSE(calculateCellVol(X, Y, Z), porv[3] / poro[3], 1e-2);
|
||||
}
|
||||
|
||||
auto nnc_list = grid1.get_nnc_ijk();
|
||||
|
||||
BOOST_CHECK_EQUAL(nnc_list.size(), nnc_global_ref.size());
|
||||
|
||||
for (size_t n=0; n< nnc_list.size(); n++){
|
||||
BOOST_CHECK_EQUAL(std::get<0>(nnc_list[n]), std::get<0>(nnc_global_ref[n]));
|
||||
BOOST_CHECK_EQUAL(std::get<1>(nnc_list[n]), std::get<1>(nnc_global_ref[n]));
|
||||
BOOST_CHECK_EQUAL(std::get<2>(nnc_list[n]), std::get<2>(nnc_global_ref[n]));
|
||||
BOOST_CHECK_EQUAL(std::get<3>(nnc_list[n]), std::get<3>(nnc_global_ref[n]));
|
||||
BOOST_CHECK_EQUAL(std::get<4>(nnc_list[n]), std::get<4>(nnc_global_ref[n]));
|
||||
BOOST_CHECK_EQUAL(std::get<5>(nnc_list[n]), std::get<5>(nnc_global_ref[n]));
|
||||
BOOST_REQUIRE_CLOSE(std::get<6>(nnc_list[n]), std::get<6>(nnc_global_ref[n]), 1e-3);
|
||||
}
|
||||
|
||||
// testing LGR1 - cartesian grid
|
||||
|
||||
EGrid lgr1(testEgridFile, "LGR1");
|
||||
|
||||
BOOST_CHECK_EQUAL(lgr1.is_radial(), false);
|
||||
|
||||
auto nijk_lgr1 = lgr1.dimension();
|
||||
BOOST_CHECK_EQUAL(nijk_lgr1 == ref_dim_lgr1, true);
|
||||
|
||||
BOOST_CHECK_EQUAL(lgr1.activeCells(), 128);
|
||||
BOOST_CHECK_EQUAL(lgr1.totalNumberOfCells(), 128);
|
||||
|
||||
std::array<double,8> lgr1_X_1_2_0 = { 2124.983, 2149.979, 2124.984, 2149.980, 2124.978, 2149.974, 2124.979, 2149.975 };
|
||||
std::array<double,8> lgr1_Y_1_2_0 = { 2050.000, 2050.000, 2075.000, 2075.000, 2050.000, 2050.000, 2075.000, 2075.000 };
|
||||
std::array<double,8> lgr1_Z_1_2_0 = { 2002.182, 2002.618, 2002.182, 2002.618, 2003.431, 2003.868, 2003.431, 2003.868 };
|
||||
|
||||
// cell 2,3,1 => zero based 1,2,0
|
||||
lgr1.getCellCorners({1, 2, 0}, X, Y, Z);
|
||||
|
||||
for (size_t n = 0; n < 8; n++){
|
||||
BOOST_REQUIRE_CLOSE(lgr1_X_1_2_0[n], X[n], 1e-3);
|
||||
BOOST_REQUIRE_CLOSE(lgr1_Y_1_2_0[n], Y[n], 1e-3);
|
||||
BOOST_REQUIRE_CLOSE(lgr1_Z_1_2_0[n], Z[n], 1e-3);
|
||||
}
|
||||
|
||||
{
|
||||
int act_ind = lgr1.active_index(1,2,0);
|
||||
auto porv = init1.getInitData<float>("PORV", "LGR1");
|
||||
auto poro = init1.getInitData<float>("PORO", "LGR1");
|
||||
|
||||
BOOST_REQUIRE_CLOSE(calculateCellVol(X, Y, Z), porv[act_ind] / poro[act_ind], 1e-2);
|
||||
}
|
||||
|
||||
auto nnc_list_lgr1 = lgr1.get_nnc_ijk();
|
||||
|
||||
for (size_t n=0; n< nnc_list_lgr1.size(); n++){
|
||||
BOOST_CHECK_EQUAL(std::get<0>(nnc_list_lgr1[n]), std::get<0>(nnc_lgr1_ref[n]));
|
||||
BOOST_CHECK_EQUAL(std::get<1>(nnc_list_lgr1[n]), std::get<1>(nnc_lgr1_ref[n]));
|
||||
BOOST_CHECK_EQUAL(std::get<2>(nnc_list_lgr1[n]), std::get<2>(nnc_lgr1_ref[n]));
|
||||
BOOST_CHECK_EQUAL(std::get<3>(nnc_list_lgr1[n]), std::get<3>(nnc_lgr1_ref[n]));
|
||||
BOOST_CHECK_EQUAL(std::get<4>(nnc_list_lgr1[n]), std::get<4>(nnc_lgr1_ref[n]));
|
||||
BOOST_CHECK_EQUAL(std::get<5>(nnc_list_lgr1[n]), std::get<5>(nnc_lgr1_ref[n]));
|
||||
BOOST_REQUIRE_CLOSE(std::get<6>(nnc_list_lgr1[n]), std::get<6>(nnc_lgr1_ref[n]), 1e-3);
|
||||
}
|
||||
|
||||
|
||||
// testing LGR2 - radial grid
|
||||
|
||||
EGrid lgr2(testEgridFile, "LGR2");
|
||||
|
||||
BOOST_CHECK_EQUAL(lgr2.is_radial(), true);
|
||||
|
||||
auto nijk_lgr2 = lgr2.dimension();
|
||||
BOOST_CHECK_EQUAL(nijk_lgr2 == ref_dim_lgr2, true);
|
||||
|
||||
BOOST_CHECK_EQUAL(lgr2.activeCells(), 192);
|
||||
BOOST_CHECK_EQUAL(lgr2.totalNumberOfCells(), 192);
|
||||
|
||||
std::array<double,8> lgr2_X_4_1_0 = { 8.469, 23.561, 0.000, 0.000, 8.469, 23.561, 0.000, 0.000 };
|
||||
std::array<double,8> lgr2_Y_4_1_0 = { 8.471, 23.564, 11.978, 33.322, 8.471, 23.564, 11.978, 33.322 };
|
||||
std::array<double,8> lgr2_Z_4_1_0 = { 2011.892, 2012.155, 2011.744, 2011.744, 2013.892, 2014.155, 2013.744, 2013.744 };
|
||||
|
||||
// cell 5,2,1 => zero based 4,1,0
|
||||
lgr2.getCellCorners({4, 1, 0}, X, Y, Z);
|
||||
|
||||
for (size_t n = 0; n < 8; n++){
|
||||
BOOST_CHECK_EQUAL(abs(lgr2_X_4_1_0[n]- X[n]) < 1e-3, true);
|
||||
BOOST_CHECK_EQUAL(abs(lgr2_Y_4_1_0[n]- Y[n]) < 1e-3, true);
|
||||
BOOST_CHECK_EQUAL(abs(lgr2_Z_4_1_0[n]- Z[n]) < 1e-3, true);
|
||||
}
|
||||
|
||||
auto hostcells_ijk = lgr1.hostCellsIJK();
|
||||
auto hostcells_gind = lgr1.hostCellsGlobalIndex();
|
||||
|
||||
for (size_t n = 0; n < hostcells_gind.size(); n++)
|
||||
BOOST_CHECK_EQUAL(grid1.ijk_from_global_index(hostcells_gind[n]) == hostcells_ijk[n], true);
|
||||
}
|
||||
|
||||
|
||||
Loading…
Reference in New Issue
Block a user