User guide¶

Save turbulence box¶

[7]:

# make tmp folder
os.makedirs('tmp', exist_ok=True)

Netcdf¶

The turbulence box can be saved as netcdf (one file including all attributes)

[8]:

mtf.to_netcdf(folder='tmp', filename=None)

If filename=None (default), the name is autogenerated from the properties:

[9]:

mtf.name

[9]:

'Hipersim_mann_l33.6_ae0.1000_g3.9_h0_1024xd64xd64_4.000x2.00x2.00_s0001'

HAWC2¶

The turbulence box can be saved as hawc2 (three binary files) by

[10]:

mtf.to_hawc2(folder='tmp', basename=None,htc_dict={'wind.wsp':9})

[10]:

<wetb.hawc2.htc_contents.HTCSection at 0x1ae00b0bcd0>

[11]:

os.listdir('tmp')

[11]:

['Hipersim_mann_l33.6_ae0.1000_g3.9_h0_1024xd64xd64_4.000x2.00x2.00_s0001.nc',
 'Hipersim_mann_l33.6_ae0.1000_g3.9_h0_1024xd64xd64_4.000x2.00x2.00_s0001u.turb',
 'Hipersim_mann_l33.6_ae0.1000_g3.9_h0_1024xd64xd64_4.000x2.00x2.00_s0001v.turb',
 'Hipersim_mann_l33.6_ae0.1000_g3.9_h0_1024xd64xd64_4.000x2.00x2.00_s0001w.turb']

Run HAWC2 simulation and compare with turbulence field¶

The to_hawc2 method returns the a wind section for the HAWC2 htc input file. Its contents can be overriden via the htc_dict argument

[59]:

htc_wind_section = mtf.to_hawc2(folder='tmp', basename=None,htc_dict={'wind.wsp':9, 'wind.center_pos0':[0,0,-100]})
print (htc_wind_section)

begin wind;
  density       1.225;
  wsp   9;
  tint  1;
  horizontal_input      1;      0=false, 1=true
  windfield_rotations   0 0 0;  yaw, tilt, rotation
  center_pos0   0 0 -100;       hub heigth
  shear_format  1 0;    0=none,1=constant,2=log,3=power,4=linear
  turb_format   1;      0=none, 1=mann,2=flex
  tower_shadow_method   0;      0=none, 1=potential flow, 2=jet
  begin mann;
    filename_u  tmp/Hipersim_mann_l33.6_ae0.1000_g3.9_h0_1024xd64xd64_4.000x2.00x2.00_s0001u.turb;
    filename_v  tmp/Hipersim_mann_l33.6_ae0.1000_g3.9_h0_1024xd64xd64_4.000x2.00x2.00_s0001v.turb;
    filename_w  tmp/Hipersim_mann_l33.6_ae0.1000_g3.9_h0_1024xd64xd64_4.000x2.00x2.00_s0001w.turb;
    box_dim_u   1024 4;
    box_dim_v   64 2;
    box_dim_w   64 2;
    dont_scale  1;
  end mann;
end wind;
;
;----------------------------------------------------------------------------------------------------------------------------------------------------------------
;
;

Run HAWC2 with a sensor in \((x,y,z)=(0,5,10)\) and load result

[61]:

from wetb.hawc2.htc_file import HTCFile
from h2lib import H2Lib
from wetb.gtsdf import gtsdf
htc = HTCFile(modelpath='.')
htc.wind = htc_wind_section
# HAWC2 global coordinate system: [y,x,-z]
# box center in h2_gl: [0,0,-100]
# box origo in h2_gl [-63,0,-37]
# point (0,5,10) in [-58,0,-47]
center_pos = htc.wind.center_pos0
htc.output.data_format = 'gtsdf64'
htc.output.add_sensor('wind', 'free_wind', [1,-58,0,-47])

htc.set_name('tmp')
htc.save()
with H2Lib(suppress_output=0) as h2:
    h2.init(htc.filename)
    h2.run(100)
time, data, info = gtsdf.load(htc.output_files()[0])

Plot wind speed from HAWC2 sensor and turbulence box interpolation

[58]:

h2_uvw = [data[:,1], data[:,0],-data[:,2]]
ws = htc.wind.wsp[0]
x = -time * ws # convert time to box x coordinate
turb_uvw = mtf(x, x*0+5, x*0+10).T + np.array([ws,0,0])[:,np.newaxis]

for uvw, h2, turb in zip('uvw',h2_uvw, turb_uvw):
    plt.plot(time, turb, label=f'Turb {uvw}')
    plt.plot(time, h2, '--', label=f'Hawc2 {uvw}')
    np.testing.assert_array_almost_equal(h2,turb,7)
plt.legend()

[58]:

<matplotlib.legend.Legend at 0x1ae02b98310>

Interpolation¶

The MannTurbulenceField has a buildin linear interpolation function that handles out-of-box interpolation depending on whether the axes have been doubled. For axes which have been doubled (default for y and z), the box is mirrored, while axes that have not been doubled (default for x), the box is repeated.

[12]:

ip_x, ip_z = np.meshgrid(np.linspace(-1000,5000,2000), np.linspace(-100,400,200))
ip_y = ip_x*0

[13]:

u = mtf(ip_x, ip_y, ip_z)[...,0]

[14]:

plt.figure(figsize=(16,2))
plt.contourf(ip_x,ip_z,u.reshape(ip_x.shape), cmap='RdBu_r')
x1,y1 = mtf.Nx*mtf.dx, mtf.Nz*mtf.dz
plt.plot([0,x1, x1, 0,0],[0,0, y1,y1,0],'k', label='box')
plt.legend()
plt.axis('scaled')

[14]:

(-1000.0, 5000.0, -100.0, 400.0)

The raw turbulence box is marked with the black rectangle. The remaining turbulence is mirrored (vertical direction) or repeated (longitudinal direction)

Multiple turbulence boxes (reuse spectral tensor)¶

In case, you need to generate multiple turbulence boxes where only the seed and alphaepsilon varies, the computational time can be reduced by reusing the spectral tensor. The eaiest way to do this is set the argument cache_spectral_tensor=True

[15]:

for seed, alphaepsilon in [(1,1), (1, 2), (2, 2)]:
    mtf = MannTurbulenceField.generate(alphaepsilon=alphaepsilon,
                                       L=33.6,
                                       Gamma=3.9,
                                       Nxyz=(1024,64,64),
                                       dxyz=(4,2,2),
                                       seed=seed,
                                       HighFreqComp=0,
                                       double_xyz=(False, True, True),
                                       cache_spectral_tensor=True) # reuse spectal tensor
    plt.figure(figsize=(18,2))
    mtf.to_xarray().sel(uvw='u',y=0).plot(x='x')
    plt.title(f'Seed: {seed}, alphaepsilon: {alphaepsilon}')

In this case the cached spectral tensor is stored in the file mannsqrtphi_l33.6_g3.9_h0_1024xd64xd64_4.000x2.00x2.00.npy. Note, remember to delete the cached spectral tensor files after use as they can be quite big, in this case around 300MB

[16]:

!ls mannsqrtphi_l33.6_g3.9_h0_1024xd64xd64_4.000x2.00x2.00.npy


(.py311_dynamiks_jupyter7) ]9;8;"USERNAME"\@]9;8;"COMPUTERNAME"\ c:\mmpe\programming\python\Topfarm\hipersim\doc
$ ]9;12\dir mannsqrtphi_l33.6_g3.9_h0_1024xd64xd64_4.000x2.00x2.00.npy
 Volume in drive C is Windows
 Volume Serial Number is 4637-B881

 Directory of c:\mmpe\programming\python\Topfarm\hipersim\doc

28-10-2024  14:51       302.579.840 mannsqrtphi_l33.6_g3.9_h0_1024xd64xd64_4.000x2.00x2.00.npy
               1 File(s)    302.579.840 bytes
               0 Dir(s)  677.719.523.328 bytes free

Spectra¶

The spectra of the turbulence realization can be obtained by the spectra method

[17]:

k_rea, spectra_rea = mtf.spectra()

Furthermore, the theoretical spectra of the spectral tensor can be found by integration

[18]:

k_int, spectra_int = mtf.spectra_integrated()

or from a build-in look-up table

[19]:

k_lut, spectra_lut = mtf.spectra_lookup()

[20]:

c_lst = ['#1f77b4', '#ff7f0e', '#2ca02c', '#d62728']
s_lst = ["uu","vv","ww","uw"]

[21]:

for phi_rea, phi_int, phi_lut, c, l in zip(spectra_rea, spectra_int, spectra_lut, c_lst, s_lst):
    plt.semilogx(k_rea, phi_rea*k_rea, color=c, label=l)
    plt.semilogx(k_int, phi_int*k_int, '--', color=c)
    plt.semilogx(k_lut, phi_lut*k_lut, ':', color=c)

plt.plot([], '-', color='gray', label='Realization')
plt.plot([], '--', color='gray', label='LookupTable')
plt.plot([], ':', color='gray', label='Integrated')
plt.xlabel('Wave number, $k_1$ $  [m^{-1}$]')
plt.ylabel('$k_1 S(k_1)[m^2s^{-2}]$')
plt.legend()

[21]:

<matplotlib.legend.Legend at 0x25484c52010>

Coherence¶

The coherence of the turbulence in two points separated by (dy, dz), can be obtained by the coherence method

Furthermore, the theoretical coherence of the spectral tensor can be found by integration

[22]:

dy = 2*mtf.dy # 4m
dz = 3*mtf.dz # 6m

[23]:

for comp in ['u','v','w']:
    k_int, coherence_int = mtf.coherence_integrated(dy,dz, component=comp)
    k_rea, coherence_rea = mtf.coherence(dy,dz, component=comp)
    plt.figure()
    s = r'$\Delta_y$ = %.2fm, $\Delta_z$ = %.2fm' % (dy, dz)
    plt.title(f'{comp}, {s}')
    plt.plot(k_int, coherence_int, label='Theoretical')
    plt.plot(k_rea, coherence_rea, label='Realization')
    plt.xlabel('Wave number, $k_1$ [$m^{-1}$]')
    plt.ylabel('co-coherence [-]')
    plt.legend()

Variance¶

Theoretical spectrum¶

The variance of the theoretical Mann model spectrum can be obtained by the spectrum_variance method

[24]:

mtf.spectrum_variance()

[24]:

38.559824397503185

The variance is obtained by discrete integration of the theoretical spectrum

[25]:

k = mtf.get_k()
k_lut, spectra_lut = mtf.spectra_lookup(k)

def plot_variance_integration(k,S):
    plt.semilogx(k, k*S[0],'.-')
    dk = np.diff(k[:2])
    variance = 2 * np.sum(S[0]*dk)
    plt.bar(k, S[0]*k, width=dk, alpha=.5, label=f'var={variance:.3f}')
    plt.title('uu spectrum')
    plt.xlabel('Wave number, $k_1$ $  [m^{-1}$]')
    plt.ylabel('$k_1 S(k_1)[m^2s^{-2}]$')
    plt.legend()
plot_variance_integration(k_lut,spectra_lut)

As default get_variance obtains the wave number limits from the length of the box, Nx*dx and the distance between grid points, dx, but other values can be specified. This is usefull if you want to compare the variance of a 1 hour box with the variance of e.g. 10min measurements.

[26]:

mtf.spectrum_variance(Lx=1024, dx=16)

[26]:

26.040102636402843

[27]:

k = mtf.get_k(Lx=1024, dx=16)
k_lut, spectra_lut = mtf.spectra_lookup(k)
plot_variance_integration(k_lut, spectra_lut)

Realization¶

The variance can be obtained from the realization as the mean of the variance of each grid point in the (y,z) plane:

[28]:

mtf.uvw[0].var(0).mean()

[28]:

38.490097

or by discrete integrating the spectrum:

[29]:

k_rea, spectra_rea = mtf.spectra(log10_bin_size=None) # skip binning
plot_variance_integration(k_rea, spectra_rea)

Turbulence intensity¶

[30]:

mtf_lst = [MannTurbulenceField.generate(alphaepsilon=.1,L=33.6,Gamma=3.9,
                                   Nxyz=(1024,64,64), dxyz=(4,2,2), seed=seed,
                                   HighFreqComp=0, double_xyz=(False, True, True),
                                   cache_spectral_tensor=True)
           for seed in [0,4,8]]
mtf = mtf_lst[0]
U = 10

Theoretical spectrum¶

The turbulence intensity of the theoretical Mann model spectrum can be obtained with the spectra_TI method

[31]:

mtf.spectrum_TI(U=U)

[31]:

0.1388521234938508

Wave number limits can be specified via signal length, T, and cut-off frequency, cutoff_frq, associated with the TI measure. The cut-off frequency is e.g. caused by inertia in the cups of a cup anemometer and probe volume of a lidar.

[32]:

mtf.spectrum_TI(U=U, T=600, cutoff_frq=10)

[32]:

0.14324855863543767

Realization¶

The turbulence intensity of a box realization can be obtained by

[33]:

mtf_lst[0].uvw[0].std(0).mean() / U

[33]:

0.1373198390007019

Realizations may, however have very different turbulence intensities due to uncertainty on the low frequencies

[34]:

T = mtf.dx * mtf.Nx / U # Box length [s]
cutoff_frq = U / (2 * mtf.dx) # Nyquist frequency [Hz]

# Calculate and plot theoretical spectrum and ti
k_lut, spectra_lut = mtf.spectra_lookup(k1 = mtf.get_k())
spectra_ti = mtf.spectrum_TI(U=10)
plt.semilogx(k_lut, k_lut*spectra_lut[0], '--k', label=f'Theoretical, ti: {np.round(spectra_ti,3):0.3f}')

# Calculate and plot spectrum and ti of the four realizations
for mtf in mtf_lst:
    k_rea, spectra_rea = mtf.spectra(min_bin_count=1, log10_bin_size=None)
    plt.semilogx(k_rea, k_rea*spectra_rea[0], label=f'Seed: {mtf.seed}, ti: {np.round(np.mean(mtf.uvw[0].std(0)) / U,3):0.3f}')

plt.title('uu spectrum')
plt.xlabel('Wave number, $k_1$ $  [m^{-1}$]')
plt.ylabel('$k_1 S(k_1)[m^2s^{-2}]$')
plt.legend()

[34]:

<matplotlib.legend.Legend at 0x25483e6b150>

Turbulence scaling¶

In case, a specific turbulence intensity is needed, the turbulence field can be scaled accordingly.

[35]:

for i, mtf in enumerate(mtf_lst):
    print(f"Box {i}")
    print (f'Before: Box TI={mtf.uvw[0].std(0).mean()/U:.3f}, alphaepsilon:{mtf.alphaepsilon:.3f}, theoretical spectrum TI {mtf.spectrum_TI(U):.2f}')
    mtf.scale_TI(TI=0.1, U=10)
    print (f'After: Box TI={mtf.uvw[0].std(0).mean()/U:.3f}, alphaepsilon:{mtf.alphaepsilon:.3f}, theoretical spectrum TI {mtf.spectrum_TI(U):.2f}')
    print()

Box 0
Before: Box TI=0.137, alphaepsilon:0.100, theoretical spectrum TI 0.14
After: Box TI=0.099, alphaepsilon:0.052, theoretical spectrum TI 0.10

Box 1
Before: Box TI=0.167, alphaepsilon:0.100, theoretical spectrum TI 0.14
After: Box TI=0.120, alphaepsilon:0.052, theoretical spectrum TI 0.10

Box 2
Before: Box TI=0.113, alphaepsilon:0.100, theoretical spectrum TI 0.14
After: Box TI=0.081, alphaepsilon:0.052, theoretical spectrum TI 0.10

Get alpha epsilon¶

The \(\alpha\epsilon^{2/3}\) value that corresponds to a theoretical spectrum with a specific TI can be obtained by

[36]:

mtf.get_alpha_epsilon(TI=.1, U=U)

[36]:

0.051867456121753086

It can even be calculated without generating the turbulence box:

[37]:

from hipersim.turbgen.spectral_tensor import MannTurbulenceInput
mti = MannTurbulenceInput(alphaepsilon=1, # dummy value
                          L=33.6,Gamma=3.9,
                          Nxyz=(1024,64,64), dxyz=(4,2,2))
mti.get_alpha_epsilon(TI=.1, U=U)

[37]:

0.051867456121753086

High frequency compensation¶

The HighFreqComp argument defines whether high-frequency compensation is applied. There are three options:

0 or False (default): No high-frequency compensation applied
1 or True: A fast high-Frequency compensation method is applied. This method differs from the method in Mann (1998)
2: The high-Frequency compensation method from the C++ version is applied. The method corresponds to Eq. A.6 in Mann (1998) except that the convolution is only applied in the (k2,k3), i.e. -2<=n_l<=2, l=(2,3). Note that this method is very slow

[38]:

axes = plt.subplots(2,2, figsize=(10,8))[1].flatten()

for hfc in [0,1,2]:
    t = time.time()
    mtf = MannTurbulenceField.generate(Nxyz=(256,64,64), HighFreqComp=hfc, double_xyz=(0,0,0), verbose=0)
    k,S = mtf.spectra()
    print (f"Spectra with HighFreqComp={hfc} computed in {time.time()-t:.0f}s")
    for ax, Phi, in zip(axes, S):
        ls = ['-', '--'][hfc==2]
        ax.semilogx(k, Phi*k, ls, label=f'HighFreqComp={hfc}')

k,S = mtf.spectra_lookup()
for ax, Phi, s in zip(axes, S, s_lst):
    ax.set_title(s)
    ax.semilogx(k, Phi*k, label=f'Theoretical', zorder=-32)
    ax.set_xlabel('Wave number, $k_1$ $  [m^{-1}$]')
    ax.set_ylabel('$k_1 S(k_1)[m^2s^{-2}]$')
    ax.grid()
    ax.legend()
    ax.set_xlim([0.1,5])
    ax.set_ylim([0,.5])
axes[3].set_ylim([-0.1,.02])
plt.tight_layout()

Spectra with HighFreqComp=0 computed in 0s
Spectra with HighFreqComp=1 computed in 1s
Spectra with HighFreqComp=2 computed in 22s

Contraints¶

It is possible to modify a turbulence field

[39]:

mtf = MannTurbulenceField.generate(Nxyz=(128, 16, 8), dxyz=(2, 3, 4))

to satisfy a number of constraints

[40]:

constraints = np.array([[0, 0, 4, 5, 1, 2], # constr 1, x,y,z, u,v,w
                        [50, 0, 4, 6, 2, 1],
                        [100, 0, 4, 5, 1, 2]])

via the constrain method

[41]:

axes = plt.subplots(3, 1, figsize=(10,10))[1]
for ax, c in zip(axes, 'uvw'):
    mtf.to_xarray().sel(uvw=c, y=0, z=4).plot(ax=ax, label='Raw turbulence')

mtf.constrain(constraints)


for i, (ax, c) in enumerate(zip(axes, 'uvw'), 3):
    mtf.to_xarray().sel(uvw=c, y=0, z=4).plot(ax=ax, label='Constrained turbulence')

    ax.plot([], '.k', label='Contraints')
    for con in constraints:
        ax.plot(con[0], con[i], '.k')
    ax.legend()

Computing Mann tensor / correlation arrays for constrained simulation:
Correlation computations complete
Time elapsed is 0.02411961555480957
Populating covariance matrix for the constraints:
Constraint-constraint covariance matrix has been assembled
Time elapsed is 0.0010647773742675781
Applying constraints...
Constrained simulation complete
Total time elapsed is 0.027193307876586914

User guide¶

Generate and access turbulence field¶

Save turbulence box¶

Netcdf¶

HAWC2¶

Run HAWC2 simulation and compare with turbulence field¶

Interpolation¶

Multiple turbulence boxes (reuse spectral tensor)¶

Spectra¶

Coherence¶

Variance¶

Theoretical spectrum¶

Realization¶

Turbulence intensity¶

Theoretical spectrum¶

Realization¶

Turbulence scaling¶

Get alpha epsilon¶

High frequency compensation¶

Contraints¶

Hipersim

Navigation

Related Topics