Skip to content

Commit

Permalink
Handle MAPAXES and MAPUNITS
Browse files Browse the repository at this point in the history
  • Loading branch information
@daavid00
daavid00 committed Oct 14, 2024
1 parent cfdc43a commit 68c61ad
Show file tree
Hide file tree
Showing 6 changed files with 97 additions and 23 deletions.
67 changes: 57 additions & 10 deletions opm/io/eclipse/EGrid.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -37,7 +37,7 @@ namespace Opm { namespace EclIO {

using NNCentry = std::tuple<int, int, int, int, int,int, float>;

EGrid::EGrid(const std::string& filename, const std::string& grid_name)
EGrid::EGrid(const std::string& filename, bool apply_mapaxes, const std::string& grid_name)
: EclFile(filename), inputFileName { filename }, m_grid_name {grid_name}
{
initFileName = inputFileName.parent_path() / inputFileName.stem();
Expand All @@ -55,6 +55,8 @@ EGrid::EGrid(const std::string& filename, const std::string& grid_name)
nnc2_array_index = -1;
coordsys_array_index = -1;
m_radial = false;
m_mapaxes_loaded = false;
length_factor = 1.0;

int hostnum_index = -1;

Expand All @@ -78,13 +80,29 @@ EGrid::EGrid(const std::string& filename, const std::string& grid_name)
lgrname = lgr_names[nnchead[1] - 1];
}

if (array_name[n] == "MAPUNITS"){
if ((array_name[n] == "MAPUNITS") && (apply_mapaxes == true)){
auto mapunits = this->get<std::string>(n);
m_mapunits = mapunits[0];
if (m_mapunits == "METRES")
length_factor = 1.0;
else if (m_mapunits == "FEET")
length_factor = 0.3048;
else if (m_mapunits == "CM")
length_factor = 0.01;
else{
std::string message = "Unit system " + m_mapunits + " not supported for MAPUNITS";
OPM_THROW(std::invalid_argument, message);
}
}

if (array_name[n] == "MAPAXES")
m_mapaxes = this->get<float>(n);
if ((array_name[n] == "MAPAXES") && (apply_mapaxes == true)){
auto mapaxes = this->get<float>(n);
for (size_t i = 0; i < 6; i++) {
m_mapaxes[i] = length_factor * mapaxes[i];
}
mapaxes_init();
m_mapaxes_loaded = true;
}

if (lgrname == grid_name) {
if (array_name[n] == "GRIDHEAD") {
Expand Down Expand Up @@ -322,11 +340,31 @@ std::array<int, 3> EGrid::ijk_from_global_index(int globInd) const
return result;
}

void EGrid::mapaxes_transform(double& x, double& y) const {
double tmpx = x;
x = origin[0] + tmpx * unit_x[0] + y * unit_y[0];
y = origin[1] + tmpx * unit_x[1] + y * unit_y[1];
}

void EGrid::mapaxes_init()
{
origin = {m_mapaxes[2], m_mapaxes[3]};
unit_x = {m_mapaxes[4] - m_mapaxes[2], m_mapaxes[5] - m_mapaxes[3]};
unit_y = {m_mapaxes[0] - m_mapaxes[2], m_mapaxes[1] - m_mapaxes[3]};

auto norm_x = 1.0 / std::hypot(unit_x[0], unit_x[1]);
auto norm_y = 1.0 / std::hypot(unit_y[0], unit_y[1]);

unit_x[0] *= norm_x;
unit_x[1] *= norm_x;
unit_y[0] *= norm_y;
unit_y[1] *= norm_y;
}

void EGrid::getCellCorners(const std::array<int, 3>& ijk,
std::array<double,8>& X,
std::array<double,8>& Y,
std::array<double,8>& Z)
std::array<double, 8>& X,
std::array<double, 8>& Y,
std::array<double, 8>& Z)
{
if (coord_array.empty())
load_grid_data();
Expand All @@ -351,7 +389,7 @@ void EGrid::getCellCorners(const std::array<int, 3>& ijk,
for (int n = 0; n < 4; n++)
zind.push_back(zind[n] + nijk[0]*nijk[1]*4);

for (int n = 0; n< 8; n++)
for (int n = 0; n < 8; n++)
Z[n] = zcorn_array[zind[n]];

for (int n = 0; n < 4; n++) {
Expand All @@ -366,8 +404,8 @@ void EGrid::getCellCorners(const std::array<int, 3>& ijk,
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);
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];
Expand All @@ -390,6 +428,15 @@ void EGrid::getCellCorners(const std::array<int, 3>& ijk,
}


void EGrid::getCellCornersWithMapAxes(const std::array<int, 3>& ijk, std::array<double,8>& X,
std::array<double,8>& Y, std::array<double,8>& Z)
{
this->getCellCorners(ijk, X, Y, Z );

for (int n = 0; n < 8; n++)
mapaxes_transform(X[n], Y[n]);
}


void EGrid::getCellCorners(int globindex, std::array<double,8>& X,
std::array<double,8>& Y, std::array<double,8>& Z)
Expand Down
17 changes: 13 additions & 4 deletions opm/io/eclipse/EGrid.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -32,7 +32,7 @@ namespace Opm { namespace EclIO {
class EGrid : public EclFile
{
public:
explicit EGrid(const std::string& filename, const std::string& grid_name = "global");
explicit EGrid(const std::string& filename, bool apply_mapaxes = false, const 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;
Expand All @@ -44,6 +44,7 @@ class EGrid : public EclFile

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);
void getCellCornersWithMapAxes(const std::array<int, 3>& ijk, std::array<double,8>& X, std::array<double,8>& Y, std::array<double,8>& Z);

std::vector<std::array<float, 3>> getXYZ_layer(int layer, bool bottom=false);
std::vector<std::array<float, 3>> getXYZ_layer(int layer, const std::array<int, 4>& box, bool bottom=false);
Expand All @@ -53,27 +54,35 @@ class EGrid : public EclFile

void load_grid_data();
void load_nnc_data();
void mapaxes_init();
void mapaxes_transform(double& x, double& y) const;
bool with_mapaxes() const { return m_mapaxes_loaded; }
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
// 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; }

const std::vector<float>& get_mapaxes() const { return m_mapaxes; }
const std::array<double, 6>& get_mapaxes() const { return m_mapaxes; }
const std::string& get_mapunits() const { return m_mapunits; }

private:
std::filesystem::path inputFileName, initFileName;
std::string m_grid_name;
bool m_radial;

std::vector<float> m_mapaxes;
std::array<double, 6> m_mapaxes;
std::string m_mapunits;
mutable bool m_mapaxes_loaded;
std::array<double, 4> origin;
std::array<double, 2> unit_x;
std::array<double, 2> unit_y;
double length_factor;

std::array<int, 3> nijk;
std::array<int, 3> host_nijk;
Expand Down
17 changes: 10 additions & 7 deletions python/cxx/eclipse_io.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -294,24 +294,27 @@ npArray get_erst_vector(Opm::EclIO::ERst * file_ptr, const std::string& key, siz


std::tuple<std::array<double,8>, std::array<double,8>, std::array<double,8>>
get_xyz_from_ijk(Opm::EclIO::EGrid * file_ptr,int i, int j, int k)
get_xyz_from_ijk(Opm::EclIO::EGrid * file_ptr, int i, int j, int k)
{
std::array<double,8> X = {0.0};
std::array<double,8> Y = {0.0};
std::array<double,8> Z = {0.0};

std::array<int, 3> ijk = {i, j, k };
std::array<int, 3> ijk = {i, j, k};

if (file_ptr->with_mapaxes())
file_ptr->getCellCornersWithMapAxes(ijk, X, Y, Z);
else
file_ptr->getCellCorners(ijk, X, Y, Z);

file_ptr->getCellCorners(ijk, X, Y, Z);

return std::make_tuple( X, Y, Z);
return std::make_tuple(X, Y, Z);
}

std::tuple<std::array<double,8>, std::array<double,8>, std::array<double,8>>
get_xyz_from_active_index(Opm::EclIO::EGrid * file_ptr, int actIndex)
{
std::array<int, 3> ijk = file_ptr->ijk_from_active_index(actIndex);
return get_xyz_from_ijk(file_ptr,ijk[0], ijk[1], ijk[2]);
return get_xyz_from_ijk(file_ptr, ijk[0], ijk[1], ijk[2]);
}

py::array get_cellvolumes_mask(Opm::EclIO::EGrid * file_ptr, std::vector<int> mask)
Expand Down Expand Up @@ -463,7 +466,7 @@ void python::common::export_IO(py::module& m) {
.def("units", &ESmryBind::units);

py::class_<Opm::EclIO::EGrid>(m, "EGrid")
.def(py::init<const std::string &>())
.def(py::init<const std::string &, const bool>(), py::arg("filename"), py::arg("apply_mapaxes") = true)
.def_property_readonly("active_cells", &Opm::EclIO::EGrid::activeCells)
.def_property_readonly("dimension", &Opm::EclIO::EGrid::dimension)
.def("ijk_from_global_index", &Opm::EclIO::EGrid::ijk_from_global_index)
Expand Down
2 changes: 1 addition & 1 deletion python/opm_embedded/__init__.pyi
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2480,7 +2480,7 @@ class Dimension:
def scaling(self) -> float: ...

class EGrid:
def __init__(self, arg0: str) -> None: ...
def __init__(self, arg0: str, arg1: str, apply_mapaxes: bool) -> None: ...
def active_index(self, arg0: int, arg1: int, arg2: int) -> int: ...
@overload
def cellvolumes(self) -> numpy.ndarray: ...
Expand Down
Binary file added python/tests/data/MAPAXES.EGRID
View file
Binary file not shown.
17 changes: 16 additions & 1 deletion python/tests/test_egrid.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -61,6 +61,21 @@ def test_ijk_active_and_global_indices(self):
self.assertEqual( grid1.active_index(1, 15, 3), 1000 )


def test_mapaxes(self):

Xref=[304.8, 305.8, 304.8, 305.8, 304.8, 305.8, 304.8, 305.8]
Yref=[457.2, 457.2, 456.2, 456.2, 457.2, 457.2, 456.2, 456.2]
Zref=[0., 0., 0., 0., 1., 1., 1., 1.]

grid = EGrid(test_path("data/MAPAXES.EGRID"))

X1, Y1, Z1 = grid.xyz_from_ijk(0, 0, 0)

for n in range(0,8):
self.assertAlmostEqual(Xref[n], X1[n], 2)
self.assertAlmostEqual(Yref[n], Y1[n], 2)
self.assertAlmostEqual(Zref[n], Z1[n], 2)

def test_coordinates(self):

Xref=[2899.45166015625, 2999.390869140625, 2899.45166015625, 2999.390869140625,
Expand All @@ -79,7 +94,7 @@ def test_coordinates(self):
self.assertAlmostEqual(Yref[n], Y1[n], 8)
self.assertAlmostEqual(Zref[n], Z1[n], 8)

actInd = grid1.active_index(9, 7, 0);
actInd = grid1.active_index(9, 7, 0)
X2, Y2, Z2 = grid1.xyz_from_active_index(actInd)

for n in range(0,8):
Expand Down

0 comments on commit 68c61ad

Please sign in to comment.