swfunction Add new method for height above FFL and refactoring (#300)

Breaking change: the `compute` method now return dictionary of grid property instead of grid property. The key for dictionary are: SW, HTOP, HBOT, HCENTER

Author: alifbe

docs/properties.rst

@@ -69,7 +69,7 @@ factors are constant numbers:
             invert=True,
             method=METHOD,
         )
-        sw = sw_func.compute(CALC)
+        sw = sw_func.compute(CALC)["SW"]
         sw.values[poro.values < 0.05] = 1.0  # 100% water when low porosity
         sw.to_roxar(project, GRIDNAME, SW_RESULT)
 

src/fmu/tools/properties/swfunction.py

@@ -25,9 +25,9 @@
 logger = logging.getLogger(__name__)
 
 ALLOWED_METHODS = [
-    "cell_center",
     "cell_center_above_ffl",
     "cell_corners_above_ffl",
+    "truncated_cell_corners_above_ffl",
 ]
 
 
@@ -67,8 +67,6 @@ class SwFunction:
             is called several time for the same grid.
         hbot: Optional. May speed up computation if provided, in case the function
             is called several time for the same grid.
-        debug: If True, several "check parameters" will we given, either as grid
-            properties (if working in RMS) or as ROFF files (if working outside RMS).
         tag: Optional string to identify debug parameters.
 
 
@@ -92,13 +90,12 @@ class SwFunction:
     invert: bool = False  # if the SwJ function is on "reverse" form
 
     # if None they will be computed here; otherwise they can be given explicitly:
-    hcenter: xtgeo.GridProperty = None
-    htop: xtgeo.GridProperty = None
-    hbot: xtgeo.GridProperty = None
+    hcenter: xtgeo.GridProperty | None = None
+    htop: xtgeo.GridProperty | None = None
+    hbot: xtgeo.GridProperty | None = None
 
     # debug flag for additional parameters. Given in RMS if project is given; otherwise
     # as files on your working folder
-    debug: bool = False
     tag: str = ""  # identification tag to add to debug params
 
     # derived and internal
@@ -118,12 +115,6 @@ def __post_init__(self) -> None:
         else:
             self._compute_htop_hbot()
 
-        if self.debug:
-            if self.project:
-                self.grid.to_roxar(self.project, "DEBUG_" + self.gridname)
-            else:
-                self.grid.to_file("debug_grid.roff")
-
     def _process_input(self) -> None:
         """Work with a, b, x, ie. inversing and convert from float."""
         logger.info("Process a, b, etc")
@@ -181,20 +172,6 @@ def _compute_htop_hbot(self) -> None:
         self.hcenter = hmid
         logger.info("Use method %s", self.method)
 
-        if self.debug:
-            htop.name = "TOP_SW" + self.tag
-            hbot.name = "BOT_SW" + self.tag
-            hmid.name = "CENTER_SW" + self.tag
-            if self.project:
-                logger.debug("TAG is ", self.tag)
-                htop.to_roxar(self.project, self.gridname, "HTOP" + self.tag)
-                hbot.to_roxar(self.project, self.gridname, "HBOT" + self.tag)
-                hmid.to_roxar(self.project, self.gridname, "HCENTER" + self.tag)
-            else:
-                htop.to_file(f"debug_htop{self.tag}.roff")
-                hbot.to_file(f"debug_hbot{self.tag}.roff")
-                hmid.to_file(f"debug_hcenter{self.tag}.roff")
-
     def _sw_function_direct(self) -> None:
         """Use function on form Sw = A*(M + X*h)^B; generic function!"""
         assert isinstance(self.a, xtgeo.GridProperty)  # mypy
@@ -203,6 +180,7 @@ def _sw_function_direct(self) -> None:
         assert isinstance(self.m, xtgeo.GridProperty)  # mypy
         assert isinstance(self.swira, xtgeo.GridProperty)  # mypy
         assert isinstance(self.swmax, xtgeo.GridProperty)  # mypy
+        assert isinstance(self.hcenter, xtgeo.GridProperty)
 
         # the direct function is mostly used to compare with integrated approach, as QC
         height = self.hcenter.values
@@ -237,20 +215,16 @@ def _sw_function_integrate_w_mterm(self) -> None:
         assert isinstance(self.m, xtgeo.GridProperty)  # mypy
         assert isinstance(self.swira, xtgeo.GridProperty)  # mypy
         assert isinstance(self.swmax, xtgeo.GridProperty)  # mypy
+        assert isinstance(self.htop, xtgeo.GridProperty)  # mypy
+        assert isinstance(self.hbot, xtgeo.GridProperty)  # mypy
+
         ht = (
             ((1.0 / self.a.values) ** (1.0 / self.b.values)) - self.m.values
         ) / self.x.values  # threshold height
 
         h2 = self.htop.values.copy()  # h_top or H2 in integration
         h1 = self.hbot.values.copy()  # h_bot or H1 in integration
 
-        if self.debug:
-            tmp = xtgeo.GridProperty(self.grid, values=ht, name="HT" + self.tag)
-            if self.project:
-                tmp.to_roxar(self.project, self.gridname, "THRESHOLD HEIGHT" + self.tag)
-            else:
-                tmp.to_file(f"debug_ht{self.tag}.roff")
-
         water = h2 * 0.0
         water = np.ma.where(h2 < ht, 1.0, water)
         water = np.ma.where(h2 <= 0.0, 1.0, water)  # may occur for negative ht
@@ -317,7 +291,9 @@ def _compute_direct(self) -> None:
         if self._sw.values.min() < 0.0:
             raise RuntimeError(f"SW min out of range: {self._sw.values.min()}")
 
-    def compute(self, compute_method: str = "integrated") -> xtgeo.GridProperty:
+    def compute(
+        self, compute_method: str = "integrated"
+    ) -> dict[str, xtgeo.GridProperty]:
         """Common compute function for saturation, and returns the Sw property"""
         if compute_method == "integrated":
             self._compute_integrated()
@@ -331,4 +307,11 @@ def compute(self, compute_method: str = "integrated") -> xtgeo.GridProperty:
             "Bug: Grid mask and Sw mask are not equal, contact developer!"
         )
 
-        return self._sw
+        output_props = [self._sw, self.htop, self.hbot, self.hcenter]
+        output_props_name = ["SW", "HTOP", "HBOT", "HCENTER"]
+
+        for prop, prop_name in zip(output_props, output_props_name):
+            assert isinstance(prop, xtgeo.GridProperty)
+            prop.name = self.tag + prop_name
+
+        return dict(zip(output_props_name, output_props))

tests/properties/test_swfunction.py

@@ -18,6 +18,8 @@
         (0.3, -5, 12.5, True, "cell_center_above_ffl", 0.547251),
         (0.3, -5, 12.5, False, "cell_center_above_ffl", 0.549613),
         (0.3, -5, 12.5, False, "cell_corners_above_ffl", 0.549613),
+        (1, -2, 13, True, "truncated_cell_corners_above_ffl", 0.408742),
+        (1, -2, 13, False, "truncated_cell_corners_above_ffl", 0.412536),
     ],
 )
 def test_swj_simple(avalue, bvalue, ffl, direct, cellmethod, expected_mean):
@@ -38,7 +40,7 @@ def test_swj_simple(avalue, bvalue, ffl, direct, cellmethod, expected_mean):
         invert=True,
         method=cellmethod,
     )
-    sw = sw_obj.compute("direct" if direct else "integrated")
+    sw = sw_obj.compute("direct" if direct else "integrated")["SW"]
     assert sw.values.mean() == pytest.approx(expected_mean, rel=0.01)
 
 
@@ -63,7 +65,8 @@ def test_swj_simple_x_zero():
         invert=True,
         method="cell_center_above_ffl",
     )
-    sw = sw_obj.compute("integrated")
+
+    sw = sw_obj.compute("integrated")["SW"]
     assert not sw.values.mask[0, 0, 0]
 
     # now make a grid with masked cells
@@ -86,7 +89,7 @@ def test_swj_simple_x_zero():
         invert=True,
         method="cell_center_above_ffl",
     )
-    sw = sw_obj.compute("integrated")
+    sw = sw_obj.compute("integrated")["SW"]
     assert sw.values.mask[0, 0, 0]
 
 
@@ -119,7 +122,7 @@ def test_swj_simple_threshold_2grids():
         invert=True,
         method=cellmethod,
     )
-    sw = sw_obj.compute("integrated")
+    sw = sw_obj.compute("integrated")["SW"]
     assert sw.values.mean() == pytest.approx(0.9689, rel=0.01)
 
     sw_obj = SwFunction(
@@ -131,7 +134,7 @@ def test_swj_simple_threshold_2grids():
         invert=True,
         method=cellmethod,
     )
-    sw = sw_obj.compute("integrated")
+    sw = sw_obj.compute("integrated")["SW"]
 
     assert sw.values.mean() == pytest.approx(0.9689, rel=0.01)
     assert float(sw.values[0, 0, 70]) == pytest.approx(1, rel=0.0001)
@@ -143,6 +146,8 @@ def test_swj_simple_threshold_2grids():
         (True, "cell_center_above_ffl", 0.7057, 0.046719),  # n/a vs RMS
         (False, "cell_center_above_ffl", 0.70736, 0.046724),  # n/a vs RMS
         (False, "cell_corners_above_ffl", 0.674485, 0.046791),  # n/a vs RMS
+        (True, "truncated_cell_corners_above_ffl", 0.67737, 0.046719),  # n/a vs RMS
+        (False, "truncated_cell_corners_above_ffl", 0.67931, 0.046724),  # n/a vs RMS
     ],
 )
 def test_swj_simple_reek(direct, cellmethod, expected_mean, exp_cell1, testdata_path):
@@ -174,7 +179,7 @@ def test_swj_simple_reek(direct, cellmethod, expected_mean, exp_cell1, testdata_
         method=cellmethod,
     )
 
-    sw = sw_obj.compute("direct" if direct else "integrated")
+    sw = sw_obj.compute("direct" if direct else "integrated")["SW"]
 
     assert sw.values.mean() == pytest.approx(expected_mean, rel=0.01)
 
@@ -252,17 +257,16 @@ def test_sw_bvw():
         x=xvalue,
         ffl=ffl,
         method="cell_center_above_ffl",
-        debug=False,
     )
-    sw = sw_obj.compute("direct")
+    sw = sw_obj.compute("direct")["SW"]
 
     sw10 = float(sw.values[:, :, 20])  # 10 meter above contact
     assert sw10 == pytest.approx(manual10)
 
     sw20 = float(sw.values[:, :, 10])  # 20 meter above contact
     assert sw20 == pytest.approx(0.13755086)
 
-    sw = sw_obj.compute("integrated")
+    sw = sw_obj.compute("integrated")["SW"]
     sw10_i = float(sw.values[:, :, 20])
     assert sw10_i == pytest.approx(sw10, abs=0.0001)
 
@@ -332,9 +336,8 @@ def test_sw_brooks_corey():
         x=xvalue,
         ffl=ffl,
         method="cell_center_above_ffl",
-        debug=False,
     )
-    sw = sw_obj.compute("direct")
+    sw = sw_obj.compute("direct")["SW"]
 
     sw10 = float(sw.values[:, :, 20])  # 10 meter above contact
     assert sw10 == pytest.approx(0.40303321)
@@ -344,6 +347,6 @@ def test_sw_brooks_corey():
     sw20 = float(sw.values[:, :, 10])  # 20 meter above contact
     assert sw20 == pytest.approx(0.28731234)
 
-    sw = sw_obj.compute("integrated")
+    sw = sw_obj.compute("integrated")["SW"]
     sw10_i = float(sw.values[:, :, 20])
     assert sw10_i == pytest.approx(sw10, abs=0.001)