Wind Farm Simulation

In this example, a simple wind farm case is presented where the different capabilities of PyWake such as calculating AEP, extracting power values and plotting time series are shown.

Install PyWake if needed

[1]:

# Install PyWake if needed
try:
    import py_wake
except ModuleNotFoundError:
    !pip install git+https://gitlab.windenergy.dtu.dk/TOPFARM/PyWake.git

First we import some basic Python elements

[2]:

import numpy as np
import matplotlib.pyplot as plt
import xarray as xr

Now we import the site, wind turbine and wake deficit model to use in the simulation

[3]:

# import and setup site and windTurbines
%load_ext py_wake.utils.notebook_extensions

from tqdm.notebook import tqdm
from py_wake.utils import layouts
from py_wake.utils.profiling import timeit, profileit
from py_wake.examples.data.hornsrev1 import Hornsrev1Site, V80
from py_wake.literature.gaussian_models import Bastankhah_PorteAgel_2014
from py_wake.utils.plotting import setup_plot
xr.set_options(display_expand_data=False)

site = Hornsrev1Site()
x, y = site.initial_position.T
windTurbines = V80()

wf_model = Bastankhah_PorteAgel_2014(site, windTurbines, k=0.0324555)

#this allows you to see what type of engineering models you are simulating
print(wf_model)

Bastankhah_PorteAgel_2014(PropagateDownwind, BastankhahGaussianDeficit-wake, LinearSum-superposition)

Simple simulation - all wind directions and wind speeds

To run the wind farm simulation, we must call the WindFarmModel (wf_model) element that was previously created. As default, the model will run for all wind directions and wind speeds defined. The default properties are:

site.default_wd: 0-360\(^\circ\) in bins of 1\(^\circ\)
site.default_ws: 3-25 m/s in bins of 1 m/s

The default values for wind speeds and wind directions can be overwritten by either a scalar or an array.

[4]:

sim_res = wf_model(x, y,     # wind turbine positions
                   h=None,   # wind turbine heights (defaults to the heights defined in windTurbines)
                   type=0,   # Wind turbine types
                   wd=None,  # Wind direction
                   ws=None,  # Wind speed
                  )

Chunkification and Parallelization

PyWake makes it easy to split the wind farm simulation computation into smaller subproblems. This allows:

Simulation of large wind farms or time series with less memory usage
Parallel execution for faster simulation

1) Chunkfication

To split the simulation into smaller sub-tasks, just specify the desired number of wd_chunks and ws_chunks.

[19]:

# split problem into 4x2 subtasks and simulate sequentially on one CPU
sim_res = wf_model(x, y,
                   wd_chunks=4,
                   ws_chunks=2)

Time series of (wd,ws)-flow cases is split into chunks, by specifying either ws_chunks or wd_chunks

[20]:

sim_res_time = wf_model(x, y, # wind turbine positions
                        wd=wd, # Wind direction time series
                        ws=ws, # Wind speed time series
                        time=time_stamp, # time stamps
                        wd_chunks=4,
                  )

Wind speed or wind direction chunks

[21]:

%%skip
n_lst = [1,2,4,8,16,23,32]
res = {'ws_chunks': np.array([(n,) + profileit(wf_model)(x, y, wd_chunks=1, ws_chunks=n)[1:] for n in tqdm(n_lst)]),
       'wd_chunks': np.array([(n,) + profileit(wf_model)(x, y, wd_chunks=n, ws_chunks=1)[1:] for n in tqdm(n_lst)])}

ax1,ax2 = plt.subplots(1,2,figsize=(12,4))[1]

for k, d in res.items():
    n,t,m = d.T
    ax1.plot(n,t, label=k)
    ax2.plot(n,m, label=k)
setup_plot(ax=ax1, ylabel='Time [s]', xlabel='Number of chunks')
setup_plot(ax=ax2, ylabel='Memory usage [MB]', xlabel='chunks')
plt.savefig('RunWindFarmSimuation_wdws_chunks.svg')

Cell skipped. Precomputed result shown below. Remove '%%skip' to force run.

Result computed on the Sophia HPC cluster. Note that the number of wind speed chunks is limited to the number of wind speeds, in this case 23.

It is clearly seen that chunkification of wind directions are more efficient than chunkification of wind speeds with respect to both time and memory usage.

It is also seen that chunkification reduces the memory usage at the cost of computation time. With 32 wind direction chunks, the computational time is more than doubled, i.e. the speedup from running in parallel on a system with 32 CPUs is expected to be less than 16 times.

2) Parallelization

Running a wind farm simulation in parallel is just as easy - simply specify the number of CPUs to use to the input argument n_cpu or None to use all available CPUs.

As seen above, wind directions chunks are more efficient. Hence, the problem is as default split into n_cpu wind direction chunks, where n_cpu is the numbers of CPUs to use.

As default, n_cpu=1 (sequential execution), so to run in parallel, you need to specify a number of CPUs, e.g. n_cpu=4. Alternatively, n_cpu=None will use all available CPUs.

[22]:

sim_res = wf_model(x, y,
                   n_cpu=None # run wind directions in parallel on all available CPUs
                  )

CPU utilization

The plot below shows the time used to calculate the AEP on 1-32 CPUs as a function of number of wind turbines.

[23]:

%%skip
from py_wake.utils import layouts
from py_wake.utils.profiling import timeit
from tqdm.notebook import tqdm

n_lst = np.arange(100,600,100)

def run(n, n_cpu):
    x,y = layouts.rectangle(n,20,5*windTurbines.diameter())
    return (n, n_cpu, np.mean(timeit(wf_model, min_runs=5)(x,y,n_cpu=n_cpu)[1]))

res = {f'{n_cpu} CPUs': np.array([run(n, n_cpu=n_cpu) for n in tqdm(n_lst)]) for n_cpu in [1, 4, 16, 32]}

ax1,ax2 = plt.subplots(1,2, figsize=(12,4))[1]
for k,v in res.items():
    n,n_cpu,t = v.T
    ax1.plot(n, t, label=k)
    ax2.plot(n, res['1 CPUs'][:,2]/n_cpu/t*100, label=k)
setup_plot(ax=ax1,xlabel='No. wind turbines',ylabel='Time [s]')
setup_plot(ax=ax2,xlabel='No. wind turbines',ylabel='CPU utilization [%]')
plt.savefig('images/RunWindFarmSimulation_time_cpuwt.svg')

Cell skipped. Precomputed result shown below. Remove '%%skip' to force run.

Result precomputed on the Sophia HPC cluster on a node with 32 CPUs.

The plot clearly shows that the time reduces with more CPUs. In this case, however, the gain from 16 to 32 CPUs is very limited. Note, this may highly dependend on the system architecture.

Flow map

Finally, SimulationResult has a flow_map method which returns a FlowMap object. This allows the user to visualize the wake behind each turbine’s rotor in the wind farm and represents the velocity deficit that occurs due to wake interactions in the wind farm.

[24]:

from py_wake.examples.data.iea37 import IEA37Site, IEA37_WindTurbines
from py_wake.literature.gaussian_models import Bastankhah_PorteAgel_2014

site = IEA37Site(16)
x, y = site.initial_position.T
windTurbines = IEA37_WindTurbines()

wf_model = Bastankhah_PorteAgel_2014(site, windTurbines, k=0.0324555)
sim_res = wf_model(x,y)

[25]:

#change the wind speed and wind direction to visualize different flow cases
wsp = 9.8
wdir = 270
flow_map = sim_res.flow_map(grid=None, # defaults to HorizontalGrid(resolution=500, extend=0.2), see below
                            wd=wdir,
                            ws=wsp)

To plot the wake map we call the plot_wake_map property.

You can change the values of the wind speed and wind direction to see how different flow maps look like.

[26]:

plt.figure()
flow_map.plot_wake_map()
plt.xlabel('x [m]')
plt.ylabel('y [m]')
plt.title('Wake map for'+ f' {wdir} deg and {wsp} m/s')

[26]:

Text(0.5, 1.0, 'Wake map for 270 deg and 9.8 m/s')

../_images/notebooks_RunWindFarmSimulation_59_1.png

The wake map can also be customized in terms of size and contour levels.

[27]:

flow_map.plot_wake_map(levels=10, # contourf levels (int or list of levels)
                       cmap='jet', # color map
                       plot_colorbar=False,
                       plot_windturbines=False,
                       ax=None)
plt.axis('equal')
plt.xlabel('x [m]')
plt.ylabel('y [m]')

[27]:

Text(0, 0.5, 'y [m]')

../_images/notebooks_RunWindFarmSimulation_61_1.png

There is also a Grid argument that can be set up to further customize the way the flow map is shown in different planes

The grid argument should be either

a HorizontalGrid(same as XYGrid), YZGrid or
a tuple(X, Y, x, y, h) where X, Y is the meshgrid for visualizing the data and x, y, h are the flattened grid points

HorizontalGrid (XYGrid)

[28]:

from py_wake import HorizontalGrid

for grid in [None,                                                     # defaults to HorizontalGrid(resolution=500, extend=0.2)
             HorizontalGrid(x=None, y=None, resolution=100, extend=1), # custom resolution and extend
             HorizontalGrid(x = np.arange(-1000,1000,100),
                            y = np.arange(-500,500,100))               # custom x and y
            ]:
    plt.figure()
    sim_res.flow_map(grid=grid, wd=270, ws=[9.8]).plot_wake_map()
    plt.xlabel('x [m]')
    plt.ylabel('y [m]')

../_images/notebooks_RunWindFarmSimulation_64_0.png

../_images/notebooks_RunWindFarmSimulation_64_1.png

../_images/notebooks_RunWindFarmSimulation_64_2.png

YZGrid

Plotting the flow map in the vertical YZ plane through the red dashed line can be done using the YZGrid

[29]:

from py_wake import YZGrid
plt.figure(figsize=(12,4))
sim_res.flow_map(YZGrid(x=-100, y=None, resolution=100), wd=270, ws=None).plot_wake_map()
plt.xlabel('y [m]')
plt.ylabel('z [m]')

[29]:

Text(0, 0.5, 'z [m]')

../_images/notebooks_RunWindFarmSimulation_66_1.png

3) XZGrid

Plotting the flow map in the vertical XZ plane through the red dashed line can be done using the XZGrid

[30]:

from py_wake import XZGrid
plt.figure(figsize=(12,4))
sim_res.flow_map(grid=XZGrid(y=0, resolution=1000), wd=270, ws=None).plot_wake_map()
plt.xlabel('x [m]')
plt.ylabel('z [m]')

[30]:

Text(0, 0.5, 'z [m]')

../_images/notebooks_RunWindFarmSimulation_68_1.png

Wind Farm Simulation

Simple simulation - all wind directions and wind speeds

Simulation Results

Selecting data

Plotting data

AEP calculation

Time series

Chunkification and Parallelization

Flow map

HorizontalGrid (XYGrid)

YZGrid