Add wake deflection to the TurbOPark wake deficit model (#439)

* add deflection to the turbopark model

* Add yawed reg test

* Add warning when using wake deflection with the TurbOPark model

* simplifying check for yaw angles in turbopark deflection

Co-authored-by: Rafael M Mudafort <rafmudaf@gmail.com>

floris/simulation/solver.py

@@ -703,6 +703,7 @@ def turbopark_solver(farm: Farm, flow_field: FlowField, grid: TurbineGrid, model
     w_wake = np.zeros_like(flow_field.w_initial_sorted)
     shape = (farm.n_turbines,) + np.shape(flow_field.u_initial_sorted)
     velocity_deficit = np.zeros(shape)
+    deflection_field = np.zeros_like(flow_field.u_initial_sorted)
 
     turbine_turbulence_intensity = flow_field.turbulence_intensity * np.ones((flow_field.n_wind_directions, flow_field.n_wind_speeds, farm.n_turbines, 1, 1))
     ambient_turbulence_intensity = flow_field.turbulence_intensity
@@ -769,15 +770,37 @@ def turbopark_solver(farm: Farm, flow_field: FlowField, grid: TurbineGrid, model
 
         # Model calculations
         # NOTE: exponential
-        deflection_field = model_manager.deflection_model.function(
-            x_i,
-            y_i,
-            effective_yaw_i,
-            turbulence_intensity_i,
-            ct_i,
-            rotor_diameter_i,
-            **deflection_model_args
-        )
+        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 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]
+
+                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,
+                    fCt=farm.turbine_fCts,
+                    turbine_type_map=farm.turbine_type_map_sorted,
+                    ix_filter=[ii]
+                )
+                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,
+                    y_ii,
+                    yaw_ii,
+                    turbulence_intensity_ii,
+                    ct_ii,
+                    rotor_diameter_ii,
+                    **deflection_model_args
+                )
+
+                deflection_field[:,:,ii:ii+1,:,:] = deflection_field_ii[:,:,i:i+1,:,:]
 
         if model_manager.enable_transverse_velocities:
             v_wake, w_wake = calculate_transverse_velocity(
@@ -816,6 +839,7 @@ def turbopark_solver(farm: Farm, flow_field: FlowField, grid: TurbineGrid, model
             rotor_diameter_i,
             farm.rotor_diameters_sorted[:, :, :, None, None],
             i,
+            deflection_field,
             **deficit_model_args
         )
 

floris/simulation/wake_velocity/turbopark.py

@@ -72,6 +72,7 @@ def function(
         rotor_diameter_i: np.ndarray,
         rotor_diameters: np.ndarray,
         i: int,
+        deflection_field: np.ndarray,
         # enforces the use of the below as keyword arguments and adherence to the
         # unpacking of the results from prepare_function()
         *,
@@ -89,8 +90,8 @@ def function(
         x_dist = (x_i - x) * downstream_mask / rotor_diameters
 
         # Radial distance between turbine i and the centerlines of wakes from all real/image turbines
-        r_dist = np.sqrt((y_i - y) ** 2 + (z_i - z) ** 2)
-        r_dist_image = np.sqrt((y_i - y) ** 2 + (z_i - (-z)) ** 2)
+        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
 

floris/tools/floris_interface.py

@@ -118,12 +118,7 @@ def calculate_wake(
         # )
 
         # TODO decide where to handle this sign issue
-        if (yaw_angles is not None) and not (np.all(yaw_angles == 0.0)):
-            if self.floris.wake.model_strings["velocity_model"] == "turbopark":
-                # TODO: Implement wake steering for the TurbOPark model
-                raise ValueError(
-                    "Non-zero yaw angles given and for TurbOPark model; wake steering with this model is not yet implemented."
-                )
+        if (yaw_angles is not None) and not (np.all(yaw_angles==0.)):
             self.floris.farm.yaw_angles = yaw_angles
 
         # Initialize solution space
@@ -149,12 +144,7 @@ def calculate_no_wake(
         """
 
         # TODO decide where to handle this sign issue
-        if (yaw_angles is not None) and not (np.all(yaw_angles == 0.0)):
-            if self.floris.wake.model_strings["velocity_model"] == "turbopark":
-                # TODO: Implement wake steering for the TurbOPark model
-                raise ValueError(
-                    "Non-zero yaw angles given and for TurbOPark model; wake steering with this model is not yet implemented."
-                )
+        if (yaw_angles is not None) and not (np.all(yaw_angles==0.)):
             self.floris.farm.yaw_angles = yaw_angles
 
         # Initialize solution space

tests/reg_tests/cumulative_curl_regression_test.py

@@ -18,7 +18,7 @@
 from floris.simulation import Ct, power, axial_induction, average_velocity
 from tests.conftest import N_TURBINES, N_WIND_DIRECTIONS, N_WIND_SPEEDS, print_test_values, assert_results_arrays
 
-DEBUG = True
+DEBUG = False
 VELOCITY_MODEL = "cc"
 DEFLECTION_MODEL = "gauss"
 

tests/reg_tests/turbopark_regression_test.py

@@ -18,7 +18,7 @@
 from floris.simulation import Ct, power, axial_induction, average_velocity
 from tests.conftest import N_TURBINES, N_WIND_DIRECTIONS, N_WIND_SPEEDS, print_test_values, assert_results_arrays
 
-DEBUG = True
+DEBUG = False
 VELOCITY_MODEL = "turbopark"
 DEFLECTION_MODEL = "gauss"
 COMBINATION_MODEL = "fls"
@@ -53,6 +53,35 @@
 )
 
 
+yawed_baseline = np.array(
+    [
+        # 8 m/s
+        [
+            [7.9803783, 0.7605249, 1683956.5765064, 0.2548147],
+            [5.9926862, 0.8357973, 704869.1763857, 0.2973903],
+            [5.3145419, 0.8725432, 479691.1339821, 0.3214945],
+        ],
+        # 9 m/s
+        [
+            [8.9779256, 0.7596713, 2397236.5542849, 0.2543815],
+            [6.7429885, 0.8025994, 1023094.6963579, 0.2778511],
+            [5.9836502, 0.8362597, 701734.6626599, 0.2976758],
+        ],
+        # 10 m/s
+        [
+            [9.9754729, 0.7499157, 3283591.8023665, 0.2494847],
+            [7.5085974, 0.7756254, 1412651.4697014, 0.2631590],
+            [6.6781823, 0.8052071, 992930.8979929, 0.2793232],
+        ],
+        # 11 m/s
+        [
+            [10.9730201, 0.7276532, 4344222.0129382, 0.2386508],
+            [8.3071319, 0.7620861, 1916104.8725891, 0.2561179],
+            [7.3875052, 0.7795398, 1347926.7384587, 0.2652341],
+        ],
+    ]
+)
+
 # 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.
@@ -199,6 +228,72 @@ def test_regression_rotation(sample_inputs_fixture):
     assert np.allclose(t3_270, t2_360)
 
 
+def test_regression_yaw(sample_inputs_fixture):
+    """
+    Tandem turbines with the upstream turbine yawed
+    """
+    sample_inputs_fixture.floris["wake"]["model_strings"]["velocity_model"] = VELOCITY_MODEL
+    sample_inputs_fixture.floris["wake"]["model_strings"]["deflection_model"] = DEFLECTION_MODEL
+
+    floris = Floris.from_dict(sample_inputs_fixture.floris)
+
+    yaw_angles = np.zeros((N_WIND_DIRECTIONS, N_WIND_SPEEDS, N_TURBINES))
+    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
+
+    velocities = floris.flow_field.u
+    yaw_angles = floris.farm.yaw_angles
+    test_results = np.zeros((n_wind_directions, n_wind_speeds, n_turbines, 4))
+
+    farm_avg_velocities = average_velocity(
+        velocities,
+    )
+    farm_cts = Ct(
+        velocities,
+        yaw_angles,
+        floris.farm.turbine_fCts,
+        floris.farm.turbine_type_map,
+    )
+    farm_powers = power(
+        floris.flow_field.air_density,
+        velocities,
+        yaw_angles,
+        floris.farm.pPs,
+        floris.farm.turbine_power_interps,
+        floris.farm.turbine_type_map,
+    )
+    farm_axial_inductions = axial_induction(
+        velocities,
+        yaw_angles,
+        floris.farm.turbine_fCts,
+        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]
+
+    if DEBUG:
+        print_test_values(
+            farm_avg_velocities,
+            farm_cts,
+            farm_powers,
+            farm_axial_inductions,
+        )
+
+    assert_results_arrays(test_results[0], yawed_baseline)
+
+
 def test_regression_small_grid_rotation(sample_inputs_fixture):
     """
     Where wake models are masked based on the x-location of a turbine, numerical precision