diff --git a/floris/simulation/solver.py b/floris/simulation/solver.py index 942f427a1..0291b6945 100644 --- a/floris/simulation/solver.py +++ b/floris/simulation/solver.py @@ -1188,28 +1188,33 @@ def empirical_gauss_solver( v_wake = np.zeros_like(flow_field.v_initial_sorted) w_wake = np.zeros_like(flow_field.w_initial_sorted) - x_locs = np.mean(grid.x_sorted, axis=(3, 4))[:,:,:,None] - downstream_distance_D = x_locs - np.transpose(x_locs, axes=(0,1,3,2)) + x_locs = np.mean(grid.x_sorted, axis=(2, 3))[:,:,None] + downstream_distance_D = x_locs - np.transpose(x_locs, axes=(0,2,1)) downstream_distance_D = downstream_distance_D / \ - np.repeat(farm.rotor_diameters_sorted[:,:,:,None], grid.n_turbines, axis=-1) + np.repeat(farm.rotor_diameters_sorted[:,:,None], grid.n_turbines, axis=-1) downstream_distance_D = np.maximum(downstream_distance_D, 0.1) # For ease - mixing_factor = np.zeros_like(downstream_distance_D) - mixing_factor[:,:,:,:] = model_manager.turbulence_model.atmospheric_ti_gain*\ - flow_field.turbulence_intensity*np.eye(grid.n_turbines) + # Initialize the mixing factor model using TI if specified + initial_mixing_factor = model_manager.turbulence_model.atmospheric_ti_gain*\ + flow_field.turbulence_intensity*np.eye(grid.n_turbines) + mixing_factor = np.repeat( + initial_mixing_factor[None,:,:], + flow_field.n_findex, + axis=0 + ) # 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] - flow_field.u_sorted[:, :, i:i+1] - flow_field.v_sorted[:, :, i:i+1] + flow_field.u_sorted[:, i:i+1] + flow_field.v_sorted[:, i:i+1] ct_i = Ct( velocities=flow_field.u_sorted, @@ -1226,7 +1231,7 @@ def empirical_gauss_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, @@ -1242,13 +1247,14 @@ def empirical_gauss_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] - 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] + axial_induction_i = axial_induction_i[:, 0:1, None, None] + 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] - effective_yaw_i = np.zeros_like(yaw_angle_i) - effective_yaw_i += yaw_angle_i + # Secondary steering not currently implemented in EmGauss model + # effective_yaw_i = np.zeros_like(yaw_angle_i) + # effective_yaw_i += yaw_angle_i average_velocities = average_velocity( flow_field.u_sorted, @@ -1256,7 +1262,7 @@ def empirical_gauss_solver( cubature_weights=grid.cubature_weights ) tilt_angle_i = farm.calculate_tilt_for_eff_velocities(average_velocities) - tilt_angle_i = tilt_angle_i[:, :, i:i+1, None, None] + tilt_angle_i = tilt_angle_i[:, i:i+1, None, None] if model_manager.enable_secondary_steering: raise NotImplementedError( @@ -1268,7 +1274,7 @@ def empirical_gauss_solver( if model_manager.enable_yaw_added_recovery: # Influence of yawing on turbine's own wake - mixing_factor[:, :, i:i+1, i] += \ + mixing_factor[:, i:i+1, i] += \ yaw_added_wake_mixing( axial_induction_i, yaw_angle_i, @@ -1278,7 +1284,7 @@ def empirical_gauss_solver( # Extract total wake induced mixing for turbine i mixing_i = np.linalg.norm( - mixing_factor[:, :, i:i+1, :, None], + mixing_factor[:, i:i+1, :, None], ord=2, axis=3, keepdims=True ) @@ -1287,7 +1293,7 @@ def empirical_gauss_solver( deflection_field_y, deflection_field_z = model_manager.deflection_model.function( x_i, y_i, - effective_yaw_i, + yaw_angle_i, tilt_angle_i, mixing_i, ct_i, @@ -1318,6 +1324,7 @@ def empirical_gauss_solver( ) # Calculate wake overlap for wake-added turbulence (WAT) + import ipdb; ipdb.set_trace() # AT THIS POINT area_overlap = np.sum(velocity_deficit * flow_field.u_initial_sorted > 0.05, axis=(3, 4))\ / (grid.grid_resolution * grid.grid_resolution) diff --git a/floris/simulation/wake_velocity/empirical_gauss.py b/floris/simulation/wake_velocity/empirical_gauss.py index fe2bed0e0..db2308d22 100644 --- a/floris/simulation/wake_velocity/empirical_gauss.py +++ b/floris/simulation/wake_velocity/empirical_gauss.py @@ -170,7 +170,7 @@ def function( self.mixing_gain_velocity * mixing_i, ) sigma_y[upstream_mask] = \ - np.tile(sigma_y0, np.shape(sigma_y)[2:])[upstream_mask] + np.tile(sigma_y0, np.shape(sigma_y)[1:])[upstream_mask] sigma_z = empirical_gauss_model_wake_width( x - x_i, @@ -181,7 +181,7 @@ def function( self.mixing_gain_velocity * mixing_i, ) sigma_z[upstream_mask] = \ - np.tile(sigma_z0, np.shape(sigma_z)[2:])[upstream_mask] + np.tile(sigma_z0, np.shape(sigma_z)[1:])[upstream_mask] # 'Standard' wake component r, C = rCalt(