notebook.community

Edit and run 





In [1]:



    


%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import h5py
from importlib import reload
import PIVutils
import PODutils
import modred as mr

saveFolder = '/Users/Owen/Dropbox/Python Codes/ASIIT/Data/'
saveFile = 'RNV45-thumbs.hdf5'
imgFolder = saveFolder[:-5] + 'RNV45-PODthumbs'

import os
if not os.path.exists(imgFolder):
    os.makedirs(imgFolder)



    











In [2]:



    


import sys
sys.executable



    


















    
Out[2]:








'/Users/Owen/anaconda/bin/python'

















In [3]:



    


#PIVutils = reload(PIVutils)
#X, Y, U, V, Swirl = PIVutils.loadDataset('/Users/Owen/Dropbox/Data/ABL/Heat Flux Data/Processed Results/N/Neutral45_2.mat',['X','Y','U','V','Swirl'],[])
#X, Y, U, V, Swirl, Cond, Prof = PIVutils.loadDataset('/Users/Owen/Dropbox/Data/ABL/Heat Flux Data/Processed Results/N/Neutral45.mat',['X','Y','U','V','Swirl'],['Cond','Prof'])
#X, Y, U, V, Swirl, Cond, Prof = PIVutils.loadDataset('/Users/Owen/Dropbox/Data/ABL/SBL PIV data/RNV45-RI2.mat',\
#                                                     ['X','Y','U','V','Swirl'],['Cond','Prof'])



    











In [57]:



    


X,Y,x,y,U,V,S,missVecs_Pro,Yind_Pro,Ypos_Pro,Ur,Vr,Sr,missVecs_Ret,Yind_Ret,Ypos_Ret,Yvec,Cond,Prof,ThumbParams = \
                        PIVutils.loadDataset(saveFolder + saveFile,
                         ['X','Y','x','y','U','V','S','missVecs','Yind_Pro','Ypos_Pro','Ur','Vr','Sr','missVecs_Ret','Yind_Ret','Ypos_Ret','Yvec'], 
                         ['Cond','Prof','ThumbParams'])
BoxSize = ThumbParams['BoxSize']
uSize = U.shape



    


















    






['Cond', 'Prof', 'S', 'Sr', 'ThumbParams', 'U', 'Ur', 'V', 'Vr', 'X', 'Y', 'Yind_Pro', 'Yind_Ret', 'Ypos_Pro', 'Ypos_Ret', 'Yvec', 'missVecs', 'missVecs_Ret', 'x', 'y']

Subtract mean profile





In [5]:



    


#Calculate mean velocity profiles
Umean = Prof['U']

#Now get the thumbnails
thumb = 0
for i in range(uSize[2]):
    ind = Yind_Pro[i]
    #U[:,:,thumb] = U[:,:,thumb]-np.transpose(np.tile(Umean[ind-BoxSize:ind+BoxSize+1],(1,uSize[0])),(1,0))
    U[:,:,thumb] = U[:,:,thumb]-np.tile(Umean[ind-BoxSize:ind+BoxSize+1],(1,uSize[0]))#Now get the thumbnails
    thumb+=1
#del Umean



    


















    






/Users/Owen/anaconda/lib/python3.5/site-packages/ipykernel/__main__.py:10: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future

Now do POD on the resulting thumbnails





In [6]:



    


#Set all NaNs to zero
U[np.isnan(U)] = 0
V[np.isnan(V)] = 0

uSize = U.shape



    











In [7]:



    


uAll = np.concatenate((U.reshape(uSize[0]*uSize[1],uSize[2]), V.reshape(uSize[0]*uSize[1],uSize[2])), axis = 0)
#uAll.shape



    











In [8]:



    


num_modes = 50;
modes, eig_vals = mr.compute_POD_matrices_snaps_method(uAll, list(range(num_modes)))



    











In [9]:



    


menergy = eig_vals/np.sum(eig_vals)
menergy_sum = np.zeros(len(menergy))
for i in range(len(menergy)):
    menergy_sum[i] = np.sum(menergy[:i]);



    











In [10]:



    


reload(PODutils)
Umodes, Vmodes= PODutils.reconstructPODmodes(modes,uSize,num_modes,2)



    











In [11]:



    


#Calculate the mode coefficients
C = modes.transpose()*uAll
C = np.asarray(C)

Plot modal energy and contribution total





In [12]:



    


ind = np.arange(num_modes)  # the x locations for the groups
width = 1       # the width of the bars

f = plt.figure()
ax = plt.gca()
ax2 = plt.twinx()
rect = ax.bar(ind,menergy[:num_modes], width, color='gray')
line = ax2.plot(ind,menergy_sum[:num_modes],'--r')

ax.set_xlabel("Mode Number",fontsize=14)
ax.set_ylabel("Scaled Mode Energy",fontsize=14)
ax2.set_ylabel("Integrated Energy",fontsize=14,color='red')

f.savefig(saveFolder + '/ModeEnergy-Pro-10.pdf', transparent=True, bbox_inches='tight', pad_inches=0)

Plot some modes





In [13]:



    


import PODutils
PODutils.plotPODmodes2D(X[:2*BoxSize+1,:2*BoxSize+1],Y[:2*BoxSize+1,:2*BoxSize+1],Umodes,Vmodes,list(range(25)),saveFolder)



    


















    






/Users/Owen/anaconda/lib/python3.5/site-packages/ipykernel/__main__.py:2: VisibleDeprecationWarning: converting an array with ndim > 0 to an index will result in an error in the future
  from ipykernel import kernelapp as app
/Users/Owen/anaconda/lib/python3.5/site-packages/matplotlib/pyplot.py:524: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).
  max_open_warning, RuntimeWarning)

Plot the variation of the coefficients





In [14]:



    


reload(PODutils)
PODutils.plotPODcoeff(C,list(range(6)),100,8,1,saveFolder)



    


















    
























In [15]:



    


#ind = 0
#ypos =[]
#for i in range(len(cent_Pro)):
#    for j in range(len(cent_Pro[i])):
#        ypos = np.concatenate([ypos, [Y[cent_Pro[i][j][0],1]]])
        
        #[ind] = cent_Pro[i][j][0]*scale+Y[1,1]
        #ind+=1



    











In [16]:



    


# Calcualte pdf of wall-normal locations
xedges = Yvec
xedges = np.concatenate([[xedges[0]-(xedges[1]-xedges[0])],xedges, [xedges[-1]+xedges[-1]-xedges[-2]]])
p, x = np.histogram(Ypos_Pro, xedges,normed=True)
x = x[:-1] + (x[1] - x[0])/2 
plt.plot(x, p)
plt.xlabel('$y/\delta$', fontsize=15)
plt.ylabel('pdf of thumbnail locations')

plt.savefig(saveFolder + '/PDFswirl-Pro-10.pdf', transparent=True, bbox_inches='tight', pad_inches=0)



    


















    
























In [22]:



    


reload(PODutils)
C2 = C.copy()
C2 = np.concatenate([Ypos_Pro.reshape((-1, 1)).T, C2], axis=0)
PODutils.plotYposPODcoeff(Yvec,C2,list(range(10)),50,10,1,saveFolder)



    


















    
























In [23]:



    


del C,C2,modes,Umodes,Vmodes

Do POD on retrograde swirl blobs (Box = 21x21)





In [58]:



    


#Set all NaNs to zero
Ur[np.isnan(Ur)] = 0
Vr[np.isnan(Vr)] = 0

uSize = Ur.shape
print(uSize)



    


















    






(21, 21, 3489)

















In [59]:



    


thumb = 0

for i in range(uSize[2]):
    ind = Yind_Ret[i]
    #U[:,:,thumb] = U[:,:,thumb]-np.transpose(np.tile(Umean[ind-BoxSize:ind+BoxSize+1],(1,uSize[0])),(1,0))
    Ur[:,:,thumb] = Ur[:,:,thumb]-np.tile(Umean[ind-BoxSize:ind+BoxSize+1],(1,uSize[0]))
    thumb+=1
#del Umean



    


















    






/Users/Owen/anaconda/lib/python3.5/site-packages/ipykernel/__main__.py:6: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future

















In [60]:



    


uAll = np.concatenate((Ur.reshape(uSize[0]*uSize[1],uSize[2]), Vr.reshape(uSize[0]*uSize[1],uSize[2])), axis = 0)
#uAll.shape



    











In [61]:



    


import modred as mr

num_modes = 50;
modes, eig_vals = mr.compute_POD_matrices_snaps_method(uAll, list(range(num_modes)))



    











In [62]:



    


menergy = eig_vals/np.sum(eig_vals)
menergy_sum = np.zeros(len(menergy))
for i in range(len(menergy)):
    menergy_sum[i] = np.sum(menergy[:i]);



    











In [63]:



    


reload(PODutils)
Umodes, Vmodes= PODutils.reconstructPODmodes(modes,uSize,num_modes,2)



    











In [64]:



    


#Calculate the mode coefficients
C = modes.transpose()*uAll
C = np.asarray(C)

Plot modal energy and contribution total





In [65]:



    


ind = np.arange(num_modes)  # the x locations for the groups
width = 1       # the width of the bars

f = plt.figure()
ax = plt.gca()
ax2 = plt.twinx()
rect = ax.bar(ind,menergy[:num_modes], width, color='gray')
line = ax2.plot(ind,menergy_sum[:num_modes],'--r')

ax.set_xlabel("Mode Number",fontsize=14)
ax.set_ylabel("Scaled Mode Energy",fontsize=14)
ax2.set_ylabel("Integrated Energy",fontsize=14,color='red')

f.savefig(saveFolder + '/ModeEnergy-Ret-10.pdf', transparent=True, bbox_inches='tight', pad_inches=0)

Plot some modes





In [66]:



    


import PODutils
PODutils.plotPODmodes2D(X[:2*BoxSize+1,:2*BoxSize+1],Y[:2*BoxSize+1,:2*BoxSize+1],Umodes,Vmodes,list(range(25)))



    


















    






/Users/Owen/anaconda/lib/python3.5/site-packages/ipykernel/__main__.py:2: VisibleDeprecationWarning: converting an array with ndim > 0 to an index will result in an error in the future
  from ipykernel import kernelapp as app
/Users/Owen/anaconda/lib/python3.5/site-packages/matplotlib/pyplot.py:524: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).
  max_open_warning, RuntimeWarning)

Plot the variation of the coefficients





In [67]:



    


reload(PODutils)
PODutils.plotPODcoeff(C,list(range(6)),75,8,1)



    


















    
























In [68]:



    


# Calcualte pdf of wall-normal locations
xedges = Yvec
xedges = np.concatenate([[xedges[0]-(xedges[1]-xedges[0])],xedges, [xedges[-1]+xedges[-1]-xedges[-2]]])
p, x = np.histogram(Ypos_Ret, xedges,normed=True)
x = x[:-1] + (x[1] - x[0])/2 
plt.plot(x, p)
plt.xlabel('$y/\delta$', fontsize=15)
plt.ylabel('pdf of thumbnail locations')

plt.savefig(saveFolder + '/PDFswirl-Ret-10.pdf', transparent=True, bbox_inches='tight', pad_inches=0)



    


















    
























In [69]:



    


reload(PODutils)
C2 = C.copy()
C2 = np.concatenate([Ypos_Ret.reshape((-1, 1)).T, C2], axis=0)
PODutils.plotYposPODcoeff(Yvec,C2,list(range(10)),50,10,1)

Do POD on random set of data





In [24]:



    


reload(PIVutils)
BoxSize = 10

Ut, Vt, St = PIVutils.getRandomThumbnails2D(U,Vf,Swirl,num_features_Pro,BoxSize)



    


















    






---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-24-25784f891fef> in <module>()
      2 BoxSize = 10
      3 
----> 4 Ut, Vt, St = PIVutils.getRandomThumbnails2D(Uf,Vf,Swirl,num_features_Pro,BoxSize)

NameError: name 'Uf' is not defined
















In [ ]:



    


#Set all NaNs to zero
Ut[np.isnan(Ut)] = 0
Vt[np.isnan(Vt)] = 0

uSize = Ut.shape



    











In [ ]:



    


uAll = np.concatenate((Ut.reshape(uSize[0]*uSize[1],uSize[2]), Vt.reshape(uSize[0]*uSize[1],uSize[2])), axis = 0)
#uAll.shape



    











In [ ]:



    


import modred as mr

num_modes = 50;
modes, eig_vals = mr.compute_POD_matrices_snaps_method(uAll, list(range(num_modes)))



    











In [ ]:



    


menergy = eig_vals/np.sum(eig_vals)
menergy_sum = np.zeros(len(menergy))
for i in range(len(menergy)):
    menergy_sum[i] = np.sum(menergy[:i]);



    











In [ ]:



    


reload(PODutils)
Umodes, Vmodes= PODutils.reconstructPODmodes(modes,uSize,num_modes,2)



    











In [ ]:



    


#Calculate the mode coefficients
C = modes.transpose()*uAll
C = np.asarray(C)

Plot modal energy and contribution total





In [ ]:



    


ind = np.arange(num_modes)  # the x locations for the groups
width = 1       # the width of the bars

f = plt.figure()
ax = plt.gca()
ax2 = plt.twinx()
rect = ax.bar(ind,menergy[:num_modes], width, color='gray')
line = ax2.plot(ind,menergy_sum[:num_modes],'--r')

ax.set_xlabel("Mode Number",fontsize=14)
ax.set_ylabel("Scaled Mode Energy",fontsize=14)
ax2.set_ylabel("Integrated Energy",fontsize=14,color='red')

f.savefig(saveFolder + '/ModeEnergy-Pro-10.pdf', transparent=True, bbox_inches='tight', pad_inches=0)

Plot some modes





In [ ]:



    


import PODutils
PODutils.plotPODmodes2D(X[:2*BoxSize+1,:2*BoxSize+1],Y[:2*BoxSize+1,:2*BoxSize+1],Umodes,Vmodes,list(range(25)),saveFolder)

Plot the variation of the coefficients





In [ ]:



    


reload(PODutils)
PODutils.plotPODcoeff(C,list(range(6)),100,8,1,saveFolder)



    











In [ ]:



    


ind = 0
ypos =[]
for i in range(len(cent_Pro)):
    for j in range(len(cent_Pro[i])):
        ypos = np.concatenate([ypos, [Y[cent_Pro[i][j][0],1]]])
        
        #[ind] = cent_Pro[i][j][0]*scale+Y[1,1]
        #ind+=1



    











In [ ]:



    


# Calcualte pdf of wall-normal locations
xedges = Y[:,0]/Cond["delta"]
xedges = np.concatenate([[xedges[0]-(xedges[1]-xedges[0])],xedges, [xedges[-1]+xedges[-1]-xedges[-2]]])
p, x = np.histogram(ypos/Cond["delta"], xedges,normed=True)
x = x[:-1] + (x[1] - x[0])/2 
plt.plot(x, p)
plt.xlabel('$y/\delta$', fontsize=15)
plt.ylabel('pdf of thumbnail locations')

plt.savefig(saveFolder + '/PDFswirl-Pro-10.pdf', transparent=True, bbox_inches='tight', pad_inches=0)



    











In [ ]:



    


reload(PODutils)
C2 = C.copy()
C2 = np.concatenate([ypos.reshape((-1, 1)).T/Cond["delta"], C2], axis=0)
PODutils.plotYposPODcoeff(Y[:,0]/Cond["delta"],C2,list(range(10)),50,10,1,saveFolder)



    











In [ ]:



    






    











In [ ]:

Content source: owenjhwilliams/ASIIT

Similar notebooks:

notebook.community | gallery | about