Readding aeroacoustic examples

openfast_toolbox/aeroacoustics/examples/_plot_directivity.py

@@ -0,0 +1,86 @@
+import numpy as np
+import pandas as pd
+import matplotlib.pyplot as plt
+from parse import *
+import re, os, platform
+from openfast_toolbox.io.fast_output_file import FASTOutputFile
+
+#########################################################################################################################################
+## User inputs
+# Save plot and/or data?
+save_fig = False
+save_data = False
+fig_ext = '.png'
+
+# Number of revolutions (n) to average spectra
+n = 1
+
+#########################################################################################################################################
+## Paths to files
+if platform.system() == 'Windows':
+ FAST_directory = os.path.dirname( os.path.dirname( os.path.dirname( os.path.dirname( os.path.realpath(__file__) ) ) ) ) + os.sep + 'reg_tests' + os.sep + 'r-tests' + os.sep + 'glue-codes' + os.sep + 'openfast' + os.sep + 'IEA_LB_RWT-AeroAcoustics'
+else:
+ FAST_directory = os.path.dirname( os.path.dirname( os.path.dirname( os.path.dirname( os.path.realpath(__file__) ) ) ) ) + os.sep + 'openfast' + os.sep + 'reg_tests' + os.sep + 'r-tests' + os.sep + 'glue-codes' + os.sep + 'openfast' + os.sep + 'IEA_LB_RWT-AeroAcoustics'
+AAfilename = FAST_directory + os.sep + 'IEA_LB_RWT-AeroAcoustics_1.out'
+OFfilename = FAST_directory + os.sep + 'IEA_LB_RWT-AeroAcoustics.out'
+locfilename = FAST_directory + os.sep + 'AA_ObserverLocations_Map.dat'
+output_dir = os.path.dirname( os.path.realpath(__file__) )
+outputfilename = output_dir + os.sep + "data_output1"
+
+#########################################################################################################################################
+## Read in data, manipulate it, and plot it
+# reads in file data
+AA_1 = FASTOutputFile(AAfilename).toDataFrame()
+OF = FASTOutputFile(OFfilename).toDataFrame()
+location = pd.read_csv(locfilename,delimiter='\s+',skiprows=[0,1],names=['x','y','z'])
+
+# determine number of observers
+num_obs = AA_1.shape[1]-1
+
+# calculate sample time for n revolutions
+rpm = OF[["RotSpeed_[rpm]"]].mean()[0]
+yaw = OF[["YawPzn_[deg]"]].mean()[0] / 180. * np.pi
+time_revs = n*60/rpm
+tot_time = AA_1["Time_[s]"].max()
+if time_revs < tot_time:
+ sample_time = tot_time - time_revs
+else:
+ print("Error: Time for number of revolutions exceeds simulation time. Reduce n.")
+ raise SystemExit('')
+
+# slice AA dataframe for t > sample_time
+AA_1 = AA_1[AA_1["Time_[s]"] > sample_time]
+AA_1=AA_1.drop("Time_[s]",axis=1)
+
+# average P over rotor revolution
+AA_1 = AA_1.mean()
+
+# merge location info with SPL info
+AA_1=AA_1.reset_index()
+AA_1=AA_1.drop("index",axis=1)
+AA_1=pd.merge(location,AA_1,left_index=True,right_index=True)
+AA_1=AA_1.rename(index=str,columns={0:"SPL"})
+
+# contour plot of SPL for each location
+if num_obs < 3:
+ print("Error: Need at least 3 observers to generate contour.")
+else:
+ x=AA_1['x'];
+ y=AA_1['y'];
+ z=AA_1['SPL'];
+ fs = 10
+ fig,ax=plt.subplots()
+ ax.set_aspect('equal')
+ ax.set_xlabel('x [m]', fontsize=fs+2, fontweight='bold')
+ ax.set_ylabel('y [m]', fontsize=fs+2, fontweight='bold')
+ tcf=ax.tricontourf(x,y,z, range(58, 84, 1))
+ fig.colorbar(tcf,orientation="vertical").set_label(label = 'Overall SPL [dB]', fontsize=fs+2,weight='bold')
+ if save_fig == True:
+ fig_name = 'directivity_map'
+ fig.savefig(output_dir + os.sep + fig_name + fig_ext)
+ plt.show()
+
+# export to csv
+if save_data == True:
+ AA_1.to_csv(r'{}-data.csv'.format(outputfilename))
+

openfast_toolbox/aeroacoustics/examples/_plot_mechanisms.py

@@ -0,0 +1,112 @@
+import numpy as np
+import pandas as pd
+import matplotlib.pyplot as plt
+from parse import *
+import re, os, platform
+from openfast_toolbox.io.fast_output_file import FASTOutputFile
+
+#########################################################################################################################################
+## User inputs
+# Save plot and/or data?
+save_fig = False
+save_data = False
+fig_ext = '.png'
+
+# Number of revolutions (n) to average spectra
+n = 1
+
+#########################################################################################################################################
+## Paths to files
+if platform.system() == 'Windows':
+ FAST_directory = os.path.dirname( os.path.dirname( os.path.dirname( os.path.dirname( os.path.realpath(__file__) ) ) ) ) + os.sep + 'reg_tests' + os.sep + 'r-tests' + os.sep + 'glue-codes' + os.sep + 'openfast' + os.sep + 'IEA_LB_RWT-AeroAcoustics'
+else:
+ FAST_directory = os.path.dirname( os.path.dirname( os.path.dirname( os.path.dirname( os.path.realpath(__file__) ) ) ) ) + os.sep + 'openfast' + os.sep + 'reg_tests' + os.sep + 'r-tests' + os.sep + 'glue-codes' + os.sep + 'openfast' + os.sep + 'IEA_LB_RWT-AeroAcoustics'
+AAfilename = FAST_directory + os.sep + 'IEA_LB_RWT-AeroAcoustics_3.out'
+OFfilename = FAST_directory + os.sep + 'IEA_LB_RWT-AeroAcoustics.out'
+output_dir = os.path.dirname( os.path.realpath(__file__) )
+outputfilename = output_dir + os.sep + "data_output2"
+
+#########################################################################################################################################
+## Read in data, manipulate it, and plot it
+
+# Read in file data
+AA_3 = FASTOutputFile(AAfilename).toDataFrame()
+OF = FASTOutputFile(OFfilename).toDataFrame()
+
+# Determine number of observers
+num_obs = int((AA_3.shape[1]-1)/(7*34))
+
+# Calculate sample time for n revolutions
+rpm = OF[["RotSpeed_[rpm]"]].mean()[0]
+time_revs = n*60/rpm
+tot_time = AA_3["Time_[s]"].max()
+if time_revs < tot_time:
+ sample_time = tot_time - time_revs
+else:
+ print("Error: Time for number of revolutions exceeds simulation time. Reduce n.")
+ raise SystemExit('')
+
+# Slice AA dataframe for t > sample_time
+AA_3 = AA_3[AA_3["Time_[s]"] > sample_time]
+AA_3=AA_3.drop("Time_[s]",axis=1)
+
+# Average SPL for each observer
+AA_3 = AA_3.mean()
+
+# Manipulate PD dataframes
+# convert to dataframe with appropriate columns
+cols = ['Observer','Mechanism','Frequency (Hz)','SPL (dB)']
+aa_3 = pd.DataFrame(columns=cols)
+for i in AA_3.index:
+ nums = re.findall(r"[-+]?\d*\.\d+|\d+",i)
+ aa_3.loc[len(aa_3)] = [nums[0],nums[2],nums[1],AA_3[i]]
+
+AA_3 = aa_3
+
+# rename mechanism for legend
+for i in range(0,AA_3.last_valid_index()+1):
+ if AA_3.loc[i,"Mechanism"]=='1':
+ AA_3.loc[i,"Mechanism"]="LBL-VS"
+ if AA_3.loc[i,"Mechanism"]=='2':
+ AA_3.loc[i,"Mechanism"]="TBL-TE-PS"
+ if AA_3.loc[i,"Mechanism"]=='3':
+ AA_3.loc[i,"Mechanism"]="TBL-TE-SS"
+ if AA_3.loc[i,"Mechanism"]=='4':
+ AA_3.loc[i,"Mechanism"]="TBL-TE-AoA"
+ if AA_3.loc[i,"Mechanism"]=='5':
+ AA_3.loc[i,"Mechanism"]="TE Bluntness"
+ if AA_3.loc[i,"Mechanism"]=='6':
+ AA_3.loc[i,"Mechanism"]="Tip Vortex"
+ if AA_3.loc[i,"Mechanism"]=='7':
+ AA_3.loc[i,"Mechanism"]="TI"
+
+AA_3["Observer"]=AA_3["Observer"].apply(pd.to_numeric)
+AA_3["Frequency (Hz)"]=AA_3["Frequency (Hz)"].apply(pd.to_numeric)
+AA_3["SPL (dB)"]=AA_3["SPL (dB)"].apply(pd.to_numeric)
+
+
+# Plot spectra
+fs = 10
+for j in range(num_obs):
+ fig,ax=plt.subplots()
+ plt.xscale('log')
+ ax.set_xlabel('Frequency (Hz)', fontsize=fs+2, fontweight='bold')
+ ax.set_ylabel('SPL (dB)', fontsize=fs+2, fontweight='bold')
+ for i in range(7):
+ plt.plot(AA_3["Frequency (Hz)"][j*34*7 + i : j*34*7 + i + 34 * 7:7], AA_3["SPL (dB)"][j*34*7 + i : j*34*7 + i + 34 * 7:7], label = AA_3.loc[i,"Mechanism"])
+ ax.set_title('Observer ' + str(j), fontsize=fs+2, fontweight='bold')
+ plt.grid(color=[0.8,0.8,0.8], linestyle='--')
+ ax.set_ylim(0,)
+ ax.legend()
+ if save_fig == True:
+ fig_name = 'spectra_Obs' + str(j) + fig_ext
+ fig.savefig(output_dir + os.sep + fig_name)
+ plt.show()
+
+# Export to csv
+if save_data == True:
+ AA_3.to_csv(r'{}-data.csv'.format(outputfilename))
+
+
+
+

openfast_toolbox/aeroacoustics/examples/_plot_rotor_map.py

@@ -0,0 +1,92 @@
+import numpy as np
+import pandas as pd
+import matplotlib.pyplot as plt
+from parse import *
+import re, os, platform
+import matplotlib.colors
+from openfast_toolbox.io.fast_output_file import FASTOutputFile
+
+#########################################################################################################################################
+## User inputs
+# Save plot and/or data?
+save_fig = False
+save_data = False
+fig_ext = '.png'
+R = 65.
+R_hub = 2.
+
+# Number of revolutions (n) to average spectra
+n = 1
+
+#########################################################################################################################################
+## Paths to files
+if platform.system() == 'Windows':
+ FAST_directory = os.path.dirname( os.path.dirname( os.path.dirname( os.path.dirname( os.path.realpath(__file__) ) ) ) ) + os.sep + 'reg_tests' + os.sep + 'r-tests' + os.sep + 'glue-codes' + os.sep + 'openfast' + os.sep + 'IEA_LB_RWT-AeroAcoustics'
+else:
+ FAST_directory = os.path.dirname( os.path.dirname( os.path.dirname( os.path.dirname( os.path.realpath(__file__) ) ) ) ) + os.sep + 'openfast' + os.sep + 'reg_tests' + os.sep + 'r-tests' + os.sep + 'glue-codes' + os.sep + 'openfast' + os.sep + 'IEA_LB_RWT-AeroAcoustics'
+AAfilename = FAST_directory + os.sep + 'IEA_LB_RWT-AeroAcoustics_4.out'
+OFfilename = FAST_directory + os.sep + 'IEA_LB_RWT-AeroAcoustics.out'
+locfilename = FAST_directory + os.sep + 'AA_ObserverLocations.dat'
+output_dir = os.path.dirname( os.path.realpath(__file__) )
+outputfilename = output_dir + os.sep + "data_output4"
+
+#########################################################################################################################################
+
+
+location = pd.read_csv(locfilename,delimiter='\s+',skiprows=[0,1],names=['x','y','z'])
+
+
+AA_1 = FASTOutputFile(AAfilename).toDataFrame()
+OF = FASTOutputFile(OFfilename).toDataFrame()
+
+
+with open(AAfilename, 'r') as f:
+ f.readline()
+ f.readline()
+ f.readline()
+ n_obs = int(f.readline().split()[-1])
+ n_blades = int(f.readline().split()[-1])
+ n_nodes = int(f.readline().split()[-1])
+f.close()
+
+k = np.ones(n_obs)
+for i in range(n_obs):
+ if location['x'][i] < 0:
+ k[i] = -1
+
+phi = OF['Azimuth_[deg]'] / 180. * np.pi
+phi_interp = np.interp(AA_1['Time_[s]'], OF['Time_[s]'], phi)
+index = []
+for i in range(1, len(phi_interp)):
+ if phi_interp[i] < phi_interp[i-1]:
+ index.append(i)
+
+y_b = np.linspace(R_hub, R, n_nodes)
+x_b = np.zeros_like(y_b)
+
+n_pts = index[-1] - index[-2]
+
+for j in range(n_obs):
+ x = np.zeros((n_pts,n_nodes))
+ y = np.zeros((n_pts,n_nodes))
+
+ for i in range(n_pts):
+ x[i,:] = x_b * np.cos(k[j]*phi_interp[i + index[-2]]) - y_b * np.sin(k[j]*phi_interp[i + index[-2]])
+ y[i,:] = x_b * np.sin(k[j]*phi_interp[i + index[-2]]) + y_b * np.cos(k[j]*phi_interp[i + index[-2]])
+
+ z = np.array(AA_1)[index[-2]:index[-1], 1+j:1 + 30*n_obs + j:n_obs]
+ fs = 10
+ fig,ax=plt.subplots()
+ ax.set_aspect('equal')
+ ax.set_xlabel('y [m]', fontsize=fs+2, fontweight='bold')
+ ax.set_ylabel('z [m]', fontsize=fs+2, fontweight='bold')
+ ax.set_title('Observer ' + str(j), fontsize=fs+2, fontweight='bold')
+ tcf=ax.tricontourf(x.flatten(),y.flatten(),z.flatten(), range(20,75))
+ fig.colorbar(tcf,orientation="vertical").set_label(label = 'Overall SPL [dB]', fontsize=fs+2,weight='bold')
+ if save_fig == True:
+ fig_name = 'rotor_map_Obs' + str(j) + fig_ext
+ fig.savefig(output_dir + os.sep + fig_name)
+ plt.show()
+
+
+

openfast_toolbox/aeroacoustics/examples/_plot_spectra.py

@@ -0,0 +1,87 @@
+import numpy as np
+import pandas as pd
+import matplotlib.pyplot as plt
+from parse import *
+import re, os, platform
+from openfast_toolbox.io.fast_output_file import FASTOutputFile
+
+#########################################################################################################################################
+## User inputs
+# Save plot and/or data?
+save_fig = False
+save_data = False
+fig_ext = '.png'
+
+# Number of revolutions (n) to average spectra
+n = 1
+
+#########################################################################################################################################
+## Paths to files
+if platform.system() == 'Windows':
+ FAST_directory = os.path.dirname( os.path.dirname( os.path.dirname( os.path.dirname( os.path.realpath(__file__) ) ) ) ) + os.sep + 'reg_tests' + os.sep + 'r-tests' + os.sep + 'glue-codes' + os.sep + 'openfast' + os.sep + 'IEA_LB_RWT-AeroAcoustics'
+else:
+ FAST_directory = os.path.dirname( os.path.dirname( os.path.dirname( os.path.dirname( os.path.realpath(__file__) ) ) ) ) + os.sep + 'openfast' + os.sep + 'reg_tests' + os.sep + 'r-tests' + os.sep + 'glue-codes' + os.sep + 'openfast' + os.sep + 'IEA_LB_RWT-AeroAcoustics'
+AAfilename = FAST_directory + os.sep + 'IEA_LB_RWT-AeroAcoustics_2.out'
+OFfilename = FAST_directory + os.sep + 'IEA_LB_RWT-AeroAcoustics.out'
+output_dir = os.path.dirname( os.path.realpath(__file__) )
+outputfilename = output_dir + os.sep + "data_output2"
+
+#########################################################################################################################################
+## Read in data, manipulate it, and plot it
+
+# Read in file data
+AA_2 = FASTOutputFile(AAfilename).toDataFrame()
+OF = FASTOutputFile(OFfilename).toDataFrame()
+
+# Determine number of observers
+num_obs = int((AA_2.shape[1]-1)/34)
+
+# Calculate sample time for n revolutions
+rpm = OF[["RotSpeed_[rpm]"]].mean()[0]
+time_revs = n*60/rpm
+tot_time = AA_2["Time_[s]"].max()
+if time_revs < tot_time:
+ sample_time = tot_time - time_revs
+else:
+ print("Error: Time for number of revolutions exceeds simulation time. Reduce n.")
+ raise SystemExit('')
+
+# Slice AA dataframe for t > sample_time
+AA_2 = AA_2[AA_2["Time_[s]"] > sample_time]
+AA_2=AA_2.drop("Time_[s]",axis=1)
+
+# Average SPL for each observer
+AA_2 = AA_2.mean()
+
+# Manipulate PD dataframes
+cols = ['Observer','Frequency (Hz)','SPL (dB)'] 
+aa_2 = pd.DataFrame(columns=cols)
+for i in AA_2.index:
+ nums = re.findall(r"[-+]?\d*\.\d+|\d+",i)
+ aa_2.loc[len(aa_2)] = [nums[0],nums[1],AA_2[i]]
+AA_2 = aa_2
+AA_2["Frequency (Hz)"]=AA_2["Frequency (Hz)"].apply(pd.to_numeric)
+AA_2["SPL (dB)"]=AA_2["SPL (dB)"].apply(pd.to_numeric)
+
+# Plot spectra
+fs = 10
+fig,ax=plt.subplots()
+plt.xscale('log')
+ax.set_xlabel('Frequency (Hz)', fontsize=fs+2, fontweight='bold')
+ax.set_ylabel('SPL (dB)', fontsize=fs+2, fontweight='bold')
+for i in range(num_obs):
+ plt.plot(AA_2["Frequency (Hz)"][i*34:i*34 + 34], AA_2["SPL (dB)"][i*34:i*34 + 34], label = 'Observer ' + str(i))
+plt.grid(color=[0.8,0.8,0.8], linestyle='--')
+ax.legend()
+if save_fig == True:
+ fig_name = 'spectra' + fig_ext
+ fig.savefig(output_dir + os.sep + fig_name)
+plt.show()
+
+# Export to csv
+if save_data == True:
+ AA_2.to_csv(r'{}-data.csv'.format(outputfilename))
+
+
+
+