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 16, 2024
1 parent cfdc43a commit 3084f81
Show file tree
Hide file tree
Showing 7 changed files with 150 additions and 69 deletions.
73 changes: 60 additions & 13 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,11 +428,20 @@ 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)
void EGrid::getCellCorners(int globindex, std::array<double, 8>& X,
std::array<double, 8>& Y, std::array<double, 8>& Z)
{
return getCellCorners(ijk_from_global_index(globindex),X,Y,Z);
return getCellCorners(ijk_from_global_index(globindex), X, Y, Z);
}


Expand Down
23 changes: 16 additions & 7 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 @@ -42,8 +42,9 @@ class EGrid : public EclFile
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);
void getCellCorners(const std::array<int, 3>& ijk, std::array<double,8>& X, std::array<double,8>& Y, std::array<double,8>& Z);
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 Expand Up @@ -107,7 +116,7 @@ class EGrid : public EclFile
std::vector<float> get_zcorn_from_disk(int layer, bool bottom);

void getCellCorners(const std::array<int, 3>& ijk, const std::vector<float>& zcorn_layer,
std::array<double,4>& X, std::array<double,4>& Y, std::array<double,4>& Z);
std::array<double, 4>& X, std::array<double, 4>& Y, std::array<double, 4>& Z);

};

Expand Down
18 changes: 11 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,13 +466,14 @@ 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)
.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)
Expand Down
3 changes: 2 additions & 1 deletion python/opm_embedded/__init__.pyi
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2480,12 +2480,13 @@ 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: ...
@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]]: ...
Expand Down
Binary file added python/tests/data/MAPAXES.EGRID
View file
Binary file not shown.
Loading

0 comments on commit 3084f81

Please sign in to comment.