Improved model generation

CHANGELOG.md

@@ -5,6 +5,7 @@ This project adheres to [Semantic Versioning](https://semver.org/).
 ## Unreleased
 
 ### Added
+- [#298](https://github.com/equinor/flownet/pull/298) Connections between (well)nodes that go through non-reservoir are now removed. Angle threshold export to user.
 - [#284](https://github.com/equinor/flownet/pull/284) Added the option to specify cumulative phase rates as observations (WOPT, WWPT, WGPT, WGIT, WWIT)
 
 ### Fixes

src/flownet/config_parser/_config_parser.py

@@ -424,6 +424,18 @@ def _to_abs_path(path: Optional[str]) -> str:
                         MK.Type: types.String,
                         MK.Default: "RATE",
                     },
+                    "angle_threshold": {
+                        MK.Type: types.Number,
+                        MK.Default: 150,
+                        MK.Description: "Angle threshold used, after Delaunay triagulation,"
+                        "to remove sides/tubes opposite angles larger than the supplied treshold.",
+                    },
+                    "n_non_reservoir_evaluation": {
+                        MK.Type: types.Number,
+                        MK.Default: 10,
+                        MK.Description: "Number of points along a tube to check whether they are in"
+                        "non reservoir for removal purposes.",
+                    },
                     "hyperopt": {
                         MK.Type: types.NamedDict,
                         MK.Content: {

src/flownet/network_model/_generate_connections.py

@@ -4,66 +4,72 @@
 from operator import itemgetter
 
 import numpy as np
+from numpy.core.function_base import linspace
 import pandas as pd
 from scipy.spatial import Delaunay, distance  # pylint: disable=no-name-in-module
 
 from ._mitchell import mitchell_best_candidate_modified_3d
+from ._hull import check_in_hull
 from ..utils.types import Coordinate
 
 
-def _create_record(
-    well_pairs: np.ndarray, dist_matrix: np.ndarray
-) -> List[Tuple[Any, ...]]:
-    """
-    Creates a record of every well_pair triangle for which one of the angles is too large
-    (currently set to be above 150 degrees)
+def _is_angle_too_large(
+    angle_threshold: float, side_a: float, side_b: float, side_c: float
+) -> bool:
+    """Function checks if there is an angle larger than a specified angle in
+    degrees for a triangle with given side lengths
 
     Args:
-        well_pairs: All well pairs created through Delaunay triangulation
-        dist_matrix: Euclidean distance between well coordinates
+        angle_threshold: threshold angle in degrees
+        side_a: Length of side a
+        side_b: Length of side b
+        side_c: Length of side c
 
     Returns:
-         A record of connections (from_well, to_well) of undesirables.
+        True if an angle larger than the specified angle
 
     """
-    record: List[Tuple[Any, ...]] = []
+    calculate_angle = np.rad2deg(
+        math.acos(
+            (math.pow(side_a, 2) + math.pow(side_b, 2) - math.pow(side_c, 2))
+            / (2 * side_a * side_b)
+        )
+    )
 
-    def is_angle_too_large(side_a: float, side_b: float, side_c: float) -> bool:
-        """Function checks if there is an angle larger than 150 degrees for a triangle with given side lengths
+    return calculate_angle > angle_threshold
 
-        Args:
-            side_a: Length of side a
-            side_b: Length of side b
-            side_c: Length of side c
 
-        Returns:
-            True if an angle larger than 150 degrees exists
+def _create_record(
+    connection_pairs: np.ndarray, dist_matrix: np.ndarray, angle_threshold: float
+) -> List[Tuple[Any, ...]]:
+    """
+    Creates a record of every connection_pairs triangle for which one of the angles is too small/large
 
-        """
-        return (
-            np.rad2deg(
-                math.acos(
-                    (math.pow(side_a, 2) + math.pow(side_b, 2) - math.pow(side_c, 2))
-                    / (2 * side_a * side_b)
-                )
-            )
-            > 150
-        )
+    Args:
+        connection_pairs: All connection_pairs created through Delaunay triangulation
+        dist_matrix: Euclidean distance between coordinate pairs
+        angle_threshold: angle in Delaunay triangles for which to ignore connections
+
+    Returns:
+         A record of undesirable connection_pairs.
+
+    """
+    record: List[Tuple[Any, ...]] = []
 
     for vertex_a in range(0, 2):
         for vertex_b in range(vertex_a + 1, 3):
             for vertex_c in range(vertex_b + 1, 4):
                 edge_a = dist_matrix[vertex_a, vertex_b]
                 edge_b = dist_matrix[vertex_a, vertex_c]
                 edge_c = dist_matrix[vertex_b, vertex_c]
-                if is_angle_too_large(edge_a, edge_b, edge_c):
-                    record.append(tuple(well_pairs[vertex_b, vertex_c]))
+                if _is_angle_too_large(angle_threshold, edge_a, edge_b, edge_c):
+                    record.append(tuple(connection_pairs[vertex_b, vertex_c]))
 
-                if is_angle_too_large(edge_c, edge_b, edge_a):
-                    record.append(tuple(well_pairs[vertex_a, vertex_b]))
+                if _is_angle_too_large(angle_threshold, edge_c, edge_b, edge_a):
+                    record.append(tuple(connection_pairs[vertex_a, vertex_b]))
 
-                if is_angle_too_large(edge_a, edge_c, edge_b):
-                    record.append(tuple(well_pairs[vertex_a, vertex_c]))
+                if _is_angle_too_large(angle_threshold, edge_a, edge_c, edge_b):
+                    record.append(tuple(connection_pairs[vertex_a, vertex_c]))
     return record
 
 
@@ -187,9 +193,13 @@ def are_points_from_same_existing_entity(point1: Coordinate, point2: Coordinate)
                 conn_matrix[flow_node_a, flow_node_b] = conn_matrix[
                     flow_node_b, flow_node_a
                 ] = 1
-        conn_matrix = _remove_recorded_connections(
-            _create_record(well_pairs, dist_matrix), conn_matrix
-        )
+        if configuration.flownet.angle_threshold:
+            conn_matrix = _remove_recorded_connections(
+                _create_record(
+                    well_pairs, dist_matrix, configuration.flownet.angle_threshold
+                ),
+                conn_matrix,
+            )
     print("done.")
 
     # Are there nodes in the connection matrix that has no connections? Definitely a problem
@@ -300,7 +310,7 @@ def __get_entity_str(df_coordinates: pd.DataFrame, coordinate: Coordinate) -> st
     return entity
 
 
-# pylint: disable=too-many-arguments
+# pylint: disable=too-many-arguments,too-many-locals
 def _create_entity_connection_matrix(
     df_coordinates: pd.DataFrame,
     starts: List[Coordinate],
@@ -309,6 +319,8 @@ def _create_entity_connection_matrix(
     aquifer_ends: List[Coordinate],
     max_distance_fraction: float,
     max_distance: float,
+    concave_hull_bounding_boxes: Optional[np.ndarray] = None,
+    n_non_reservoir_evaluation: Optional[int] = 10,
 ) -> pd.DataFrame:
     """
     Converts the the coordinates given for starts and ends to the desired DataFrame format for simulation input.
@@ -321,6 +333,8 @@ def _create_entity_connection_matrix(
         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
+        concave_hull_bounding_boxes: Numpy array with x, y, z min/max boundingboxes for each grid block
+        n_non_reservoir_evaluation: Number of equally spaced points along a connection to check fornon-reservoir.
 
     Returns:
         Connection coordinate DataFrame on Flow desired format.
@@ -343,6 +357,19 @@ def _create_entity_connection_matrix(
         str_start_entity = __get_entity_str(df_coordinates, start)
         str_end_entity = __get_entity_str(df_coordinates, end)
 
+        if concave_hull_bounding_boxes is not None:
+            tube_coordinates = linspace(
+                start=start,
+                stop=end,
+                num=n_non_reservoir_evaluation,
+                endpoint=False,
+                dtype=float,
+                axis=1,
+            ).T
+
+            if not all(check_in_hull(concave_hull_bounding_boxes, tube_coordinates)):
+                continue
+
         df_out = df_out.append(
             {
                 "xstart": start[0],
@@ -539,4 +566,6 @@ def create_connections(
         aquifer_ends,
         configuration.flownet.max_distance_fraction,
         configuration.flownet.max_distance,
+        concave_hull_bounding_boxes=concave_hull_bounding_boxes,
+        n_non_reservoir_evaluation=configuration.flownet.n_non_reservoir_evaluation,
     ).reset_index()

src/flownet/network_model/_hull.py

@@ -0,0 +1,45 @@
+from typing import List
+
+import numpy as np
+
+
+def check_in_hull(
+    concave_hull_bounding_boxes: np.ndarray, coordinates: np.ndarray
+) -> List[bool]:
+    """Checks if all coordinates are inside the concave hull bounding boxes.
+
+    Args:
+        concave_hull_bounding_boxes (np.ndarray): Bounding boxes for the concave hull
+        connections (np.ndarray): All coordianates to check
+
+    Returns:
+        List[bool]: List
+    """
+
+    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]
+
+    in_hull = np.asarray([False] * coordinates.shape[0])
+
+    for c_index, coordinate in enumerate(coordinates):
+        if not in_hull[c_index]:
+            in_hull[c_index] = (
+                (
+                    (coordinate[0] >= xmin_grid_cells)
+                    & (coordinate[0] <= xmax_grid_cells)
+                )
+                & (
+                    (coordinate[1] >= ymin_grid_cells)
+                    & (coordinate[1] <= ymax_grid_cells)
+                )
+                & (
+                    (coordinate[2] >= zmin_grid_cells)
+                    & (coordinate[2] <= zmax_grid_cells)
+                )
+            ).any()
+
+    return in_hull

src/flownet/network_model/_mitchell.py

@@ -5,7 +5,7 @@
 import numpy as np
 
 from ..utils.types import Coordinate
-
+from ._hull import check_in_hull
 
 # pylint: disable=too-many-branches,too-many-statements
 def mitchell_best_candidate_modified_3d(
@@ -110,14 +110,6 @@ 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
@@ -135,23 +127,7 @@ def mitchell_best_candidate_modified_3d(
             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()
-
+                in_hull = check_in_hull(concave_hull_bounding_boxes, candidates)
             else:
                 # Test whether all points are inside the convex hull of the perforations
                 if np.all(z == z[0]):

src/flownet/realization/_schedule.py

@@ -162,17 +162,26 @@ def _calculate_welspecs(self):
                 (self._df_production_data["WELL_NAME"] == well_name)
             ]["PHASE"].values[0]
 
-            self.append(
-                WELSPECS(
-                    date=date,
-                    well_name=well_name,
-                    group_name="WELLS",  # Wells directly connected to FIELD gives
-                    # segmentation fault in Flow
-                    i=self.get_compdat(well_name)[0].i,
-                    j=self.get_compdat(well_name)[0].j,
-                    phase=phase,
+            try:
+                self.append(
+                    WELSPECS(
+                        date=date,
+                        well_name=well_name,
+                        group_name="WELLS",  # Wells directly connected to FIELD gives
+                        # segmentation fault in Flow
+                        i=self.get_compdat(well_name)[0].i,
+                        j=self.get_compdat(well_name)[0].j,
+                        phase=phase,
+                    )
                 )
-            )
+            except IndexError:
+                print(
+                    f"""The schedule could not be created for well '{well_name}'.\n
+                This most likely is a result of this well not having any connections.\n
+                Try adding more additional nodes, relax angle constraint for the Delaunay triangles,\n
+                maximum distance and/or convex hull."""
+                )
+                raise
 
     def _calculate_wconhist(self):
         """