Skip to content

Commit

Permalink
Merge pull request #4257 from daavid00/mapaxes
Browse files Browse the repository at this point in the history
Handle MAPAXES and MAPUNITS
  • Loading branch information
bska authored Dec 11, 2024
2 parents 268f5b8 + 6dba4cb commit 9768aa7
Show file tree
Hide file tree
Showing 6 changed files with 173 additions and 66 deletions.
62 changes: 49 additions & 13 deletions opm/io/eclipse/EGrid.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -35,7 +35,7 @@

namespace Opm { namespace EclIO {

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

EGrid::EGrid(const std::string& filename, const std::string& grid_name)
: EclFile(filename), inputFileName { filename }, m_grid_name {grid_name}
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;
double length_factor = 1.0;

int hostnum_index = -1;

Expand All @@ -81,10 +83,25 @@ EGrid::EGrid(const std::string& filename, const std::string& grid_name)
if (array_name[n] == "MAPUNITS"){
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"){
const auto& mapAx = this->get<float>(n);
std::transform(mapAx.begin(), mapAx.end(), this->m_mapaxes.begin(),
[length_factor](const float elm) { return elm * length_factor; });
mapaxes_init();
m_mapaxes_loaded = true;
}

if (lgrname == grid_name) {
if (array_name[n] == "GRIDHEAD") {
Expand Down Expand Up @@ -322,11 +339,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 +388,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 +403,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 +427,10 @@ void EGrid::getCellCorners(const std::array<int, 3>& ijk,
}



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
20 changes: 14 additions & 6 deletions opm/io/eclipse/EGrid.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -42,8 +42,8 @@ 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);

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 +53,33 @@ class EGrid : public EclFile

void load_grid_data();
void load_nnc_data();
bool with_mapaxes() const { return m_mapaxes_loaded; }
void mapaxes_transform(double& x, double& y) const;
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;
bool m_mapaxes_loaded;
std::array<double, 4> origin;
std::array<double, 2> unit_x;
std::array<double, 2> unit_y;

std::array<int, 3> nijk;
std::array<int, 3> host_nijk;
Expand Down Expand Up @@ -107,8 +113,10 @@ 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);

void mapaxes_init();

};

}} // namespace Opm::EclIO
Expand Down
43 changes: 36 additions & 7 deletions python/cxx/eclipse_io.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -156,7 +156,7 @@ class EclOutputBind {

public:

EclOutputBind(const std::string& filename,const bool formatted, const bool append)
EclOutputBind(const std::string& filename, const bool formatted, const bool append)
{
if (append == true)
m_output = std::make_unique<Opm::EclIO::EclOutput>(filename, formatted, std::ios::app);
Expand Down Expand Up @@ -292,26 +292,51 @@ npArray get_erst_vector(Opm::EclIO::ERst * file_ptr, const std::string& key, siz
return get_erst_by_index(file_ptr, array_index, rstep);
}


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};

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_ijk_mapaxes(Opm::EclIO::EGrid * file_ptr, int i, int j, int k, bool mapaxes)
{
auto xyz = get_xyz_from_ijk(file_ptr, i, j, k);

if (file_ptr->with_mapaxes() && mapaxes){
for (int n = 0; n < 8; n++)
file_ptr->mapaxes_transform(std::get<0>(xyz)[n], std::get<1>(xyz)[n]);
}

return xyz;
}

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]);
}

std::tuple<std::array<double,8>, std::array<double,8>, std::array<double,8>>
get_xyz_from_active_index_mapaxes(Opm::EclIO::EGrid * file_ptr, int actIndex, bool mapaxes)
{
auto xyz = get_xyz_from_active_index(file_ptr, actIndex);

if (file_ptr->with_mapaxes() && mapaxes){
for (int n = 0; n < 8; n++)
file_ptr->mapaxes_transform(std::get<0>(xyz)[n], std::get<1>(xyz)[n]);
}

return xyz;
}

py::array get_cellvolumes_mask(Opm::EclIO::EGrid * file_ptr, std::vector<int> mask)
Expand Down Expand Up @@ -463,15 +488,19 @@ 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 std::string &>(), py::arg("filename"),
py::arg("grid_name") = "global")
.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_ijk", &get_xyz_from_ijk_mapaxes)
.def("xyz_from_active_index", &get_xyz_from_active_index)
.def("xyz_from_active_index", &get_xyz_from_active_index_mapaxes)
.def("cellvolumes", &get_cellvolumes)
.def("cellvolumes", &get_cellvolumes_mask);

Expand Down
9 changes: 8 additions & 1 deletion python/opm_embedded/__init__.pyi
Original file line number Diff line number Diff line change
Expand Up @@ -2480,17 +2480,24 @@ class Dimension:
def scaling(self) -> float: ...

class EGrid:
def __init__(self, arg0: str) -> None: ...
def __init__(self, filename: str, grid_name: str) -> 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]]: ...
@overload
def xyz_from_active_index(self, arg0: int) -> Tuple[List[float[8]], List[float[8]], List[float[8]]]: ...
@overload
def xyz_from_active_index(self, arg0: int, apply_mapaxes: bool) -> Tuple[List[float[8]], List[float[8]], List[float[8]]]: ...
@overload
def xyz_from_ijk(self, arg0: int, arg1: int, arg2: int) -> Tuple[List[float[8]], List[float[8]], List[float[8]]]: ...
@overload
def xyz_from_ijk(self, arg0: int, arg1: int, arg2: int, apply_mapaxes: bool) -> Tuple[List[float[8]], List[float[8]], List[float[8]]]: ...
@property
def active_cells(self) -> int: ...
@property
Expand Down
Binary file added python/tests/data/MAPAXES.EGRID
Binary file not shown.
Loading

0 comments on commit 9768aa7

Please sign in to comment.