diff --git a/docs/references.bib b/docs/references.bib index da217dd4f..a58dbed39 100644 --- a/docs/references.bib +++ b/docs/references.bib @@ -273,6 +273,19 @@ @Article{bay_2022 DOI = {10.5194/wes-2022-17} } +@article{gunn2019evmodel, + title = {Improvements to the Eddy Viscosity Wind Turbine Wake Model}, + volume = {1222}, + issn = {1742-6596}, + url = {https://dx.doi.org/10.1088/1742-6596/1222/1/012003}, + doi = {10.1088/1742-6596/1222/1/012003}, + pages = {012003}, + number = {1}, + journaltitle = {Journal of Physics: Conference Series}, + author = {Gunn, K.}, + note = {Publisher: {IOP} Publishing}, +} + @article{Pedersen_2022_turbopark2, url = {https://dx.doi.org/10.1088/1742-6596/2265/2/022063}, year = {2022}, @@ -299,3 +312,16 @@ @article{SinnerFleming2024grs title = {Robust wind farm layout optimization}, journal = {Journal of Physics: Conference Series}, } + +@proceedings{Kuo2014soed, + author = {Kuo, Jim Y. J. and Romero, David A. and Amon, Cristina H.}, + title = "{A Novel Wake Interaction Model for Wind Farm Layout Optimization}", + volume = {Volume 6B: Energy}, + series = {ASME International Mechanical Engineering Congress and Exposition}, + pages = {V06BT07A074}, + year = {2014}, + month = {11}, + doi = {10.1115/IMECE2014-39073}, + url = {https://doi.org/10.1115/IMECE2014-39073}, + eprint = {https://asmedigitalcollection.asme.org/IMECE/proceedings-pdf/IMECE2014/46521/V06BT07A074/4267418/v06bt07a074-imece2014-39073.pdf}, +} diff --git a/docs/wake_models.ipynb b/docs/wake_models.ipynb index 78919191a..eacbd57b7 100644 --- a/docs/wake_models.ipynb +++ b/docs/wake_models.ipynb @@ -202,6 +202,26 @@ "model_plot(\"../examples/inputs/turboparkgauss_cubature.yaml\", include_wake_deflection=False)" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Eddy Viscosity\n", + "\n", + "The eddy viscosity model is a wake model originally presented by {cite:t}`ainslie1988calculating` and more recently advanced by {cite:t}`gunn2019evmodel`. The definition of the wake model requires solving an ordinary differential equation (ODE) to determine the wake centerline deficit at various downstream distances. Following {cite:t}`gunn2019evmodel` and others, the FLORIS implementation makes the self-similarity assumption to simplify the wake centerline ODE. The wake profile is Gaussian.\n", + "\n", + "Further, the wake width for individual turbine wakes is \"expanded\" upon encountering another turbine, and the eddy viscosity model is designed to be run with the \"sum of energy deficit\" (SOED) wake combination model." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "model_plot(\"../examples/inputs/ev.yaml\", include_wake_deflection=False)" + ] + }, { "attachments": {}, "cell_type": "markdown", @@ -360,6 +380,15 @@ "combination_plot(\"sosfs\")" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Sum of Energy Deficit (SOED)\n", + "\n", + "This model combines the wakes by conserving momentum, as described in {cite:t}`Kuo2014soed`." + ] + }, { "attachments": {}, "cell_type": "markdown", diff --git a/examples/examples_eddy_viscosity/001_demonstrate_eddy_viscosity_model.py b/examples/examples_eddy_viscosity/001_demonstrate_eddy_viscosity_model.py new file mode 100644 index 000000000..c2f71742d --- /dev/null +++ b/examples/examples_eddy_viscosity/001_demonstrate_eddy_viscosity_model.py @@ -0,0 +1,124 @@ +"""Example: Eddy viscosity model +This example demonstrates the wake profile of the eddy viscosity wake model, +presented originally by Ainslie (1988) and updated by Gunn (2019). +Links: +- Ainslie (1988): https://doi.org/10.1016/0167-6105(88)90037-2 +- Gunn (2019): https://dx.doi.org/10.1088/1742-6596/1222/1/012003 +""" + +import matplotlib.pyplot as plt +import numpy as np + +import floris.flow_visualization as flowviz +import floris.layout_visualization as layoutviz +from floris import FlorisModel + + +# Instantiate FLORIS model +fmodel = FlorisModel("../inputs/ev.yaml") + +## Plot the flow velocity profiles +# Set up a two-turbine farm +D = 126 +fmodel.set(layout_x=[0, 500, 1000], layout_y=[0, 0, 0]) + +fig, ax = plt.subplots(1, 1) +fig.set_size_inches(10, 4) +ax.scatter(fmodel.layout_x, fmodel.layout_y, color="red", label="Turbine location") + +# Set a single wind condition of westerly wind +wd_array = np.array([270]) +ws_array = 8.0 * np.ones_like(wd_array) +ti_array = 0.06 * np.ones_like(wd_array) +fmodel.set(wind_directions=wd_array, wind_speeds=ws_array, turbulence_intensities=ti_array) + +# Create points to sample the flow in the turbines' wakes, both at the centerline and +# across the wake's width +x_locs = np.linspace(0, 2000, 400) +y_locs = np.linspace(-100, 100, 41) +points_x, points_y = np.meshgrid(x_locs, y_locs) +points_x = points_x.flatten() +points_y = points_y.flatten() +points_z = 90 * np.ones_like(points_x) + +# Collect the points +u_at_points = fmodel.sample_flow_at_points(points_x, points_y, points_z) +u_at_points = u_at_points.reshape((len(y_locs), len(x_locs), 1)) + +# Plot the flow velocities +for y_idx, y in enumerate(y_locs): + a = 1-np.abs(y/100) + if y_idx == (len(y_locs)-1)/2: + label = r"$r = 0D$ (centerline)" + elif y_idx == 3*(len(y_locs)-1)/4: + label = r"$r \approx 0.4D$" + elif y_idx == len(y_locs)-1: + label = r"$r \approx 0.8D$" + else: + label=None + ax.plot(x_locs, u_at_points[y_idx, :, 0].flatten(), color="black", alpha=a, label=label) +ax.grid() +ax.legend() +ax.set_xlabel("x location [m]") +ax.set_ylabel("Wind Speed [m/s]") + +## Visualize the flow in aligned and slightly misaligned conditions for a 9-turbine farm +fig, ax = plt.subplots(2,2) +fig.set_size_inches(10, 10) +D = 126.0 +x_locs = np.array([0, 5*D, 10*D]) +y_locs = x_locs +points_x, points_y = np.meshgrid(x_locs, y_locs) +fmodel.set( + layout_x = points_x.flatten(), + layout_y = points_y.flatten() +) + +# Aligned +layoutviz.plot_turbine_rotors(fmodel, ax=ax[0,0]) +layoutviz.plot_turbine_labels(fmodel, ax=ax[0,0]) + +fmodel.set( + wind_speeds=[8.0, 8.0], + wind_directions=[270.0, 270.0+15.0], + turbulence_intensities=[0.06, 0.06] +) +fmodel.run() +cut_plane = fmodel.calculate_horizontal_plane(height=90, findex_for_viz=0) + +flowviz.visualize_cut_plane(cut_plane, ax=ax[0,0]) +ax[0,0].set_title("Aligned flow") + +# Plot and print the power output of the turbines in each row +np.set_printoptions(formatter={"float": "{0:0.3f}".format}) +print("Aligned case:") +for i in range(3): + idxs = [3*i, 3*i+1, 3*i+2] + ax[1,0].scatter([0,1,2], fmodel.get_turbine_powers()[0, idxs]/1e6, label="Column {0}".format(i)) + print(idxs, " -- ", fmodel.get_turbine_powers()[0, idxs]/1e6) +ax[1,0].grid() +ax[1,0].legend() +ax[1,0].set_xlabel("Turbine in column") +ax[1,0].set_ylabel("Power [MW]") +ax[1,0].set_ylim([0.5, 1.8]) + +# Misaligned +layoutviz.plot_turbine_rotors(fmodel, yaw_angles=(-15.0)*np.ones(9), ax=ax[0,1]) +layoutviz.plot_turbine_labels(fmodel, ax[0,1]) + +cut_plane = fmodel.calculate_horizontal_plane(height=90, findex_for_viz=1) +flowviz.visualize_cut_plane(cut_plane, ax=ax[0,1]) +ax[0,1].set_title("Misaligned flow") + +print("\nMisaligned case:") +for i in range(3): + idxs = [3*i, 3*i+1, 3*i+2] + ax[1,1].scatter([0,1,2], fmodel.get_turbine_powers()[1, idxs]/1e6, label="Column {0}".format(i)) + print(idxs, " -- ", fmodel.get_turbine_powers()[1, idxs]/1e6) +ax[1,1].grid() +ax[1,1].legend() +ax[1,1].set_xlabel("Turbine in column") +ax[1,1].set_ylabel("Power [MW]") +ax[1,1].set_ylim([0.5, 1.8]) + +plt.show() diff --git a/examples/examples_eddy_viscosity/002_reproduce_published_results.py b/examples/examples_eddy_viscosity/002_reproduce_published_results.py new file mode 100644 index 000000000..00095385f --- /dev/null +++ b/examples/examples_eddy_viscosity/002_reproduce_published_results.py @@ -0,0 +1,153 @@ +"""Example: Reproduce published eddy viscosity results +This example attempts to reproduce the results of Ainslie (1988) and Gunn (2019) +using the FLORIS implementation of the eddy viscosity model. + +Links: +- Ainslie (1988): https://doi.org/10.1016/0167-6105(88)90037-2 +- Gunn (2019): https://dx.doi.org/10.1088/1742-6596/1222/1/012003 +""" + +import matplotlib.pyplot as plt +import numpy as np + +from floris import FlorisModel +from floris.turbine_library import build_cosine_loss_turbine_dict + + +# Build a constant CT turbine model for use in comparisons (not realistic) +u_0 = 8.0 # wind speed [m/s] + +# Load the EV model +fmodel = FlorisModel("../inputs/ev.yaml") + +## First, attempt to reproduce Figs. 3 and 4 from Ainslie (1988). + +# It is not clear exactly how the parametrization of Gunn, which is the one +# that is implemented in the EddyViscosityVelocity model, can be matched to +# the parameters given by Ainslie. Rather than try to do this, we simply +# generate plots similar to Figs. 3 and 4. using the default parameters following +# Gunn (2019). + +# Generate figure to plot on +fig, ax = plt.subplots(1, 2, sharex=True, sharey=True) +fig.set_size_inches(8, 4) + +HH = 50 +D = 50 + +for i, wd_std_ev in enumerate([0.0, np.rad2deg(0.1)]): + fmodel.set_param(["wake", "wake_velocity_parameters", "eddy_viscosity", "wd_std_ev"], wd_std_ev) + for C_T, ls in zip([0.5, 0.7, 0.9], ["-", "--", ":"]): + const_CT_turb = build_cosine_loss_turbine_dict( + turbine_data_dict={ + "wind_speed":[0.0, 30.0], + "power":[0.0, 1.0], # Not realistic but won't be used here + "thrust_coefficient":[C_T, C_T] + }, + turbine_name="ConstantCT", + rotor_diameter=D, + hub_height=HH, + ref_tilt=0.0, + ) + + # Load the EV model and set a constant CT turbine + fmodel.set( + layout_x=[0], + layout_y=[0], + turbine_type=[const_CT_turb], + wind_speeds=[u_0], + wind_directions=[270], + turbulence_intensities=[0.14], # As specified by Ainslie + wind_shear=0.0, + reference_wind_height=HH + ) + + points_x = np.linspace(2*D, 10*D, 1000) + points_y = np.zeros_like(points_x) + points_z = HH * np.ones_like(points_x) + + u_at_points = fmodel.sample_flow_at_points(points_x, points_y, points_z) + + # Plot results (not different axis scales in Ainslie) + ax[i].plot( + points_x/D, 1-u_at_points[0, :]/u_0, color="k", linestyle=ls, label=rf"$C_T$ = {C_T}" + ) + ax[i].set_title(r"WD std. dev. $\sigma_{\theta} = $"+"{:.1f} deg".format(wd_std_ev)) + +ax[0].set_xlabel("Downstream distance [D]") +ax[1].set_xlabel("Downstream distance [D]") +ax[0].set_ylabel("Centerline velocity deficit [-]") +ax[0].set_ylim([0, 1]) +ax[0].legend() +ax[0].grid() +ax[1].legend() +ax[1].grid() + + +## Second, reproduce Figure 1b (right) from Gunn (2019) + +# Reset to the default model parameters +fmodel = FlorisModel("../inputs/ev.yaml") + +# Adjustments +C_T2 = 0.34 +D2 = 1.52 + +y_offset = -20 +HH = 0.8 + +# Match depends on ambient turbulence intensity. 7.5% appears close. Also, switch off meandering. +TI = 0.075 +fmodel.set_param(["wake", "wake_velocity_parameters", "eddy_viscosity", "wd_std_ev"], 0.0) + +turb_1 = build_cosine_loss_turbine_dict( + turbine_data_dict={ + "wind_speed":[0.0, 30.0], + "power":[0.0, 1.0], # Not realistic but won't be used here + "thrust_coefficient":[0.8, 0.8] + }, + turbine_name="ConstantCT1", + rotor_diameter=1, + hub_height=HH, + ref_tilt=0.0, +) +turb_2 = build_cosine_loss_turbine_dict( + turbine_data_dict={ + "wind_speed":[0.0, 30.0], + "power":[0.0, 1.0], # Not realistic but won't be used here + "thrust_coefficient":[C_T2, C_T2] + }, + turbine_name="ConstantCT2", + rotor_diameter=D2, + hub_height=HH, + ref_tilt=0.0, +) + +fmodel.set( + layout_x=[0, 1.5], + layout_y=[0, y_offset], + turbine_type=[turb_1, turb_2], + wind_speeds=[u_0], + wind_directions=[270.0], + turbulence_intensities=[TI], + wind_shear=0.0, + reference_wind_height=HH +) + +n_pts = 1000 +points_x = np.concatenate((np.linspace(0, 10, n_pts), np.linspace(1.5, 11.5, n_pts))) +points_y = np.concatenate((np.zeros(n_pts), y_offset*np.ones(n_pts))) +points_z = HH * np.ones_like(points_x) +u_at_points = fmodel.sample_flow_at_points(points_x, points_y, points_z) + +fig, ax = plt.subplots(1,1) +ax.plot([0.0, 0.0], [0.3, 1.1], color="black") +ax.plot([1.5, 1.5], [0.3, 1.1], color="black") +ax.plot(points_x[:n_pts], u_at_points[0, :n_pts]/u_0, color="C0") +ax.plot(points_x[n_pts:], u_at_points[0, n_pts:]/u_0, color="C2", linestyle="--") +ax.grid() +ax.set_xlabel(r"$X$") +ax.set_ylabel(r"$\tilde{U}_c$") +ax.set_ylim([0.3, 1.1]) + +plt.show() diff --git a/examples/inputs/ev.yaml b/examples/inputs/ev.yaml new file mode 100644 index 000000000..dd2bd40d5 --- /dev/null +++ b/examples/inputs/ev.yaml @@ -0,0 +1,119 @@ + +name: Eddy viscosity +description: Three turbines using the eddy viscosity model +floris_version: v4.x + +logging: + console: + enable: true + level: WARNING + file: + enable: false + level: WARNING + +solver: + type: turbine_grid + turbine_grid_points: 3 + +farm: + layout_x: + - 0.0 + - 630.0 + - 1260.0 + layout_y: + - 0.0 + - 0.0 + - 0.0 + turbine_type: + - nrel_5MW + +flow_field: + air_density: 1.225 + reference_wind_height: -1 # -1 is code for use the hub height + turbulence_intensities: + - 0.06 + wind_directions: + - 270.0 + wind_shear: 0.12 + wind_speeds: + - 8.0 + wind_veer: 0.0 + +wake: + model_strings: + combination_model: soed + deflection_model: none + turbulence_model: crespo_hernandez + velocity_model: eddy_viscosity + + enable_secondary_steering: false + enable_yaw_added_recovery: false + enable_transverse_velocities: false + enable_active_wake_mixing: false + + wake_deflection_parameters: + gauss: + ad: 0.0 + alpha: 0.58 + bd: 0.0 + beta: 0.077 + dm: 1.0 + ka: 0.38 + kb: 0.004 + jimenez: + ad: 0.0 + bd: 0.0 + kd: 0.05 + empirical_gauss: + horizontal_deflection_gain_D: 3.0 + vertical_deflection_gain_D: -1 + deflection_rate: 30 + mixing_gain_deflection: 0.0 + yaw_added_mixing_gain: 0.0 + + wake_velocity_parameters: + cc: + a_s: 0.179367259 + b_s: 0.0118889215 + c_s1: 0.0563691592 + c_s2: 0.13290157 + a_f: 3.11 + b_f: -0.68 + c_f: 2.41 + alpha_mod: 1.0 + gauss: + alpha: 0.58 + beta: 0.077 + ka: 0.38 + kb: 0.004 + jensen: + we: 0.05 + empirical_gauss: + wake_expansion_rates: + - 0.023 + - 0.008 + breakpoints_D: + - 10 + sigma_0_D: 0.28 + smoothing_length_D: 2.0 + mixing_gain_velocity: 2.0 + eddy_viscosity: + i_const_1: 0.05 + i_const_2: 16 + i_const_3: 0.5 + i_const_4: 10 + c_0: 2.0 + c_1: 1.5 + k_l: 0.0283 + k_a: 0.5 + von_Karman_constant: 0.4 + wd_std_ev: 3.0 + + wake_turbulence_parameters: + crespo_hernandez: + initial: 0.1 + constant: 0.5 + ai: 0.8 + downstream: -0.32 + wake_induced_mixing: + atmospheric_ti_gain: 0.0 diff --git a/floris/core/__init__.py b/floris/core/__init__.py index e37f9c113..01183860f 100644 --- a/floris/core/__init__.py +++ b/floris/core/__init__.py @@ -47,10 +47,12 @@ from .wake import WakeModelManager from .solver import ( cc_solver, + streamtube_expansion_solver, empirical_gauss_solver, full_flow_cc_solver, full_flow_empirical_gauss_solver, full_flow_sequential_solver, + full_flow_streamtube_expansion_solver, full_flow_turbopark_solver, sequential_solver, turbopark_solver, diff --git a/floris/core/core.py b/floris/core/core.py index d12005ba5..b2be53655 100644 --- a/floris/core/core.py +++ b/floris/core/core.py @@ -20,11 +20,13 @@ full_flow_cc_solver, full_flow_empirical_gauss_solver, full_flow_sequential_solver, + full_flow_streamtube_expansion_solver, full_flow_turbopark_solver, Grid, PointsGrid, sequential_solver, State, + streamtube_expansion_solver, TurbineCubatureGrid, TurbineGrid, turbopark_solver, @@ -198,6 +200,13 @@ def steady_state_atmospheric_condition(self): self.grid, self.wake ) + elif vel_model=="eddy_viscosity": + streamtube_expansion_solver( + self.farm, + self.flow_field, + self.grid, + self.wake + ) else: sequential_solver( self.farm, @@ -225,6 +234,8 @@ def solve_for_viz(self): full_flow_turbopark_solver(self.farm, self.flow_field, self.grid, self.wake) elif vel_model=="empirical_gauss": full_flow_empirical_gauss_solver(self.farm, self.flow_field, self.grid, self.wake) + elif vel_model=="eddy_viscosity": + full_flow_streamtube_expansion_solver(self.farm, self.flow_field, self.grid, self.wake) else: full_flow_sequential_solver(self.farm, self.flow_field, self.grid, self.wake) @@ -255,10 +266,12 @@ def solve_for_points(self, x, y, z): if vel_model == "cc" or vel_model == "turbopark": raise NotImplementedError( "solve_for_points is currently only available with the "+\ - "gauss, jensen, and empirical_gauss models." + "gauss, jensen, empirical_gauss, and eddy_viscosity models." ) elif vel_model == "empirical_gauss": full_flow_empirical_gauss_solver(self.farm, self.flow_field, field_grid, self.wake) + elif vel_model=="eddy_viscosity": + full_flow_streamtube_expansion_solver(self.farm, self.flow_field, field_grid, self.wake) else: full_flow_sequential_solver(self.farm, self.flow_field, field_grid, self.wake) diff --git a/floris/core/solver.py b/floris/core/solver.py index a7c3d8796..ea65a0e1a 100644 --- a/floris/core/solver.py +++ b/floris/core/solver.py @@ -1526,3 +1526,368 @@ def full_flow_empirical_gauss_solver( flow_field.u_sorted = flow_field.u_initial_sorted - wake_field flow_field.v_sorted += v_wake flow_field.w_sorted += w_wake + + +# @profile +def streamtube_expansion_solver( + farm: Farm, + flow_field: FlowField, + grid: TurbineGrid, + model_manager: WakeModelManager +) -> None: + # Algorithm + # For each turbine, calculate its effect on every downstream turbine. + # Also calculates the expansion in the streamtube of upstream turbines. + + # <> + deflection_model_args = model_manager.deflection_model.prepare_function(grid, flow_field) + deficit_model_args = model_manager.velocity_model.prepare_function(grid, flow_field) + #combination_model_args = model_manager.combination_model.prepare_function(grid, flow_field) + + + # Ambient turbulent intensity should be a copy of n_findex-long turbulence_intensity + # with dimensions expanded for (n_turbines, grid, grid) + ambient_turbulence_intensities = flow_field.turbulence_intensities.copy() + ambient_turbulence_intensities = ambient_turbulence_intensities[:, None, None, None] + # TODO: should TI be updated at each turbine? Seems not? + turbine_turbulence_intensity = flow_field.turbulence_intensities[:, None, None, None] + turbine_turbulence_intensity = np.repeat(turbine_turbulence_intensity, farm.n_turbines, axis=1) + + # Declare storage for centerline velocities, wake_widths, and thrust coefficients + centerline_velocities = np.zeros((flow_field.n_findex, farm.n_turbines, farm.n_turbines)) + wake_widths_squared = np.zeros((flow_field.n_findex, farm.n_turbines, farm.n_turbines)) + thrust_coefficients = np.zeros((flow_field.n_findex, farm.n_turbines)) + + # Calculate the velocity deficit sequentially from upstream to downstream turbines + for tindex in range(grid.n_turbines): + + # Get the current turbine quantities + x_i = np.mean(grid.x_sorted[:, tindex:tindex+1], axis=(2, 3)) + x_i = x_i[:, :, None, None] + y_i = np.mean(grid.y_sorted[:, tindex:tindex+1], axis=(2, 3)) + y_i = y_i[:, :, None, None] + z_i = np.mean(grid.z_sorted[:, tindex:tindex+1], axis=(2, 3)) + z_i = z_i[:, :, None, None] + + ct_i = thrust_coefficient( + velocities=flow_field.u_sorted, + turbulence_intensities=flow_field.turbulence_intensity_field_sorted, + air_density=flow_field.air_density, + yaw_angles=farm.yaw_angles_sorted, + tilt_angles=farm.tilt_angles_sorted, + power_setpoints=farm.power_setpoints_sorted, + awc_modes=farm.awc_modes_sorted, + awc_amplitudes=farm.awc_amplitudes_sorted, + thrust_coefficient_functions=farm.turbine_thrust_coefficient_functions, + tilt_interps=farm.turbine_tilt_interps, + correct_cp_ct_for_tilt=farm.correct_cp_ct_for_tilt_sorted, + turbine_type_map=farm.turbine_type_map_sorted, + turbine_power_thrust_tables=farm.turbine_power_thrust_tables, + ix_filter=[tindex], + average_method=grid.average_method, + cubature_weights=grid.cubature_weights, + multidim_condition=flow_field.multidim_conditions + ) + thrust_coefficients[:, tindex] = ct_i[:, 0] + # Since we are filtering for the i'th turbine in the thrust coefficient function, + # get the first index here (0:1) + axial_induction_i = axial_induction( + velocities=flow_field.u_sorted, + turbulence_intensities=flow_field.turbulence_intensity_field_sorted, + air_density=flow_field.air_density, + yaw_angles=farm.yaw_angles_sorted, + tilt_angles=farm.tilt_angles_sorted, + power_setpoints=farm.power_setpoints_sorted, + awc_modes=farm.awc_modes_sorted, + awc_amplitudes=farm.awc_amplitudes_sorted, + axial_induction_functions=farm.turbine_axial_induction_functions, + tilt_interps=farm.turbine_tilt_interps, + correct_cp_ct_for_tilt=farm.correct_cp_ct_for_tilt_sorted, + turbine_type_map=farm.turbine_type_map_sorted, + turbine_power_thrust_tables=farm.turbine_power_thrust_tables, + ix_filter=[tindex], + average_method=grid.average_method, + cubature_weights=grid.cubature_weights, + multidim_condition=flow_field.multidim_conditions + ) + # Since we are filtering for the i'th turbine in the axial induction function, + # get the first index here (0:1) + + # Can all of this be avoided? Don't need these at all grid points anyway? + ct_i = ct_i[:, 0:1, None, None] + axial_induction_i = axial_induction_i[:, 0:1, None, None] + turbulence_intensity_i = turbine_turbulence_intensity[:, tindex:tindex+1] + yaw_angle_i = farm.yaw_angles_sorted[:, tindex:tindex+1, None, None] + hub_height_i = farm.hub_heights_sorted[:, tindex:tindex+1, None, None] + rotor_diameter_i = farm.rotor_diameters_sorted[:, tindex:tindex+1, None, None] + + effective_yaw_i = np.zeros_like(yaw_angle_i) + effective_yaw_i += yaw_angle_i + + if (model_manager.enable_secondary_steering + or model_manager.enable_transverse_velocities + or model_manager.enable_yaw_added_recovery + ): + raise NotImplementedError( + "Secondary effects not available for this model." + ) + + if np.any(farm.yaw_angles_sorted): + raise NotImplementedError( + "WARNING: Deflection with the eddy viscosity model has not been implemented." + " Set yaw_angles to zero to run the eddy viscosity model." + ) + + _ = model_manager.deflection_model.function( + x_i, + y_i, + effective_yaw_i, + turbulence_intensity_i, + ct_i, + rotor_diameter_i, + **deflection_model_args, + ) + + _, wake_width_squared_i = model_manager.velocity_model.function( + x_i[:,:,0,0], + turbulence_intensity_i[:,:,0,0], + ct_i[:,:,0,0], + hub_height_i[:,:,0,0], + rotor_diameter_i[:,:,0,0], + **deficit_model_args, + ) + + # Store the centerline velocities and wake widths for each turbine--turbine pair + # centerline_velocities[:, tindex, :] = centerline_velocity_i + wake_widths_squared[:, tindex, :] = wake_width_squared_i + + # Now, apply the streamtube expansion. + (centerline_velocities, + wake_widths_squared) = model_manager.velocity_model.streamtube_expansion( + x_i[:,:,0,0], + (grid.x_sorted.mean(axis=(2,3)) - x_i[:,:,0,0]), + (grid.y_sorted.mean(axis=(2,3)) - y_i[:,:,0,0]), + (farm.hub_heights_sorted - z_i[:,:,0,0]), + thrust_coefficients, + axial_induction_i[:,:,0,0], + wake_widths_squared, + rotor_diameter_i[:,:,0,0], + **deficit_model_args, + ) + + # Now, can compute the actual velocities at each turbine location + velocities = model_manager.velocity_model.evaluate_velocities( + centerline_velocities, + thrust_coefficients, + farm.rotor_diameters_sorted, + grid.y_sorted.mean(axis=(2,3)), + farm.hub_heights_sorted, + **deficit_model_args, + ) + + + # Combine + velocity_field = model_manager.combination_model.function(velocities) + + # Compute absolute velocity field based on all turbines up to i + flow_field.u_sorted = velocity_field * flow_field.u_initial_sorted + +def full_flow_streamtube_expansion_solver( + farm: Farm, + flow_field: FlowField, + flow_field_grid: FlowFieldGrid, + model_manager: WakeModelManager +) -> None: + + # Get the flow quantities and turbine performance + turbine_grid_farm = copy.deepcopy(farm) + turbine_grid_flow_field = copy.deepcopy(flow_field) + + turbine_grid_farm.construct_turbine_map() + turbine_grid_farm.construct_turbine_thrust_coefficient_functions() + turbine_grid_farm.construct_turbine_axial_induction_functions() + turbine_grid_farm.construct_turbine_power_functions() + turbine_grid_farm.construct_hub_heights() + turbine_grid_farm.construct_rotor_diameters() + turbine_grid_farm.construct_turbine_TSRs() + turbine_grid_farm.construct_turbine_ref_tilts() + turbine_grid_farm.construct_turbine_tilt_interps() + turbine_grid_farm.construct_turbine_correct_cp_ct_for_tilt() + turbine_grid_farm.set_tilt_to_ref_tilt(flow_field.n_findex) + + turbine_grid = TurbineGrid( + turbine_coordinates=turbine_grid_farm.coordinates, + turbine_diameters=turbine_grid_farm.rotor_diameters, + wind_directions=turbine_grid_flow_field.wind_directions, + grid_resolution=3, + ) + turbine_grid_farm.expand_farm_properties( + turbine_grid_flow_field.n_findex, + turbine_grid.sorted_coord_indices + ) + turbine_grid_flow_field.initialize_velocity_field(turbine_grid) + turbine_grid_farm.initialize(turbine_grid.sorted_indices) + + # Perform turbine grid calculation to compute inflows + streamtube_expansion_solver( + turbine_grid_farm, + turbine_grid_flow_field, + turbine_grid, + model_manager + ) + + ### Referring to the quantities from above, calculate the wake in the full grid + + # Use full flow_field here to use the full grid in the wake models + deflection_model_args = model_manager.deflection_model.prepare_function( + flow_field_grid, flow_field + ) + deficit_model_args = model_manager.velocity_model.prepare_function(flow_field_grid, flow_field) + + + ambient_turbulence_intensities = flow_field.turbulence_intensities.copy() + ambient_turbulence_intensities = ambient_turbulence_intensities[:, None, None, None] + # TODO: should TI be updated at each turbine? Seems not? + turbine_turbulence_intensity = flow_field.turbulence_intensities[:, None, None, None] + turbine_turbulence_intensity = np.repeat(turbine_turbulence_intensity, farm.n_turbines, axis=1) + + # Declare storage for centerline velocities, wake_widths, and thrust coefficients + # TODO: Should these be expanded to all grid points, possibly? + n_x = flow_field_grid.x_sorted.shape[1] + centerline_velocities = np.zeros((flow_field.n_findex, farm.n_turbines, n_x)) + wake_widths_squared = np.zeros((flow_field.n_findex, farm.n_turbines, n_x)) + thrust_coefficients = np.zeros((flow_field.n_findex, farm.n_turbines)) + + # Calculate the velocity deficit sequentially from upstream to downstream turbines + for tindex in range(flow_field_grid.n_turbines): + + # Get the current turbine quantities + x_i = np.mean(turbine_grid.x_sorted[:, tindex:tindex+1], axis=(2,3)) + x_i = x_i[:, :, None, None] + y_i = np.mean(turbine_grid.y_sorted[:, tindex:tindex+1], axis=(2,3)) + y_i = y_i[:, :, None, None] + z_i = np.mean(turbine_grid.z_sorted[:, tindex:tindex+1], axis=(2,3)) + z_i = z_i[:, :, None, None] + + # Thrust, axial induction need to be recomputed as not stored in turbine solve. + ct_i = thrust_coefficient( + velocities=turbine_grid_flow_field.u_sorted, + turbulence_intensities=turbine_grid_flow_field.turbulence_intensity_field_sorted, + air_density=turbine_grid_flow_field.air_density, + yaw_angles=turbine_grid_farm.yaw_angles_sorted, + tilt_angles=turbine_grid_farm.tilt_angles_sorted, + power_setpoints=turbine_grid_farm.power_setpoints_sorted, + awc_modes=turbine_grid_farm.awc_modes_sorted, + awc_amplitudes=turbine_grid_farm.awc_amplitudes_sorted, + thrust_coefficient_functions=turbine_grid_farm.turbine_thrust_coefficient_functions, + tilt_interps=turbine_grid_farm.turbine_tilt_interps, + correct_cp_ct_for_tilt=turbine_grid_farm.correct_cp_ct_for_tilt_sorted, + turbine_type_map=turbine_grid_farm.turbine_type_map_sorted, + turbine_power_thrust_tables=turbine_grid_farm.turbine_power_thrust_tables, + ix_filter=[tindex], + average_method=turbine_grid.average_method, + cubature_weights=turbine_grid.cubature_weights, + multidim_condition=turbine_grid_flow_field.multidim_conditions, + ) + thrust_coefficients[:, tindex] = ct_i[:, 0] + # Since we are filtering for the i'th turbine in the thrust coefficient function, + # get the first index here (0:1) + axial_induction_i = axial_induction( + velocities=turbine_grid_flow_field.u_sorted, + turbulence_intensities=turbine_grid_flow_field.turbulence_intensity_field_sorted, + air_density=turbine_grid_flow_field.air_density, + yaw_angles=turbine_grid_farm.yaw_angles_sorted, + tilt_angles=turbine_grid_farm.tilt_angles_sorted, + power_setpoints=turbine_grid_farm.power_setpoints_sorted, + awc_modes=turbine_grid_farm.awc_modes_sorted, + awc_amplitudes=turbine_grid_farm.awc_amplitudes_sorted, + axial_induction_functions=turbine_grid_farm.turbine_axial_induction_functions, + tilt_interps=turbine_grid_farm.turbine_tilt_interps, + correct_cp_ct_for_tilt=turbine_grid_farm.correct_cp_ct_for_tilt_sorted, + turbine_type_map=turbine_grid_farm.turbine_type_map_sorted, + turbine_power_thrust_tables=turbine_grid_farm.turbine_power_thrust_tables, + ix_filter=[tindex], + average_method=turbine_grid.average_method, + cubature_weights=turbine_grid.cubature_weights, + multidim_condition=turbine_grid_flow_field.multidim_conditions, + ) + # Since we are filtering for the i'th turbine in the axial induction function, + # get the first index here (0:1) + ct_i = ct_i[:, 0:1, None, None] + axial_induction_i = axial_induction_i[:, 0:1, None, None] + turbulence_intensity_i = turbine_turbulence_intensity[:, tindex:tindex+1] + yaw_angle_i = turbine_grid_farm.yaw_angles_sorted[:, tindex:tindex+1, None, None] + hub_height_i = turbine_grid_farm.hub_heights_sorted[:, tindex:tindex+1, None, None] + rotor_diameter_i = turbine_grid_farm.rotor_diameters_sorted[:, tindex:tindex+1, None, None] + + effective_yaw_i = np.zeros_like(yaw_angle_i) + effective_yaw_i += yaw_angle_i + + if model_manager.enable_secondary_steering: + raise NotImplementedError( + "Secondary steering not available for this model.") + + if model_manager.enable_transverse_velocities: + raise NotImplementedError( + "Transverse velocities not used in this model.") + + if np.any(turbine_grid_farm.yaw_angles_sorted): + raise NotImplementedError( + "WARNING: Deflection with the eddy viscosity model has not been implemented." + " Set yaw_angles to zero to run the eddy viscosity model." + ) + + _ = model_manager.deflection_model.function( + x_i, + y_i, + effective_yaw_i, + turbulence_intensity_i, + ct_i, + rotor_diameter_i, + **deflection_model_args, + ) + + # NOTE: exponential + _, wake_width_squared_i = model_manager.velocity_model.function( + x_i[:,:,0,0], + turbulence_intensity_i[:,:,0,0], + ct_i[:,:,0,0], + hub_height_i[:,:,0,0], + rotor_diameter_i[:,:,0,0], + **deficit_model_args, + ) + + # Store the centerline velocities and wake widths for each turbine--turbine pair + # centerline_velocities[:, tindex, :] = centerline_velocity_i + wake_widths_squared[:, tindex, :] = wake_width_squared_i + + # Apply the streamtube expansion. + (centerline_velocities, + wake_widths_squared) = model_manager.velocity_model.streamtube_expansion( + x_i[:,:,0,0], + (turbine_grid.x_sorted.mean(axis=(2,3)) - x_i[:,:,0,0]), + (turbine_grid.y_sorted.mean(axis=(2,3)) - y_i[:,:,0,0]), + (turbine_grid_farm.hub_heights_sorted - z_i[:,:,0,0]), + thrust_coefficients, + axial_induction_i[:,:,0,0], + wake_widths_squared, + rotor_diameter_i[:,:,0,0], + # + **deficit_model_args, + ) + + # Now, can compute the actual velocities at each TODO grid point + velocities = model_manager.velocity_model.evaluate_velocities( + centerline_velocities, + thrust_coefficients, + turbine_grid_farm.rotor_diameters_sorted, + turbine_grid.y_sorted.mean(axis=(2,3)), + turbine_grid_farm.hub_heights_sorted, + **deficit_model_args, + ) + + # Combine + velocity_field = model_manager.combination_model.function(velocities) + + # Compute absolute velocity field based on all turbines up to i + flow_field.u_sorted = velocity_field * flow_field.u_initial_sorted diff --git a/floris/core/wake.py b/floris/core/wake.py index e58a85cb4..ce72e228f 100644 --- a/floris/core/wake.py +++ b/floris/core/wake.py @@ -6,6 +6,7 @@ from floris.core.wake_combination import ( FLS, MAX, + SOED, SOSFS, ) from floris.core.wake_deflection import ( @@ -21,6 +22,7 @@ ) from floris.core.wake_velocity import ( CumulativeGaussCurlVelocityDeficit, + EddyViscosityVelocity, EmpiricalGaussVelocityDeficit, GaussVelocityDeficit, JensenVelocityDeficit, @@ -34,7 +36,8 @@ "combination_model": { "fls": FLS, "max": MAX, - "sosfs": SOSFS + "sosfs": SOSFS, + "soed": SOED, }, "deflection_model": { "jimenez": JimenezVelocityDeflection, @@ -54,6 +57,7 @@ "jensen": JensenVelocityDeficit, "turbopark": TurbOParkVelocityDeficit, "empirical_gauss": EmpiricalGaussVelocityDeficit, + "eddy_viscosity": EddyViscosityVelocity, "turboparkgauss": TurboparkgaussVelocityDeficit, }, } diff --git a/floris/core/wake_combination/__init__.py b/floris/core/wake_combination/__init__.py index 246aab65c..7617fdcf9 100644 --- a/floris/core/wake_combination/__init__.py +++ b/floris/core/wake_combination/__init__.py @@ -1,4 +1,5 @@ from floris.core.wake_combination.fls import FLS from floris.core.wake_combination.max import MAX +from floris.core.wake_combination.soed import SOED from floris.core.wake_combination.sosfs import SOSFS diff --git a/floris/core/wake_combination/soed.py b/floris/core/wake_combination/soed.py new file mode 100644 index 000000000..7bfec7f38 --- /dev/null +++ b/floris/core/wake_combination/soed.py @@ -0,0 +1,31 @@ +import logging +from typing import Any, Dict + +import numpy as np +from attrs import define + +from floris.core import BaseModel, FlowField + + +logger = logging.getLogger(name="floris") + +@define +class SOED(BaseModel): + """ + Sum of energy deficit, as described in Kuo et al, is used by the eddy vicosity model. + + Kuo et al, 2014 + https://mechanicaldesign.asmedigitalcollection.asme.org/IMECE/proceedings/IMECE2014/46521/V06BT07A074/263017 + """ + + def function( + self, + U_tilde_field: np.ndarray, + ): + n_turbines = U_tilde_field.shape[1] + + U_tilde_combined = np.sqrt( + np.maximum(1 - n_turbines + np.sum(U_tilde_field**2, axis=1), 0.0) + ) + + return U_tilde_combined diff --git a/floris/core/wake_velocity/__init__.py b/floris/core/wake_velocity/__init__.py index 07762a5cd..3010cdc2c 100644 --- a/floris/core/wake_velocity/__init__.py +++ b/floris/core/wake_velocity/__init__.py @@ -1,5 +1,6 @@ from floris.core.wake_velocity.cumulative_gauss_curl import CumulativeGaussCurlVelocityDeficit +from floris.core.wake_velocity.eddy_viscosity import EddyViscosityVelocity from floris.core.wake_velocity.empirical_gauss import EmpiricalGaussVelocityDeficit from floris.core.wake_velocity.gauss import GaussVelocityDeficit from floris.core.wake_velocity.jensen import JensenVelocityDeficit diff --git a/floris/core/wake_velocity/eddy_viscosity.py b/floris/core/wake_velocity/eddy_viscosity.py new file mode 100644 index 000000000..bba503136 --- /dev/null +++ b/floris/core/wake_velocity/eddy_viscosity.py @@ -0,0 +1,410 @@ +import logging +from typing import Any, Dict + +import matplotlib.pyplot as plt +import numexpr as ne +import numpy as np +from attrs import define, field +from scipy.integrate import solve_ivp + +from floris.core import ( + BaseModel, + Farm, + FlowField, + Grid, + Turbine, +) +from floris.utilities import ( + cosd, + sind, + tand, +) + + +logger = logging.getLogger(name="floris") + +@define +class EddyViscosityVelocity(BaseModel): + + k_l: float = field(default=0.015*np.sqrt(3.56)) + k_a: float = field(default=0.5) + von_Karman_constant: float = field(default=0.41) + + i_const_1: float = field(default=0.05) + i_const_2: float = field(default=16) + i_const_3: float = field(default=0.5) + i_const_4: float = field(default=10) + + # Below are likely not needed [or, I'll need to think more about it] + filter_const_1: float = field(default=0.65) + filter_const_2: float = field(default=4.5) + filter_const_3: float = field(default=23.32) + filter_const_4: float = field(default=1/3) + filter_cutoff_D: float = field(default=0.0) + + c_0: float = field(default=2.0) + c_1: float = field(default=1.5) + + wd_std_ev: float = field(default=3.0) # Also try with 0.0 for no meandering + + def prepare_function( + self, + grid: Grid, + flow_field: FlowField, + ) -> Dict[str, Any]: + + kwargs = { + "x": grid.x_sorted, + "y": grid.y_sorted, + "z": grid.z_sorted, + "u_initial": flow_field.u_initial_sorted, + "wind_veer": flow_field.wind_veer + } + return kwargs + + # @profile + def function( + self, + x_i: np.ndarray, + turbulence_intensity_i: np.ndarray, + ct_i: np.ndarray, + hub_height_i: float, + rotor_diameter_i: np.ndarray, + # enforces the use of the below as keyword arguments and adherence to the + # unpacking of the results from prepare_function() + *, + x: np.ndarray, + y: np.ndarray, + z: np.ndarray, + u_initial: np.ndarray, + wind_veer: np.ndarray, + ) -> np.ndarray: + + # Non-dimensionalize and center distances + x_tilde = ((x.mean(axis=(2,3)) - x_i) / rotor_diameter_i) + + # Compute centerline velocities + # TODO: This is using an "updated" TI. Is that appropriate? + U_tilde_c_initial = initial_centerline_velocity( + ct_i, + turbulence_intensity_i, + self.i_const_1, + self.i_const_2, + self.i_const_3, + self.i_const_4 + ) + + # Solve ODE to find centerline velocities at each x + U_tilde_c = np.zeros_like(x_tilde) + for findex in range(x_tilde.shape[0]): + x_tilde_unique, unique_ind = np.unique(x_tilde[findex, :], return_inverse=True) + sorting_indices = np.argsort(x_tilde_unique) + x_tilde_sorted = x_tilde_unique[sorting_indices] + valid_indices = x_tilde_sorted >= 2 + x_tilde_eval = x_tilde_sorted[valid_indices] + if len(x_tilde_eval) == 0: # No downstream locations to fill + U_tilde_c[findex, :] = U_tilde_c_initial[findex] + continue + sol = solve_ivp( + fun=centerline_ode, + t_span=[2, x_tilde_eval[-1]], + y0=U_tilde_c_initial[findex,:], + method='RK45', + t_eval=x_tilde_eval, + args=( + turbulence_intensity_i[findex,0], + ct_i[findex,0], + hub_height_i[findex,0], + rotor_diameter_i[findex,0], + self.k_a, + self.k_l, + self.von_Karman_constant, + self.filter_cutoff_D, + self.filter_const_1, + self.filter_const_2, + self.filter_const_3, + self.filter_const_4 + ) + ) + + if sol.status == -1: + raise RuntimeError( + f"Eddy viscosity ODE solver failed to converge for findex={findex}.\n" + + "t_span: " + np.array2string(np.array([2, x_tilde_eval[-1]])) + "\n" + + "y0: " + np.array2string(U_tilde_c_initial[findex,:]) + "\n" + + "t_eval: " + np.array2string(x_tilde_eval) + "\n\n" + + "This may be caused by an initial condition for the ODE that is " + + "greater than 1." + ) + + # Extract the solution + if (sol.t != x_tilde_eval).any(): + raise ValueError("ODE solver did not return requested values") + U_tilde_c_eval = sol.y.flatten() + + U_tilde_c_fill = np.full_like( + x_tilde_sorted[x_tilde_sorted < 2], + U_tilde_c_initial[findex,:] + ) + U_tilde_c_sorted = np.concatenate((U_tilde_c_fill, U_tilde_c_eval)) + + # "Unsort", and broadcast back to shape of x_tilde + U_tilde_c_unique = U_tilde_c_sorted[np.argsort(sorting_indices)] + U_tilde_c_findex = U_tilde_c_unique[unique_ind] + U_tilde_c[findex, :] = U_tilde_c_findex + + + # Compute wake width + w_tilde_sq = wake_width_squared(ct_i, U_tilde_c) + + # Correct for wake meandering + U_tilde_c_meandering = wake_meandering_centerline_correction( + U_tilde_c, w_tilde_sq, x_tilde, self.wd_std_ev + ) + + # Recompute wake width + w_tilde_sq_meandering = wake_width_squared(ct_i, U_tilde_c_meandering) + + # Set all upstream values (including current turbine's position) to no wake + U_tilde_c_meandering[x_tilde < 0.1] = 1 + w_tilde_sq_meandering[x_tilde < 0.1] = 0 + + # Return velocities NOT as deficits + return U_tilde_c_meandering, w_tilde_sq_meandering + + def streamtube_expansion( + self, + x_i, + x_from_i, + y_from_i, + z_from_i, + ct_all, + axial_induction_i, + w_tilde_sq_tt, + rotor_diameter_i, + *, + x, + y, + z, + u_initial, + wind_veer, + ): + # Non-dimensionalize and center distances + x_tilde_points = (x.mean(axis=(2,3)) - x_i) / rotor_diameter_i + x_tilde = x_from_i / rotor_diameter_i + y_tilde = y_from_i / rotor_diameter_i + z_tilde = z_from_i / rotor_diameter_i + + # Compute wake width + e_tilde = wake_width_streamtube_correction_term( + axial_induction_i, + x_tilde_points, + x_tilde, + y_tilde, + z_tilde, + self.c_0, + self.c_1 + ) + + w_tilde_sq_tt = expanded_wake_width_squared(w_tilde_sq_tt, e_tilde) + + U_tilde_c_tt = expanded_wake_centerline_velocity(ct_all[:,:,None], w_tilde_sq_tt) + + return U_tilde_c_tt, w_tilde_sq_tt + + def evaluate_velocities( + self, + U_tilde_c_tt, + ct_all, + rotor_diameters, + y_turbines, + z_turbines, + *, + x, + y, + z, + u_initial, + wind_veer, + ): + # Non-dimensionalize and center distances + y_tilde_rel = ( + (y[:,None,:,:,:] - y_turbines[:,:,None,None,None]) + / rotor_diameters[:,:,None,None,None] + ) + z_tilde_rel = ( + (z[:,None,:,:,:] - z_turbines[:,:,None,None,None]) + / rotor_diameters[:,:,None,None,None] + ) + # TODO: Check working as expected with correct D, hh being applied + # when there are multiple turbine types + + # Compute radial positions + r_tilde_sq = y_tilde_rel**2 + z_tilde_rel**2 + + U_tilde_r_tt = compute_off_center_velocities( + U_tilde_c_tt, + ct_all[:,:,None], + r_tilde_sq + ) + + return U_tilde_r_tt + + +def compute_off_center_velocities(U_tilde_c, Ct, r_tilde_sq): + """ + Compute the off-centerline velocities using the eddy viscosity model + y_, z_ supposed to be defined from the center of the rotor. + """ + w_tilde_sq = wake_width_squared(Ct, U_tilde_c) + U_tilde_c_mask = U_tilde_c == 1 + # As long as U_tilde_c is 1, w_tilde_sq won't affect result, but this + # silences a division by zero warning + w_tilde_sq[U_tilde_c_mask] = 1 + U_tilde = ( + 1 + - (1 - U_tilde_c[:,:,:,None,None]) + * np.exp(-r_tilde_sq/w_tilde_sq[:,:,:,None,None]) + ) + + return U_tilde + +def wake_width_squared(Ct, U_tilde_c): + """ + Compute the wake width squared using the eddy viscosity model + """ + U_tilde_c_mask = U_tilde_c < 1 + + w_tilde_sq = np.zeros_like(U_tilde_c) + Ct = _resize_Ct(Ct, U_tilde_c) + w_tilde_sq[U_tilde_c_mask] = ( + Ct[U_tilde_c_mask] / (4 * (1 - U_tilde_c[U_tilde_c_mask]) * (1 + U_tilde_c[U_tilde_c_mask])) + ) + w_tilde_sq.reshape(U_tilde_c.shape) + return w_tilde_sq + +def centerline_ode( + x_tilde, + U_tilde_c, + ambient_ti, + Ct, + hh, + D, + k_a, + k_l, + von_Karman_constant, + filter_cutoff_D, + filter_const_1, + filter_const_2, + filter_const_3, + filter_const_4 + ): + """ + Define the ODE for the centerline velocities + """ + + # Local component, nondimensionalized by U_inf*D (compared to Gunn 2019's K_l) + K_l_tilde = k_l * np.sqrt(wake_width_squared(Ct, U_tilde_c)) * (1 - U_tilde_c) + + # Ambient component, nondimensionalized by U_inf*D (compared to Gunn 2019's K_a, eq. (9)) + K_a_tilde = k_a * ambient_ti * von_Karman_constant * (hh/D) + + F = filter_function( + x_tilde, filter_cutoff_D, filter_const_1, filter_const_2, filter_const_3, filter_const_4 + ) + eddy_viscosity_tilde = F*(K_l_tilde + K_a_tilde) + + dU_tilde_c_dx_tilde = ( + 16 * eddy_viscosity_tilde + * (U_tilde_c**3 - U_tilde_c**2 - U_tilde_c + 1) + / (U_tilde_c * Ct) + ) + + return [dU_tilde_c_dx_tilde] + +def filter_function( + x_tilde, + filter_cutoff_D, + filter_const_1, + filter_const_2, + filter_const_3, + filter_const_4 + ): + """ + Gunn uses 1 for all x, whereas Ainslie uses the filter function. + To reproduce Gunn, set filter_cutoff_D to 0. + """ + + # The form given by Ainslie appears to produce complex numbers for x_tilde < filter_const_2. + # Here, we take only the real component. + if filter_cutoff_D == 0: + return 1 + else: # TEMP + logger.warning("Cannot solve ODE with filter function. Returning 1.") + return 1 + F = np.where(x_tilde < filter_cutoff_D, + np.real(filter_const_1 + ((x_tilde - filter_const_2) / filter_const_3)**filter_const_4), + 1 + ) + F = np.clip(F, 0, 1) + + return F + +def initial_centerline_velocity(Ct, ambient_ti, i_const_1, i_const_2, i_const_3, i_const_4): + + # The below are from Ainslie (1988) + initial_velocity_deficit = ( + Ct + - i_const_1 + - (i_const_2 * Ct - i_const_3) * ambient_ti / i_const_4 + ) + + if (initial_velocity_deficit < 0).any(): + logger.warning( + "Initial velocity deficit is negative. Setting to 0." + ) + initial_velocity_deficit[initial_velocity_deficit < 0] = 0 + + U_c0_ = 1 - initial_velocity_deficit + + return U_c0_ + +def wake_meandering_centerline_correction(U_tilde_c, w_tilde_sq, x_tilde, wd_std_ev): + wd_std_rad = np.deg2rad(wd_std_ev) + + m = np.sqrt(1 + 2*wd_std_rad**2 * x_tilde**2/w_tilde_sq) + + U_tilde_c_meandering = 1/m * U_tilde_c + (m-1)/m + + return U_tilde_c_meandering + + +def wake_width_streamtube_correction_term(ai_i, x_pts, x_ji_, y_ji_, z_ji_, c_0, c_1): + e_i_ = np.sqrt(1-ai_i) * (1/np.sqrt(1-2*ai_i) - 1) + + e_ji_ = c_0 * e_i_ * np.exp(-(y_ji_**2 + z_ji_**2) / c_1**2) + e_ji_[x_ji_ >= 0] = 0 # Does not affect wakes of self or downstream turbines + downstream_mask = (x_pts > 0).astype(int) + e_ji_ = e_ji_[:,:,None] * downstream_mask[:,None,:] # Affects only downstream locations + + return e_ji_ + +def expanded_wake_width_squared(w_tilde_sq, e_tilde): + return (np.sqrt(w_tilde_sq) + e_tilde)**2 + +def expanded_wake_centerline_velocity(Ct, w_tilde_sq): + w_tilde_sq_mask = w_tilde_sq > 1e-6 + expanded_U_tilde_c = np.ones_like(w_tilde_sq) + Ct = _resize_Ct(Ct, w_tilde_sq) + expanded_U_tilde_c[w_tilde_sq_mask] = np.sqrt( + 1 - Ct[w_tilde_sq_mask]/(4*w_tilde_sq[w_tilde_sq_mask]) + ) + expanded_U_tilde_c.reshape(w_tilde_sq.shape) + + return expanded_U_tilde_c + +def _resize_Ct(Ct, resize_like): + if type(Ct) == np.ndarray: + Ct = np.repeat(Ct, resize_like.shape[-1], axis=-1) + else: + Ct = Ct * np.ones_like(resize_like) + return Ct diff --git a/floris/core/wake_velocity/recreate_Fig1.py b/floris/core/wake_velocity/recreate_Fig1.py new file mode 100644 index 000000000..c91a8e813 --- /dev/null +++ b/floris/core/wake_velocity/recreate_Fig1.py @@ -0,0 +1,58 @@ +import matplotlib.pyplot as plt +import numpy as np +from eddy_viscosity import ( + compute_centerline_velocities, + compute_off_center_velocities, + wake_width_squared, +) + + +D_T1 = 1 +Ct_T1 = 0.8 +x_0_T1 = 0 + +# Suggested settings +D_T2 = 1.52 +Ct_T2 = 0.34 +x_0_T2 = 1.5 + +# Retuned settings +Ct_T3 = 0.39 +D_T3 = 2.0 + + +hh = 1.0 +ambient_ti = 0.06 +U_inf = 8.0 + +x_test_T1 = np.linspace(2*D_T1, 10*D_T1, 100) +x_test_T2 = np.linspace(2*D_T2, 10*D_T2, 100) + +U_c_T1, x_out_T1 = compute_centerline_velocities(x_test_T1, U_inf, ambient_ti, Ct_T1, hh, D_T1) +U_c_T2, x_out_T2 = compute_centerline_velocities(x_test_T2, U_inf, ambient_ti, Ct_T2, hh, D_T2) +U_c_T3, x_out_T3 = compute_centerline_velocities(x_test_T2, U_inf, ambient_ti, Ct_T3, hh, D_T3) + +fig, ax = plt.subplots(2,1) +ax[0].plot(x_out_T1+x_0_T1, U_c_T1, color="C0", linestyle="solid", label="T1") +ax[0].plot(x_out_T2+x_0_T2, U_c_T2, color="C2", linestyle="dashed", label="T2, suggested") +ax[0].plot(x_out_T3+x_0_T2, U_c_T3, color="C3", linestyle="dashed", label="T2, retuned") +ax[0].set_xlabel("x_ [D]") +ax[0].set_ylabel("U_c_ [-]") +ax[0].grid() +ax[0].set_xlim([0, 12]) +ax[0].legend() + +w_T1 = np.sqrt(wake_width_squared(Ct_T1, U_c_T1)) +w_T2 = np.sqrt(wake_width_squared(Ct_T2, U_c_T2)) +w_T3 = np.sqrt(wake_width_squared(Ct_T3, U_c_T3)) + +ax[1].plot(x_out_T1+x_0_T1, w_T1, color="C0", linestyle="solid", label="T1") +ax[1].plot(x_out_T2+x_0_T2, w_T2*D_T2, color="C2", linestyle="dashed", label="T2, suggested") +ax[1].plot(x_out_T3+x_0_T2, w_T3*D_T3, color="C3", linestyle="dashed", label="T2, retuned") +ax[1].set_xlabel("x_ [D]") +ax[1].set_ylabel("w [m]") +ax[1].grid() +ax[1].set_xlim([0, 12]) +#ax[1].legend() + +plt.show() diff --git a/tests/conftest.py b/tests/conftest.py index cde943043..1a6d5486d 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -503,6 +503,17 @@ def __init__(self): "awc_wake_exp": 1.2, "awc_wake_denominator": 400 }, + "eddy_viscosity": { + "i_const_1": 0.05, + "i_const_2": 16, + "i_const_3": 0.5, + "i_const_4": 10, + "c_0": 2.0, + "c_1": 1.5, + "k_l": 0.0283, + "k_a": 0.5, + "von_Karman_constant": 0.4 + } }, "wake_turbulence_parameters": { "crespo_hernandez": { diff --git a/tests/reg_tests/eddy_viscosity_regression_test.py b/tests/reg_tests/eddy_viscosity_regression_test.py new file mode 100644 index 000000000..62cd85130 --- /dev/null +++ b/tests/reg_tests/eddy_viscosity_regression_test.py @@ -0,0 +1,483 @@ + +import numpy as np +import pytest + +from floris.core import ( + average_velocity, + axial_induction, + Core, + power, + rotor_effective_velocity, + thrust_coefficient, +) +from tests.conftest import ( + assert_results_arrays, + N_FINDEX, + N_TURBINES, + print_test_values, +) + + +DEBUG = False +VELOCITY_MODEL = "eddy_viscosity" +DEFLECTION_MODEL = "gauss" +COMBINATION_MODEL = "soed" + +baseline = np.array( + [ + # 8 m/s + [ + [7.9736858, 0.7871515, 1753954.4591792, 0.2693224], + [6.6666422, 0.8305317, 1037303.6026684, 0.2941674], + [6.3776334, 0.8436754, 907368.5661344, 0.3023105], + ], + # 9 m/s + [ + [8.9703965, 0.7858774, 2496427.8618358, 0.2686331], + [7.5017741, 0.7976252, 1461746.8459320, 0.2750696], + [7.1929915, 0.8081969, 1288784.8178808, 0.2810234], + ], + # 10 m/s + [ + [9.9671073, 0.7838789, 3417797.0050916, 0.2675559], + [8.3384475, 0.7866918, 2024117.2497181, 0.2690735], + [8.0002374, 0.7871277, 1771343.4133243, 0.2693096], + ], + # 11 m/s + [ + [10.9638180, 0.7565157, 4519404.3072862, 0.2532794], + [9.2200413, 0.7853932, 2723153.7974906, 0.2683716], + [8.8125257, 0.7860809, 2378437.2315890, 0.2687430], + ], + ] +) + +# yawed_baseline = np.array( +# [ +# # 8 m/s +# [ +# [7.9736858, 0.7841561, 1741508.6722008, 0.2671213], +# [6.0816475, 0.8571363, 774296.7271893, 0.3110134], +# [5.5272875, 0.8877222, 579850.4298177, 0.3324606], +# ], +# # 9 m/s +# [ +# [8.9703965, 0.7828869, 2480428.8963141, 0.2664440], +# [6.8472506, 0.8223180, 1118503.0309148, 0.2892383], +# [6.3747452, 0.8438067, 906070.0511419, 0.3023935], +# ], +# # 10 m/s +# [ +# [9.9671073, 0.7808960, 3395681.0032992, 0.2653854], +# [7.6174285, 0.7940006, 1530191.8035935, 0.2730642], +# [7.2119500, 0.8075204, 1299067.3876318, 0.2806375], +# ], +# # 11 m/s +# [ +# [10.9638180, 0.7536370, 4488242.9153943, 0.2513413], +# [8.5159500, 0.7864631, 2156780.3499849, 0.2689497], +# [8.0047998, 0.7871218, 1774753.2988553, 0.2693064], +# ], +# ] +# ) + +full_flow_baseline = np.array( + [ + [ + [ + [7.88772361, 8.0, 8.10178821], + [7.88772361, 8.0, 8.10178821], + [7.88772361, 8.0, 8.10178821], + [7.88772361, 8.0, 8.10178821], + [7.88772361, 8.0, 8.10178821], + ], + [ + [7.88758793, 7.99986028, 8.10164885], + [7.69657254, 7.80314836, 7.9054495 ], + [5.74483323, 5.79320367, 5.90074201], + [7.69657254, 7.80314836, 7.9054495 ], + [7.88758793, 7.99986028, 8.10164885], + ], + [ + [7.8151773 , 7.92629295, 8.02727307], + [7.48541427, 7.58840385, 7.68856061], + [5.67767815, 5.73066043, 5.83176442], + [7.48541427, 7.58840385, 7.68856061], + [7.8151773 , 7.92629295, 8.02727307], + ], + [ + [7.75903294, 7.86931452, 7.969605 ], + [7.38563303, 7.48715077, 7.58607141], + [5.79942256, 5.85940559, 5.95681284], + [7.38563303, 7.48715077, 7.58607141], + [7.75903294, 7.86931452, 7.969605 ], + ], + [ + [7.74436376, 7.85435448, 7.95453772], + [7.37530096, 7.47827906, 7.57545893], + [6.89021982, 6.98370227, 7.07721321], + [7.37530096, 7.47827906, 7.57545893], + [7.74436376, 7.85435448, 7.95453772], + ] + ] + ] +) + +# Note: compare the yawed vs non-yawed results. The upstream turbine +# power should be lower in the yawed case. The following turbine +# powers should higher in the yawed case. + + +def test_regression_tandem(sample_inputs_fixture): + """ + Tandem turbines + """ + sample_inputs_fixture.core["wake"]["model_strings"]["velocity_model"] = VELOCITY_MODEL + sample_inputs_fixture.core["wake"]["model_strings"]["deflection_model"] = DEFLECTION_MODEL + sample_inputs_fixture.core["wake"]["model_strings"]["combination_model"] = COMBINATION_MODEL + + + floris = Core.from_dict(sample_inputs_fixture.core) + floris.initialize_domain() + floris.steady_state_atmospheric_condition() + + n_turbines = floris.farm.n_turbines + n_findex = floris.flow_field.n_findex + + velocities = floris.flow_field.u + turbulence_intensities = floris.flow_field.turbulence_intensity_field + air_density = floris.flow_field.air_density + yaw_angles = floris.farm.yaw_angles + tilt_angles = floris.farm.tilt_angles + power_setpoints = floris.farm.power_setpoints + awc_modes = floris.farm.awc_modes + awc_amplitudes = floris.farm.awc_amplitudes + test_results = np.zeros((n_findex, n_turbines, 4)) + + farm_avg_velocities = average_velocity( + velocities, + ) + farm_cts = thrust_coefficient( + velocities, + turbulence_intensities, + air_density, + yaw_angles, + tilt_angles, + power_setpoints, + awc_modes, + awc_amplitudes, + floris.farm.turbine_thrust_coefficient_functions, + floris.farm.turbine_tilt_interps, + floris.farm.correct_cp_ct_for_tilt, + floris.farm.turbine_type_map, + floris.farm.turbine_power_thrust_tables, + ) + farm_powers = power( + velocities, + turbulence_intensities, + air_density, + floris.farm.turbine_power_functions, + yaw_angles, + tilt_angles, + power_setpoints, + awc_modes, + awc_amplitudes, + floris.farm.turbine_tilt_interps, + floris.farm.turbine_type_map, + floris.farm.turbine_power_thrust_tables, + ) + farm_axial_inductions = axial_induction( + velocities, + turbulence_intensities, + air_density, + yaw_angles, + tilt_angles, + power_setpoints, + awc_modes, + awc_amplitudes, + floris.farm.turbine_axial_induction_functions, + floris.farm.turbine_tilt_interps, + floris.farm.correct_cp_ct_for_tilt, + floris.farm.turbine_type_map, + floris.farm.turbine_power_thrust_tables, + ) + for i in range(n_findex): + for j in range(n_turbines): + test_results[i, j, 0] = farm_avg_velocities[i, j] + test_results[i, j, 1] = farm_cts[i, j] + test_results[i, j, 2] = farm_powers[i, j] + test_results[i, j, 3] = farm_axial_inductions[i, j] + + if DEBUG: + print_test_values( + farm_avg_velocities, + farm_cts, + farm_powers, + farm_axial_inductions, + max_findex_print=4 + ) + + assert_results_arrays(test_results[0:4], baseline) + + +def test_regression_rotation(sample_inputs_fixture): + """ + Turbines in tandem and rotated. + The result from 270 degrees should match the results from 360 degrees. + + Wind from the West (Left) + + ^ + | + y + + 1|1 3 + | + | + | + 0|0 2 + |----------| + 0 1 x-> + + + Wind from the North (Top), rotated + + ^ + | + y + + 1|3 2 + | + | + | + 0|1 0 + |----------| + 0 1 x-> + + In 270, turbines 2 and 3 are waked. In 360, turbines 0 and 2 are waked. + The test compares turbines 2 and 3 with 0 and 2 from 270 and 360. + """ + TURBINE_DIAMETER = 126.0 + + sample_inputs_fixture.core["wake"]["model_strings"]["velocity_model"] = VELOCITY_MODEL + sample_inputs_fixture.core["wake"]["model_strings"]["deflection_model"] = DEFLECTION_MODEL + sample_inputs_fixture.core["wake"]["model_strings"]["combination_model"] = COMBINATION_MODEL + sample_inputs_fixture.core["farm"]["layout_x"] = [ + 0.0, + 0.0, + 5 * TURBINE_DIAMETER, + 5 * TURBINE_DIAMETER, + ] + sample_inputs_fixture.core["farm"]["layout_y"] = [ + 0.0, + 5 * TURBINE_DIAMETER, + 0.0, + 5 * TURBINE_DIAMETER + ] + sample_inputs_fixture.core["flow_field"]["wind_directions"] = [270.0, 360.0] + sample_inputs_fixture.core["flow_field"]["wind_speeds"] = [8.0, 8.0] + sample_inputs_fixture.core["flow_field"]["turbulence_intensities"] = [0.1, 0.1] + + floris = Core.from_dict(sample_inputs_fixture.core) + floris.initialize_domain() + floris.steady_state_atmospheric_condition() + + farm_avg_velocities = average_velocity(floris.flow_field.u) + + t0_270 = farm_avg_velocities[0, 0] # upstream + t1_270 = farm_avg_velocities[0, 1] # upstream + t2_270 = farm_avg_velocities[0, 2] # waked + t3_270 = farm_avg_velocities[0, 3] # waked + + t0_360 = farm_avg_velocities[1, 0] # waked + t1_360 = farm_avg_velocities[1, 1] # upstream + t2_360 = farm_avg_velocities[1, 2] # waked + t3_360 = farm_avg_velocities[1, 3] # upstream + + assert np.allclose(t0_270, t1_360) + assert np.allclose(t1_270, t3_360) + assert np.allclose(t2_270, t0_360) + assert np.allclose(t3_270, t2_360) + + +def test_regression_yaw(sample_inputs_fixture): + """ + Tandem turbines with the upstream turbine yawed + """ + sample_inputs_fixture.core["wake"]["model_strings"]["velocity_model"] = VELOCITY_MODEL + sample_inputs_fixture.core["wake"]["model_strings"]["deflection_model"] = DEFLECTION_MODEL + sample_inputs_fixture.core["wake"]["model_strings"]["combination_model"] = COMBINATION_MODEL + + floris = Core.from_dict(sample_inputs_fixture.core) + + yaw_angles = np.zeros((N_FINDEX, N_TURBINES)) + yaw_angles[:,0] = 5.0 + floris.farm.yaw_angles = yaw_angles + + floris.initialize_domain() + with pytest.raises(NotImplementedError): + floris.steady_state_atmospheric_condition() + # Once implemented, copy code from test_regression_yaw on jensen_jimenez_regression_test.py and + # update yawed_baseline values + +def test_symmetry(sample_inputs_fixture): + """ + This utilizes a 5x5 wind farm with the layout in a regular grid oriented along the cardinal + directions. The wind direction in this test is from 270 degrees, directly aligned with the + columns of the farm. The objective of this test is to check that the solve is symmetric. + """ + + sample_inputs_fixture.core["wake"]["model_strings"]["velocity_model"] = VELOCITY_MODEL + sample_inputs_fixture.core["wake"]["model_strings"]["deflection_model"] = DEFLECTION_MODEL + sample_inputs_fixture.core["wake"]["model_strings"]["combination_model"] = COMBINATION_MODEL + X, Y = np.meshgrid( + 6.0 * 126.0 * np.arange(0, 5, 1), + 6.0 * 126.0 * np.arange(0, 5, 1) + ) + X = X.flatten() + Y = Y.flatten() + + sample_inputs_fixture.core["farm"]["layout_x"] = X + sample_inputs_fixture.core["farm"]["layout_y"] = Y + + floris = Core.from_dict(sample_inputs_fixture.core) + floris.initialize_domain() + floris.steady_state_atmospheric_condition() + + # farm_avg_velocities = average_velocity(floris.flow_field.u) + velocities = floris.flow_field.u + turbulence_intensities = floris.flow_field.turbulence_intensity_field + air_density = floris.flow_field.air_density + yaw_angles = floris.farm.yaw_angles + tilt_angles = floris.farm.tilt_angles + power_setpoints = floris.farm.power_setpoints + awc_modes = floris.farm.awc_modes + awc_amplitudes = floris.farm.awc_amplitudes + + farm_powers = power( + velocities, + turbulence_intensities, + air_density, + floris.farm.turbine_power_functions, + yaw_angles, + tilt_angles, + power_setpoints, + awc_modes, + awc_amplitudes, + floris.farm.turbine_tilt_interps, + floris.farm.turbine_type_map, + floris.farm.turbine_power_thrust_tables, + ) + + rtol = 1e-7 + assert np.allclose(farm_powers[0,0:5], farm_powers[0,20:25], rtol=rtol) # Outer columns + assert np.allclose(farm_powers[0,5:10], farm_powers[0,15:20], rtol=rtol) # Inner columns + +def test_regression_small_grid_rotation(sample_inputs_fixture): + """ + This utilizes a 5x5 wind farm with the layout in a regular grid oriented along the cardinal + directions. The wind direction in this test is from 285 degrees which is slightly north of + west. The objective of this test is to create a case with a very slight rotation of the wind + farm to target the rotation and masking routines. + + Where wake models are masked based on the x-location of a turbine, numerical precision + can cause masking to fail unexpectedly. For example, in the configuration here one of + the turbines has these delta x values; + + [[4.54747351e-13 4.54747351e-13 4.54747351e-13 4.54747351e-13 4.54747351e-13] + [4.54747351e-13 4.54747351e-13 4.54747351e-13 4.54747351e-13 4.54747351e-13] + [4.54747351e-13 4.54747351e-13 4.54747351e-13 4.54747351e-13 4.54747351e-13] + [4.54747351e-13 4.54747351e-13 4.54747351e-13 4.54747351e-13 4.54747351e-13] + [4.54747351e-13 4.54747351e-13 4.54747351e-13 4.54747351e-13 4.54747351e-13]] + + and therefore the masking statement is False when it should be True. This causes the current + turbine to be affected by its own wake. This test requires that at least in this particular + configuration the masking correctly filters grid points. + """ + sample_inputs_fixture.core["wake"]["model_strings"]["velocity_model"] = VELOCITY_MODEL + sample_inputs_fixture.core["wake"]["model_strings"]["deflection_model"] = DEFLECTION_MODEL + sample_inputs_fixture.core["wake"]["model_strings"]["combination_model"] = COMBINATION_MODEL + X, Y = np.meshgrid( + 6.0 * 126.0 * np.arange(0, 5, 1), + 6.0 * 126.0 * np.arange(0, 5, 1) + ) + X = X.flatten() + Y = Y.flatten() + + sample_inputs_fixture.core["farm"]["layout_x"] = X + sample_inputs_fixture.core["farm"]["layout_y"] = Y + + floris = Core.from_dict(sample_inputs_fixture.core) + floris.initialize_domain() + floris.steady_state_atmospheric_condition() + + # farm_avg_velocities = average_velocity(floris.flow_field.u) + velocities = floris.flow_field.u + turbulence_intensities = floris.flow_field.turbulence_intensity_field + air_density = floris.flow_field.air_density + yaw_angles = floris.farm.yaw_angles + tilt_angles = floris.farm.tilt_angles + power_setpoints = floris.farm.power_setpoints + awc_modes = floris.farm.awc_modes + awc_amplitudes = floris.farm.awc_amplitudes + + farm_powers = power( + velocities, + turbulence_intensities, + air_density, + floris.farm.turbine_power_functions, + yaw_angles, + tilt_angles, + power_setpoints, + awc_modes, + awc_amplitudes, + floris.farm.turbine_tilt_interps, + floris.farm.turbine_type_map, + floris.farm.turbine_power_thrust_tables, + ) + + # A "column" is oriented parallel to the wind direction + # Columns 1 - 4 should have the same power profile + # Column 5 leading turbine is completely unwaked + # and the rest of the turbines have a partial wake from their immediate upstream turbine + rtol = 1e-2 # Fails for default rtol=1e-5 + assert np.allclose(farm_powers[8,0:5], farm_powers[8,5:10], rtol=rtol) + assert np.allclose(farm_powers[8,0:5], farm_powers[8,10:15], rtol=rtol) + # The following fails, but it's not clear that it should pass. Setting rtol=1e-2 makes it pass. + # assert np.allclose(farm_powers[8,0:5], farm_powers[8,15:20], rtol=rtol) + assert np.allclose(farm_powers[8,20], farm_powers[8,0], rtol=rtol) + assert np.allclose(farm_powers[8,21], farm_powers[8,21:25], rtol=rtol) + + +def test_full_flow_solver(sample_inputs_fixture): + """ + Full flow solver test with the flow field planar grid. + This requires one wind condition, and the grid is deliberately coarse to allow for + visually comparing results, as needed. + The u-component of velocity is compared, and the array has the shape + (n_findex, n_turbines, n grid points in x, n grid points in y, 3 grid points in z). + """ + + sample_inputs_fixture.core["wake"]["model_strings"]["velocity_model"] = VELOCITY_MODEL + sample_inputs_fixture.core["wake"]["model_strings"]["deflection_model"] = DEFLECTION_MODEL + sample_inputs_fixture.core["wake"]["model_strings"]["combination_model"] = COMBINATION_MODEL + sample_inputs_fixture.core["solver"] = { + "type": "flow_field_planar_grid", + "normal_vector": "z", + "planar_coordinate": sample_inputs_fixture.core["farm"]["turbine_type"][0]["hub_height"], + "flow_field_grid_points": [5, 5], + "flow_field_bounds": [None, None], + } + sample_inputs_fixture.core["flow_field"]["wind_directions"] = [270.0] + sample_inputs_fixture.core["flow_field"]["wind_speeds"] = [8.0] + sample_inputs_fixture.core["flow_field"]["turbulence_intensities"] = [0.1] + + floris = Core.from_dict(sample_inputs_fixture.core) + floris.solve_for_viz() + + velocities = floris.flow_field.u_sorted + + if DEBUG: + print(velocities) + + assert_results_arrays(velocities, full_flow_baseline)