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

added: EGrid class for reading GRID data from ECL files

Browse Source
This commit is contained in:
Torbjørn Skille 2019-03-28 14:33:30 +01:00 committed by Arne Morten Kvarving
parent 46f2542487
commit 181eeffbdb
6 changed files with 433 additions and 1 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
CMakeLists.txt
View File

@ -164,6 +164,7 @@ if(ENABLE_ECL_INPUT)
examples/test_util/EclFile.cpp
examples/test_util/EclOutput.cpp
examples/test_util/EclUtil.cpp
examples/test_util/EGrid.cpp
examples/test_util/summaryComparator.cpp
examples/test_util/summaryIntegrationTest.cpp
examples/test_util/summaryRegressionTest.cpp)
@ -174,7 +175,8 @@ if(ENABLE_ECL_INPUT)
# Add the tests
set(_libs testutil opmcommon
${Boost_UNIT_TEST_FRAMEWORK_LIBRARY})
foreach(test test_compareSummary test_EclFilesComparator test_EclIO)
foreach(test test_compareSummary test_EclFilesComparator test_EclIO
test_EGrid)
opm_add_test(${test} CONDITION ENABLE_ECL_INPUT
LIBRARIES ${_libs}
WORKING_DIRECTORY ${PROJECT_BINARY_DIR}/tests)

1
CMakeLists_files.cmake
View File

@ -321,6 +321,7 @@ if(ENABLE_ECL_INPUT)
list (APPEND TEST_DATA_FILES
tests/ECLFILE.INIT
tests/ECLFILE.FINIT
tests/SPE1CASE1.EGRID
)
list (APPEND EXAMPLE_SOURCE_FILES
examples/opmi.cpp

186
examples/test_util/EGrid.cpp Normal file
View File

@ -0,0 +1,186 @@
/*
Copyright 2019 Equinor ASA.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include "EGrid.hpp"
#include <opm/common/ErrorMacros.hpp>
#include <string>
#include <string.h>
#include <sstream>
#include <iterator>
#include <iomanip>
#include <algorithm>
#include <numeric>
EGrid::EGrid(const std::string &filename) : EclFile(filename)
{
EclFile file(filename);
auto gridhead = get<int>("GRIDHEAD");
nijk[0] = gridhead[1];
nijk[1] = gridhead[2];
nijk[2] = gridhead[3];
if (file.hasKey("ACTNUM")) {
auto actnum = get<int>("ACTNUM");
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);
}
}
} 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");
}
int EGrid::global_index(int i, int j, int k) const
{
if (i < 0 || i >= nijk[0] || j < 0 || j >= nijk[1] || k < 0 || k >= nijk[2]) {
OPM_THROW(std::invalid_argument, "i, j or/and k out of range");
}
return i + j * nijk[0] + k * nijk[0] * nijk[1];
}
int EGrid::active_index(int i, int j, int k) const
{
int n = i + j * nijk[0] + k * nijk[0] * nijk[1];
if (i < 0 || i >= nijk[0] || j < 0 || j >= nijk[1] || k < 0 || k >= nijk[2]) {
OPM_THROW(std::invalid_argument, "i, j or/and k out of range");
}
return act_index[n];
}
std::array<int, 3> EGrid::ijk_from_active_index(int actInd) const
{
if (actInd < 0 || actInd >= nactive) {
OPM_THROW(std::invalid_argument, "active index out of range");
}
int _glob = glob_index[actInd];
std::array<int, 3> result;
result[2] = _glob / (nijk[0] * nijk[1]);
int rest = _glob % (nijk[0] * nijk[1]);
result[1] = rest / nijk[0];
result[0] = rest % nijk[0];
return result;
}
std::array<int, 3> EGrid::ijk_from_global_index(int globInd) const
{
if (globInd < 0 || globInd >= nijk[0] * nijk[1] * nijk[2]) {
OPM_THROW(std::invalid_argument, "global index out of range");
}
std::array<int, 3> result;
result[2] = globInd / (nijk[0] * nijk[1]);
int rest = globInd % (nijk[0] * nijk[1]);
result[1] = rest / nijk[0];
result[0] = rest % nijk[0];
return result;
}
void EGrid::getCellCorners(const std::array<int, 3>& ijk,
std::vector<double>& X,
std::vector<double>& Y,
std::vector<double>& Z) const
{
std::vector<int> zind;
std::vector<int> pind;
if (X.size() < 8 || Y.size() < 8 || Z.size() < 8) {
OPM_THROW(std::invalid_argument,
"In routine cellConrner. X, Y and Z should be a vector of size 8");
}
if (ijk[0] < 0 || ijk[0] >= nijk[0] || ijk[1] < 0 || ijk[1] >= nijk[1] || ijk[2] < 0 || ijk[2] >= nijk[2]) {
OPM_THROW(std::invalid_argument, "i, j and/or k out of range");
}
// 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);
// 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< 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];
double xb = coord_array[pind[n] + 3];
double yb = coord_array[pind[n] + 4];
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]);
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::vector<double>& X,
std::vector<double>& Y, std::vector<double>& Z) const
{
return getCellCorners(ijk_from_global_index(globindex),X,Y,Z);
}

62
examples/test_util/EGrid.hpp Normal file
View File

@ -0,0 +1,62 @@
/*
Copyright 2019 Equinor ASA.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef EGRID_HPP
#define EGRID_HPP
#include "EclFile.hpp"
#include <array>
#include <iostream>
#include <string>
#include <fstream>
#include <vector>
#include <ctime>
#include <map>
class EGrid : public EclFile
{
public:
explicit EGrid(const std::string& filename);
int global_index(int i, int j, int k) const;
int active_index(int i, int j, int k) const;
const std::array<int, 3>& dimension() const { return nijk; }
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::vector<double>& X, std::vector<double>& Y, std::vector<double>& Z) const;
void getCellCorners(const std::array<int, 3>& ijk, std::vector<double>& X, std::vector<double>& Y, std::vector<double>& Z) const;
int activeCells() const { return nactive; }
int totalNumberOfCells() const { return nijk[0] * nijk[1] * nijk[2]; }
private:
std::array<int, 3> nijk;
int nNNC, nactive;
std::vector<int> act_index;
std::vector<int> glob_index;
std::vector<float> coord_array;
std::vector<float> zcorn_array;
};
#endif

BIN
tests/SPE1CASE1.EGRID Normal file
View File

Binary file not shown.

181
tests/test_EGrid.cpp Normal file
View File

@ -0,0 +1,181 @@
/*
+ Copyright 2019 Equinor ASA.
+
+ This file is part of the Open Porous Media project (OPM).
+
+ OPM is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ OPM is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with OPM. If not, see <http://www.gnu.org/licenses/>.
+ */
#include "config.h"
#include <iostream>
#include <stdio.h>
#include <iomanip>
#include <math.h>
#include <algorithm>
#include <tuple>
#include <examples/test_util/EGrid.hpp>
#define BOOST_TEST_MODULE Test EGrid
#include <boost/test/unit_test.hpp>
template<typename InputIterator1, typename InputIterator2>
bool
range_equal(InputIterator1 first1, InputIterator1 last1,
InputIterator2 first2, InputIterator2 last2)
{
while(first1 != last1 && first2 != last2)
{
if(*first1 != *first2) return false;
++first1;
++first2;
}
return (first1 == last1) && (first2 == last2);
}
bool compare_files(const std::string& filename1, const std::string& filename2)
{
std::ifstream file1(filename1);
std::ifstream file2(filename2);
std::istreambuf_iterator<char> begin1(file1);
std::istreambuf_iterator<char> begin2(file2);
std::istreambuf_iterator<char> end;
return range_equal(begin1, end, begin2, end);
}
template <typename T>
bool operator==(const std::vector<T> & t1, const std::vector<T> & t2)
{
return std::equal(t1.begin(), t1.end(), t2.begin(), t2.end());
}
// 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
// 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
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);
auto dim = grid1.dimension();
BOOST_CHECK_EQUAL(dim[0],10);
BOOST_CHECK_EQUAL(dim[1],10);
BOOST_CHECK_EQUAL(dim[2],3);
int globInd = grid1.global_index(3, 2, 1);
BOOST_CHECK_EQUAL(globInd, 123); // 10*10*1 + 10*2 + 3 = 100+20+3 = 123
BOOST_CHECK_EQUAL(grid1.global_index(0, 0, 0), 0);
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);
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);
BOOST_CHECK_EQUAL(actInd, 117); // global index 123, - 6 inactive
BOOST_CHECK_EQUAL(grid1.active_index(0, 0, 0), 0);
BOOST_CHECK_EQUAL(grid1.active_index(dim[0] - 1, dim[1] - 1, dim[2] - 1), nAct - 1);
// check active_index valid range, should throw exception if outside
BOOST_CHECK_THROW(grid1.active_index(-1, 0, 0), std::invalid_argument);
BOOST_CHECK_THROW(grid1.active_index(-1, -1, 0), std::invalid_argument);
BOOST_CHECK_THROW(grid1.active_index(-1, -1, -1), std::invalid_argument);
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);
BOOST_CHECK_THROW(grid1.ijk_from_global_index(-1), std::invalid_argument);
BOOST_CHECK_THROW(grid1.ijk_from_global_index(nTot), std::invalid_argument);
}
BOOST_AUTO_TEST_CASE(getCellCorners) {
std::string testFile="SPE1CASE1.EGRID";
std::vector<double> ref_X = {3000,4000,3000,4000,3000,4000,3000,4000};
std::vector<double> ref_Y = {2000,2000,3000,3000,2000,2000,3000,3000};
std::vector<double> ref_Z = {8345,8345,8345,8345,8375,8375,8375,8375};
EGrid grid1(testFile);
std::vector<double> X(8,0.0);
std::vector<double> Y(8,0.0);
std::vector<double> Z(8,0.0);
// cell 4,3,2 => zero based 3,2,1
grid1.getCellCorners({3, 2, 1}, X, Y, Z);
BOOST_CHECK_EQUAL(X == ref_X, true);
BOOST_CHECK_EQUAL(Y == ref_Y, true);
BOOST_CHECK_EQUAL(Z == ref_Z, true);
X.assign(8,0.0);
Y.assign(8,0.0);
Z.assign(8,0.0);
int globInd=grid1.global_index(3,2,1);
grid1.getCellCorners(globInd,X,Y,Z);
BOOST_CHECK_EQUAL(X == ref_X, true);
BOOST_CHECK_EQUAL(Y == ref_Y, true);
BOOST_CHECK_EQUAL(Z == ref_Z, true);
}