Feature: Identify impacting turbines geometrically and reduce FLORIS accordingly

examples_artificial_data/00_setup_floris_model/04_geometric_dependencies.py

@@ -0,0 +1,57 @@
+import matplotlib.pyplot as plt
+import numpy as np
+from floris.tools import FlorisInterface
+
+from flasc.floris_tools import get_all_impacting_turbines_geometrical
+
+# Demonstrate the get_all_impacting_turbines_geometrical
+# function in floris_tools
+
+# Load a large FLORIS object
+fi = FlorisInterface("../floris_input_artificial/gch.yaml")
+D = 126.0
+X, Y = np.meshgrid(7.0 * D * np.arange(20), 5.0 * D * np.arange(20))
+fi.reinitialize(layout_x=X.flatten(), layout_y=Y.flatten())
+
+# Specify which turbines are of interest
+turbine_weights = np.zeros(len(X.flatten()), dtype=float)
+turbine_weights[np.hstack([a + range(10) for a in np.arange(50, 231, 20)])] = 1.0
+
+# Get all impacting turbines for each wind direction using simple geometry rules
+df_impacting = get_all_impacting_turbines_geometrical(
+    fi=fi, turbine_weights=turbine_weights, wd_array=np.arange(0.0, 360.0, 30.0)
+)
+
+# Produce plots showcasing which turbines are estimated to be impacting
+for ii in range(df_impacting.shape[0]):
+    wd = df_impacting.loc[ii, "wd"]
+
+    fig, ax = plt.subplots()
+    ax.plot(fi.layout_x, fi.layout_y, "o", color="lightgray", label="All turbines")
+
+    ids = df_impacting.loc[ii, "impacting_turbines"]
+    no_turbines_total = len(fi.layout_x)
+    no_turbines_reduced = len(ids)
+    ax.plot(fi.layout_x[ids], fi.layout_y[ids], "o", color="black", label="Impacting turbines")
+
+    ids = np.where(turbine_weights > 0.001)[0]
+    ax.plot(fi.layout_x[ids], fi.layout_y[ids], "o", color="red", label="Turbines of interest")
+
+    ax.set_xlabel("X location (m)")
+    ax.set_ylabel("Y location (m)")
+    ax.axis("equal")
+    ax.grid(True)
+    ax.legend()
+    percentage = 100.0 * no_turbines_reduced / no_turbines_total
+    ax.set_title(
+        f"Wind direction: {wd:.1f} deg. Turbines modelled: "
+        f"{no_turbines_reduced:d}/{no_turbines_total} ({percentage:.1f}%)."
+    )
+
+    # Make a statement on number of wake-steered turbines vs. total farm size
+    print(
+        f"wd={wd:.1f} deg. Reduced from {no_turbines_total:d} "
+        f"to {no_turbines_reduced} ({percentage:.1f}%)."
+    )
+
+plt.show()

flasc/floris_tools.py

@@ -91,6 +91,44 @@ def merge_floris_objects(fi_list, reference_wind_height=None):
     return fi_merged
 
 
+def reduce_floris_object(fi, turbine_list, copy=False):
+    """Reduce a large FLORIS object to a subset selection of wind turbines.
+
+    Args:
+        fi (FlorisInterface): FLORIS object.
+        turbine_list (list, array-like): List of turbine indices which should be maintained.
+
+    Returns:
+        fi_reduced (FlorisInterface): The reduced FlorisInterface object.
+    """
+
+    # Copy, if necessary
+    if copy:
+        fi_reduced = fi.copy()
+    else:
+        fi_reduced = fi
+
+    # Get the turbine locations from the floris object
+    x = np.array(fi.layout_x, dtype=float, copy=True)
+    y = np.array(fi.layout_y, dtype=float, copy=True)
+
+    # Get turbine definitions from floris object
+    fi_turbine_type = fi.floris.farm.turbine_type
+    if len(fi_turbine_type) == 1:
+        fi_turbine_type = fi_turbine_type * len(fi.layout_x)
+    elif not len(fi_turbine_type) == len(fi.layout_x):
+        raise UserWarning("Incompatible format of turbine_type in FlorisInterface.")
+
+    # Construct the merged FLORIS model based on the first entry in fi_list
+    fi_reduced.reinitialize(
+        layout_x=x[turbine_list],
+        layout_y=y[turbine_list],
+        turbine_type=list(np.array(fi_turbine_type)[turbine_list]),
+    )
+
+    return fi_reduced
+
+
 def interpolate_floris_from_df_approx(
     df,
     df_approx,
@@ -599,6 +637,109 @@ def get_turbs_in_radius(
     return turbs_within_radius
 
 
+def get_all_impacting_turbines_geometrical(
+    fi, turbine_weights, wd_array=np.arange(0.0, 360.0, 3.0), wake_slope=0.30
+):
+    """Determine which turbines affect the turbines of interest
+    (i.e., those with a turbine_weights > 0.00001). This function
+    uses very simplified geometric functions to very quickly
+    derive which turbines are supposedly waking at least one
+    turbine in the farm of interest.
+
+    Args:
+        fi ([floris object]): FLORIS object of the farm of interest.
+        turbine_weights [list]: List of with turbine weights with length
+        equal to the number of wind turbines, and typically filled with
+        0s (neighbouring farms) and 1s (farm of interest).
+        wd_step (float, optional): Wind direction discretization step.
+        Defaults to 3.0.
+        wake_slope (float, optional): linear slope of the wake (dy/dx)
+        plot_lines (bool, optional): Enable plotting wakes/turbines.
+        Defaults to False.
+
+    Returns:
+        df_impacting_simple ([pd.DataFrame]): A Pandas Dataframe in which each row
+        contains a wind direction and a list of turbine numbers. The turbine
+        numbers are those turbines that should be modelled to accurately
+        capture the wake losses for the wind farm of interest. Turbine numbers
+        that are not in the 'impacting_turbines' can safely be removed from
+        the simulation without affecting any of the turbines that have a nonzero
+        turbine weight.
+    """
+
+    # Get farm layout
+    x = fi.layout_x
+    y = fi.layout_y
+    n_turbs = len(x)
+    D = [t["rotor_diameter"] for t in fi.floris.farm.turbine_definitions]
+    D = np.array(D, dtype=float)
+
+    # Rotate farm and determine freestream/waked turbines
+    is_impacting_list = []
+    for wd in wd_array:
+        is_impacting = [None for _ in range(n_turbs)]
+
+        # Rotate according to freestream wind direction
+        x_rot = np.cos((wd - 270.0) * np.pi / 180.0) * x - np.sin((wd - 270.0) * np.pi / 180.0) * y
+        y_rot = np.sin((wd - 270.0) * np.pi / 180.0) * x + np.cos((wd - 270.0) * np.pi / 180.0) * y
+
+        # Get turbine indices of the farm turbines of interest, and find its most downstream location
+        turb_ids_of_interest = np.where(turbine_weights > 0.0001)[0]
+        x_rot_most_downstream_of_interest = np.max(x_rot[turb_ids_of_interest])
+
+        # Check for each turbine
+        for ii in range(n_turbs):
+            # Check easy skips: turbine is in farm of interest
+            if ii in turb_ids_of_interest:
+                is_impacting[ii] = True
+                continue
+
+            # Check easy skips: further downstream than last turbine
+            if x_rot[ii] >= x_rot_most_downstream_of_interest:
+                is_impacting[ii] = False
+                continue
+
+            x0 = x_rot[ii]
+            y0 = y_rot[ii]
+
+            def yw_upper(x):
+                y = (y0 + D[ii]) + (x - x0) * wake_slope
+                if isinstance(y, (float, np.float64, np.float32)):
+                    if x < (x0 + 0.01):
+                        y = -np.Inf
+                else:
+                    y[x < x0 + 0.01] = -np.Inf
+                return y
+
+            def yw_lower(x):
+                y = (y0 - D[ii]) - (x - x0) * wake_slope
+                if isinstance(y, (float, np.float64, np.float32)):
+                    if x < (x0 + 0.01):
+                        y = -np.Inf
+                else:
+                    y[x < x0 + 0.01] = -np.Inf
+                return y
+
+            def is_in_wake(xt, yt):
+                return (yt < yw_upper(xt)) & (yt > yw_lower(xt))
+
+            is_impacting[ii] = any(
+                is_in_wake(x_rot[turb_ids_of_interest], y_rot[turb_ids_of_interest])
+            )
+
+        is_impacting_list.append(np.where(is_impacting)[0])
+
+    n_turbines_reduced = [len(ids) for ids in is_impacting_list]
+    df_impacting_simple = pd.DataFrame(
+        {
+            "wd": wd_array,
+            "impacting_turbines": is_impacting_list,
+            "n_turbines_reduced": n_turbines_reduced,
+        }
+    )
+    return df_impacting_simple
+
+
 def get_upstream_turbs_floris(fi, wd_step=0.1, wake_slope=0.10, plot_lines=False):
     """Determine which turbines are operating in freestream (unwaked)
     flow, for the entire wind rose. This function will return a data-