From 12629086eebdc44c71f2af2676ab896eafc36335 Mon Sep 17 00:00:00 2001 From: bayc Date: Thu, 14 Dec 2023 09:06:15 -0700 Subject: [PATCH 1/4] update turbopark solver for 4D --- floris/simulation/solver.py | 65 +++++++++++++++++++------------------ 1 file changed, 33 insertions(+), 32 deletions(-) diff --git a/floris/simulation/solver.py b/floris/simulation/solver.py index c47b247d0..ab0502160 100644 --- a/floris/simulation/solver.py +++ b/floris/simulation/solver.py @@ -876,19 +876,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] @@ -921,7 +921,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, @@ -937,23 +937,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, @@ -967,18 +968,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, @@ -992,8 +993,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, @@ -1005,7 +1006,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( @@ -1042,9 +1043,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, @@ -1068,10 +1069,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 @@ -1096,7 +1097,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) ) From 128ee8231529e4d059ff1bbe0982c8007fe5ecdb Mon Sep 17 00:00:00 2001 From: bayc Date: Thu, 14 Dec 2023 09:06:30 -0700 Subject: [PATCH 2/4] update turbopark model for 4D --- floris/simulation/wake_velocity/turbopark.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/floris/simulation/wake_velocity/turbopark.py b/floris/simulation/wake_velocity/turbopark.py index cf0443347..0b52c0476 100644 --- a/floris/simulation/wake_velocity/turbopark.py +++ b/floris/simulation/wake_velocity/turbopark.py @@ -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) @@ -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 From b094607b4b98e3749c5c0a0211b43e97841b8386 Mon Sep 17 00:00:00 2001 From: bayc Date: Thu, 14 Dec 2023 09:08:01 -0700 Subject: [PATCH 3/4] update turbopark regression test for 4D --- tests/reg_tests/turbopark_regression_test.py | 24 +++++++++----------- 1 file changed, 11 insertions(+), 13 deletions(-) diff --git a/tests/reg_tests/turbopark_regression_test.py b/tests/reg_tests/turbopark_regression_test.py index 6e34c5096..d5f9c8f41 100644 --- a/tests/reg_tests/turbopark_regression_test.py +++ b/tests/reg_tests/turbopark_regression_test.py @@ -24,8 +24,8 @@ ) from tests.conftest import ( assert_results_arrays, - N_TURBINES, N_FINDEX, + N_TURBINES, print_test_values, ) @@ -112,17 +112,16 @@ 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)) + np.ones((n_findex, 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, @@ -166,13 +165,12 @@ def test_regression_tandem(sample_inputs_fixture): 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( @@ -182,7 +180,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): From 2ac9fecacc23300a922f8deb7c79f24e17c5f69b Mon Sep 17 00:00:00 2001 From: Rafael M Mudafort Date: Thu, 14 Dec 2023 11:48:33 -0600 Subject: [PATCH 4/4] Update regression test API --- tests/reg_tests/turbopark_regression_test.py | 72 ++++++++------------ 1 file changed, 30 insertions(+), 42 deletions(-) diff --git a/tests/reg_tests/turbopark_regression_test.py b/tests/reg_tests/turbopark_regression_test.py index 19b7f9003..5d138cdc3 100644 --- a/tests/reg_tests/turbopark_regression_test.py +++ b/tests/reg_tests/turbopark_regression_test.py @@ -117,10 +117,6 @@ def test_regression_tandem(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((n_findex, n_turbines)) - * floris.farm.ref_tilt_cp_cts - ) test_results = np.zeros((n_findex, n_turbines, 4)) farm_avg_velocities = average_velocity( @@ -132,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, @@ -143,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, @@ -159,7 +155,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, @@ -246,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) @@ -272,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, @@ -300,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, @@ -311,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, @@ -327,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( @@ -349,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): """ @@ -389,7 +378,6 @@ 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, @@ -397,7 +385,7 @@ def test_regression_small_grid_rotation(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, @@ -415,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])