diff --git a/examples_artificial_data/00_setup_floris_model/04_geometric_dependencies.py b/examples_artificial_data/00_setup_floris_model/04_geometric_dependencies.py new file mode 100644 index 00000000..871dd5d8 --- /dev/null +++ b/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() diff --git a/flasc/floris_tools.py b/flasc/floris_tools.py index 6fcc5405..9b3d9bda 100644 --- a/flasc/floris_tools.py +++ b/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-