diff --git a/opm/io/eclipse/EGrid.cpp b/opm/io/eclipse/EGrid.cpp index 926d0183d58..03685ecddc6 100644 --- a/opm/io/eclipse/EGrid.cpp +++ b/opm/io/eclipse/EGrid.cpp @@ -37,7 +37,7 @@ namespace Opm { namespace EclIO { using NNCentry = std::tuple; -EGrid::EGrid(const std::string& filename, const std::string& grid_name) +EGrid::EGrid(const std::string& filename, const std::string& grid_name, bool apply_mapaxes) : EclFile(filename), inputFileName { filename }, m_grid_name {grid_name} { initFileName = inputFileName.parent_path() / inputFileName.stem(); @@ -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; @@ -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(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(n); + if ((array_name[n] == "MAPAXES") && (apply_mapaxes == true)){ + auto mapaxes = this->get(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") { @@ -322,11 +340,31 @@ std::array 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& ijk, - std::array& X, - std::array& Y, - std::array& Z) + std::array& X, + std::array& Y, + std::array& Z) { if (coord_array.empty()) load_grid_data(); @@ -351,7 +389,7 @@ void EGrid::getCellCorners(const std::array& 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++) { @@ -366,8 +404,8 @@ void EGrid::getCellCorners(const std::array& 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]; @@ -390,11 +428,20 @@ void EGrid::getCellCorners(const std::array& ijk, } +void EGrid::getCellCornersWithMapAxes(const std::array& ijk, std::array& X, + std::array& Y, std::array& 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& X, - std::array& Y, std::array& Z) +void EGrid::getCellCorners(int globindex, std::array& X, + std::array& Y, std::array& Z) { - return getCellCorners(ijk_from_global_index(globindex),X,Y,Z); + return getCellCorners(ijk_from_global_index(globindex), X, Y, Z); } diff --git a/opm/io/eclipse/EGrid.hpp b/opm/io/eclipse/EGrid.hpp index 693237e9129..5fbde1d6458 100644 --- a/opm/io/eclipse/EGrid.hpp +++ b/opm/io/eclipse/EGrid.hpp @@ -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, const std::string& grid_name = "global", bool apply_mapaxes = false); int global_index(int i, int j, int k) const; int active_index(int i, int j, int k) const; @@ -42,8 +42,9 @@ class EGrid : public EclFile std::array ijk_from_active_index(int actInd) const; std::array ijk_from_global_index(int globInd) const; - void getCellCorners(int globindex, std::array& X, std::array& Y, std::array& Z); - void getCellCorners(const std::array& ijk, std::array& X, std::array& Y, std::array& Z); + void getCellCorners(int globindex, std::array& X, std::array& Y, std::array& Z); + void getCellCorners(const std::array& ijk, std::array& X, std::array& Y, std::array& Z); + void getCellCornersWithMapAxes(const std::array& ijk, std::array& X, std::array& Y, std::array& Z); std::vector> getXYZ_layer(int layer, bool bottom=false); std::vector> getXYZ_layer(int layer, const std::array& box, bool bottom=false); @@ -53,18 +54,21 @@ 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& hostCellsGlobalIndex() const { return host_cells; } std::vector> hostCellsIJK(); - // zero based: i1, j1,k1, i2,j2,k2, transmisibility + // zero based: i1,j1,k1, i2,j2,k2, transmisibility using NNCentry = std::tuple; std::vector get_nnc_ijk(); const std::vector& list_of_lgrs() const { return lgr_names; } - const std::vector& get_mapaxes() const { return m_mapaxes; } + const std::array& get_mapaxes() const { return m_mapaxes; } const std::string& get_mapunits() const { return m_mapunits; } private: @@ -72,8 +76,13 @@ class EGrid : public EclFile std::string m_grid_name; bool m_radial; - std::vector m_mapaxes; + std::array m_mapaxes; std::string m_mapunits; + mutable bool m_mapaxes_loaded; + std::array origin; + std::array unit_x; + std::array unit_y; + double length_factor; std::array nijk; std::array host_nijk; @@ -107,7 +116,7 @@ class EGrid : public EclFile std::vector get_zcorn_from_disk(int layer, bool bottom); void getCellCorners(const std::array& ijk, const std::vector& zcorn_layer, - std::array& X, std::array& Y, std::array& Z); + std::array& X, std::array& Y, std::array& Z); }; diff --git a/python/cxx/eclipse_io.cpp b/python/cxx/eclipse_io.cpp index 8c724dcf31a..8d47faf88bd 100644 --- a/python/cxx/eclipse_io.cpp +++ b/python/cxx/eclipse_io.cpp @@ -294,24 +294,27 @@ npArray get_erst_vector(Opm::EclIO::ERst * file_ptr, const std::string& key, siz std::tuple, std::array, std::array> -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 X = {0.0}; std::array Y = {0.0}; std::array Z = {0.0}; - std::array ijk = {i, j, k }; + std::array 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, std::array> get_xyz_from_active_index(Opm::EclIO::EGrid * file_ptr, int actIndex) { std::array 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 mask) @@ -463,13 +466,15 @@ void python::common::export_IO(py::module& m) { .def("units", &ESmryBind::units); py::class_(m, "EGrid") - .def(py::init()) + .def(py::init(), py::arg("filename"), + py::arg("grid_name") = "global", 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) .def("ijk_from_active_index", &Opm::EclIO::EGrid::ijk_from_active_index) .def("active_index", &Opm::EclIO::EGrid::active_index) .def("global_index", &Opm::EclIO::EGrid::global_index) + .def("export_mapaxes", &Opm::EclIO::EGrid::get_mapaxes) .def("xyz_from_ijk", &get_xyz_from_ijk) .def("xyz_from_active_index", &get_xyz_from_active_index) .def("cellvolumes", &get_cellvolumes) diff --git a/python/opm_embedded/__init__.pyi b/python/opm_embedded/__init__.pyi index 913628e6a93..b98665f4a00 100644 --- a/python/opm_embedded/__init__.pyi +++ b/python/opm_embedded/__init__.pyi @@ -2480,12 +2480,13 @@ class Dimension: def scaling(self) -> float: ... class EGrid: - def __init__(self, arg0: str) -> None: ... + def __init__(self, filename: str, grid_name: str, apply_mapaxes: bool) -> None: ... def active_index(self, arg0: int, arg1: int, arg2: int) -> int: ... @overload def cellvolumes(self) -> numpy.ndarray: ... @overload def cellvolumes(self, arg0: List[int]) -> numpy.ndarray: ... + def export_mapaxes(self) -> numpy.ndarray: ... def global_index(self, arg0: int, arg1: int, arg2: int) -> int: ... def ijk_from_active_index(self, arg0: int) -> List[int[3]]: ... def ijk_from_global_index(self, arg0: int) -> List[int[3]]: ... diff --git a/python/tests/data/MAPAXES.EGRID b/python/tests/data/MAPAXES.EGRID new file mode 100644 index 00000000000..53d152b2d6c Binary files /dev/null and b/python/tests/data/MAPAXES.EGRID differ diff --git a/python/tests/test_egrid.py b/python/tests/test_egrid.py index d0e09e0e1c7..f18520b15d7 100755 --- a/python/tests/test_egrid.py +++ b/python/tests/test_egrid.py @@ -1,7 +1,5 @@ import unittest -import sys import numpy as np -import datetime from opm.io.ecl import EGrid try: @@ -22,67 +20,98 @@ def test_ijk_active_and_global_indices(self): self.assertEqual(grid1.active_cells, 2794) - nI,nJ,nK = grid1.dimension - tot_ant_cells = nI*nJ*nK + nI, nJ, nK = grid1.dimension + tot_ant_cells = nI * nJ * nK - self.assertEqual( nI, 13) - self.assertEqual( nJ, 22) - self.assertEqual( nK, 11) + self.assertEqual(nI, 13) + self.assertEqual(nJ, 22) + self.assertEqual(nK, 11) - i,j,k = grid1.ijk_from_global_index(0) + i, j, k = grid1.ijk_from_global_index(0) - self.assertEqual( i, 0) - self.assertEqual( j, 0) - self.assertEqual( k, 0) + self.assertEqual(i, 0) + self.assertEqual(j, 0) + self.assertEqual(k, 0) - i,j,k = grid1.ijk_from_global_index(1000) + i, j, k = grid1.ijk_from_global_index(1000) - self.assertEqual( i, 12) - self.assertEqual( j, 10) - self.assertEqual( k, 3) + self.assertEqual(i, 12) + self.assertEqual(j, 10) + self.assertEqual(k, 3) - self.assertEqual( grid1.global_index(12, 10, 3), 1000 ) + self.assertEqual(grid1.global_index(12, 10, 3), 1000) - i,j,k = grid1.ijk_from_global_index(tot_ant_cells - 1) + i, j, k = grid1.ijk_from_global_index(tot_ant_cells - 1) - self.assertEqual( i, nI -1 ) - self.assertEqual( j, nJ -1 ) - self.assertEqual( k, nK -1 ) + self.assertEqual(i, nI - 1) + self.assertEqual(j, nJ - 1) + self.assertEqual(k, nK - 1) with self.assertRaises(ValueError): i,j,k = grid1.ijk_from_global_index(tot_ant_cells) - i,j,k = grid1.ijk_from_active_index( 1000 ) + i,j,k = grid1.ijk_from_active_index(1000) - self.assertEqual( i, 1 ) - self.assertEqual( j, 15 ) - self.assertEqual( k, 3 ) + self.assertEqual(i, 1) + self.assertEqual(j, 15) + self.assertEqual(k, 3) - self.assertEqual( grid1.active_index(1, 15, 3), 1000 ) + self.assertEqual( grid1.active_index(1, 15, 3), 1000) + + + def test_mapaxes(self): + + axsref = [304.8, 152.4, 304.8, 457.2, 609.6, 457.2] + X1ref = [304.8, 305.8, 304.8, 305.8, 304.8, 305.8, 304.8, 305.8] + Y1ref = [457.2, 457.2, 456.2, 456.2, 457.2, 457.2, 456.2, 456.2] + + grid1 = EGrid(test_path("data/MAPAXES.EGRID")) + + X1, Y1, _ = grid1.xyz_from_ijk(0, 0, 0) + + for n in range(0, 8): + self.assertAlmostEqual(X1ref[n], X1[n], 8) + self.assertAlmostEqual(Y1ref[n], Y1[n], 8) + + axs = grid1.export_mapaxes() + + for n in range(0, 6): + self.assertAlmostEqual(axsref[n], axs[n], 8) + + X2ref = [0., 1., 0., 1., 0., 1., 0., 1.] + Y2ref = [0., 0., 1., 1., 0., 0., 1., 1.] + + grid2 = EGrid(test_path("data/MAPAXES.EGRID"), apply_mapaxes = False) + + X2, Y2, _ = grid2.xyz_from_ijk(0, 0, 0) + + for n in range(0, 8): + self.assertAlmostEqual(X2ref[n], X2[n], 8) + self.assertAlmostEqual(Y2ref[n], Y2[n], 8) def test_coordinates(self): - Xref=[2899.45166015625, 2999.390869140625, 2899.45166015625, 2999.390869140625, - 2899.4176237656716, 2999.3568089317187, 2899.417623015281, 2999.356808099622] - Yref=[2699.973388671875, 2699.973388671875, 2799.969482421875, 2799.969482421875, - 2699.9818918149376, 2699.9818918149376, 2799.978009571257, 2799.9780095915808] - Zref=[2565.301025390625, 2568.791015625, 2564.42822265625, 2567.918212890625, - 2575.29443359375, 2578.784423828125, 2574.421875, 2577.911865234375] + Xref = [2899.45166015625, 2999.390869140625, 2899.45166015625, 2999.390869140625, + 2899.4176237656716, 2999.3568089317187, 2899.417623015281, 2999.356808099622] + Yref = [2699.973388671875, 2699.973388671875, 2799.969482421875, 2799.969482421875, + 2699.9818918149376, 2699.9818918149376, 2799.978009571257, 2799.9780095915808] + Zref = [2565.301025390625, 2568.791015625, 2564.42822265625, 2567.918212890625, + 2575.29443359375, 2578.784423828125, 2574.421875, 2577.911865234375] grid1 = EGrid(test_path("data/9_EDITNNC.EGRID")) X1, Y1, Z1 = grid1.xyz_from_ijk(9, 7, 0) - for n in range(0,8): + for n in range(0, 8): self.assertAlmostEqual(Xref[n], X1[n], 8) 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): + for n in range(0, 8): self.assertAlmostEqual(Xref[n], X2[n], 8) self.assertAlmostEqual(Yref[n], Y2[n], 8) self.assertAlmostEqual(Zref[n], Z2[n], 8) @@ -92,8 +121,8 @@ def test_cell_volume(self): grid1 = EGrid(test_path("data/9_EDITNNC.EGRID")) - nI,nJ,nK = grid1.dimension - tot_ant_cells = nI*nJ*nK + nI, nJ, nK = grid1.dimension + tot_ant_cells = nI * nJ * nK celVolAll = grid1.cellvolumes() @@ -102,16 +131,16 @@ def test_cell_volume(self): self.assertTrue(min(celVolAll) > 0.0) - mask =[0]*tot_ant_cells + mask =[0] * tot_ant_cells for k in range(nK): - globInd=grid1.global_index(0,0,k) + globInd = grid1.global_index(0, 0, k) mask[globInd] = 1 celVol = grid1.cellvolumes(mask) self.assertTrue(min(celVol) == 0.0) - self.assertEqual(np.count_nonzero(celVol) , nK) + self.assertEqual(np.count_nonzero(celVol), nK) if __name__ == "__main__":