Read vicondar output and create input for the CS framework

python_scripts/from_vicondar_to_cs_files.py

+# -*- 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 
+
+"""
+
+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 
+
+
+# For each file in the folder 
+PATH_VICONDAR_OUTPUT = './'
+
+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)*(z_ref/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)
+        t_lidar = test_mat['time'][0].reshape(-1)
+        
+        # Gets coordinates of the turbulence box
+        dy, dz, Ny, Nz = 6.5, 6.5, 32, 32
+        ly, lz = dy*Ny, dz*Nz
+        Yvalues = dy/2 + np.arange(0,(ny),1)*dy
+        Zvalues = dz/2 + np.arange(0,(nz),1)*dz
+              
+        # Gets coordinates of the lidar scanning patterns (arbitrary example)         
+        LidarSys = 'circ7p' # Name of the lidar 
+        
+        Rrotor = 89.5 
+        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]))]
+        y_lidar,z_lidar = Yvalues[idx_y]-cc,Zvalues[idx_z]-cc
+        
+        # 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(z_ref,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)))  
+        filename.write("\nInput['Vdatafile']=" + "'turb%s_av_setbv.bin'" %(tu_seed))
+        filename.write("\nInput['Wdatafile']=" + "'turb%s_av_setbw.bin'" %(tu_seed))
+        filename.write("\nInput['LidarSys']=" + "'%s'" %LidarSys)
+        filename.write("\nInput['CrossSpectrumFilenames']=" + "'./MannCrossSpectrum_Gamma'")
+        filename.write("\nInput['CalculateSpectrum']=" + '%d' %2)
+        filename.write("\nInput['SpectrumSteps']=" + '%d' %128)
+        filename.write("\nInput['UseBoxCovariance']=" + '%d' %0)
+        filename.write("\nInput['Gamma']=" + '%.2f' %3.9) # Mann spectral velocity tensor parameters
+        filename.write("\nInput['L']=" + '%.4f' %29.4) # 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)
+        filename.write("\nInput['PrintDetailedOutput']=" + '%d' %0)
+        filename.write("\nInput['CrossSpectrumComponent1']=" + '%d' %1)
+        filename.write("\nInput['CrossSpectrumComponent2']=" + '%d' %1)
+        filename.write("\nInput['OutputFileName']=" + "'./turb%s_ti%s_av_setb_%s_c_vol%du.bin'" %(tu_seed,int(TI_amb),LidarSys,vsize))
+        filename.write("\nInput['Constraints']= [")  
+        filename.write(", \n".join(constraint_cells))
+        filename.write("]\n")
+        filename.write("\nprint('Start applying constraints:')")
+        filename.write("\nConstrainTurbulenceBox_U_S(Input)")
+        
+        # Close the file after writing to it.
+        filename.close()
+       
\ No newline at end of file