Performance Improvements and Code Simplification (#494)

* add working example to test speedups as size grows in three dimensions

* update pre-commit hooks

* update example

* udpate to run for just a single set of parameters

* update for base timing results

* add small positive changes

* update runtime, boolean operations, and gamma

* compartmentalize the vortex calculation

* apply vortex calculation to transverse velocity

* remove comments

* small improvements

* reformat

* small improvements to crespo hernandez

* add more numexpression where it counts most

* add more numexpression

* remove np.array around boolean arrays

* add control run time for develop

* remove @Profile

* refactor vortex calculation

* refactor characteristic wake width

* refactor precalculate

* slick minor speedup

* jensen speedup by calculating in one fell swoop

* fix type annotation

* repalce filtering with where statements

* remove commented out old code

* switch numpy add reduce for multiple additions

* unblacken solver for most changes

* unblacken cumulative gauss for most changes

* unblacken gauss for most changes

* make exponents consistent

* Maintain profiling decorators

These are super useful for both memory and performance (runtime) profiling

* Correct a merge conflict

* update ruff settings

* Fix formatting

* Fix formatting

* Fix formatting

* Expand in-line comments for Jensen model

Also removes the implementation of the model that doesn't use numexpr. This was helpful in initial development of v3, but it's no longer practical.

* Revert code structure but maintain speed ups

This change should be revisited. There’s good changes, especially putting the vortex calculations in functions, but it should more directly relate to the model papers

* Remove redundant np.array type-casts

* Remove speed test example

---------

Co-authored-by: Rafael M Mudafort <rafael.mudafort@nrel.gov>
Co-authored-by: Rafael M Mudafort <rafmudaf@gmail.com>

floris/simulation/solver.py

@@ -130,9 +130,9 @@ def sequential_solver(
         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]
+        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
@@ -148,7 +148,7 @@ def sequential_solver(
                 hub_height_i,
                 ct_i,
                 TSR_i,
-                axial_induction_i
+                axial_induction_i,
             )
             effective_yaw_i += added_yaw
 
@@ -161,7 +161,7 @@ def sequential_solver(
             turbulence_intensity_i,
             ct_i,
             rotor_diameter_i,
-            **deflection_model_args
+            **deflection_model_args,
         )
 
         if model_manager.enable_transverse_velocities:
@@ -177,7 +177,7 @@ def sequential_solver(
                 yaw_angle_i,
                 ct_i,
                 TSR_i,
-                axial_induction_i
+                axial_induction_i,
             )
 
         if model_manager.enable_yaw_added_recovery:
@@ -204,7 +204,7 @@ def sequential_solver(
             ct_i,
             hub_height_i,
             rotor_diameter_i,
-            **deficit_model_args
+            **deficit_model_args,
         )
 
         wake_field = model_manager.combination_model.function(
@@ -217,7 +217,7 @@ def sequential_solver(
             grid.x_sorted,
             x_i,
             rotor_diameter_i,
-            axial_induction_i
+            axial_induction_i,
         )
 
         # Calculate wake overlap for wake-added turbulence (WAT)
@@ -232,9 +232,9 @@ def sequential_solver(
         ti_added = (
             area_overlap
             * np.nan_to_num(wake_added_turbulence_intensity, posinf=0.0)
-            * np.array(grid.x_sorted > x_i)
-            * np.array(np.abs(y_i - grid.y_sorted) < 2 * rotor_diameter_i)
-            * np.array(grid.x_sorted <= downstream_influence_length + x_i)
+            * (grid.x_sorted > x_i)
+            * (np.abs(y_i - grid.y_sorted) < 2 * rotor_diameter_i)
+            * (grid.x_sorted <= downstream_influence_length + x_i)
         )
 
         # Combine turbine TIs with WAT
@@ -291,7 +291,7 @@ def full_flow_sequential_solver(
     turbine_grid_farm.expand_farm_properties(
         turbine_grid_flow_field.n_wind_directions,
         turbine_grid_flow_field.n_wind_speeds,
-        turbine_grid.sorted_coord_indices
+        turbine_grid.sorted_coord_indices,
     )
     turbine_grid_flow_field.initialize_velocity_field(turbine_grid)
     turbine_grid_farm.initialize(turbine_grid.sorted_indices)
@@ -376,7 +376,7 @@ def full_flow_sequential_solver(
                 hub_height_i,
                 ct_i,
                 TSR_i,
-                axial_induction_i
+                axial_induction_i,
             )
             effective_yaw_i += added_yaw
 
@@ -389,7 +389,7 @@ def full_flow_sequential_solver(
             turbulence_intensity_i,
             ct_i,
             rotor_diameter_i,
-            **deflection_model_args
+            **deflection_model_args,
         )
 
         if model_manager.enable_transverse_velocities:
@@ -405,7 +405,7 @@ def full_flow_sequential_solver(
                 yaw_angle_i,
                 ct_i,
                 TSR_i,
-                axial_induction_i
+                axial_induction_i,
             )
 
         # NOTE: exponential
@@ -420,7 +420,7 @@ def full_flow_sequential_solver(
             ct_i,
             hub_height_i,
             rotor_diameter_i,
-            **deficit_model_args
+            **deficit_model_args,
         )
 
         wake_field = model_manager.combination_model.function(
@@ -476,10 +476,10 @@ def cc_solver(
         rotor_diameter_i = farm.rotor_diameters_sorted[: ,:, i:i+1, None, None]
 
         mask2 = (
-            np.array(grid.x_sorted < x_i + 0.01)
-            * np.array(grid.x_sorted > x_i - 0.01)
-            * np.array(grid.y_sorted < y_i + 0.51 * rotor_diameter_i)
-            * np.array(grid.y_sorted > y_i - 0.51 * rotor_diameter_i)
+            (grid.x_sorted < x_i + 0.01)
+            * (grid.x_sorted > x_i - 0.01)
+            * (grid.y_sorted < y_i + 0.51 * rotor_diameter_i)
+            * (grid.y_sorted > y_i - 0.51 * rotor_diameter_i)
         )
         turb_inflow_field = (
             turb_inflow_field * ~mask2
@@ -536,8 +536,8 @@ def cc_solver(
 
         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]
-        TSR_i = farm.TSRs_sorted[: ,:, i:i+1, None, None]
+        hub_height_i = farm.hub_heights_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
@@ -567,7 +567,7 @@ def cc_solver(
             turbulence_intensity_i,
             turb_Cts[:, :, i:i+1],
             rotor_diameter_i,
-            **deflection_model_args
+            **deflection_model_args,
         )
 
         if model_manager.enable_transverse_velocities:
@@ -584,7 +584,7 @@ def cc_solver(
                 turb_Cts[:, :, i:i+1],
                 TSR_i,
                 axial_induction_i,
-                scale=2.0
+                scale=2.0,
             )
 
         if model_manager.enable_yaw_added_recovery:
@@ -612,7 +612,7 @@ def cc_solver(
             farm.rotor_diameters_sorted[:, :, :, None, None],
             turb_u_wake,
             Ctmp,
-            **deficit_model_args
+            **deficit_model_args,
         )
 
         wake_added_turbulence_intensity = model_manager.turbulence_model.function(
@@ -635,16 +635,17 @@ def cc_solver(
         ti_added = (
             area_overlap
             * np.nan_to_num(wake_added_turbulence_intensity, posinf=0.0)
-            * np.array(grid.x_sorted > x_i)
-            * np.array(np.abs(y_i - grid.y_sorted) < 2 * rotor_diameter_i)
-            * np.array(grid.x_sorted <= downstream_influence_length + x_i)
+            * (grid.x_sorted > x_i)
+            * (np.abs(y_i - grid.y_sorted) < 2 * rotor_diameter_i)
+            * (grid.x_sorted <= downstream_influence_length + x_i)
         )
 
         # Combine turbine TIs with WAT
         turbine_turbulence_intensity = np.maximum(
-            np.sqrt( ti_added ** 2 + ambient_turbulence_intensity ** 2 ),
+            np.sqrt(ti_added ** 2 + ambient_turbulence_intensity ** 2),
             turbine_turbulence_intensity
         )
+
         flow_field.v_sorted += v_wake
         flow_field.w_sorted += w_wake
     flow_field.u_sorted = turb_inflow_field
@@ -660,7 +661,7 @@ def full_flow_cc_solver(
     farm: Farm,
     flow_field: FlowField,
     flow_field_grid: FlowFieldGrid,
-    model_manager: WakeModelManager
+    model_manager: WakeModelManager,
 ) -> None:
     # Get the flow quantities and turbine performance
     turbine_grid_farm = copy.deepcopy(farm)
@@ -692,7 +693,7 @@ def full_flow_cc_solver(
     turbine_grid_farm.expand_farm_properties(
         turbine_grid_flow_field.n_wind_directions,
         turbine_grid_flow_field.n_wind_speeds,
-        turbine_grid.sorted_coord_indices
+        turbine_grid.sorted_coord_indices,
     )
     turbine_grid_flow_field.initialize_velocity_field(turbine_grid)
     turbine_grid_farm.initialize(turbine_grid.sorted_indices)
@@ -764,9 +765,9 @@ def full_flow_cc_solver(
         turbulence_intensity_i = \
             turbine_grid_flow_field.turbulence_intensity_field_sorted_avg[:, :, i:i+1]
         yaw_angle_i = turbine_grid_farm.yaw_angles_sorted[:, :, i:i+1, None, None]
-        hub_height_i = turbine_grid_farm.hub_heights_sorted[: ,:, i:i+1, None, None]
-        rotor_diameter_i = turbine_grid_farm.rotor_diameters_sorted[: ,:, i:i+1, None, None]
-        TSR_i = turbine_grid_farm.TSRs_sorted[: ,:, i:i+1, None, None]
+        hub_height_i = turbine_grid_farm.hub_heights_sorted[:, :, i:i+1, None, None]
+        rotor_diameter_i = turbine_grid_farm.rotor_diameters_sorted[:, :, i:i+1, None, None]
+        TSR_i = turbine_grid_farm.TSRs_sorted[:, :, i:i+1, None, None]
 
         effective_yaw_i = np.zeros_like(yaw_angle_i)
         effective_yaw_i += yaw_angle_i
@@ -783,7 +784,7 @@ def full_flow_cc_solver(
                 turb_Cts[:, :, i:i+1],
                 TSR_i,
                 axial_induction_i,
-                scale=2.0
+                scale=2.0,
             )
             effective_yaw_i += added_yaw
 
@@ -796,7 +797,7 @@ def full_flow_cc_solver(
             turbulence_intensity_i,
             turb_Cts[:, :, i:i+1],
             rotor_diameter_i,
-            **deflection_model_args
+            **deflection_model_args,
         )
 
         if model_manager.enable_transverse_velocities:
@@ -813,7 +814,7 @@ def full_flow_cc_solver(
                 turb_Cts[:, :, i:i+1],
                 TSR_i,
                 axial_induction_i,
-                scale=2.0
+                scale=2.0,
             )
 
         # NOTE: exponential
@@ -830,7 +831,7 @@ def full_flow_cc_solver(
             turbine_grid_farm.rotor_diameters_sorted[:, :, :, None, None],
             turb_u_wake,
             Ctmp,
-            **deficit_model_args
+            **deficit_model_args,
         )
 
         flow_field.v_sorted += v_wake
@@ -928,9 +929,9 @@ def turbopark_solver(
         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]
+        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
@@ -946,7 +947,7 @@ def turbopark_solver(
                 hub_height_i,
                 ct_i,
                 TSR_i,
-                axial_induction_i
+                axial_induction_i,
             )
             effective_yaw_i += added_yaw
 
@@ -980,7 +981,7 @@ def turbopark_solver(
                     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]
+                rotor_diameter_ii = farm.rotor_diameters_sorted[:, :, ii:ii+1, None, None]
 
                 deflection_field_ii = model_manager.deflection_model.function(
                     x_ii,
@@ -989,10 +990,10 @@ def turbopark_solver(
                     turbulence_intensity_ii,
                     ct_ii,
                     rotor_diameter_ii,
-                    **deflection_model_args
+                    **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(
@@ -1007,7 +1008,7 @@ def turbopark_solver(
                 yaw_angle_i,
                 ct_i,
                 TSR_i,
-                axial_induction_i
+                axial_induction_i,
             )
 
         if model_manager.enable_yaw_added_recovery:
@@ -1033,7 +1034,7 @@ def turbopark_solver(
             farm.rotor_diameters_sorted[:, :, :, None, None],
             i,
             deflection_field,
-            **deficit_model_args
+            **deficit_model_args,
         )
 
         wake_field = model_manager.combination_model.function(
@@ -1064,9 +1065,9 @@ def turbopark_solver(
         ti_added = (
             area_overlap
             * np.nan_to_num(wake_added_turbulence_intensity, posinf=0.0)
-            * np.array(grid.x_sorted > x_i)
-            * np.array(np.abs(y_i - grid.y_sorted) < 2 * rotor_diameter_i)
-            * np.array(grid.x_sorted <= downstream_influence_length + x_i)
+            * (grid.x_sorted > x_i)
+            * (np.abs(y_i - grid.y_sorted) < 2 * rotor_diameter_i)
+            * (grid.x_sorted <= downstream_influence_length + x_i)
         )
 
         # Combine turbine TIs with WAT
@@ -1111,7 +1112,6 @@ def full_flow_turbopark_solver(
     # turbine_grid_farm.construc_turbine_pPs()
     # turbine_grid_farm.construct_coordinates()
 
-
     # turbine_grid = TurbineGrid(
     #     turbine_coordinates=turbine_grid_farm.coordinates,
     #     reference_turbine_diameter=turbine_grid_farm.rotor_diameters,

floris/simulation/turbine.py

@@ -142,7 +142,7 @@ def compute_tilt_angles_for_floating_turbines(
         else:
             tilt_angles += (
                 tilt_interp[turb_type](rotor_effective_velocities)
-                * np.array(turbine_type_map == turb_type)
+                * (turbine_type_map == turb_type)
             )
 
     # TODO: Not sure if this is the best way to do this? Basically replaces the initialized
@@ -267,7 +267,7 @@ def power(
         # type to the main thrust coefficient array
         p += (
             power_interp[turb_type](rotor_effective_velocities)
-            * np.array(turbine_type_map == turb_type)
+            * (turbine_type_map == turb_type)
         )
 
     return p * ref_density_cp_ct
@@ -355,7 +355,7 @@ def Ct(
         # type to the main thrust coefficient array
         thrust_coefficient += (
             fCt[turb_type](average_velocities)
-            * np.array(turbine_type_map == turb_type)
+            * (turbine_type_map == turb_type)
         )
     thrust_coefficient = np.clip(thrust_coefficient, 0.0001, 0.9999)
     effective_thrust = thrust_coefficient * cosd(yaw_angle) * cosd(tilt_angle - ref_tilt_cp_ct)

floris/simulation/wake_deflection/empirical_gauss.py

@@ -126,7 +126,7 @@ def function(
         A_z = (deflection_gain_z * ct_i * tilt_r) / (1 + self.mixing_gain_deflection * mixing_i)
 
         # Apply downstream mask in the process
-        x_normalized = (x - x_i) * np.array(x > x_i + 0.1) / rotor_diameter_i
+        x_normalized = (x - x_i) * (x > x_i + 0.1) / rotor_diameter_i
 
         log_term = np.log(
             (x_normalized - self.deflection_rate) / (x_normalized + self.deflection_rate)