From 97932d65c8c7cef27ce31e849621ecc5a3805941 Mon Sep 17 00:00:00 2001 From: Paul Date: Wed, 13 Mar 2024 09:31:36 -0600 Subject: [PATCH] Add example testing TI in emg --- examples/xx_add_ti_to_emg.py | 243 +++++++++++++++++++++++++++++++++++ 1 file changed, 243 insertions(+) create mode 100644 examples/xx_add_ti_to_emg.py diff --git a/examples/xx_add_ti_to_emg.py b/examples/xx_add_ti_to_emg.py new file mode 100644 index 000000000..a79a65bac --- /dev/null +++ b/examples/xx_add_ti_to_emg.py @@ -0,0 +1,243 @@ +""" Script for reconsidering TI recovery in the EMG model + + +""" + + +from typing import ( + Any, + Dict, + List, + Optional, +) + +import matplotlib.pyplot as plt +import numpy as np +import pandas as pd + +from floris import FlorisModel, WindRose +from floris.layout_visualization import ( + plot_turbine_labels, + plot_turbine_points, + plot_waking_directions, +) + + +# Tuning parameters for new ti model +atmospheric_ti_gain = 0.0 # note default: 0.0 +wake_expansion_rates = [0.023, 0.008] # note default: [0.023, 0.008] +breakpoints_D = [10] # Node default: [10] + +# Layout parameters +n_col = 6 # Number of columns +n_t = 10 # Number of turbines per column +dist_c = 7.0 # Distance between columns +dist_t = 3.0 # Distance between turbines + +# Atmospheric parameters +wind_directions = np.arange(250, 291, 1.0) +wind_speeds = np.array([8.0, 9.0]) + +# Parameters +D = 126.0 + + +def nested_get(dic: Dict[str, Any], keys: List[str]) -> Any: + """Get a value from a nested dictionary using a list of keys. + Based on: stackoverflow.com/questions/14692690/access-nested-dictionary-items-via-a-list-of-keys + + Args: + dic (Dict[str, Any]): The dictionary to get the value from. + keys (List[str]): A list of keys to traverse the dictionary. + + Returns: + Any: The value at the end of the key traversal. + """ + for key in keys: + dic = dic[key] + return dic + + +def nested_set(dic: Dict[str, Any], keys: List[str], value: Any, idx: Optional[int] = None) -> None: + """Set a value in a nested dictionary using a list of keys. + Based on: stackoverflow.com/questions/14692690/access-nested-dictionary-items-via-a-list-of-keys + + Args: + dic (Dict[str, Any]): The dictionary to set the value in. + keys (List[str]): A list of keys to traverse the dictionary. + value (Any): The value to set. + idx (Optional[int], optional): If the value is an list, the index to change. + Defaults to None. + """ + dic_in = dic.copy() + + for key in keys[:-1]: + dic = dic.setdefault(key, {}) + if idx is None: + # Parameter is a scaler, set directly + dic[keys[-1]] = value + else: + # Parameter is a list, need to first get the list, change the values at idx + + # # Get the underlying list + par_list = nested_get(dic_in, keys) + par_list[idx] = value + dic[keys[-1]] = par_list + + +def set_fi_param( + fm_in: FlorisModel, param: List[str], value: Any, param_idx: Optional[int] = None +) -> FlorisModel: + """Set a parameter in a FlorisInterface object. + + Args: + fi_in (FlorisInterface): The FlorisInterface object to modify. + param (List[str]): A list of keys to traverse the FlorisInterface dictionary. + value (Any): The value to set. + idx (Optional[int], optional): The index to set the value at. Defaults to None. + + Returns: + FlorisInterface: The modified FlorisInterface object. + """ + fm_dict_mod = fm_in.core.as_dict() + nested_set(fm_dict_mod, param, value, param_idx) + return FlorisModel(fm_dict_mod) + + +# Initialize FLORIS with the given input file. +# The Floris class is the entry point for most usage. +fmodel_emg = FlorisModel("inputs/emgauss.yaml") + +# Get a copy for the ti-enabled case +fmodel_emg_ti = fmodel_emg.copy() + +# Get the dictionary +# fm_dict = fmodel_emg.core.as_dict() +# print(fm_dict['wake']['wake_velocity_parameters']['empirical_gauss'].keys()) +# Update the parameters +fmodel_emg_ti = set_fi_param( + fm_in=fmodel_emg_ti, + param=["wake", "wake_velocity_parameters", "empirical_gauss", "wake_expansion_rates"], + value=wake_expansion_rates[0], + param_idx=0, +) +fmodel_emg_ti = set_fi_param( + fm_in=fmodel_emg_ti, + param=["wake", "wake_velocity_parameters", "empirical_gauss", "wake_expansion_rates"], + value=wake_expansion_rates[1], + param_idx=1, +) +fmodel_emg_ti = set_fi_param( + fm_in=fmodel_emg_ti, + param=["wake", "wake_velocity_parameters", "empirical_gauss", "breakpoints_D"], + value=breakpoints_D[0], + param_idx=0, +) +fmodel_emg_ti = set_fi_param( + fm_in=fmodel_emg_ti, + param=["wake", "wake_turbulence_parameters", "wake_induced_mixing", "atmospheric_ti_gain"], + value=atmospheric_ti_gain, +) + + +# Use a nested list comprehension to create the layout_x and layout_y arrays +layout_x = [dist_c * D * i for i in range(n_col) for j in range(n_t)] +layout_y = [dist_t * D * j for i in range(n_col) for j in range(n_t)] + +# Set the wind farm layout using the set method +fmodel_emg.set(layout_x=layout_x, layout_y=layout_y) +fmodel_emg_ti.set(layout_x=layout_x, layout_y=layout_y) + +# Show the layout of the farm +fig, ax = plt.subplots() +plot_turbine_points(fmodel_emg, ax) +plot_turbine_labels(fmodel_emg, ax) +plot_waking_directions(fmodel_emg, ax, limit_dist_D=dist_c * 1.01) + +# Set up a wind rose input +wind_rose = WindRose( + wind_directions=wind_directions, + wind_speeds=wind_speeds, + ti_table=0.06, +) +fmodel_emg.set(wind_data=wind_rose) +fmodel_emg_ti.set(wind_data=wind_rose) + +# Run the FLORIS model +fmodel_emg.run() +fmodel_emg_ti.run() + +# Get the power from the turbines +turbine_power_emg = fmodel_emg.get_turbine_powers() +turbine_power_emg_ti = fmodel_emg_ti.get_turbine_powers() + +# Currently the matrices have n_t * n_col columns, each corresponding to a turbine +# average the columns according to the column, this means averaging every n_t columns +turbine_power_emg = turbine_power_emg.reshape(turbine_power_emg.shape[0], -1, n_t) +turbine_power_emg = np.mean(turbine_power_emg, axis=2) +turbine_power_emg_ti = turbine_power_emg_ti.reshape(turbine_power_emg_ti.shape[0], -1, n_t) +turbine_power_emg_ti = np.mean(turbine_power_emg_ti, axis=2) + + +# Now determine the gridded wind speeds and directions +wd_grid, ws_grid = np.meshgrid(wind_directions, wind_speeds, indexing="ij") + +wd_flat = wd_grid.flatten() +ws_flat = ws_grid.flatten() + +# Now build up a dataframe +df_power_emg = pd.DataFrame( + data=turbine_power_emg, +) + +# Add additional columns +df_power_emg["wind_speed"] = ws_flat +df_power_emg["wind_direction"] = wd_flat +df_power_emg["model"] = "emg" + +# Repeat for the ti-enabled case +df_power_emg_ti = pd.DataFrame( + data=turbine_power_emg_ti, +) + +# Add additional columns +df_power_emg_ti["wind_speed"] = ws_flat +df_power_emg_ti["wind_direction"] = wd_flat +df_power_emg_ti["model"] = "emg_ti" + +# Concatenate the two dataframes +df = pd.concat([df_power_emg, df_power_emg_ti], axis=0) + +# Now make a figure comparing the two models +num_ws = len(wind_speeds) +fig, axarr = plt.subplots(num_ws, n_col - 1, figsize=(10, 10), sharex=True, sharey=True) + +for i, ws in enumerate(wind_speeds): + for j in range(n_col - 1): + ax = axarr[i, j] + # Get the data for the current wind speed + df_ws = df[df["wind_speed"] == ws] + + # Get the data for the current column + df_ws["ratio"] = df_ws[j + 1] / df_ws[0] + + # Now compare the ratios against wind direction for the emg and emg_ti models + df_ws_emg = df_ws[df_ws["model"] == "emg"] + df_ws_emg_ti = df_ws[df_ws["model"] == "emg_ti"] + + ax.plot(df_ws_emg["wind_direction"], df_ws_emg["ratio"], label="emg") + ax.plot(df_ws_emg_ti["wind_direction"], df_ws_emg_ti["ratio"], label="emg_ti") + + ax.set_xlabel("Wind direction") + ax.grid(True) + + if j == 0: + ax.set_ylabel(f"Power ratio for wind speed {ws}") + + if i == 0: + ax.set_title(f"Column {j + 1} of {n_col}") + + if i == 0 and j == 0: + ax.legend() + +plt.show()