Newer
Older
# -*- coding: utf-8 -*-
"""
Created on Sat Mar 14 09:49:17 2020
@author: davcon
This script reads output files from the VICONDAR framework and
create inputs to the constrained simulation framework.
CS: stands for Constrained simulation framework
Note that as a convention the output filenames from vicondar
have wind field characteristics incorporated, e.g.,:
"uval_av_sens_dist_Sh20_SD01_V06_TI08_s1000_circ7p_vol30.mat"
Sh (Shear) = 0.2 [-]
V (Wind speed) = 6 m/s
TI (Turbulence intensity) = 0.08 [-]
s (seed) = 1000
circ7p = lidar scanning configuration's name
vol (lidar's probe volume) = 30 [m]
The user can easily change this convention and adjust this script accordingly.
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
"""
import numpy as np
import scipy.io
import os
from scipy.interpolate import interp1d
import numpy as np
import matplotlib.pyplot as plt
import os
import scipy.io as sio
import scipy
import scipy.optimize as optimization
from scipy.interpolate import interp1d,interp2d
from scipy.interpolate import griddata
import pickle
from scipy.optimize import minimize
import fnmatch
from datetime import datetime, timedelta
import re
from mann_parameters import *
from mann_parameters import var2ae
import pandas as pd
# User-defined inputs
PATH_VICONDAR_OUTPUT = './'
# Turbulence box inputs
dy, dz, Ny, Nz, Nx = 6.5, 6.5, 32, 32, 8192
ly, lz = dy*Ny, dz*Nz
Yvalues = dy/2 + np.arange(0,(Ny),1)*dy
Zvalues = dz/2 + np.arange(0,(Nz),1)*dz
Tsim = 600 # Simulation time in seconds
# Arbritary lidar scanning patterns
LidarSys = 'circ7p' # Name of the lidar
Rrotor = 89.5 # Turbine parameters (rotor)
hh = 119 # Turbine parameters (hub-height)
frac = Rrotor* 0.85
phi = np.deg2rad(60)
lidar_7pcirc = np.array([-frac,-frac*np.cos(phi),0,frac*np.cos(phi),frac,frac*np.cos(phi),-frac*np.cos(phi)]),\
np.array([0,-frac*np.sin(phi),0,frac*np.sin(phi),0,-frac*np.sin(phi),frac*np.sin(phi)])
lidar_7pcirc = lidar_7pcirc + Yvalues[int(len(Yvalues)/2)]
idx_y = [np.argmin(np.abs(lidar_7pcirc[0][i]- Yvalues)) for i in range(len(lidar_7pcirc[0]))]
idx_z = [np.argmin(np.abs(lidar_7pcirc[1][i]- Yvalues)) for i in range(len(lidar_7pcirc[0]))]
cc = Yvalues[int(len(Yvalues)/2)] # center of the box
y_lidar,z_lidar = Yvalues[idx_y]-cc,Zvalues[idx_z]-cc
# Mann parameters
L_mann = 29.4
Gamma_Mann = 3.9
for filename in os.listdir(PATH_VICONDAR_OUTPUT):
if fnmatch.fnmatch(filename, '*.mat'):
print(filename)
# Get wind field statitistcs and turbulence seed number etc. from the filename
temp=re.findall(r'\d+', filename)
# Shear
alpha = temp[0]
# Mean wind speed at hub-height
U_amb = int(temp[2])
# Ambient turbulence
TI_amb = temp[3]
# Turbulence seed
seed = temp[4]
tu_seed = int(seed)
# Probe volume size of the lidar
vsize = int(temp[6])
# Vertical wind profile
Uz = int(U_amb)*(Zvalues/hh)**(int(alpha)/100)
# Load VICONDAR outputs
test_mat=scipy.io.loadmat(PATH_VICONDAR_OUTPUT + '/%s' %filename)
u_lidar = test_mat['uval'][0].reshape(-1)
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
# Center the lidar pattern to the turbulence box
y_lidar = Yvalues[int(len(Yvalues)/2)] + y_lidar[~np.isnan(y_lidar)]
z_lidar = Zvalues[int(len(Zvalues)/2)] + z_lidar[~np.isnan(z_lidar)]
# Number of points per scan
ppt_scans = np.shape(u_lidar)[0]
# Number of scan within the time period (e.g., 10-min period)
nr_scans = np.shape(u_lidar[0].reshape(-1))[0]
# Initialize variables
u_lidar_tot = np.zeros(shape=(ppt_scans,nr_scans))
t_lidar_tot = np.zeros_like(u_lidar_tot); y_lidar_tot = np.zeros_like(u_lidar_tot)
z_lidar_tot = np.zeros_like(u_lidar_tot)
# Reshape vectors according to VICONDAR notation
for i in range(ppt_scans):
u_lidar_tot[i,:]= u_lidar[i].reshape(-1)
t_lidar_tot[i,:]= t_lidar[i].reshape(-1)
y_lidar_tot[i,:]= z_lidar[i].reshape(-1)
z_lidar_tot[i,:]= y_lidar[i].reshape(-1)
# Interpolate the vertical ambient wind speed at the location of the constraints
u_amb_z = interp1d(Zvalues,Uz)(z_lidar)
# Compute wind speed residuals to be input to the constrained simulation framework
u_lidar_res = np.asarray([(u_lidar_tot[i_coord,:] -u_amb_z[i_coord]) for i_coord in range(ppt_scans)])
# Reshape vectors according to the CS convention
t_cnst_array = (t_lidar_tot.T).reshape(-1)
u_res_cnst_array = (u_lidar_res.T).reshape(-1)
z_cnst_array = (z_lidar_tot.T).reshape(-1)
y_cnst_array = (y_lidar_tot.T).reshape(-1)
# Create a single array with all info to be input to the CS framework
constraint_array =np.zeros(shape=(len(t_cnst_array),4))
for iC in range(0,len(t_cnst_array)):
constraint_array[iC,:] = np.array([t_cnst_array[iC], y_cnst_array[iC],z_cnst_array[iC],u_res_cnst_array[iC]])
# Fit ae (of the Mann model) from TI statistics
var_u = ((int(TI_amb)/100)*int(U_amb))**2
ae = var2ae(var_u, L=29.4, G=3.9)
idSim_idx = 0
# Name of the 'wake/simulation' case
iWake = 'av'
# ==============================================
# Create the CS file for simulations
constraint_cells=[]
constraint_cells=['[%.4f, %.4f, %.4f, %.4f]' %( t_cnst_array[iC], y_cnst_array[iC],z_cnst_array[iC],u_res_cnst_array[iC]) for iC in range(0,len(t_cnst_array))]
title = "./cs_seed%s_%s_%s_vol%d" %(seed,iWake,LidarSys,vsize)
# Add .py to the end of the script.
title = title + '.py'
# Convert all letters to lower case.
title = title.lower()
# Remove spaces from the title.
title = title.replace(' ', '_')
# Create a file that can be written to.
filename = open(title, 'w')
# Write the data to the file.
filename.write('\n#Script for imposing constraints in a turbulence box, seed no.%s' %tu_seed)
filename.write('\n#INPUT TO CONSTRAINING PROCEDURE\n')
filename.write('\nimport pickle')
filename.write('\nimport pandas as pd')
filename.write('\nimport numpy as np')
filename.write('\nimport matplotlib.pyplot as plt')
filename.write('\nimport os')
filename.write('\nimport scipy.io')
filename.write('\nfrom scipy import stats')
filename.write('\nimport scipy.io as sio')
filename.write('\nfrom MannTensor_lib import MannTensor,TrapezoidalSum2D')
filename.write('\nimport cmath')
filename.write('\nfrom scipy.interpolate import interp1d')
filename.write('\nfrom ConstrainTurbulenceBox_U_S import ConstrainTurbulenceBox_U_S\n')
filename.write('\nfrom ConstrainTurbulenceBox_U_S import ConstrainTurbulenceBox_U_S\n')
filename.write('\n#Initialization of variables:')
filename.write('\nInput={} \n')
filename.write("\nInput['Udatafile']=" + "'./turb%s_ti%s_av_setbu.bin'" %(tu_seed,int(TI_amb))) # Input name of the unconstrained turb. box
filename.write("\nInput['Vdatafile']=" + "'turb%s_av_setbv.bin'" %(tu_seed)) # Input name of the unconstrained turb. box
filename.write("\nInput['Wdatafile']=" + "'turb%s_av_setbw.bin'" %(tu_seed)) # Input name of the unconstrained turb. box
filename.write("\nInput['LidarSys']=" + "'%s'" %LidarSys) # Lidar configuration
filename.write("\nInput['CrossSpectrumFilenames']=" + "'./MannCrossSpectrum_Gamma'")
filename.write("\nInput['CalculateSpectrum']=" + '%d' %1) # Inputs required for the constrained simulation algorithm
filename.write("\nInput['SpectrumSteps']=" + '%d' %128) # Inputs required for the constrained simulation algorithm
filename.write("\nInput['UseBoxCovariance']=" + '%d' %0) # Inputs required for the constrained simulation algorithm
filename.write("\nInput['Gamma']=" + '%.2f' %Gamma_Mann) # Mann spectral velocity tensor parameters
filename.write("\nInput['L']=" + '%.4f' %L_mann) # Mann spectral velocity tensor parameters
filename.write("\nInput['AlphaEpsilon']=" + '%.4f' %ae) # Mann spectral velocity tensor parameters
filename.write("\nInput['Nx']=" + '%d' %Nx) # Mann turbulence box dimensions
filename.write("\nInput['Ny']=" + '%d' %Ny) # Mann turbulence box dimensions
filename.write("\nInput['Nz']=" + '%d' %Nz) # Mann turbulence box dimensions
filename.write("\nInput['T']=" + '%d' %Tsim) # Simulation time in seconds
filename.write("\nInput['Umean']=" + '%d' %U_amb) # Ambient wind speed at hub-height
filename.write("\nInput['BoxWidth']=" + '%d' %(Ny*dy)) # in meters
filename.write("\nInput['BoxHeight']=" + '%d' %(Nz*dz)) # in meters
filename.write("\nInput['SaveToFile']=" + '%d' %1)# Flag for saving constrained turb. fields
filename.write("\nInput['CrossSpectrumComponent1']=" + '%d' %1) # Flag for calculating the uu-spectrum
filename.write("\nInput['CrossSpectrumComponent2']=" + '%d' %1) # Flag for calculating the uu-spectrum
filename.write("\nInput['OutputFileName']=" + "'./turb%s_ti%s_av_setb_%s_c_vol%du.bin'" %(tu_seed,int(TI_amb),LidarSys,vsize)) # Output name of the constrained turb. box
filename.write(", \n".join(constraint_cells)) # Print the constraints and their locations