Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

4D update for turbopark #63

Merged
merged 5 commits into from
Dec 14, 2023
Merged

4D update for turbopark #63

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
1262908
update turbopark solver for 4D
bayc Dec 14, 2023
128ee82
update turbopark model for 4D
bayc Dec 14, 2023
b094607
update turbopark regression test for 4D
bayc Dec 14, 2023
df02c92
Merge branch '4d' into pr/bayc/63
rafmudaf Dec 14, 2023
2ac9fec
Update regression test API
rafmudaf Dec 14, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
65 changes: 33 additions & 32 deletions floris/simulation/solver.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -874,19 +874,19 @@ def turbopark_solver(

turbine_turbulence_intensity = (
flow_field.turbulence_intensity
* np.ones((flow_field.n_wind_directions, flow_field.n_wind_speeds, farm.n_turbines, 1, 1))
* np.ones((flow_field.n_findex, farm.n_turbines, 1, 1))
)
ambient_turbulence_intensity = flow_field.turbulence_intensity

# Calculate the velocity deficit sequentially from upstream to downstream turbines
for i in range(grid.n_turbines):
# Get the current turbine quantities
x_i = np.mean(grid.x_sorted[:, :, i:i+1], axis=(3, 4))
x_i = x_i[:, :, :, None, None]
y_i = np.mean(grid.y_sorted[:, :, i:i+1], axis=(3, 4))
y_i = y_i[:, :, :, None, None]
z_i = np.mean(grid.z_sorted[:, :, i:i+1], axis=(3, 4))
z_i = z_i[:, :, :, None, None]
x_i = np.mean(grid.x_sorted[:, i:i+1], axis=(2, 3))
x_i = x_i[:, :, None, None]
y_i = np.mean(grid.y_sorted[:, i:i+1], axis=(2, 3))
y_i = y_i[:, :, None, None]
z_i = np.mean(grid.z_sorted[:, i:i+1], axis=(2, 3))
z_i = z_i[:, :, None, None]

u_i = flow_field.u_sorted[:, :, i:i+1]
v_i = flow_field.v_sorted[:, :, i:i+1]
Expand Down Expand Up @@ -919,7 +919,7 @@ def turbopark_solver(
)
# Since we are filtering for the i'th turbine in the Ct function,
# get the first index here (0:1)
ct_i = ct_i[:, :, 0:1, None, None]
ct_i = ct_i[:, 0:1, None, None]
axial_induction_i = axial_induction(
velocities=flow_field.u_sorted,
yaw_angle=farm.yaw_angles_sorted,
Expand All @@ -935,23 +935,24 @@ def turbopark_solver(
)
# Since we are filtering for the i'th turbine in the axial induction function,
# get the first index here (0:1)
axial_induction_i = axial_induction_i[:, :, 0:1, None, None]
turbulence_intensity_i = turbine_turbulence_intensity[:, :, i:i+1]
yaw_angle_i = farm.yaw_angles_sorted[:, :, i:i+1, None, None]
hub_height_i = farm.hub_heights_sorted[:, :, i:i+1, None, None]
rotor_diameter_i = farm.rotor_diameters_sorted[:, :, i:i+1, None, None]
TSR_i = farm.TSRs_sorted[:, :, i:i+1, None, None]
axial_induction_i = axial_induction_i[:, 0:1, None, None]
turbulence_intensity_i = turbine_turbulence_intensity[:, i:i+1]
yaw_angle_i = farm.yaw_angles_sorted[:, i:i+1, None, None]
hub_height_i = farm.hub_heights_sorted[:, i:i+1, None, None]
rotor_diameter_i = farm.rotor_diameters_sorted[:, i:i+1, None, None]
TSR_i = farm.TSRs_sorted[:, i:i+1, None, None]

effective_yaw_i = np.zeros_like(yaw_angle_i)
effective_yaw_i += yaw_angle_i


if model_manager.enable_secondary_steering:
added_yaw = wake_added_yaw(
u_i,
v_i,
flow_field.u_initial_sorted,
grid.y_sorted[:, :, i:i+1] - y_i,
grid.z_sorted[:, :, i:i+1],
grid.y_sorted[:, i:i+1] - y_i,
grid.z_sorted[:, i:i+1],
rotor_diameter_i,
hub_height_i,
ct_i,
Expand All @@ -965,18 +966,18 @@ def turbopark_solver(
# NOTE: exponential
if not np.all(farm.yaw_angles_sorted):
model_manager.deflection_model.logger.warning(
"WARNING: Deflection with the TurbOPark model has not been fully validated."
"This is an initial implementation, and we advise you use at your own risk"
"WARNING: Deflection with the TurbOPark model has not been fully validated. "
"This is an initial implementation, and we advise you use at your own risk "
"and perform a thorough examination of the results."
)
for ii in range(i):
x_ii = np.mean(grid.x_sorted[:, :, ii:ii+1], axis=(3, 4))
x_ii = x_ii[:, :, :, None, None]
y_ii = np.mean(grid.y_sorted[:, :, ii:ii+1], axis=(3, 4))
y_ii = y_ii[:, :, :, None, None]
x_ii = np.mean(grid.x_sorted[:, ii:ii+1], axis=(2, 3))
x_ii = x_ii[:, :, None, None]
y_ii = np.mean(grid.y_sorted[:, ii:ii+1], axis=(2, 3))
y_ii = y_ii[:, :, None, None]

yaw_ii = farm.yaw_angles_sorted[:, :, ii:ii+1, None, None]
turbulence_intensity_ii = turbine_turbulence_intensity[:, :, ii:ii+1]
yaw_ii = farm.yaw_angles_sorted[:, ii:ii+1, None, None]
turbulence_intensity_ii = turbine_turbulence_intensity[:, ii:ii+1]
ct_ii = Ct(
velocities=flow_field.u_sorted,
yaw_angle=farm.yaw_angles_sorted,
Expand All @@ -990,8 +991,8 @@ def turbopark_solver(
average_method=grid.average_method,
cubature_weights=grid.cubature_weights
)
ct_ii = ct_ii[:, :, 0:1, None, None]
rotor_diameter_ii = farm.rotor_diameters_sorted[:, :, ii:ii+1, None, None]
ct_ii = ct_ii[:, 0:1, None, None]
rotor_diameter_ii = farm.rotor_diameters_sorted[:, ii:ii+1, None, None]

deflection_field_ii = model_manager.deflection_model.function(
x_ii,
Expand All @@ -1003,7 +1004,7 @@ def turbopark_solver(
**deflection_model_args,
)

deflection_field[:, :, ii:ii+1, :, :] = deflection_field_ii[:, :, i:i+1, :, :]
deflection_field[:, ii:ii+1, :, :] = deflection_field_ii[:, i:i+1, :, :]

if model_manager.enable_transverse_velocities:
v_wake, w_wake = calculate_transverse_velocity(
Expand Down Expand Up @@ -1040,9 +1041,9 @@ def turbopark_solver(
y_i,
z_i,
turbine_turbulence_intensity,
Cts[:, :, :, None, None],
Cts[:, :, None, None],
rotor_diameter_i,
farm.rotor_diameters_sorted[:, :, :, None, None],
farm.rotor_diameters_sorted[:, :, None, None],
i,
deflection_field,
**deficit_model_args,
Expand All @@ -1066,10 +1067,10 @@ def turbopark_solver(
# turbines; could use WAT_upstream
# Calculate wake overlap for wake-added turbulence (WAT)
area_overlap = (
np.sum(velocity_deficit * flow_field.u_initial_sorted > 0.05, axis=(3, 4))
np.sum(velocity_deficit * flow_field.u_initial_sorted > 0.05, axis=(2, 3))
/ (grid.grid_resolution * grid.grid_resolution)
)
area_overlap = area_overlap[:, :, :, None, None]
area_overlap = area_overlap[:, :, None, None]

# Modify wake added turbulence by wake area overlap
downstream_influence_length = 15 * rotor_diameter_i
Expand All @@ -1094,7 +1095,7 @@ def turbopark_solver(
flow_field.turbulence_intensity_field_sorted = turbine_turbulence_intensity
flow_field.turbulence_intensity_field_sorted_avg = np.mean(
turbine_turbulence_intensity,
axis=(3,4)
axis=(2, 3)
)


Expand Down
6 changes: 3 additions & 3 deletions floris/simulation/wake_velocity/turbopark.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -109,7 +109,7 @@ def function(
r_dist = np.sqrt((y_i - (y + deflection_field)) ** 2 + (z_i - z) ** 2)
r_dist_image = np.sqrt((y_i - (y + deflection_field)) ** 2 + (z_i - (-z)) ** 2)

Cts[:, :, i:, :, :] = 0.00001
Cts[:, i:, :, :] = 0.00001

# Characteristic wake widths from all turbines relative to turbine i
dw = characteristic_wake_width(x_dist, ambient_turbulence_intensity, Cts, self.A)
Expand Down Expand Up @@ -137,9 +137,9 @@ def function(
delta_image = C * wtg_overlapping * self.overlap_gauss_interp(
(r_dist_image / sigma, rotor_diameter_i / 2 / sigma)
)
delta = np.concatenate((delta_real, delta_image), axis=2)
delta = np.concatenate((delta_real, delta_image), axis=1)

delta_total[:, :, i, :, :] = np.sqrt(np.sum(np.nan_to_num(delta) ** 2, axis=2))
delta_total[:, i, :, :] = np.sqrt(np.sum(np.nan_to_num(delta) ** 2, axis=1))

return delta_total

Expand Down
92 changes: 39 additions & 53 deletions tests/reg_tests/turbopark_regression_test.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -112,17 +112,12 @@ def test_regression_tandem(sample_inputs_fixture):
floris.steady_state_atmospheric_condition()

n_turbines = floris.farm.n_turbines
n_wind_speeds = floris.flow_field.n_wind_speeds
n_wind_directions = floris.flow_field.n_wind_directions
n_findex = floris.flow_field.n_findex

velocities = floris.flow_field.u
yaw_angles = floris.farm.yaw_angles
tilt_angles = floris.farm.tilt_angles
ref_tilt_cp_cts = (
np.ones((n_wind_directions, n_wind_speeds, n_turbines))
* floris.farm.ref_tilt_cp_cts
)
test_results = np.zeros((n_wind_directions, n_wind_speeds, n_turbines, 4))
test_results = np.zeros((n_findex, n_turbines, 4))

farm_avg_velocities = average_velocity(
velocities,
Expand All @@ -133,7 +128,7 @@ def test_regression_tandem(sample_inputs_fixture):
velocities,
yaw_angles,
tilt_angles,
ref_tilt_cp_cts,
floris.farm.ref_tilt_cp_cts,
floris.farm.pPs,
floris.farm.pTs,
floris.farm.turbine_tilt_interps,
Expand All @@ -144,7 +139,7 @@ def test_regression_tandem(sample_inputs_fixture):
velocities,
yaw_angles,
tilt_angles,
ref_tilt_cp_cts,
floris.farm.ref_tilt_cp_cts,
floris.farm.turbine_fCts,
floris.farm.turbine_tilt_interps,
floris.farm.correct_cp_ct_for_tilt,
Expand All @@ -160,19 +155,18 @@ def test_regression_tandem(sample_inputs_fixture):
velocities,
yaw_angles,
tilt_angles,
ref_tilt_cp_cts,
floris.farm.ref_tilt_cp_cts,
floris.farm.turbine_fCts,
floris.farm.turbine_tilt_interps,
floris.farm.correct_cp_ct_for_tilt,
floris.farm.turbine_type_map,
)
for i in range(n_wind_directions):
for j in range(n_wind_speeds):
for k in range(n_turbines):
test_results[i, j, k, 0] = farm_avg_velocities[i, j, k]
test_results[i, j, k, 1] = farm_cts[i, j, k]
test_results[i, j, k, 2] = farm_powers[i, j, k]
test_results[i, j, k, 3] = farm_axial_inductions[i, j, k]
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(
Expand All @@ -182,7 +176,7 @@ def test_regression_tandem(sample_inputs_fixture):
farm_axial_inductions,
)

assert_results_arrays(test_results[0], baseline)
assert_results_arrays(test_results[0:4], baseline)


def test_regression_rotation(sample_inputs_fixture):
Expand Down Expand Up @@ -248,15 +242,15 @@ def test_regression_rotation(sample_inputs_fixture):

farm_avg_velocities = average_velocity(floris.flow_field.u)

t0_270 = farm_avg_velocities[0, 0, 0] # upstream
t1_270 = farm_avg_velocities[0, 0, 1] # upstream
t2_270 = farm_avg_velocities[0, 0, 2] # waked
t3_270 = farm_avg_velocities[0, 0, 3] # waked
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, 0] # waked
t1_360 = farm_avg_velocities[1, 0, 1] # upstream
t2_360 = farm_avg_velocities[1, 0, 2] # waked
t3_360 = farm_avg_velocities[1, 0, 3] # upstream
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)
Expand All @@ -274,24 +268,19 @@ def test_regression_yaw(sample_inputs_fixture):
floris = Floris.from_dict(sample_inputs_fixture.floris)

yaw_angles = np.zeros((N_FINDEX, N_TURBINES))
yaw_angles[:,:,0] = 5.0
yaw_angles[:,0] = 5.0
floris.farm.yaw_angles = yaw_angles

floris.initialize_domain()
floris.steady_state_atmospheric_condition()

n_turbines = floris.farm.n_turbines
n_wind_speeds = floris.flow_field.n_wind_speeds
n_wind_directions = floris.flow_field.n_wind_directions
n_findex = floris.flow_field.n_findex

velocities = floris.flow_field.u
yaw_angles = floris.farm.yaw_angles
tilt_angles = floris.farm.tilt_angles
ref_tilt_cp_cts = (
np.ones((n_wind_directions, n_wind_speeds, n_turbines))
* floris.farm.ref_tilt_cp_cts
)
test_results = np.zeros((n_wind_directions, n_wind_speeds, n_turbines, 4))
test_results = np.zeros((n_findex, n_turbines, 4))

farm_avg_velocities = average_velocity(
velocities,
Expand All @@ -302,7 +291,7 @@ def test_regression_yaw(sample_inputs_fixture):
velocities,
yaw_angles,
tilt_angles,
ref_tilt_cp_cts,
floris.farm.ref_tilt_cp_cts,
floris.farm.pPs,
floris.farm.pTs,
floris.farm.turbine_tilt_interps,
Expand All @@ -313,7 +302,7 @@ def test_regression_yaw(sample_inputs_fixture):
velocities,
yaw_angles,
tilt_angles,
ref_tilt_cp_cts,
floris.farm.ref_tilt_cp_cts,
floris.farm.turbine_fCts,
floris.farm.turbine_tilt_interps,
floris.farm.correct_cp_ct_for_tilt,
Expand All @@ -329,19 +318,18 @@ def test_regression_yaw(sample_inputs_fixture):
velocities,
yaw_angles,
tilt_angles,
ref_tilt_cp_cts,
floris.farm.ref_tilt_cp_cts,
floris.farm.turbine_fCts,
floris.farm.turbine_tilt_interps,
floris.farm.correct_cp_ct_for_tilt,
floris.farm.turbine_type_map,
)
for i in range(n_wind_directions):
for j in range(n_wind_speeds):
for k in range(n_turbines):
test_results[i, j, k, 0] = farm_avg_velocities[i, j, k]
test_results[i, j, k, 1] = farm_cts[i, j, k]
test_results[i, j, k, 2] = farm_powers[i, j, k]
test_results[i, j, k, 3] = farm_axial_inductions[i, j, k]
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(
Expand All @@ -351,8 +339,7 @@ def test_regression_yaw(sample_inputs_fixture):
farm_axial_inductions,
)

assert_results_arrays(test_results[0], yawed_baseline)

assert_results_arrays(test_results[0:4], yawed_baseline)

def test_regression_small_grid_rotation(sample_inputs_fixture):
"""
Expand Down Expand Up @@ -391,15 +378,14 @@ def test_regression_small_grid_rotation(sample_inputs_fixture):
velocities = floris.flow_field.u
yaw_angles = floris.farm.yaw_angles
tilt_angles = floris.farm.tilt_angles
ref_tilt_cp_cts = np.ones((1, 1, len(X))) * floris.farm.ref_tilt_cp_cts

farm_eff_velocities = rotor_effective_velocity(
floris.flow_field.air_density,
floris.farm.ref_density_cp_cts,
velocities,
yaw_angles,
tilt_angles,
ref_tilt_cp_cts,
floris.farm.ref_tilt_cp_cts,
floris.farm.pPs,
floris.farm.pTs,
floris.farm.turbine_tilt_interps,
Expand All @@ -417,8 +403,8 @@ def test_regression_small_grid_rotation(sample_inputs_fixture):
# 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
assert np.allclose(farm_powers[2,0,0:5], farm_powers[2,0,5:10])
assert np.allclose(farm_powers[2,0,0:5], farm_powers[2,0,10:15])
assert np.allclose(farm_powers[2,0,0:5], farm_powers[2,0,15:20])
assert np.allclose(farm_powers[2,0,20], farm_powers[2,0,0])
assert np.allclose(farm_powers[2,0,21], farm_powers[2,0,21:25])
assert np.allclose(farm_powers[8,0:5], farm_powers[8,5:10])
assert np.allclose(farm_powers[8,0:5], farm_powers[8,10:15])
assert np.allclose(farm_powers[8,0:5], farm_powers[8,15:20])
assert np.allclose(farm_powers[8,20], farm_powers[8,0])
assert np.allclose(farm_powers[8,21], farm_powers[8,21:25])
Loading