From 0af8d2fd35c2b4b82db46ed73af5251539aeb23e Mon Sep 17 00:00:00 2001 From: mferrera Date: Sat, 13 Jan 2024 13:02:16 +0100 Subject: [PATCH] CLN: Create and return bulkvol array from C++ --- .gitignore | 1 + src/lib/CMakeLists.txt | 9 ++++---- src/lib/include/xtgeo/xtgeo.h | 1 - src/lib/src/grid3d/cellvol.cpp | 40 +++++++++++++++++----------------- src/xtgeo/grid3d/_grid_etc1.py | 7 ++---- tests/test_grid3d/test_grid.py | 12 +++++++--- 6 files changed, 37 insertions(+), 33 deletions(-) diff --git a/.gitignore b/.gitignore index 07a7287db..39b54e4c8 100644 --- a/.gitignore +++ b/.gitignore @@ -83,3 +83,4 @@ pip-wheel-metadata/ .nfs* tmp/ .DS_Store +compile_commands.json diff --git a/src/lib/CMakeLists.txt b/src/lib/CMakeLists.txt index bf4a96294..5e245e5eb 100644 --- a/src/lib/CMakeLists.txt +++ b/src/lib/CMakeLists.txt @@ -6,6 +6,11 @@ find_package(SWIG 3.0.1 COMPONENTS REQUIRED) include(UseSWIG) find_package(pybind11 REQUIRED) +message(STATUS "Python executable : ${Python_EXECUTABLE}") +message(STATUS "Python include dirs : ${Python_INCLUDE_DIRS}") +message(STATUS "numpy include path : ${Python_NumPy_INCLUDE_DIRS}") +message(STATUS "pybind11 include path : ${pybind11_INCLUDE_DIRS}") + # TODO: replace globbing with unique list, as globbing is bad practice FILE(GLOB SOURCES ${SRC}/*.c) add_library(xtgeo STATIC ${SOURCES}) @@ -47,10 +52,6 @@ target_link_libraries(${SWIG_TARGET} Python::Module Python::NumPy) -message(STATUS "Python executable : ${Python_EXECUTABLE}") -message(STATUS "Python include dirs: ${Python_INCLUDE_DIRS}") -message(STATUS "numpy include path : ${Python_NumPy_INCLUDE_DIRS}") - # scikit-build-core docs recommend this if(WIN32) set_property(TARGET ${SWIG_TARGET} PROPERTY SUFFIX ".${Python_SOABI}.pyd") diff --git a/src/lib/include/xtgeo/xtgeo.h b/src/lib/include/xtgeo/xtgeo.h index f6f4f8530..a0bfae631 100644 --- a/src/lib/include/xtgeo/xtgeo.h +++ b/src/lib/include/xtgeo/xtgeo.h @@ -16,7 +16,6 @@ #define _GNU_SOURCE 1 #include -#include #include #define PI 3.14159265358979323846264338327950288 diff --git a/src/lib/src/grid3d/cellvol.cpp b/src/lib/src/grid3d/cellvol.cpp index 0ccd84e48..eb7f63d20 100644 --- a/src/lib/src/grid3d/cellvol.cpp +++ b/src/lib/src/grid3d/cellvol.cpp @@ -7,43 +7,43 @@ namespace py = pybind11; -void -grid3d_cellvol(size_t ncol, - size_t nrow, - size_t nlay, - py::array_t coordsv, - py::array_t zcornsv, - py::array_t actnumsv, - py::array_t cellvolsv, - int precision, - bool asmasked = false) +py::array_t +grid3d_cellvol(const size_t ncol, + const size_t nrow, + const size_t nlay, + const py::array_t coordsv, + const py::array_t zcornsv, + const py::array_t actnumsv, + const int precision, + const bool asmasked = false) { - py::buffer_info coordsv_buf = coordsv.request(); - py::buffer_info zcornsv_buf = zcornsv.request(); - py::buffer_info actnumsv_buf = actnumsv.request(); - py::buffer_info cellvolsv_buf = cellvolsv.request(); + auto coordsv_buf = coordsv.request(); + auto zcornsv_buf = zcornsv.request(); + auto actnumsv_buf = actnumsv.request(); double *coordsv_ = static_cast(coordsv_buf.ptr); float *zcornsv_ = static_cast(zcornsv_buf.ptr); int *actnumsv_ = static_cast(actnumsv_buf.ptr); - double *cellvolsv_ = static_cast(cellvolsv_buf.ptr); - double corners[24]{}; + pybind11::array_t cellvols({ ncol, nrow, nlay }); + double *cellvolsv_ = static_cast(cellvols.request().ptr); + double corners[24]{}; for (size_t i = 0; i < ncol; i++) { for (size_t j = 0; j < nrow; j++) { for (size_t k = 0; k < nlay; k++) { - size_t ic = i * nrow * nlay + j * nlay + k; - if (asmasked && actnumsv_[ic] == 0) { - cellvolsv_[ic] = UNDEF; + size_t idx = i * nrow * nlay + j * nlay + k; + if (asmasked && actnumsv_[idx] == 0) { + cellvolsv_[idx] = UNDEF; continue; } grdcp3d_corners(i, j, k, ncol, nrow, nlay, coordsv_, coordsv_buf.size, zcornsv_, zcornsv_buf.size, corners); - cellvolsv_[ic] = x_hexahedron_volume(corners, 24, precision); + cellvolsv_[idx] = x_hexahedron_volume(corners, 24, precision); } } } + return cellvols; } PYBIND11_MODULE(_internal, m) diff --git a/src/xtgeo/grid3d/_grid_etc1.py b/src/xtgeo/grid3d/_grid_etc1.py index bc30fc806..455f883b9 100644 --- a/src/xtgeo/grid3d/_grid_etc1.py +++ b/src/xtgeo/grid3d/_grid_etc1.py @@ -217,16 +217,13 @@ def get_bulk_volume( discrete=False, ) - bval = np.zeros(bulk.dimensions) - - xtgeo._internal.grid3d_cellvol( + bval = xtgeo._internal.grid3d_cellvol( self._ncol, self._nrow, self._nlay, self._coordsv.ravel(), self._zcornsv.ravel(), self._actnumsv.ravel(), - bval, precision, asmasked, ) @@ -234,7 +231,7 @@ def get_bulk_volume( if asmasked: bval = np.ma.masked_greater(bval, UNDEF_LIMIT) - bulk.values = bval # type: ignore + bulk.values = bval return bulk diff --git a/tests/test_grid3d/test_grid.py b/tests/test_grid3d/test_grid.py index ed1e50cf3..728edfc9c 100644 --- a/tests/test_grid3d/test_grid.py +++ b/tests/test_grid3d/test_grid.py @@ -513,10 +513,16 @@ def test_bulkvol(): """Test cell bulk volume calculation.""" grd = xtgeo.grid_from_file(GRIDQC1) cellvol_rms = xtgeo.gridproperty_from_file(GRIDQC1_CELLVOL) + bulkvol = grd.get_bulk_volume() - bulk = grd.get_bulk_volume() - logger.info("Sum this: %s", bulk.values.sum()) - logger.info("Sum RMS: %s", cellvol_rms.values.sum()) + rms_sum = np.sum(cellvol_rms.values) + bulkvol_sum = np.sum(bulkvol.values) + + print(f"RMS sum: {rms_sum}") + print(f"bulkvol sum: {bulkvol_sum}") + + assert grd.dimensions == bulkvol.dimensions + assert np.allclose(cellvol_rms.values, bulkvol.values) @pytest.mark.benchmark(group="bulkvol")