Add layer_id to each well connection and update _grid_cell_bounding_b…

…oxes() (#279)

* layers added to input yaml and parsed by the config parser

* config layers variable passed to from_flow

* Update _grid_cell_bounding_boxes to allow layering filters

* Minor fixes

* Added value errors for layers input

* Added layer id to well connections

* fix pylint issues

* Add tests + fixes

* Made test model agnostic

* Remove assert True statement

* Change path

* Typo

* Modify path

* Fix version conflict

* Allow multiple paths to test data

* Typo

* Add leading underscore to dunction

* synchronize changes after merging branch i276-grid-cell-bounding-box

* merge updated branch grid_cell_bounding_box

* merge updated branch grid_cell_bounding_box

* Revert to no changes in perforation strategy

* change layer_id from None to 0 when no layering is defined

* Apply suggestions from code review

Accepted all suggestions. Only layering to layers still has to be changed.

Co-authored-by: Wouter J. de Bruin <9119793+wouterjdb@users.noreply.github.com>

* layering changed to layers

* forgot to change one layering to layers

Co-authored-by: Ubuntu <LonnekevB@Flownet.gwfddzla5snuvfeg5bcr4imt2c.ax.internal.cloudapp.net>
Co-authored-by: Wouter J. de Bruin <wouterjdb@hotmail.com>
Co-authored-by: Wouter J. de Bruin <9119793+wouterjdb@users.noreply.github.com>

src/flownet/ahm/_run_ahm.py

@@ -399,6 +399,7 @@ def run_flownet_history_matching(
     # Load production and well coordinate data
     field_data = FlowData(
         config.flownet.data_source.simulation.input_case,
+        layers=config.flownet.data_source.simulation.layers,
         perforation_handling_strategy=config.flownet.perforation_handling_strategy,
     )
     df_production_data: pd.DataFrame = field_data.production

src/flownet/config_parser/_config_parser.py

@@ -274,6 +274,20 @@ def _to_abs_path(path: Optional[str]) -> str:
                                             },
                                         },
                                     },
+                                    "layers": {
+                                        MK.Type: types.List,
+                                        MK.Content: {
+                                            MK.Item: {
+                                                MK.Type: types.List,
+                                                MK.Content: {
+                                                    MK.Item: {
+                                                        MK.Type: types.Integer,
+                                                        MK.AllowNone: True,
+                                                    }
+                                                },
+                                            },
+                                        },
+                                    },
                                 },
                             },
                             "resampling": {
@@ -1569,6 +1583,26 @@ def parse_config(
 
     config = suite.snapshot
 
+    layers = config.flownet.data_source.simulation.layers
+    if len(layers) > 0:
+        if not all(
+            d == 1
+            for d in [layers[i][0] - layers[i - 1][-1] for i in range(1, len(layers))]
+        ):
+            raise ValueError(
+                f"The layering definition "
+                f"'{layers}' is not valid.\n"
+                f"Layer groups should be adjacent, e.g. ((start, end),(end+1, ..)) ."
+            )
+
+        if not all(l == 2 for l in [len(layer) for layer in layers]):
+            raise ValueError(
+                f"The layering definition "
+                f"'{layers}' is not valid.\n"
+                f"Valid way to define the layering is as nested list with "
+                f"layer intervals. e.g. ((1, 5),(6, 9),(10, 15))."
+            )
+
     req_relp_parameters: List[str] = []
     if (
         config.model_parameters.equil.scheme != "regions_from_sim"

src/flownet/data/from_flow.py

@@ -1,6 +1,7 @@
 import warnings
+import operator
 from pathlib import Path
-from typing import Union, List
+from typing import Union, List, Optional, Tuple
 
 import numpy as np
 import pandas as pd
@@ -20,6 +21,7 @@ class FlowData(FromSource):
 
     Args:
          input_case: Full path to eclipse case to load data from
+         layers: List with definition of isolated layers, if present.
          perforation_handling_strategy: How to deal with perforations per well.
                                                  ('bottom_point', 'top_point', 'multiple')
 
@@ -28,6 +30,7 @@ class FlowData(FromSource):
     def __init__(
         self,
         input_case: Union[Path, str],
+        layers: Tuple = (),
         perforation_handling_strategy: str = "bottom_point",
     ):
         super().__init__()
@@ -39,25 +42,40 @@ def __init__(
         self._restart = EclFile(str(self._input_case.with_suffix(".UNRST")))
         self._init = EclInitFile(self._grid, str(self._input_case.with_suffix(".INIT")))
         self._wells = compdat.df(EclFiles(str(self._input_case)))
+        self._layers = layers
 
         self._perforation_handling_strategy: str = perforation_handling_strategy
 
     # pylint: disable=too-many-branches
     def _well_connections(self) -> pd.DataFrame:
         """
-        Function to extract well connection coordinates from an Flow simulation including their
+        Function to extract well connection coordinates from a Flow simulation including their
         opening and closure time. The output of this function will be filtered based on the
         configured perforation strategy.
 
         Returns:
-            columns: WELL_NAME, X, Y, Z, DATE, OPEN
+            columns: WELL_NAME, X, Y, Z, DATE, OPEN, LAYER_ID
 
         """
+        if len(self._layers) > 0 and self._grid.nz is not self._layers[-1][-1]:
+            raise ValueError(
+                f"Number of layers from config ({self._layers[-1][-1]}) is not equal to "
+                f"number of layers from flow simulation ({self._grid.nz})."
+            )
+
         new_items = []
         for _, row in self._wells.iterrows():
             X, Y, Z = self._grid.get_xyz(
                 ijk=(row["I"] - 1, row["J"] - 1, row["K1"] - 1)
             )
+            if len(self._layers) > 0:
+                for count, (i, j) in enumerate(self._layers):
+                    if row["K1"] in range(i, j + 1):
+                        layer_id = count
+                        break
+            else:
+                layer_id = 0
+
             new_row = {
                 "WELL_NAME": row["WELL"],
                 "IJK": (
@@ -70,11 +88,13 @@ def _well_connections(self) -> pd.DataFrame:
                 "Z": Z,
                 "DATE": row["DATE"],
                 "OPEN": bool(row["OP/SH"] == "OPEN"),
+                "LAYER_ID": layer_id,
             }
             new_items.append(new_row)
 
         df = pd.DataFrame(
-            new_items, columns=["WELL_NAME", "IJK", "X", "Y", "Z", "DATE", "OPEN"]
+            new_items,
+            columns=["WELL_NAME", "IJK", "X", "Y", "Z", "DATE", "OPEN", "LAYER_ID"],
         )
         df["DATE"] = pd.to_datetime(df["DATE"], format="%Y-%m-%d").dt.date
 
@@ -297,19 +317,30 @@ def _faults(self) -> pd.DataFrame:
 
         return df_faults.drop(["I", "J", "K"], axis=1)
 
-    def _grid_cell_bounding_boxes(self) -> np.ndarray:
+    def _grid_cell_bounding_boxes(self, layer_id: Optional[int] = None) -> np.ndarray:
         """
         Function to get the bounding box (x, y and z min + max) for all grid cells
 
+        Args:
+            layer_id: The FlowNet layer id to be used to create the bounding box.
+
         Returns:
             A (active grid cells x 6) numpy array with columns [ xmin, xmax, ymin, ymax, zmin, zmax ]
+            filtered on layer_id if not None.
         """
-        xyz = np.empty((8 * self._grid.get_num_active(), 3))
-        for active_index in range(self._grid.get_num_active()):
-            for corner in range(0, 8):
-                xyz[active_index * 8 + corner, :] = self._grid.get_cell_corner(
-                    corner, active_index=active_index
-                )
+        if layer_id is not None:
+            (k_min, k_max) = tuple(map(operator.sub, self._layers[layer_id], (1, 1)))
+        else:
+            (k_min, k_max) = (0, self._grid.nz)
+
+        cells = [
+            cell for cell in self._grid.cells(active=True) if (k_min <= cell.k <= k_max)
+        ]
+        xyz = np.empty((8 * len(cells), 3))
+
+        for n_cell, cell in enumerate(cells):
+            for n_corner, corner in enumerate(cell.corners):
+                xyz[n_cell * 8 + n_corner, :] = corner
 
         xmin = xyz[:, 0].reshape(-1, 8).min(axis=1)
         xmax = xyz[:, 0].reshape(-1, 8).max(axis=1)
@@ -360,3 +391,8 @@ def well_logs(self) -> pd.DataFrame:
     def grid(self) -> EclGrid:
         """the simulation grid with properties"""
         return self._grid
+
+    @property
+    def layers(self) -> Union[Tuple[Tuple[int, int]], Tuple]:
+        """Get the list of top and bottom k-indeces of a the orignal model that represents a FlowNet layer"""
+        return self._layers

src/flownet/data/perforation_strategy.py

@@ -81,7 +81,7 @@ def multiple(df: pd.DataFrame) -> pd.DataFrame:
         df: Dataframe with all well connections, through time, including state.
 
     Returns:
-        DataFrame will all connections
+        DataFrame with all connections
 
     """
     df = df[["WELL_NAME", "X", "Y", "Z", "DATE", "OPEN"]].sort_values(
@@ -122,7 +122,7 @@ def multiple_based_on_workovers(df: pd.DataFrame) -> pd.DataFrame:
         df: Dataframe with all well connections, through time, including state.
 
     Returns:
-        Dataframe with 1 or more connections per will depending on the historic straddles / plugs.
+        Dataframe with 1 or more connections per well depending on the historic straddles / plugs.
 
     """
     df = multiple(df)

tests/test_data_flow.py

@@ -0,0 +1,85 @@
+import glob
+import operator
+from pathlib import Path
+
+from numpy.testing import assert_almost_equal
+from ecl.grid import EclRegion
+
+from flownet.data.from_flow import FlowData
+
+
+def _locate_test_case() -> Path:
+    """
+    This function will try to find the test data. On the CI/CD this will be
+    any case. Locally, the Norne case will be used. If not found, a
+    FileNotFoundError will be raised.
+
+    Raises:
+        FileNotFoundError: The test data was not found.
+
+    Returns:
+        Path to the .DATA file.
+    """
+    cicd_str = "./flownet-testdata/*/input_model/*/*.DATA"
+    home_path = Path.home() / Path(
+        "flownet-testdata/norne/input_model/norne/NORNE_ATW2013.DATA"
+    )
+    flow_path = (
+        Path.cwd() / "../flownet-testdata/norne/input_model/norne/NORNE_ATW2013.DATA"
+    )
+
+    if len(glob.glob(cicd_str)) > 0:
+        # CI/CD
+        data_file = Path(glob.glob(cicd_str)[0])
+    elif home_path.exists():
+        # Test-data repository located in home directory
+        data_file = home_path
+    elif flow_path.exists():
+        # Test-data repositry located in the same location as the flownet repository
+        data_file = flow_path
+    else:
+        raise FileNotFoundError(
+            "To be able to run the tests one needs to clone the flownet-testdata from \
+            https://github.com/equinor/flownet-testdata. The test-data needs to be located \
+            in your home folder OR in the folder containing the flownet repository."
+        )
+
+    return data_file
+
+
+# pylint: disable=protected-access
+def test_grid_cell_bounding_boxes() -> None:
+    layers = ()
+    flowdata = FlowData(
+        _locate_test_case(),
+        layers,
+        "multiple_based_on_workovers",
+    )
+
+    # Test no argument and entire field being equal
+    flowdata._layers = ((1, flowdata.grid.nz),)
+    assert_almost_equal(
+        flowdata._grid_cell_bounding_boxes(), flowdata._grid_cell_bounding_boxes(0)
+    )
+
+    # Test zero'th layer id
+    flowdata._layers = ((1, 2), (3, 4))
+    result = flowdata._grid_cell_bounding_boxes(0)
+    active_cells = EclRegion(flowdata.grid, True)
+    active_cells.select_kslice(
+        *tuple(map(operator.sub, flowdata._layers[0], (1, 1))), intersect=True
+    )
+    active_cells.select_active(intersect=True)
+    assert result.shape[0] == active_cells.active_size()
+    assert result.shape[1] == 6
+
+    # Test last layer id
+    flowdata._layers = ((1, 2), (3, 4))
+    result = flowdata._grid_cell_bounding_boxes(1)
+    active_cells = EclRegion(flowdata.grid, True)
+    active_cells.select_kslice(
+        *tuple(map(operator.sub, flowdata._layers[-1], (1, 1))), intersect=True
+    )
+    active_cells.select_active(intersect=True)
+    assert result.shape[0] == active_cells.active_size()
+    assert result.shape[1] == 6