Profiling #293

rafmudaf · 2022-01-25T22:27:38Z

rafmudaf
Jan 25, 2022
Maintainer

Speed / Time

I've used cProfile in the past, but I found line_profiler to be much more user friendly. It handles pretty much everything for you and only requires adding a @profile decorator to the function to target.

The package is called line_profiler but the executable is kernprof, and it is run with this command:

kernprof -v -l script_to_profile.py

Another potentially relevant tool for profiling is scalene - a CPU, GPU and memory profiler.

Memory

Memory profiling has proved to be much more difficult. I've used memory_profiler due to its similar API to line_profiler - you simply add the same @profile decorator to a function. However, the same code differs in memory measurement from run to run, so getting accurate results seems elusive.

This package is run with:

python -m memory_profiler script.py

rafmudaf · 2022-07-26T18:40:03Z

rafmudaf
Jul 26, 2022
Maintainer Author

Here are results of profiling using a modified version the included profiling script. The script used here is below.

    sample_inputs = SampleInputs()

    sample_inputs.floris["wake"]["model_strings"]["velocity_model"] = "gauss"
    sample_inputs.floris["wake"]["model_strings"]["deflection_model"] = "gauss"
    sample_inputs.floris["wake"]["enable_secondary_steering"] = True
    sample_inputs.floris["wake"]["enable_yaw_added_recovery"] = True
    sample_inputs.floris["wake"]["enable_transverse_velocities"] = True

    factor = 100
    TURBINE_DIAMETER = sample_inputs.floris["farm"]["turbine_type"][0]["rotor_diameter"]
    sample_inputs.floris["farm"]["layout_x"] = [5 * TURBINE_DIAMETER * i for i in range(factor)]
    sample_inputs.floris["farm"]["layout_y"] = [0.0 for i in range(factor)]

    factor = 10
    sample_inputs.floris["flow_field"]["wind_directions"] = factor * [270.0]
    sample_inputs.floris["flow_field"]["wind_speeds"] = factor * [8.0]

    floris = Floris.from_dict(sample_inputs.floris)
    floris.initialize_domain()
    floris.steady_state_atmospheric_condition()

First, I compared the three high-level functions required to initialize and run the simulation. As expected, the calculation step is the most expensive.

    27         1      67484.0  67484.0      0.6      floris = Floris.from_dict(sample_inputs.floris)
    28         1       8861.0   8861.0      0.1      floris.initialize_domain()
    29         1   10880043.0 10880043.0     99.3      floris.steady_state_atmospheric_condition()

Next, I profiled the sequential_solver function and removed the lines with less than 10% of the overhead. This is denoted by ... in the line number column.

File: /Users/rmudafor/Development/floris/floris/simulation/solver.py
Function: sequential_solver at line 46

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
    46                                           @profile
    47                                           def sequential_solver(farm: Farm, flow_field: FlowField, grid: TurbineGrid, model_manager: WakeModelManager) -> None:
    48                                               # Algorithm
    49                                               # For each turbine, calculate its effect on every downstream turbine.
    50                                               # For the current turbine, we are calculating the deficit that it adds to downstream turbines.
    51                                               # Integrate this into the main data structure.
    52                                               # Move on to the next turbine.
    53                                           
    66                                               # Calculate the velocity deficit sequentially from upstream to downstream turbines
    67       101        146.0      1.4      0.0      for i in range(grid.n_turbines):
    ...
   120                                                   # Model calculations
   121                                                   # NOTE: exponential
   122       300    2383378.0   7944.6     20.8          deflection_field = model_manager.deflection_model.function(
   123       100         96.0      1.0      0.0              x_i,
   124       100         95.0      0.9      0.0              y_i,
   125       100         93.0      0.9      0.0              effective_yaw_i,
   126       100         94.0      0.9      0.0              turbulence_intensity_i,
   127       100         90.0      0.9      0.0              ct_i,
   128       100         91.0      0.9      0.0              rotor_diameter_i,
   129       100         90.0      0.9      0.0              **deflection_model_args
   130                                                   )
   131                                           
   132       100        382.0      3.8      0.0          if model_manager.enable_transverse_velocities:
   133       200    4419457.0  22097.3     38.6              v_wake, w_wake = calculate_transverse_velocity(
   134       100        102.0      1.0      0.0                  u_i,
   135       100        148.0      1.5      0.0                  flow_field.u_initial_sorted,
   136       100      25526.0    255.3      0.2                  grid.x_sorted - x_i,
   137       100      33274.0    332.7      0.3                  grid.y_sorted - y_i,
   138       100        159.0      1.6      0.0                  grid.z_sorted,
   139       100         89.0      0.9      0.0                  rotor_diameter_i,
   140       100         92.0      0.9      0.0                  hub_height_i,
   141       100         85.0      0.8      0.0                  yaw_angle_i,
   142       100         85.0      0.8      0.0                  ct_i,
   143       100         88.0      0.9      0.0                  TSR_i,
   144       100         92.0      0.9      0.0                  axial_induction_i
   145                                                       )
   ...
   160       300    2068141.0   6893.8     18.1          velocity_deficit = model_manager.velocity_model.function(
   161       100         91.0      0.9      0.0              x_i,
   162       100        117.0      1.2      0.0              y_i,
   163       100         89.0      0.9      0.0              z_i,
   164       100         91.0      0.9      0.0              axial_induction_i,
   165       100         92.0      0.9      0.0              deflection_field,
   166       100         89.0      0.9      0.0              yaw_angle_i,
   167       100         88.0      0.9      0.0              turbulence_intensity_i,
   168       100         89.0      0.9      0.0              ct_i,
   169       100         88.0      0.9      0.0              hub_height_i,
   170       100         83.0      0.8      0.0              rotor_diameter_i,
   171       100         89.0      0.9      0.0              **deficit_model_args
   172                                                   )
   ...
   179       200    1138874.0   5694.4      9.9          wake_added_turbulence_intensity = model_manager.turbulence_model.function(
   180       100        112.0      1.1      0.0              ambient_turbulence_intensity,
   181       100        155.0      1.6      0.0              grid.x_sorted,
   182       100         92.0      0.9      0.0              x_i,
   183       100         95.0      0.9      0.0              rotor_diameter_i,
   184       100         94.0      0.9      0.0              axial_induction_i
   185                                                   )

The largest expense here is calculate_transverse_velocities in the Gauss models. Below is a profile of that function. It is much more evenly distributed as there are a few large vector operations with exponents.

File: /Users/rmudafor/Development/floris/floris/simulation/wake_deflection/gauss.py
Function: calculate_transverse_velocity at line 341

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
   341                                           @profile
   342                                           def calculate_transverse_velocity(
   343                                               u_i,
   344                                               u_initial,
   345                                               delta_x,
   346                                               delta_y,
   347                                               z,
   348                                               rotor_diameter,
   349                                               hub_height,
   350                                               yaw,
   351                                               ct_i,
   352                                               tsr_i,
   353                                               axial_induction_i,
   354                                               scale=1.0
   355                                           ):
   356                                               """
   357                                               Calculate transverse velocity components for all downstream turbines
   358                                               given the vortices at the current turbine.
   359                                               """
   360                                           
   361                                               # turbine parameters
   362       100        290.0      2.9      0.0      D = rotor_diameter
   363       100        140.0      1.4      0.0      HH = hub_height
   364       100        106.0      1.1      0.0      Ct = ct_i
   365       100        103.0      1.0      0.0      TSR = tsr_i
   366       100        102.0      1.0      0.0      aI = axial_induction_i
   367                                           
   368                                               # flow parameters
   369       100      24636.0    246.4      0.6      Uinf = np.mean(u_initial, axis=(2,3,4))
   370       100        454.0      4.5      0.0      Uinf = Uinf[:,:,None,None,None]
   371                                           
   372       100        107.0      1.1      0.0      eps_gain = 0.2
   373       100        977.0      9.8      0.0      eps = eps_gain * D  # Use set value
   374                                           
   375                                               # TODO: wind sheer is hard-coded here but should be connected to the input
   376       100       3444.0     34.4      0.1      vel_top = ((HH + D / 2) / HH) ** 0.12 * np.ones((1, 1, 1, 1, 1))
   377       200       2955.0     14.8      0.1      Gamma_top = sind(yaw) * cosd(yaw) * gamma(
   378       100        115.0      1.1      0.0          D,
   379       100        106.0      1.1      0.0          vel_top,
   380       100        100.0      1.0      0.0          Uinf,
   381       100         97.0      1.0      0.0          Ct,
   382       100        110.0      1.1      0.0          scale,
   383                                               )
   384                                           
   385       100       2248.0     22.5      0.1      vel_bottom = ((HH - D / 2) / HH) ** 0.12 * np.ones((1, 1, 1, 1, 1))
   386       200       2478.0     12.4      0.1      Gamma_bottom = -1 * sind(yaw) * cosd(yaw) * gamma(
   387       100        123.0      1.2      0.0          D,
   388       100         98.0      1.0      0.0          vel_bottom,
   389       100         97.0      1.0      0.0          Uinf,
   390       100         99.0      1.0      0.0          Ct,
   391       100         99.0      1.0      0.0          scale,
   392                                               )
   393                                           
   394       100      11369.0    113.7      0.3      turbine_average_velocity = np.cbrt(np.mean(u_i ** 3, axis=(3,4)))
   395       100        315.0      3.1      0.0      turbine_average_velocity = turbine_average_velocity[:,:,:,None,None]
   396       100       1722.0     17.2      0.0      Gamma_wake_rotation = 0.25 * 2 * np.pi * D * (aI - aI ** 2) * turbine_average_velocity / TSR
   397                                           
   398                                               ### compute the spanwise and vertical velocities induced by yaw
   399                                           
   400                                               # decay the vortices as they move downstream - using mixing length
   401       100        567.0      5.7      0.0      lmda = D / 8
   402       100        111.0      1.1      0.0      kappa = 0.41
   403       100     122688.0   1226.9      3.1      lm = kappa * z / (1 + kappa * z / lmda)
   404                                               # TODO: get this from the z input?
   405       100     131905.0   1319.0      3.3      z_basis = np.linspace(np.min(z), np.max(z), np.shape(u_initial)[4])
   406       100     198955.0   1989.5      5.0      dudz_initial = np.gradient(u_initial, z_basis, axis=4)
   407       100      70040.0    700.4      1.8      nu = lm ** 2 * np.abs(dudz_initial)
   408                                           
   409       100     108350.0   1083.5      2.7      decay = eps ** 2 / (4 * nu * delta_x / Uinf + eps ** 2)   # This is the decay downstream
   410       100      19306.0    193.1      0.5      yLocs = delta_y + BaseModel.NUM_EPS
   411                                           
   412                                               # top vortex
   413       100      43138.0    431.4      1.1      zT = z - (HH + D / 2) + BaseModel.NUM_EPS
   414       100      72547.0    725.5      1.8      rT = yLocs ** 2 + zT ** 2  # TODO: This is - in the paper
   415       100     206206.0   2062.1      5.2      core_shape = 1 - np.exp(-rT / (eps ** 2))  # This looks like spanwise decay - it defines the vortex profile in the spanwise directions
   416       100      92553.0    925.5      2.3      V1 = (Gamma_top * zT) / (2 * np.pi * rT) * core_shape * decay
   417       100      96130.0    961.3      2.4      W1 = (-1 * Gamma_top * yLocs) / (2 * np.pi * rT) * core_shape * decay
   418                                           
   419                                               # bottom vortex
   420       100      32136.0    321.4      0.8      zB = z - (HH - D / 2) + BaseModel.NUM_EPS
   421       100      73730.0    737.3      1.9      rB = yLocs ** 2 + zB ** 2
   422       100     210411.0   2104.1      5.3      core_shape = 1 - np.exp(-rB / (eps ** 2))
   423       100      83104.0    831.0      2.1      V2 = (Gamma_bottom * zB) / (2 * np.pi * rB) * core_shape * decay
   424       100     100663.0   1006.6      2.5      W2 = (-1 * Gamma_bottom * yLocs) / (2 * np.pi * rB) * core_shape * decay
   425                                           
   426                                               # wake rotation vortex
   427       100      31215.0    312.1      0.8      zC = z - HH + BaseModel.NUM_EPS
   428       100      87740.0    877.4      2.2      rC = yLocs ** 2 + zC ** 2
   429       100     208554.0   2085.5      5.2      core_shape = 1 - np.exp(-rC / (eps ** 2))
   430       100      79924.0    799.2      2.0      V5 = (Gamma_wake_rotation * zC) / (2 * np.pi * rC) * core_shape * decay
   431       100      98107.0    981.1      2.5      W5 = (-1 * Gamma_wake_rotation * yLocs) / (2 * np.pi * rC) * core_shape * decay
   432                                           
   433                                           
   434                                               ### Boundary condition - ground mirror vortex
   435                                           
   436                                               # top vortex - ground
   437       100      31539.0    315.4      0.8      zTb = z + (HH + D / 2) + BaseModel.NUM_EPS
   438       100      84368.0    843.7      2.1      rTb = yLocs ** 2 + zTb ** 2
   439       100     210313.0   2103.1      5.3      core_shape = 1 - np.exp(-rTb / (eps ** 2))  # This looks like spanwise decay - it defines the vortex profile in the spanwise directions
   440       100      79925.0    799.2      2.0      V3 = (-1 * Gamma_top * zTb) / (2 * np.pi * rTb) * core_shape * decay
   441       100      96886.0    968.9      2.4      W3 = (Gamma_top * yLocs) / (2 * np.pi * rTb) * core_shape * decay
   442                                           
   443                                               # bottom vortex - ground
   444       100      31520.0    315.2      0.8      zBb = z + (HH - D / 2) + BaseModel.NUM_EPS
   445       100      85021.0    850.2      2.1      rBb = yLocs ** 2 + zBb ** 2
   446       100     202589.0   2025.9      5.1      core_shape = 1 - np.exp(-rBb / (eps ** 2))
   447       100      76861.0    768.6      1.9      V4 = (-1 * Gamma_bottom * zBb) / (2 * np.pi * rBb) * core_shape * decay
   448       100      96910.0    969.1      2.4      W4 = (Gamma_bottom * yLocs) / (2 * np.pi * rBb) * core_shape * decay
   449                                           
   450                                               # wake rotation vortex - ground effect
   451       100      30213.0    302.1      0.8      zCb = z + HH + BaseModel.NUM_EPS
   452       100      86789.0    867.9      2.2      rCb = yLocs ** 2 + zCb ** 2
   453       100     205271.0   2052.7      5.2      core_shape = 1 - np.exp(-rCb / (eps ** 2))
   454       100      77055.0    770.5      1.9      V6 = (-1 * Gamma_wake_rotation * zCb) / (2 * np.pi * rCb) * core_shape * decay
   455       100      99461.0    994.6      2.5      W6 = (Gamma_wake_rotation * yLocs) / (2 * np.pi * rCb) * core_shape * decay
   456                                           
   457                                               # total spanwise velocity
   458       100      85799.0    858.0      2.2      V = V1 + V2 + V3 + V4 + V5 + V6
   459       100     105680.0   1056.8      2.7      W = W1 + W2 + W3 + W4 + W5 + W6
   460                                           
   461                                               # no spanwise and vertical velocity upstream of the turbine
   462                                               # V[delta_x < -1] = 0.0  # Subtract by 1 to avoid numerical issues on rotation
   463                                               # W[delta_x < -1] = 0.0  # Subtract by 1 to avoid numerical issues on rotation
   464                                               # TODO Should this be <= ? Shouldn't be adding V and W on the current turbine?
   465       100      30596.0    306.0      0.8      V[delta_x < 0.0] = 0.0  # Subtract by 1 to avoid numerical issues on rotation
   466       100      23739.0    237.4      0.6      W[delta_x < 0.0] = 0.0  # Subtract by 1 to avoid numerical issues on rotation
   467                                           
   468                                               # TODO: Why would the say W cannot be negative?
   469       100      19904.0    199.0      0.5      W[W < 0] = 0
   470                                           
   471       100        166.0      1.7      0.0      return V, W

0 replies

bayc · 2022-07-26T19:07:34Z

bayc
Jul 26, 2022
Maintainer

This is really interesting to see @rafmudaf ! Do you think you could generate similar timing plots with GCH turned off as well? Just wondering if there will be any significant changes in the other calls of the code. There shouldn't be, but curious to see if there is some other impact from memory or something else.

4 replies

paulf81 Jul 26, 2022
Maintainer

really agree with this

rafmudaf Jul 26, 2022
Maintainer Author

Good idea, that is in the comment below. You can't compare the percentages directly since they're relative to their own calculations, but you can compare the "Time" column.

rafmudaf Jul 26, 2022
Maintainer Author

Here are the relevant lines.

Yes GCH

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
   122       300    2383378.0   7944.6     20.8          deflection_field = model_manager.deflection_model.function(
   133       200    4419457.0  22097.3     38.6              v_wake, w_wake = calculate_transverse_velocity(
   160       300    2068141.0   6893.8     18.1          velocity_deficit = model_manager.velocity_model.function(
   179       200    1138874.0   5694.4      9.9          wake_added_turbulence_intensity = model_manager.turbulence_model.function(

No GCH

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
   122       300    2320347.0   7734.5     36.7          deflection_field = model_manager.deflection_model.function(
   160       300    1777625.0   5925.4     28.1          velocity_deficit = model_manager.velocity_model.function(
   179       200    1039832.0   5199.2     16.5          wake_added_turbulence_intensity = model_manager.turbulence_model.function(

rafmudaf Jul 26, 2022
Maintainer Author

Line	No GCH	Yes GCH	% diff
122	2320347	2383378	2.71644715
160	1777625	2068141	16.3429295
179	1039832	1138874	9.52480785

rafmudaf · 2022-07-26T19:18:18Z

rafmudaf
Jul 26, 2022
Maintainer Author

Here's the timing of the sequential solve function with GCH off. The script is below and then the timing.

    sample_inputs = SampleInputs()

    sample_inputs.floris["wake"]["model_strings"]["velocity_model"] = "gauss"
    sample_inputs.floris["wake"]["model_strings"]["deflection_model"] = "gauss"
    sample_inputs.floris["wake"]["enable_secondary_steering"] = False
    sample_inputs.floris["wake"]["enable_yaw_added_recovery"] = False
    sample_inputs.floris["wake"]["enable_transverse_velocities"] = False

    factor = 100
    TURBINE_DIAMETER = sample_inputs.floris["farm"]["turbine_type"][0]["rotor_diameter"]
    sample_inputs.floris["farm"]["layout_x"] = [5 * TURBINE_DIAMETER * i for i in range(factor)]
    sample_inputs.floris["farm"]["layout_y"] = [0.0 for i in range(factor)]

    factor = 10
    sample_inputs.floris["flow_field"]["wind_directions"] = factor * [270.0]
    sample_inputs.floris["flow_field"]["wind_speeds"] = factor * [8.0]

    floris = Floris.from_dict(sample_inputs.floris)
    floris.initialize_domain()
    floris.steady_state_atmospheric_condition()

File: /Users/rmudafor/Development/floris/floris/simulation/solver.py
Function: sequential_solver at line 46

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
    46                                           @profile
    47                                           def sequential_solver(farm: Farm, flow_field: FlowField, grid: TurbineGrid, model_manager: WakeModelManager) -> None:
    48                                               # Algorithm
    49                                               # For each turbine, calculate its effect on every downstream turbine.
    50                                               # For the current turbine, we are calculating the deficit that it adds to downstream turbines.
    51                                               # Integrate this into the main data structure.
    52                                               # Move on to the next turbine.
   ...
    66                                               # Calculate the velocity deficit sequentially from upstream to downstream turbines
    67       101        141.0      1.4      0.0      for i in range(grid.n_turbines):
    68                                           
   ...
   120                                                   # Model calculations
   121                                                   # NOTE: exponential
   122       300    2320347.0   7734.5     36.7          deflection_field = model_manager.deflection_model.function(
   123       100         91.0      0.9      0.0              x_i,
   124       100         89.0      0.9      0.0              y_i,
   125       100         93.0      0.9      0.0              effective_yaw_i,
   126       100         85.0      0.8      0.0              turbulence_intensity_i,
   127       100         83.0      0.8      0.0              ct_i,
   128       100         81.0      0.8      0.0              rotor_diameter_i,
   129       100         83.0      0.8      0.0              **deflection_model_args
   130                                                   )
   131                                           
   132       100        353.0      3.5      0.0          if model_manager.enable_transverse_velocities:
   133                                                       v_wake, w_wake = calculate_transverse_velocity(
   134                                                           u_i,
   135                                                           flow_field.u_initial_sorted,
   136                                                           grid.x_sorted - x_i,
   137                                                           grid.y_sorted - y_i,
   138                                                           grid.z_sorted,
   139                                                           rotor_diameter_i,
   140                                                           hub_height_i,
   141                                                           yaw_angle_i,
   142                                                           ct_i,
   143                                                           TSR_i,
   144                                                           axial_induction_i
   145                                                       )
   ...
   160       300    1777625.0   5925.4     28.1          velocity_deficit = model_manager.velocity_model.function(
   161       100         84.0      0.8      0.0              x_i,
   162       100         80.0      0.8      0.0              y_i,
   163       100        101.0      1.0      0.0              z_i,
   164       100         83.0      0.8      0.0              axial_induction_i,
   165       100         77.0      0.8      0.0              deflection_field,
   166       100         82.0      0.8      0.0              yaw_angle_i,
   167       100         85.0      0.8      0.0              turbulence_intensity_i,
   168       100         86.0      0.9      0.0              ct_i,
   169       100         85.0      0.8      0.0              hub_height_i,
   170       100         82.0      0.8      0.0              rotor_diameter_i,
   171       100         95.0      0.9      0.0              **deficit_model_args
   172                                                   )
   ...
   179       200    1039832.0   5199.2     16.5          wake_added_turbulence_intensity = model_manager.turbulence_model.function(
   180       100        101.0      1.0      0.0              ambient_turbulence_intensity,
   181       100        169.0      1.7      0.0              grid.x_sorted,
   182       100         94.0      0.9      0.0              x_i,
   183       100         91.0      0.9      0.0              rotor_diameter_i,

0 replies

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Profiling #293

{{title}}

{{editor}}'s edit

{{editor}}'s edit

Replies: 3 comments 4 replies

{{title}}

{{editor}}'s edit

{{editor}}'s edit

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

Select a reply

Profiling #293

rafmudaf Jan 25, 2022 Maintainer

Speed / Time

Memory

Replies: 3 comments · 4 replies

rafmudaf Jul 26, 2022 Maintainer Author

bayc Jul 26, 2022 Maintainer

paulf81 Jul 26, 2022 Maintainer

rafmudaf Jul 26, 2022 Maintainer Author

rafmudaf Jul 26, 2022 Maintainer Author

Yes GCH

No GCH

rafmudaf Jul 26, 2022 Maintainer Author

rafmudaf Jul 26, 2022 Maintainer Author

rafmudaf
Jan 25, 2022
Maintainer

Replies: 3 comments 4 replies

rafmudaf
Jul 26, 2022
Maintainer Author

bayc
Jul 26, 2022
Maintainer

paulf81 Jul 26, 2022
Maintainer

rafmudaf Jul 26, 2022
Maintainer Author

rafmudaf Jul 26, 2022
Maintainer Author

rafmudaf Jul 26, 2022
Maintainer Author

rafmudaf
Jul 26, 2022
Maintainer Author