Update multidimensional turbine model (#65)

* Update wind condition broadcast for turbine tests

The inputs changed in conftest but this wasn’t updated

* Update multidimensional turbine module for 4D arrays

* Update multidimensional example API’s

* Unit test bug fix

* Remove a few missed extra dimensions

examples/01_opening_floris_computing_power.py

@@ -42,43 +42,38 @@
 fi.reinitialize(wind_directions=[270.0], wind_speeds=[8.0])
 
 # Set the yaw angles to 0
-yaw_angles = np.zeros([1, 2])  # 1 wind direction / speed, 2 turbines
+yaw_angles = np.zeros([1, 2])  # 1 wind direction and speed, 2 turbines
 fi.calculate_wake(yaw_angles=yaw_angles)
 
 # Get the turbine powers
 turbine_powers = fi.get_turbine_powers() / 1000.0
 
-# TODO what should we call this user/facing?
-print("The turbine power matrix should be of dimensions 1 FINDEX X 2 Turbines")
+print("The turbine power matrix should be of dimensions 1 findex X 2 Turbines")
 print(turbine_powers)
 print("Shape: ", turbine_powers.shape)
 
 # Single wind speed and multiple wind directions
 print("\n========================= Single Wind Direction and Multiple Wind Speeds ===============")
-# Note in v3 FLORIS wind directions and speeds would be expanded to all combinations
-# in v4 the assumption is that each entry wind direction and wind speed corresponds
-# to one condtions and wind directions and wind speeds arrays should be the same length
 
 wind_speeds = np.array([8.0, 9.0, 10.0])
 wind_directions = np.array([270.0, 270.0, 270.0])
 
 fi.reinitialize(wind_speeds=wind_speeds, wind_directions=wind_directions)
-yaw_angles = np.zeros([3, 2])  # 9 wind directions/ speeds, 2 turbines
+yaw_angles = np.zeros([3, 2])  # 3 wind directions/ speeds, 2 turbines
 fi.calculate_wake(yaw_angles=yaw_angles)
 turbine_powers = fi.get_turbine_powers() / 1000.0
-print("The turbine power matrix should be of dimensions 9 FINDEX X 2 Turbines")
+print("The turbine power matrix should be of dimensions 3 findex X 2 Turbines")
 print(turbine_powers)
 print("Shape: ", turbine_powers.shape)
 
 # Multiple wind speeds and multiple wind directions
 print("\n========================= Multiple Wind Directions and Multiple Wind Speeds ============")
 
-# In the case want to consider each combination this needs to be broadcast out in advance
+# To consider each combination, this needs to be broadcast out in advance
 
 wind_speeds = np.tile([8.0, 9.0, 10.0], 3)
 wind_directions = np.repeat([260.0, 270.0, 280.0], 3)
 
-
 fi.reinitialize(wind_directions=wind_directions, wind_speeds=wind_speeds)
 yaw_angles = np.zeros([9, 2])  # 9 wind directions/ speeds, 2 turbines
 fi.calculate_wake(yaw_angles=yaw_angles)

examples/30_multi_dimensional_cp_ct.py

@@ -65,40 +65,42 @@
 print('\n========================= Single Wind Direction and Wind Speed =========================')
 
 # Get the turbine powers assuming 1 wind speed and 1 wind direction
-fi.reinitialize(wind_directions=[270.], wind_speeds=[8.0])
+fi.reinitialize(wind_directions=[270.0], wind_speeds=[8.0])
 
 # Set the yaw angles to 0
-yaw_angles = np.zeros([1,1,2]) # 1 wind direction, 1 wind speed, 2 turbines
+yaw_angles = np.zeros([1, 2]) # 1 wind direction and wind speed, 2 turbines
 fi.calculate_wake(yaw_angles=yaw_angles)
 
 # Get the turbine powers
-turbine_powers = fi.get_turbine_powers_multidim()/1000.
-print('The turbine power matrix should be of dimensions 1 WD X 1 WS X 2 Turbines')
+turbine_powers = fi.get_turbine_powers_multidim() / 1000.0
+print("The turbine power matrix should be of dimensions 1 findex X 2 Turbines")
 print(turbine_powers)
 print("Shape: ",turbine_powers.shape)
 
 # Single wind speed and multiple wind directions
 print('\n========================= Single Wind Direction and Multiple Wind Speeds ===============')
 
-
 wind_speeds = np.array([8.0, 9.0, 10.0])
-fi.reinitialize(wind_speeds=wind_speeds)
-yaw_angles = np.zeros([1,3,2]) # 1 wind direction, 3 wind speeds, 2 turbines
+wind_directions = np.array([270.0, 270.0, 270.0])
+
+fi.reinitialize(wind_speeds=wind_speeds, wind_directions=wind_directions)
+yaw_angles = np.zeros([3, 2])  # 3 wind directions/ speeds, 2 turbines
 fi.calculate_wake(yaw_angles=yaw_angles)
-turbine_powers = fi.get_turbine_powers_multidim()/1000.
-print('The turbine power matrix should be of dimensions 1 WD X 3 WS X 2 Turbines')
+turbine_powers = fi.get_turbine_powers_multidim() / 1000.0
+print("The turbine power matrix should be of dimensions 3 findex X 2 Turbines")
 print(turbine_powers)
 print("Shape: ",turbine_powers.shape)
 
 # Multiple wind speeds and multiple wind directions
 print('\n========================= Multiple Wind Directions and Multiple Wind Speeds ============')
 
-wind_directions = np.array([260., 270., 280.])
-wind_speeds = np.array([8.0, 9.0, 10.0])
+wind_speeds = np.tile([8.0, 9.0, 10.0], 3)
+wind_directions = np.repeat([260.0, 270.0, 280.0], 3)
+
 fi.reinitialize(wind_directions=wind_directions, wind_speeds=wind_speeds)
-yaw_angles = np.zeros([1,3,2]) # 1 wind direction, 3 wind speeds, 2 turbines
+yaw_angles = np.zeros([9, 2])  # 9 wind directions/ speeds, 2 turbines
 fi.calculate_wake(yaw_angles=yaw_angles)
 turbine_powers = fi.get_turbine_powers_multidim()/1000.
-print('The turbine power matrix should be of dimensions 3 WD X 3 WS X 2 Turbines')
+print("The turbine power matrix should be of dimensions 9 WD/WS X 2 Turbines")
 print(turbine_powers)
 print("Shape: ",turbine_powers.shape)

examples/31_multi_dimensional_cp_ct_2Hs.py

@@ -46,9 +46,10 @@
 fi_hs_1.reinitialize(layout_x=[0., 500., 1000.], layout_y=[0., 0., 0.])
 
 # Use a sweep of wind speeds
-wind_speeds = np.arange(5,20,1.)
-fi.reinitialize(wind_directions=[270.], wind_speeds=wind_speeds)
-fi_hs_1.reinitialize(wind_directions=[270.], wind_speeds=wind_speeds)
+wind_speeds = np.arange(5, 20, 1.0)
+wind_directions = 270.0 * np.ones_like(wind_speeds)
+fi.reinitialize(wind_directions=wind_directions, wind_speeds=wind_speeds)
+fi_hs_1.reinitialize(wind_directions=wind_directions, wind_speeds=wind_speeds)
 
 # Calculate wakes with baseline yaw
 fi.calculate_wake()
@@ -63,8 +64,8 @@
 
 for t_idx in range(3):
     ax = axarr[t_idx]
-    ax.plot(wind_speeds, turbine_powers[0,:,t_idx], color='k', label='Hs=3.1 (5)')
-    ax.plot(wind_speeds, turbine_powers_hs_1[0,:,t_idx], color='r', label='Hs=1.0')
+    ax.plot(wind_speeds, turbine_powers[:,t_idx], color='k', label='Hs=3.1 (5)')
+    ax.plot(wind_speeds, turbine_powers_hs_1[:,t_idx], color='r', label='Hs=1.0')
     ax.grid(True)
     ax.set_xlabel('Wind Speed (m/s)')
     ax.set_title(f'Turbine {t_idx}')

floris/simulation/farm.py

@@ -318,40 +318,26 @@ def expand_farm_properties(self, n_findex: int, sorted_coord_indices):
             sorted_coord_indices,
             axis=1
         )
-
-        # TODO: update multidimensional turbine for 4D arrays
         if 'multi_dimensional_cp_ct' in self.turbine_definitions[0].keys() \
             and self.turbine_definitions[0]['multi_dimensional_cp_ct'] is True:
-            wd_dim = np.shape(template_shape)[0]
-            ws_dim = np.shape(template_shape)[1]
-            if wd_dim != 1 | ws_dim != 0:
-                self.turbine_fCts_sorted = np.take_along_axis(
-                    np.reshape(
-                        np.repeat(self.turbine_fCts, wd_dim * ws_dim),
-                        np.shape(template_shape)
-                    ),
-                    sorted_coord_indices,
-                    axis=2 # TODO: This should probably be 1
-                )
-                self.turbine_power_interps_sorted = np.take_along_axis(
-                    np.reshape(
-                        np.repeat(self.turbine_power_interps, wd_dim * ws_dim),
-                        np.shape(template_shape)
-                    ),
-                    sorted_coord_indices,
-                    axis=2 # TODO: This should probably be 1
-                )
-            else:
-                self.turbine_fCts_sorted = np.take_along_axis(
-                    np.reshape(self.turbine_fCts, np.shape(template_shape)),
-                    sorted_coord_indices,
-                    axis=1
-                )
-                self.turbine_power_interps_sorted = np.take_along_axis(
-                    np.reshape(self.turbine_power_interps, np.shape(template_shape)),
-                    sorted_coord_indices,
-                    axis=1
-                )
+            findex_dim = np.shape(template_shape)[0]
+
+            self.turbine_fCts_sorted = np.take_along_axis(
+                np.reshape(
+                    np.repeat(self.turbine_fCts, findex_dim),
+                    np.shape(template_shape)
+                ),
+                sorted_coord_indices,
+                axis=1
+            )
+            self.turbine_power_interps_sorted = np.take_along_axis(
+                np.reshape(
+                    np.repeat(self.turbine_power_interps, findex_dim),
+                    np.shape(template_shape)
+                ),
+                sorted_coord_indices,
+                axis=1
+            )
         self.rotor_diameters_sorted = np.take_along_axis(
             self.rotor_diameters * template_shape,
             sorted_coord_indices,

floris/simulation/solver.py

@@ -1490,24 +1490,23 @@ def sequential_multidim_solver(
     w_wake = np.zeros_like(flow_field.w_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))
+        flow_field.turbulence_intensity * 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]
+        u_i = flow_field.u_sorted[:, i:i+1]
+        v_i = flow_field.v_sorted[:, i:i+1]
 
         ct_i = Ct_multidim(
             velocities=flow_field.u_sorted,
@@ -1524,7 +1523,7 @@ def sequential_multidim_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_multidim(
             velocities=flow_field.u_sorted,
             yaw_angle=farm.yaw_angles_sorted,
@@ -1540,12 +1539,12 @@ def sequential_multidim_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
@@ -1555,8 +1554,8 @@ def sequential_multidim_solver(
                 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,
@@ -1600,12 +1599,12 @@ def sequential_multidim_solver(
                 u_i,
                 turbulence_intensity_i,
                 v_i,
-                flow_field.w_sorted[:, :, i:i+1],
-                v_wake[:, :, i:i+1],
-                w_wake[:, :, i:i+1],
+                flow_field.w_sorted[:, i:i+1],
+                v_wake[:, i:i+1],
+                w_wake[:, i:i+1],
             )
             gch_gain = 2
-            turbine_turbulence_intensity[:, :, i:i+1] = turbulence_intensity_i + gch_gain * I_mixing
+            turbine_turbulence_intensity[:, i:i+1] = turbulence_intensity_i + gch_gain * I_mixing
 
         # NOTE: exponential
         velocity_deficit = model_manager.velocity_model.function(
@@ -1637,10 +1636,10 @@ def sequential_multidim_solver(
 
         # 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
@@ -1665,5 +1664,5 @@ def sequential_multidim_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)
-    )[:, :, :, None, None]
+        axis=(2,3)
+    )[:, :, None, None]

floris/simulation/turbine.py

@@ -128,12 +128,12 @@ def rotor_effective_velocity(
 
     # Down-select inputs if ix_filter is given
     if ix_filter is not None:
-        velocities = velocities[:, :, ix_filter]
-        yaw_angle = yaw_angle[:, :, ix_filter]
-        tilt_angle = tilt_angle[:, :, ix_filter]
-        ref_tilt_cp_ct = ref_tilt_cp_ct[:, :, ix_filter]
-        pP = pP[:, :, ix_filter]
-        pT = pT[:, :, ix_filter]
+        velocities = velocities[:, ix_filter]
+        yaw_angle = yaw_angle[:, ix_filter]
+        tilt_angle = tilt_angle[:, ix_filter]
+        ref_tilt_cp_ct = ref_tilt_cp_ct[:, ix_filter]
+        pP = pP[:, ix_filter]
+        pT = pT[:, ix_filter]
         turbine_type_map = turbine_type_map[:, ix_filter]
 
     # Compute the rotor effective velocity adjusting for air density
@@ -204,7 +204,7 @@ def power(
 
     # Down-select inputs if ix_filter is given
     if ix_filter is not None:
-        rotor_effective_velocities = rotor_effective_velocities[:, :, ix_filter]
+        rotor_effective_velocities = rotor_effective_velocities[:, ix_filter]
         turbine_type_map = turbine_type_map[:, ix_filter]
 
     # Loop over each turbine type given to get power for all turbines