Node placement in volume of original model (#372)

* adding keyword place_nodes_in_reservoir_volume to config

* Bugfix hull_factor

* Update test with additional config argument

* [WIP]: Place nodes initially within the boundingbox of the layer/reservoir

* Changes from review

* Some cleaning up and clarification in config

* Error optional float hull_factor

* Added test for places nodes in original volume

* Apply suggestions from code review

Co-authored-by: Wouter J. de Bruin <[email protected]>

* Add docstring to scale concave hull perforations method

* Not re-asigning default values in mitchells and correct naming of convex hull around perforations

* Suggested change olelod hull_factor is 1.0

Co-authored-by: Ole Petter Lødøen <[email protected]>

Co-authored-by: Wouter J. de Bruin <[email protected]>
Co-authored-by: Ole Petter Lødøen <[email protected]>

src/flownet/config_parser/_config_parser.py

@@ -272,7 +272,13 @@ def _to_abs_path(path: Optional[str]) -> str:
                                 "are supported,  e.g. weekly (W), monthly (M), quarterly (Q), "
                                 "yearly (A)",
                             },
-                            "concave_hull": {MK.Type: types.Bool, MK.AllowNone: True},
+                            "concave_hull": {
+                                MK.Type: types.Bool,
+                                MK.AllowNone: True,
+                                MK.Description: "When true, the bounding boxes of the gridcells of the "
+                                "original reservoir model are used to check if the generated additional "
+                                "nodes are positioned within the reservoir volume.",
+                            },
                         },
                     },
                     "constraining": {
@@ -361,7 +367,7 @@ def _to_abs_path(path: Optional[str]) -> str:
                         MK.Type: types.List,
                         MK.Description: "List of additional flow nodes to add "
                         "for each layer or single integer which will be "
-                        "split over the layers, when they are defined.",
+                        "split over the layers, when layers are defined.",
                         MK.LayerTransformation: _integer_to_list,
                         MK.Content: {
                             MK.Item: {
@@ -376,7 +382,18 @@ def _to_abs_path(path: Optional[str]) -> str:
                         MK.Description: "Number of additional nodes to create "
                         "(using Mitchell's best candidate algorithm)",
                     },
-                    "hull_factor": {MK.Type: types.Number, MK.Default: 1.2},
+                    "place_nodes_in_volume_reservoir": {
+                        MK.Type: types.Bool,
+                        MK.AllowNone: True,
+                        MK.Description: "When true use boundary of reservoir/layer volume as "
+                        "bounding volume to place initial candidates instead of convex hull of well perforations.",
+                    },
+                    "hull_factor": {
+                        MK.Type: types.Number,
+                        MK.Default: 1.2,
+                        MK.Description: "Increase the size of the bounding volume around the well "
+                        "perforations to place additional nodes in.",
+                    },
                     "random_seed": {
                         MK.Type: types.Number,
                         MK.AllowNone: True,
@@ -1676,6 +1693,15 @@ def parse_config(
             "there is only a single layer."
         )
 
+    if (
+        config.flownet.place_nodes_in_volume_reservoir
+        and not config.flownet.data_source.concave_hull
+    ):
+        raise ValueError(
+            "concave_hull needs to be true for flownet to be able to "
+            "place candidates within the reservoir volume."
+        )
+
     req_relp_parameters: List[str] = []
     if (
         config.model_parameters.equil.scheme != "regions_from_sim"

src/flownet/network_model/_generate_connections.py

@@ -178,6 +178,7 @@ def _generate_connections(
             num_added_flow_nodes=additional_flow_nodes,
             num_candidates=configuration.flownet.additional_node_candidates,
             hull_factor=configuration.flownet.hull_factor,
+            place_nodes_in_volume_reservoir=configuration.flownet.place_nodes_in_volume_reservoir,
             concave_hull_bounding_boxes=concave_hull_bounding_boxes,
             random_seed=configuration.flownet.random_seed,
         )

src/flownet/network_model/_mitchell.py

@@ -1,4 +1,4 @@
-from typing import List, Optional
+from typing import List, Optional, Tuple
 import time
 
 from scipy.spatial import Delaunay  # pylint: disable=no-name-in-module
@@ -11,8 +11,9 @@
 def mitchell_best_candidate_modified_3d(
     perforations: List[Coordinate],
     num_added_flow_nodes: int,
-    num_candidates: int = 1000,
-    hull_factor: Optional[float] = None,
+    num_candidates: int,
+    hull_factor: float,
+    place_nodes_in_volume_reservoir: Optional[bool] = None,
     concave_hull_bounding_boxes: Optional[np.ndarray] = None,
     random_seed: Optional[int] = None,
 ) -> List[Coordinate]:
@@ -33,9 +34,11 @@ def mitchell_best_candidate_modified_3d(
             [(xr_1, yr_1, zr_1), ..., (xr_N, yr_N, zr_N)]
         num_added_flow_nodes: Number of additional flow nodes to generate
         num_candidates: Number of candidates to consider per additional flow nodes
+        place_nodes_in_volume_reservoir: When true, additional nodes will initially be placed inside
+            the bounding box of the reservoir or layer instead of the bounding box of the well perforations.
         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.
+            Default defined in config parser is 1.2.
         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.
 
@@ -55,45 +58,48 @@ def mitchell_best_candidate_modified_3d(
     # Number of real wells
     num_points = len(x)
 
-    x_moved = np.zeros(num_points)
-    y_moved = np.zeros(num_points)
-    z_moved = np.zeros(num_points)
-
-    # Determine whether the complex hull needs to be scaled
-    if hull_factor:
-        # Calculate the centroid of all real perforations
-        centroid = (sum(x) / num_points, sum(y) / num_points, sum(z) / num_points)
-
-        # Loop through all real perforation locations
-        for count, point in enumerate(perforations):
-            # Linearly scale the well perforation's location relative to the centroid
-            moved_points = np.add(hull_factor * np.subtract(point, centroid), centroid)
-            x_moved[count] = moved_points[0]
-            y_moved[count] = moved_points[1]
-            z_moved[count] = moved_points[2]
-
-    if hull_factor:
-        x_min = min(x_moved)
-        x_max = max(x_moved)
-        y_min = min(y_moved)
-        y_max = max(y_moved)
-        z_min = min(z_moved)
-        z_max = max(z_moved)
+    # Bounding box to place initial candidates in: reservoir volume or (scaled) convex hull of real perforations.
+    if place_nodes_in_volume_reservoir:
+        x_mins, x_maxs, y_mins, y_maxs, z_mins, z_maxs = np.hsplit(
+            concave_hull_bounding_boxes, 6
+        )
+
+        x_min = min(x_mins)[0]
+        x_max = max(x_maxs)[0]
+        y_min = min(y_mins)[0]
+        y_max = max(y_maxs)[0]
+        z_min = min(z_mins)[0]
+        z_max = max(z_maxs)[0]
     else:
-        x_min = min(x)
-        x_max = max(x)
-        y_min = min(y)
-        y_max = max(y)
-        z_min = min(z)
-        z_max = max(z)
-
-    # Determine the convex hull of the original or linearly scaled perforations
-    if np.all(z == z[0]):
-        # 2D cases
-        hull = Delaunay(np.column_stack([x, y]))
-    else:
-        # 3D cases
-        hull = Delaunay(np.column_stack([x, y, z]))
+        # Determine whether the convex hull needs to be scaled
+        if not np.isclose(hull_factor, 1.0):
+            x_hull, y_hull, z_hull = scale_convex_hull_perforations(
+                perforations, hull_factor
+            )
+            x_min = min(x_hull)
+            x_max = max(x_hull)
+            y_min = min(y_hull)
+            y_max = max(y_hull)
+            z_min = min(z_hull)
+            z_max = max(z_hull)
+        else:
+            x_min = min(x)
+            x_max = max(x)
+            y_min = min(y)
+            y_max = max(y)
+            z_min = min(z)
+            z_max = max(z)
+            x_hull = x
+            y_hull = y
+            z_hull = z
+
+        # Determine the convex hull of the original or linearly scaled perforations
+        if np.all(z == z[0]):
+            # 2D cases
+            perforation_hull = Delaunay(np.column_stack([x_hull, y_hull]))
+        else:
+            # 3D cases
+            perforation_hull = Delaunay(np.column_stack([x_hull, y_hull, z_hull]))
 
     # Generate all new flow nodes
     for i in range(num_points, num_points + num_added_flow_nodes):
@@ -127,13 +133,15 @@ def mitchell_best_candidate_modified_3d(
             candidates = np.vstack([x_candidate, y_candidate, z_candidate]).T
 
             if concave_hull_bounding_boxes is not None:
+                # Test whether all points are inside the bounding boxes of the gridcells of the reservoir volume.
+                # This is the always the case when place_nodes_in_volume_reservoir is True.
                 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]):
-                    in_hull = hull.find_simplex(candidates[:, (0, 1)]) >= 0
+                    in_hull = perforation_hull.find_simplex(candidates[:, (0, 1)]) >= 0
                 else:
-                    in_hull = hull.find_simplex(candidates) >= 0
+                    in_hull = perforation_hull.find_simplex(candidates) >= 0
 
         best_distance = 0
         best_candidate = 0
@@ -175,3 +183,42 @@ def mitchell_best_candidate_modified_3d(
 
     # Return the real/original and added flow node coordinates as a list of tuples.
     return [(x[i], y[i], z[i]) for i in range(len(x))]
+
+
+def scale_convex_hull_perforations(
+    perforations: List[Coordinate], hull_factor: float
+) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
+    """
+    Linear scaling of the perforation points based on the hull_factor. Factor will
+    scale the distance of each point from the centroid of all the points.
+    These scaled points are used to create a convex hull in which additional flow nodes will be placed.
+
+    Args:
+        perforations: Python list of real well coordinate tuples
+            [(xr_1, yr_1, zr_1), ..., (xr_N, yr_N, zr_N)]
+        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.
+    Returns:
+        The tuple consisting of numpy arrays of x,y,z moved points based on the hull_factor,
+        which are used further on in the code to create a convex hull around the real wells (perforations).
+
+    """
+    x, y, z = (np.asarray(t) for t in zip(*perforations))
+    num_points = len(x)
+
+    x_hull = np.zeros(num_points)
+    y_hull = np.zeros(num_points)
+    z_hull = np.zeros(num_points)
+
+    # Calculate the centroid of all real perforations
+    centroid = (sum(x) / num_points, sum(y) / num_points, sum(z) / num_points)
+
+    # Loop through all real perforation locations
+    for count, point in enumerate(perforations):
+        # Linearly scale the well perforation's location relative to the centroid
+        moved_points = np.add(hull_factor * np.subtract(point, centroid), centroid)
+        x_hull[count] = moved_points[0]
+        y_hull[count] = moved_points[1]
+        z_hull[count] = moved_points[2]
+
+    return x_hull, y_hull, z_hull

tests/test_generate_connections.py

@@ -101,6 +101,7 @@ def test_generate_connections() -> None:
     config.flownet = collections.namedtuple("flownet", "additional_flow_nodes")
     config.flownet.additional_flow_nodes = 2
     config.flownet.additional_node_candidates = 2
+    config.flownet.place_nodes_in_volume_reservoir = None
     config.flownet.hull_factor = 1
     config.flownet.random_seed = 1
     config.flownet.angle_threshold = None

tests/test_mitchell.py

@@ -1,8 +1,13 @@
 from typing import List
 
 import numpy as np
+from scipy.spatial import Delaunay  # pylint: disable=no-name-in-module
 
-from flownet.network_model._mitchell import mitchell_best_candidate_modified_3d
+from flownet.network_model._mitchell import (
+    mitchell_best_candidate_modified_3d,
+    scale_convex_hull_perforations,
+)
+from flownet.network_model._hull import check_in_hull
 from src.flownet.utils.types import Coordinate
 
 
@@ -29,7 +34,7 @@ def test_mitchells_3d() -> None:
             well_perforations_3d,
             num_added_flow_nodes=num_added_flow_nodes,
             num_candidates=100,
-            hull_factor=1.2,
+            hull_factor=1.0,
             concave_hull_bounding_boxes=None,
             random_seed=999,
         )
@@ -69,9 +74,9 @@ def test_mitchells_2d() -> None:
         tuple(elem)
         for elem in mitchell_best_candidate_modified_3d(
             well_perforations_2d,
-            num_added_flow_nodes=10,
+            num_added_flow_nodes=num_added_flow_nodes,
             num_candidates=100,
-            hull_factor=1.2,
+            hull_factor=1.0,
             concave_hull_bounding_boxes=None,
             random_seed=999,
         )
@@ -87,3 +92,88 @@ def test_mitchells_2d() -> None:
     assert np.all([y[1] >= min(y_wells) for y in coordinates])
     assert np.all([y[1] <= max(y_wells) for y in coordinates])
     assert np.all([z[2] == z_wells[0] for z in coordinates])
+
+
+def test_hull_factor_mitchell() -> None:
+    well_perforations_2d = [
+        [36.0, 452.0, 4002.0],
+        [236.0, 420.0, 4002.0],
+        [12.0, 276.0, 4002.0],
+        [212.0, 228.0, 4002.0],
+        [396.0, 276.0, 4002.0],
+        [60.0, 68.0, 4002.0],
+        [252.0, 12.0, 4002.0],
+        [452.0, 44.0, 4002.0],
+        [124.0, 340.0, 4002.0],
+        [276.0, 316.0, 4002.0],
+        [180.0, 124.0, 4002.0],
+        [340.0, 140.0, 4002.0],
+    ]
+
+    num_added_flow_nodes = 10
+    hull_factor = 2.0
+
+    coordinates: List[Coordinate] = [
+        tuple(elem)
+        for elem in mitchell_best_candidate_modified_3d(
+            well_perforations_2d,
+            num_added_flow_nodes=num_added_flow_nodes,
+            hull_factor=hull_factor,
+            num_candidates=100,
+            concave_hull_bounding_boxes=None,
+            random_seed=999,
+        )
+    ]
+
+    x_hull, y_hull, z_hull = scale_convex_hull_perforations(
+        well_perforations_2d, hull_factor
+    )
+
+    assert len(coordinates) == (len(well_perforations_2d) + num_added_flow_nodes)
+    assert np.all([x[0] >= min(x_hull) for x in coordinates])
+    assert np.all([x[0] <= max(x_hull) for x in coordinates])
+    assert np.all([y[1] >= min(y_hull) for y in coordinates])
+    assert np.all([y[1] <= max(y_hull) for y in coordinates])
+    assert np.all([z[2] == z_hull[0] for z in coordinates])
+
+
+def test_nodes_in_reservoir_volume_mitchells() -> None:
+    well_perforations_3d = [
+        [0.5, 0.5, 1.0],
+        [0.5, 2.5, 0.5],
+        [1.5, 5.5, 1.5],
+        [3.5, 0.5, 1.0],
+    ]
+
+    concave_hull_bounding_boxes = np.array(
+        [
+            [0.0, 2.0, 0.0, 2.0, 0.0, 2.0],
+            [2.0, 3.0, 0.0, 2.0, 0.0, 2.0],
+            [3.0, 5.0, 0.0, 2.0, 0.0, 2.0],
+            [0.0, 2.0, 2.0, 4.0, 0.0, 2.0],
+            [0.0, 2.0, 4.0, 6.0, 0.0, 2.0],
+            [2.0, 5.0, 4.0, 6.0, 0.0, 2.0],
+            [5.0, 7.0, 4.0, 6.0, 0.0, 2.0],
+            [5.0, 20.0, 6.0, 20.0, 0.0, 2.0],
+        ]
+    )
+
+    coordinates: List[Coordinate] = [
+        tuple(elem)
+        for elem in mitchell_best_candidate_modified_3d(
+            well_perforations_3d,
+            num_added_flow_nodes=20,
+            num_candidates=500,
+            hull_factor=1.0,
+            place_nodes_in_volume_reservoir=True,
+            concave_hull_bounding_boxes=concave_hull_bounding_boxes,
+            random_seed=999,
+        )
+    ]
+
+    perforation_hull = Delaunay(np.array(well_perforations_3d))
+    in_hull_perforations = perforation_hull.find_simplex(np.array(coordinates))
+
+    in_hull_volume = check_in_hull(concave_hull_bounding_boxes, np.array(coordinates))
+    assert in_hull_volume.all()
+    assert any(x == -1 for x in in_hull_perforations)