CLN: Create and return bulkvol array from C++

.gitignore

@@ -83,3 +83,4 @@ pip-wheel-metadata/
 .nfs*
 tmp/
 .DS_Store
+compile_commands.json

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")

src/lib/include/xtgeo/xtgeo.h

@@ -16,7 +16,6 @@
 #define _GNU_SOURCE 1
 
 #include <stdbool.h>
-#include <stdint.h>
 #include <stdio.h>
 
 #define PI 3.14159265358979323846264338327950288

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<double> coordsv,
-               py::array_t<float> zcornsv,
-               py::array_t<int> actnumsv,
-               py::array_t<double> cellvolsv,
-               int precision,
-               bool asmasked = false)
+py::array_t<double>
+grid3d_cellvol(const size_t ncol,
+               const size_t nrow,
+               const size_t nlay,
+               const py::array_t<double> coordsv,
+               const py::array_t<float> zcornsv,
+               const py::array_t<int> 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<double *>(coordsv_buf.ptr);
     float *zcornsv_ = static_cast<float *>(zcornsv_buf.ptr);
     int *actnumsv_ = static_cast<int *>(actnumsv_buf.ptr);
-    double *cellvolsv_ = static_cast<double *>(cellvolsv_buf.ptr);
 
-    double corners[24]{};
+    pybind11::array_t<double> cellvols({ ncol, nrow, nlay });
+    double *cellvolsv_ = static_cast<double *>(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)

src/xtgeo/grid3d/_grid_etc1.py

@@ -217,24 +217,21 @@ 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,
     )
 
     if asmasked:
         bval = np.ma.masked_greater(bval, UNDEF_LIMIT)
 
-    bulk.values = bval  # type: ignore
+    bulk.values = bval
 
     return bulk
 

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")