diff --git a/lib/CMakeLists.txt b/lib/CMakeLists.txt index ea45b47b1..a207cb6e9 100644 --- a/lib/CMakeLists.txt +++ b/lib/CMakeLists.txt @@ -325,6 +325,34 @@ set_target_properties( _well_info PROPERTIES LIBRARY_OUTPUT_DIRECTORY ${PROJECT_SOURCE_DIR}/python/resdata/well) +# ----------------------------------------------------------------- +# Pybind11 extension module: resdata.gravimetry._subsidence +# ----------------------------------------------------------------- +pybind11_add_module(_subsidence resdata/rd_subsidence_pybind.cpp + resdata/cwrap_pybind.cpp) +target_link_libraries(_subsidence PRIVATE resdata fmt::fmt) +target_include_directories( + _subsidence + PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/include + ${CMAKE_CURRENT_BINARY_DIR}/include + ${CMAKE_CURRENT_SOURCE_DIR}/private-include) +if(SKBUILD) + set(_subsidence_install_dir "python/resdata/gravimetry") + if(APPLE) + set_target_properties(_subsidence PROPERTIES INSTALL_RPATH + "@loader_path/../.libs") + else() + set_target_properties(_subsidence PROPERTIES INSTALL_RPATH + "$ORIGIN/../.libs") + endif() +else() + set(_subsidence_install_dir "${CMAKE_INSTALL_LIBDIR}") +endif() +install(TARGETS _subsidence LIBRARY DESTINATION ${_subsidence_install_dir}) +set_target_properties( + _subsidence PROPERTIES LIBRARY_OUTPUT_DIRECTORY + ${PROJECT_SOURCE_DIR}/python/resdata/gravimetry) + # ----------------------------------------------------------------- # Pybind11 extension module: resdata.summary._rd_sum # ----------------------------------------------------------------- diff --git a/lib/include/resdata/rd_subsidence.hpp b/lib/include/resdata/rd_subsidence.hpp index 596c5c853..88d4c9feb 100644 --- a/lib/include/resdata/rd_subsidence.hpp +++ b/lib/include/resdata/rd_subsidence.hpp @@ -1,15 +1,13 @@ -#ifndef ERT_RD_SUBSIDENCE_H -#define ERT_RD_SUBSIDENCE_H +#pragma once + +#include +#include #include #include #include #include -#ifdef __cplusplus -extern "C" { -#endif - typedef struct rd_subsidence_struct rd_subsidence_type; typedef struct rd_subsidence_survey_struct rd_subsidence_survey_type; @@ -18,32 +16,28 @@ rd_subsidence_type *rd_subsidence_alloc(rd_grid_type *rd_grid, const rd_file_type *init_file); rd_subsidence_survey_type * rd_subsidence_add_survey_PRESSURE(rd_subsidence_type *subsidence, - const char *name, + const std::string &name, const rd_file_view_type *restart_view); bool rd_subsidence_has_survey(const rd_subsidence_type *subsidence, - const char *name); + const std::string &name); double rd_subsidence_eval(const rd_subsidence_type *subsidence, - const char *base, const char *monitor, + const std::string &base, + const std::optional &monitor, rd_region_type *region, double utm_x, double utm_y, double depth, double compressibility, double poisson_ratio); double rd_subsidence_eval_geertsma(const rd_subsidence_type *subsidence, - const char *base, const char *monitor, + const std::string &base, + const std::optional &monitor, rd_region_type *region, double utm_x, double utm_y, double depth, double youngs_modulus, double poisson_ratio, double seabed); -double rd_subsidence_eval_geertsma_rporv(const rd_subsidence_type *subsidence, - const char *base, const char *monitor, - rd_region_type *region, double utm_x, - double utm_y, double depth, - double youngs_modulus, - double poisson_ratio, double seabed); - -#ifdef __cplusplus -} -#endif -#endif +double rd_subsidence_eval_geertsma_rporv( + const rd_subsidence_type *subsidence, const std::string &base, + const std::optional &monitor, rd_region_type *region, + double utm_x, double utm_y, double depth, double youngs_modulus, + double poisson_ratio, double seabed); diff --git a/lib/resdata/cwrap_pybind.cpp b/lib/resdata/cwrap_pybind.cpp index 9d07311e4..0153e99a5 100644 --- a/lib/resdata/cwrap_pybind.cpp +++ b/lib/resdata/cwrap_pybind.cpp @@ -12,6 +12,8 @@ #include #include #include +#include +#include #include #include @@ -233,3 +235,35 @@ template <> rd_file_view_type *from_cwrap(py::handle obj) { static_cast(py::repr(obj))); return cast_cwrap(obj); } +py::object ResdataSubsidence() { + static py::object cls; + if (!cls) { + cls = + py::module_::import("resdata.gravimetry").attr("ResdataSubsidence"); + } + return cls; +} + +template <> rd_subsidence_type *from_cwrap(py::handle obj) { + if (!py::isinstance(obj, ResdataSubsidence())) + throw py::type_error("Expected ResdataSubsidence, got " + + static_cast(py::repr(obj))); + + return cast_cwrap(obj); +} + +py::object ResdataRegion() { + static py::object cls; + if (!cls) { + cls = py::module_::import("resdata.grid").attr("ResdataRegion"); + } + return cls; +} + +template <> rd_region_type *from_cwrap(py::handle obj) { + if (!py::isinstance(obj, ResdataRegion())) + throw py::type_error("Expected ResdataRegion, got " + + static_cast(py::repr(obj))); + + return cast_cwrap(obj); +} diff --git a/lib/resdata/rd_subsidence.cpp b/lib/resdata/rd_subsidence.cpp index 3319e2a05..4961acecd 100644 --- a/lib/resdata/rd_subsidence.cpp +++ b/lib/resdata/rd_subsidence.cpp @@ -8,6 +8,7 @@ #include #include #include +#include #include #include @@ -59,7 +60,8 @@ struct rd_subsidence_survey_struct { std::optional> dynamic_porevolume; /* Porevolume in each grid cell at survey time */ - rd_subsidence_survey_struct(const rd_subsidence_type &sub, const char *name) + rd_subsidence_survey_struct(const rd_subsidence_type &sub, + const std::string name) : grid_cache(sub.grid_cache.get()), aquifer_cell(sub.aquifer_cell.get()), name(name), porv(sub.grid_cache->size()), pressure(sub.grid_cache->size()) {} @@ -68,7 +70,7 @@ struct rd_subsidence_survey_struct { static std::unique_ptr rd_subsidence_survey_alloc_PRESSURE(rd_subsidence_type *rd_subsidence, const rd_file_view_type *restart_view, - const char *name) { + const std::string &name) { auto survey( std::make_unique(*rd_subsidence, name)); @@ -169,17 +171,17 @@ static double rd_subsidence_survey_eval_geertsma_rporv( std::vector weight(size); if (!base_survey->dynamic_porevolume.has_value()) { - util_abort( - "%s: Keyword RPORV not defined in .UNRST file for %s. Please add " + throw std::invalid_argument(fmt::format( + "Keyword RPORV not defined in .UNRST file for {}. Please add " "RPORV keyword to output in RPTRST clause in .DATA file.\n", - __func__, base_survey->name.c_str()); + base_survey->name)); } if (monitor_survey && !monitor_survey->dynamic_porevolume.has_value()) { - util_abort( - "%s: Keyword RPORV not defined in .UNRST file for %s. Please add " - "RPORV keyword to output in RPTRST clause in .DATA file.\n", - __func__, monitor_survey->name.c_str()); + throw std::invalid_argument(fmt::format( + "Keyword RPORV not defined in .UNRST file for {}. Please add " + "RPORV keyword to output in RPTRST clause in .DATA file.", + monitor_survey->name)); } for (size_t index = 0; index < weight.size(); ++index) { @@ -216,10 +218,8 @@ rd_subsidence_type *rd_subsidence_alloc(rd_grid_type *rd_grid, rd_subsidence_survey_type * rd_subsidence_add_survey_PRESSURE(rd_subsidence_type *subsidence, - const char *name, + const std::string &name, const rd_file_view_type *restart_view) { - if (name == nullptr) - throw std::invalid_argument("Name cannot be NULL"); auto survey = rd_subsidence_survey_alloc_PRESSURE(subsidence, restart_view, name); rd_subsidence_survey_type *ret = survey.get(); @@ -228,28 +228,28 @@ rd_subsidence_add_survey_PRESSURE(rd_subsidence_type *subsidence, } bool rd_subsidence_has_survey(const rd_subsidence_type *subsidence, - const char *name) { - if (name == nullptr) - return false; + const std::string &name) { return subsidence->surveys.count(name) > 0; } double rd_subsidence_eval(const rd_subsidence_type *subsidence, - const char *base, const char *monitor, + const std::string &base, + const std::optional &monitor, rd_region_type *region, double utm_x, double utm_y, double depth, double compressibility, double poisson_ratio) { rd_subsidence_survey_type *base_survey = subsidence->surveys.at(base).get(); rd_subsidence_survey_type *monitor_survey = nullptr; if (monitor) - monitor_survey = subsidence->surveys.at(monitor).get(); + monitor_survey = subsidence->surveys.at(*monitor).get(); return rd_subsidence_survey_eval(base_survey, monitor_survey, region, utm_x, utm_y, depth, compressibility, poisson_ratio); } double rd_subsidence_eval_geertsma(const rd_subsidence_type *subsidence, - const char *base, const char *monitor, + const std::string &base, + const std::optional &monitor, rd_region_type *region, double utm_x, double utm_y, double depth, double youngs_modulus, double poisson_ratio, @@ -257,22 +257,21 @@ double rd_subsidence_eval_geertsma(const rd_subsidence_type *subsidence, rd_subsidence_survey_type *base_survey = subsidence->surveys.at(base).get(); rd_subsidence_survey_type *monitor_survey = nullptr; if (monitor) - monitor_survey = subsidence->surveys.at(monitor).get(); + monitor_survey = subsidence->surveys.at(*monitor).get(); return rd_subsidence_survey_eval_geertsma( base_survey, monitor_survey, region, utm_x, utm_y, depth, youngs_modulus, poisson_ratio, seabed); } -double rd_subsidence_eval_geertsma_rporv(const rd_subsidence_type *subsidence, - const char *base, const char *monitor, - rd_region_type *region, double utm_x, - double utm_y, double depth, - double youngs_modulus, - double poisson_ratio, double seabed) { +double rd_subsidence_eval_geertsma_rporv( + const rd_subsidence_type *subsidence, const std::string &base, + const std::optional &monitor, rd_region_type *region, + double utm_x, double utm_y, double depth, double youngs_modulus, + double poisson_ratio, double seabed) { rd_subsidence_survey_type *base_survey = subsidence->surveys.at(base).get(); rd_subsidence_survey_type *monitor_survey = nullptr; if (monitor) - monitor_survey = subsidence->surveys.at(monitor).get(); + monitor_survey = subsidence->surveys.at(*monitor).get(); return rd_subsidence_survey_eval_geertsma_rporv( base_survey, monitor_survey, region, utm_x, utm_y, depth, youngs_modulus, poisson_ratio, seabed); diff --git a/lib/resdata/rd_subsidence_pybind.cpp b/lib/resdata/rd_subsidence_pybind.cpp new file mode 100644 index 000000000..5f21c6204 --- /dev/null +++ b/lib/resdata/rd_subsidence_pybind.cpp @@ -0,0 +1,93 @@ +#include +#include +#include + +#include +#include + +#include +#include +#include +#include +#include + +#include + +namespace py = pybind11; + +PYBIND11_MODULE(_subsidence, m) { + register_exceptions(m); + m.doc() = "pybind11 bindings for rd_subsidence.cpp"; + + m.def("_alloc", [](py::handle grid, py::handle init_file) -> py::object { + auto *subsidence = + rd_subsidence_alloc(from_cwrap(grid), + from_cwrap(init_file)); + + if (subsidence == nullptr) + return py::none(); + + return py::cast(reinterpret_cast(subsidence)); + }); + + m.def("_free", [](py::handle self) { + rd_subsidence_free(from_cwrap(self)); + }); + + m.def("_add_survey_PRESSURE", [](py::handle self, std::string survey_name, + py::handle restart_file_view) { + rd_subsidence_add_survey_PRESSURE( + from_cwrap(self), survey_name, + from_cwrap(restart_file_view)); + }); + + m.def("_eval", [](py::handle self, std::string base, + std::optional monitor, py::object region, + double utm_x, double utm_y, double depth, + double compressibility, double poisson_ratio) { + rd_region_type *region_ptr = nullptr; + if (!region.is_none()) { + region_ptr = from_cwrap(region); + } + + return rd_subsidence_eval(from_cwrap(self), base, + monitor, region_ptr, utm_x, utm_y, depth, + compressibility, poisson_ratio); + }); + + m.def("_eval_geertsma", [](py::handle self, std::string base, + std::optional monitor, + py::object region, double utm_x, double utm_y, + double depth, double youngs_modulus, + double poisson_ratio, double seabed) { + rd_region_type *region_ptr = nullptr; + if (!region.is_none()) { + region_ptr = from_cwrap(region); + } + + return rd_subsidence_eval_geertsma( + from_cwrap(self), base, monitor, region_ptr, + utm_x, utm_y, depth, youngs_modulus, poisson_ratio, seabed); + }); + + m.def("_eval_geertsma_rporv", + [](py::handle self, std::string base, + std::optional monitor, py::object region, + double utm_x, double utm_y, double depth, double youngs_modulus, + double poisson_ratio, double seabed) { + rd_region_type *region_ptr = nullptr; + if (!region.is_none()) { + region_ptr = from_cwrap(region); + } + + return rd_subsidence_eval_geertsma_rporv( + from_cwrap(self), base, monitor, + region_ptr, utm_x, utm_y, depth, youngs_modulus, + poisson_ratio, seabed); + }); + + m.def("_has_survey", [](py::handle self, std::string survey_name) { + return rd_subsidence_has_survey(from_cwrap(self), + survey_name); + }); +} diff --git a/lib/tests/test_geertsma.cpp b/lib/tests/test_geertsma.cpp index 2b38c8b35..99e8fb35c 100644 --- a/lib/tests/test_geertsma.cpp +++ b/lib/tests/test_geertsma.cpp @@ -103,7 +103,7 @@ TEST_CASE_METHOD(Tmpdir, "Geertsma kernel single cell") { SECTION("receiver at surface") { double dz = rd_subsidence_eval_geertsma( - subsidence.get(), "S1", nullptr, nullptr, 1000, 1000, 0, + subsidence.get(), "S1", std::nullopt, nullptr, 1000, 1000, 0, youngs_modulus, poisson_ratio, seabed); REQUIRE_THAT(dz, WithinRel(3.944214576168326e-09, 1e-9)); } @@ -111,7 +111,7 @@ TEST_CASE_METHOD(Tmpdir, "Geertsma kernel single cell") { SECTION("receiver below surface") { double depth = topres - seabed - above; double dz = rd_subsidence_eval_geertsma( - subsidence.get(), "S1", nullptr, nullptr, 1000, 1000, depth, + subsidence.get(), "S1", std::nullopt, nullptr, 1000, 1000, depth, youngs_modulus, poisson_ratio, seabed); REQUIRE_THAT(dz, WithinRel(5.8160298201497136e-08, 1e-9)); @@ -154,13 +154,13 @@ TEST_CASE_METHOD(Tmpdir, "Geertsma kernel two source points two vintages") { const double seabed = 0.0; double dz1 = rd_subsidence_eval_geertsma( - subsidence.get(), "S1", nullptr, nullptr, 1000, 1000, 0, youngs_modulus, - poisson_ratio, seabed); + subsidence.get(), "S1", std::nullopt, nullptr, 1000, 1000, 0, + youngs_modulus, poisson_ratio, seabed); REQUIRE_THAT(dz1, WithinRel(8.65322541521704e-07, 1e-9)); double dz2 = rd_subsidence_eval_geertsma( - subsidence.get(), "S2", nullptr, nullptr, 1000, 1000, 0, youngs_modulus, - poisson_ratio, seabed); + subsidence.get(), "S2", std::nullopt, nullptr, 1000, 1000, 0, + youngs_modulus, poisson_ratio, seabed); REQUIRE_THAT(dz2, WithinRel(2.275556615015282e-06, 1e-9)); double dz = rd_subsidence_eval_geertsma( @@ -197,7 +197,7 @@ TEST_CASE_METHOD(Tmpdir, "Geertsma kernel with seabed") { double depth = topres - seabed - above; double dz = rd_subsidence_eval_geertsma( - subsidence.get(), "S1", nullptr, nullptr, 1000, 1000, depth, + subsidence.get(), "S1", std::nullopt, nullptr, 1000, 1000, depth, youngs_modulus, poisson_ratio, seabed); REQUIRE_THAT(dz, WithinRel(5.819790154474284e-08, 1e-9)); } @@ -235,11 +235,11 @@ TEST_CASE_METHOD(Tmpdir, const double seabed = 0.0; double dz1 = rd_subsidence_eval_geertsma_rporv( - subsidence.get(), "S1", nullptr, nullptr, 1000, 1000, 0, youngs_modulus, - poisson_ratio, seabed); + subsidence.get(), "S1", std::nullopt, nullptr, 1000, 1000, 0, + youngs_modulus, poisson_ratio, seabed); double dz2 = rd_subsidence_eval_geertsma_rporv( - subsidence.get(), "S2", nullptr, nullptr, 1000, 1000, 0, youngs_modulus, - poisson_ratio, seabed); + subsidence.get(), "S2", std::nullopt, nullptr, 1000, 1000, 0, + youngs_modulus, poisson_ratio, seabed); double dz = rd_subsidence_eval_geertsma_rporv( subsidence.get(), "S1", "S2", nullptr, 1000, 1000, 0, youngs_modulus, poisson_ratio, seabed); @@ -275,13 +275,13 @@ TEST_CASE_METHOD(Tmpdir, "Subsidence survey validation") { SECTION("unknown base survey throws") { REQUIRE(!rd_subsidence_has_survey(subsidence.get(), "dummy_survey")); REQUIRE_THROWS_AS(rd_subsidence_eval_geertsma_rporv( - subsidence.get(), "dummy_survey", nullptr, + subsidence.get(), "dummy_survey", std::nullopt, nullptr, 1000, 1000, 0, youngs_modulus, poisson_ratio, seabed), std::out_of_range); REQUIRE_THROWS_AS(rd_subsidence_eval(subsidence.get(), "dummy_survey", - nullptr, nullptr, 1000, 1000, 0, - compressibility, poisson_ratio), + std::nullopt, nullptr, 1000, 1000, + 0, compressibility, poisson_ratio), std::out_of_range); } diff --git a/python/resdata/gravimetry/rd_subsidence.py b/python/resdata/gravimetry/rd_subsidence.py index 3e9621102..fa3dc2ae2 100644 --- a/python/resdata/gravimetry/rd_subsidence.py +++ b/python/resdata/gravimetry/rd_subsidence.py @@ -7,7 +7,7 @@ from cwrap import BaseCClass -from resdata import ResdataPrototype +import resdata.gravimetry._subsidence as _subsidence from resdata.grid import Grid, ResdataRegion from resdata.resfile import ResdataFile, ResdataFileView from resdata.util.util import monkey_the_camel @@ -30,35 +30,18 @@ class is focused on the ECLIPSE side of things, and does not have """ TYPE_NAME = "rd_subsidence" - _alloc = ResdataPrototype( - "void* rd_subsidence_alloc( rd_grid , rd_file )", bind=False - ) - _free = ResdataPrototype("void rd_subsidence_free( rd_subsidence )") - _add_survey_PRESSURE = ResdataPrototype( - "void* rd_subsidence_add_survey_PRESSURE( rd_subsidence , char* , rd_file_view )" - ) - _eval = ResdataPrototype( - "double rd_subsidence_eval( rd_subsidence , char* , char* , rd_region , double , double , double, double, double)" - ) - _eval_geertsma = ResdataPrototype( - "double rd_subsidence_eval_geertsma( rd_subsidence , char* , char* , rd_region , double , double , double, double, double, double)" - ) - _eval_geertsma_rporv = ResdataPrototype( - "double rd_subsidence_eval_geertsma_rporv( rd_subsidence , char* , char* , rd_region , double , double , double, double, double, double)" - ) - _has_survey = ResdataPrototype( - "bool rd_subsidence_has_survey( rd_subsidence , char*)" - ) def __init__(self, grid: Grid, init_file: ResdataFile): self.init_file = init_file # Inhibit premature garbage collection of init_file - c_ptr = self._alloc(grid, init_file) + c_ptr = _subsidence._alloc(grid, init_file) + if c_ptr is None: + raise ValueError("Unable to construct ResdataSubsidence") super().__init__(c_ptr) def __contains__(self, survey_name: str) -> bool: if survey_name is None: raise TypeError("survey_name must not be None") - return self._has_survey(survey_name) + return _subsidence._has_survey(self, survey_name) def add_survey_PRESSURE( self, survey_name: str, restart_file: ResdataFileView @@ -86,7 +69,7 @@ def add_survey_PRESSURE( """ if survey_name is None: raise TypeError("survey_name must not be None") - self._add_survey_PRESSURE(survey_name, restart_file) + _subsidence._add_survey_PRESSURE(self, survey_name, restart_file) def eval_geertsma( self, @@ -98,9 +81,6 @@ def eval_geertsma( seabed, region=None, ): - if base_survey is None: - raise TypeError("base_survey must not be None") - if base_survey not in self: raise KeyError("No such survey: %s" % base_survey) @@ -108,7 +88,8 @@ def eval_geertsma( if monitor_survey not in self: raise KeyError("No such survey: %s" % monitor_survey) - return self._eval_geertsma( + return _subsidence._eval_geertsma( + self, base_survey, monitor_survey, region, @@ -130,9 +111,6 @@ def eval_geertsma_rporv( seabed, region=None, ): - if base_survey is None: - raise TypeError("base_survey must not be None") - if base_survey not in self: raise KeyError("No such survey: %s" % base_survey) @@ -140,7 +118,8 @@ def eval_geertsma_rporv( if monitor_survey not in self: raise KeyError("No such survey: %s" % monitor_survey) - return self._eval_geertsma_rporv( + return _subsidence._eval_geertsma_rporv( + self, base_survey, monitor_survey, region, @@ -183,16 +162,14 @@ def eval( The argument @compressibility is the total reservoir compressibility. """ - if base_survey is None: - raise TypeError("base_survey must not be None") - if base_survey not in self: raise KeyError("No such survey: %s" % base_survey) if monitor_survey is not None and monitor_survey not in self: raise KeyError("No such survey: %s" % monitor_survey) - return self._eval( + return _subsidence._eval( + self, base_survey, monitor_survey, region, @@ -204,7 +181,7 @@ def eval( ) def free(self): - self._free() + _subsidence._free(self) monkey_the_camel(ResdataSubsidence, "evalGeertsma", ResdataSubsidence.eval_geertsma)