Add option to use model hull to steer node/connection generation

.github/workflows/flownet.yml

@@ -83,7 +83,7 @@ jobs:
           source $VENV_PATH/bin/activate
           black --check examples/ tests/ src/ setup.py
           pylint src/ tests/ setup.py
-          mypy --ignore-missing-imports src/ tests/ setup.py
+          mypy --ignore-missing-imports src/ setup.py
 
       - name: 🤖 Run tests
         run: |

src/flownet/ahm/_run_ahm.py

@@ -119,7 +119,7 @@ def _get_distribution(
     return df
 
 
-# pylint: disable=too-many-branches
+# pylint: disable=too-many-branches,too-many-statements
 def run_flownet_history_matching(
     config: ConfigSuite.snapshot, args: argparse.Namespace
 ):
@@ -154,7 +154,13 @@ def run_flownet_history_matching(
         pd.DataFrame
     ] = field_data.faults if config.model_parameters.fault_mult else None
 
-    df_connections: pd.DataFrame = create_connections(df_coordinates, config)
+    concave_hull_bounding_boxes: Optional[np.ndarray] = None
+    if config.flownet.data_source.concave_hull:
+        concave_hull_bounding_boxes = field_data.grid_cell_bounding_boxes
+
+    df_connections: pd.DataFrame = create_connections(
+        df_coordinates, config, concave_hull_bounding_boxes=concave_hull_bounding_boxes,
+    )
 
     network = NetworkModel(
         df_connections,

src/flownet/config_parser/_config_parser.py

@@ -33,6 +33,7 @@ def create_schema(_to_abs_path) -> Dict:
                                 MK.Type: types.String,
                                 MK.Transformation: _to_abs_path,
                             },
+                            "concave_hull": {MK.Type: types.Bool, MK.AllowNone: True},
                         },
                     },
                     "pvt": {
@@ -69,6 +70,8 @@ def create_schema(_to_abs_path) -> Dict:
                     },
                     "fast_pyscal": {MK.Type: types.Bool, MK.Default: True},
                     "fault_tolerance": {MK.Type: types.Number, MK.Default: 1.0e-5},
+                    "max_distance": {MK.Type: types.Number, MK.Default: 1e12},
+                    "max_distance_fraction": {MK.Type: types.Number, MK.Default: 0},
                 },
             },
             "ert": {

src/flownet/data/from_flow.py

@@ -46,12 +46,6 @@ def _coordinates(self) -> pd.DataFrame:
         """
         Function to extract well coordinates from an Flow simulation.
 
-        Args:
-            filename: Entire path to the simulated simulation case. This
-                            case must have both and EGRID and UNRST file.
-            perforation_handling_strategy: How to deal with perforations per well.
-                                                ('bottom_point', 'top_point', 'multiple')
-
         Returns:
             columns: WELL_NAME, X, Y, Z
 
@@ -262,9 +256,37 @@ def _faults(self) -> pd.DataFrame:
 
         return df_faults.drop(["I", "J", "K"], axis=1)
 
+    def _grid_cell_bounding_boxes(self) -> np.ndarray:
+        """
+        Function to get the bounding box (x, y and z min + max) for all grid cells
+
+        Returns:
+            A (active grid cells x 6) numpy array with columns [ xmin, xmax, ymin, ymax, zmin, zmax ]
+        """
+        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
+                )
+
+        xmin = xyz[:, 0].reshape(-1, 8).min(axis=1)
+        xmax = xyz[:, 0].reshape(-1, 8).max(axis=1)
+        ymin = xyz[:, 1].reshape(-1, 8).min(axis=1)
+        ymax = xyz[:, 1].reshape(-1, 8).max(axis=1)
+        zmin = xyz[:, 2].reshape(-1, 8).min(axis=1)
+        zmax = xyz[:, 2].reshape(-1, 8).max(axis=1)
+
+        return np.vstack([xmin, xmax, ymin, ymax, zmin, zmax]).T
+
     def _get_start_date(self):
         return self._eclsum.start_date
 
+    @property
+    def grid_cell_bounding_boxes(self) -> np.ndarray:
+        """Boundingboxes for all gridcells"""
+        return self._grid_cell_bounding_boxes()
+
     @property
     def faults(self) -> pd.DataFrame:
         """dataframe with all fault data"""

src/flownet/network_model/_generate_connections.py

@@ -1,6 +1,6 @@
 import math
 import time
-from typing import List, Tuple, Any
+from typing import List, Tuple, Any, Optional
 from operator import itemgetter
 
 import numpy as np
@@ -89,7 +89,9 @@ def _remove_recorded_connections(
 
 
 def _generate_connections(
-    df_coordinates: pd.DataFrame, configuration: Any
+    df_coordinates: pd.DataFrame,
+    configuration: Any,
+    concave_hull_bounding_boxes: Optional[np.ndarray] = None,
 ) -> Tuple[List[Coordinate], List[Coordinate]]:
     """
     Uses MitchellBestCandidate and Delaunay triangulation to populate the reservoir
@@ -99,6 +101,7 @@ def _generate_connections(
     Args:
         df_coordinates: coordinates on original DataFrame format
         configuration: FlowNet configuration yaml
+        concave_hull_bounding_boxes: Numpy array with x, y, z min/max boundingboxes for each grid block
 
     Returns:
         The tuple consisting of the start_coords and end_coords for connections.
@@ -122,6 +125,7 @@ def _generate_connections(
             num_added_flow_nodes=configuration.flownet.additional_flow_nodes,
             num_candidates=configuration.flownet.additional_node_candidates,
             hull_factor=configuration.flownet.hull_factor,
+            concave_hull_bounding_boxes=concave_hull_bounding_boxes,
             random_seed=configuration.flownet.random_seed,
         )
     ]
@@ -188,7 +192,7 @@ def are_points_from_same_existing_entity(point1: Coordinate, point2: Coordinate)
         )
     print("done.")
 
-    # Are there nodes in the connection matrix that has no connections? Definately a problem
+    # Are there nodes in the connection matrix that has no connections? Definitely a problem
     # it is one of the original well perforations
     print("Checking connectivity of well-perforations...")
     _check_for_none_connectivity_amongst_entities(conn_matrix)
@@ -296,12 +300,15 @@ def __get_entity_str(df_coordinates: pd.DataFrame, coordinate: Coordinate) -> st
     return entity
 
 
+# pylint: disable=too-many-arguments
 def _create_entity_connection_matrix(
     df_coordinates: pd.DataFrame,
     starts: List[Coordinate],
     ends: List[Coordinate],
     aquifer_starts: List[Coordinate],
     aquifer_ends: List[Coordinate],
+    max_distance_fraction: float,
+    max_distance: float,
 ) -> pd.DataFrame:
     """
     Converts the the coordinates given for starts and ends to the desired DataFrame format for simulation input.
@@ -312,6 +319,8 @@ def _create_entity_connection_matrix(
         ends: List of coordinates of all the end entities
         aquifer_starts: List of coordinates for all aquifer starts
         aquifer_ends: List of coordinates of all aquifer ends
+        max_distance_fraction: Fraction of longest connection distance to be removed
+        max_distance: Maximum distance between nodes, removed otherwise
 
     Returns:
         Connection coordinate DataFrame on Flow desired format.
@@ -365,10 +374,63 @@ def _create_entity_connection_matrix(
             ignore_index=True,
         )
 
+    df_out = _remove_long_connections(df_out, max_distance_fraction, max_distance)
+
     print("done.")
     return df_out
 
 
+def _remove_long_connections(
+    df_connections: pd.DataFrame,
+    max_distance_fraction: float = 0,
+    max_distance: float = 1e12,
+) -> pd.DataFrame:
+    """
+    Helper function to remove long connections.
+
+    Args:
+        df_connections: Pandas dataframe with start point, end point and entity type
+        max_distance_fraction: Fraction of longest connections to drop. I.e., 0.1 will drop the 10% longest connections.
+        max_distance: Maximum length of a connection; connections longer than this value will be dropped.
+
+    Returns:
+        Input DataFrame without long connections
+
+    """
+    if max_distance_fraction == 0 and max_distance == 1e12:
+        return df_connections
+
+    def __calculate_distance(row: Any) -> float:
+        """
+        Helper function that calculates the distance between two nodes.
+
+        Args:
+            row: Pandas dataframe row object
+
+        Returns:
+            Distance
+
+        """
+        return (
+            (row["xend"] - row["xstart"]) ** 2
+            + (row["yend"] - row["ystart"]) ** 2
+            + (row["zend"] - row["zstart"]) ** 2
+        ) ** 0.5
+
+    df_connections["distance"] = df_connections.apply(__calculate_distance, axis=1)
+
+    max_distance = min(
+        df_connections["distance"].quantile(1 - max_distance_fraction), max_distance
+    )
+
+    df_connections = df_connections[
+        (df_connections["distance"] < max_distance)
+        | (df_connections["end_entity"] == "aquifer")
+    ]
+
+    return df_connections
+
+
 def _generate_aquifer_connections(
     starts: List[Coordinate],
     ends: List[Coordinate],
@@ -431,7 +493,9 @@ def _generate_aquifer_connections(
 
 
 def create_connections(
-    df_coordinates: pd.DataFrame, configuration: Any
+    df_coordinates: pd.DataFrame,
+    configuration: Any,
+    concave_hull_bounding_boxes: Optional[np.ndarray] = None,
 ) -> pd.DataFrame:
     """
     Creates additional flow nodes to increase complexity of field simulation structure so that history-matching can
@@ -443,12 +507,17 @@ def create_connections(
     Args:
         df_coordinates: Original structure of entity and X, Y, Z coords
         configuration: FlowNet configuration yaml as dictionary
+        concave_hull_bounding_boxes: Numpy array with x, y, z min/max boundingboxes for each grid block
 
     Returns:
         Desired restructuring of start-end coordinates into separate columns, as per Flow needs.
 
     """
-    starts, ends = _generate_connections(df_coordinates, configuration)
+    starts, ends = _generate_connections(
+        df_coordinates=df_coordinates,
+        configuration=configuration,
+        concave_hull_bounding_boxes=concave_hull_bounding_boxes,
+    )
     aquifer_starts: List[Coordinate] = []
     aquifer_ends: List[Coordinate] = []
 
@@ -463,5 +532,11 @@ def create_connections(
         )
 
     return _create_entity_connection_matrix(
-        df_coordinates, starts, ends, aquifer_starts, aquifer_ends
+        df_coordinates,
+        starts,
+        ends,
+        aquifer_starts,
+        aquifer_ends,
+        configuration.flownet.max_distance_fraction,
+        configuration.flownet.max_distance,
     )

src/flownet/network_model/_mitchell.py

@@ -7,11 +7,13 @@
 from ..utils.types import Coordinate
 
 
+# pylint: disable=too-many-branches,too-many-statements
 def mitchell_best_candidate_modified_3d(
     perforations: List[Coordinate],
     num_added_flow_nodes: int,
     num_candidates: int = 1000,
     hull_factor: Optional[float] = None,
+    concave_hull_bounding_boxes: Optional[np.ndarray] = None,
     random_seed: Optional[int] = None,
 ) -> List[Coordinate]:
 
@@ -34,6 +36,7 @@ def mitchell_best_candidate_modified_3d(
         hull_factor: Factor to linearly scale the convex hull with. Factor will
             scale the distance of each point from the centroid of all the points.
             When a None is supplied a box-shape is used.
+        concave_hull_bounding_boxes: Numpy array with x, y, z min/max boundingboxes for each grid block
         random_seed: Random seed to control the reproducibility of the FlowNet.
 
     Returns:
@@ -102,6 +105,14 @@ def mitchell_best_candidate_modified_3d(
         y_candidate = np.zeros(num_candidates)
         z_candidate = np.zeros(num_candidates)
 
+        if concave_hull_bounding_boxes is not None:
+            xmin_grid_cells = concave_hull_bounding_boxes[:, 0]
+            xmax_grid_cells = concave_hull_bounding_boxes[:, 1]
+            ymin_grid_cells = concave_hull_bounding_boxes[:, 2]
+            ymax_grid_cells = concave_hull_bounding_boxes[:, 3]
+            zmin_grid_cells = concave_hull_bounding_boxes[:, 4]
+            zmax_grid_cells = concave_hull_bounding_boxes[:, 5]
+
         # Repeat while not all random points are inside the convex hull
         while not all(in_hull):
             # Generate a set of random candidates that will be the new
@@ -116,11 +127,29 @@ def mitchell_best_candidate_modified_3d(
             np.putmask(y_candidate, np.invert(in_hull), y_candidate_tmp)
             np.putmask(z_candidate, np.invert(in_hull), z_candidate_tmp)
 
-            # Test whether all points are inside the convex hull of the perforations
-            in_hull = (
-                hull.find_simplex(np.vstack([x_candidate, y_candidate, z_candidate]).T)
-                >= 0
-            )
+            candidates = np.vstack([x_candidate, y_candidate, z_candidate]).T
+
+            if concave_hull_bounding_boxes is not None:
+                for c_index, candidate in enumerate(candidates):
+                    if not in_hull[c_index]:
+                        in_hull[c_index] = (
+                            (
+                                (candidate[0] >= xmin_grid_cells)
+                                & (candidate[0] <= xmax_grid_cells)
+                            )
+                            & (
+                                (candidate[1] >= ymin_grid_cells)
+                                & (candidate[1] <= ymax_grid_cells)
+                            )
+                            & (
+                                (candidate[2] >= zmin_grid_cells)
+                                & (candidate[2] <= zmax_grid_cells)
+                            )
+                        ).any()
+
+            else:
+                # Test whether all points are inside the convex hull of the perforations
+                in_hull = hull.find_simplex(candidates) >= 0
 
         best_distance = 0
         best_candidate = 0