diff --git a/README.md b/README.md index 8e42ea9..4c3f19f 100644 --- a/README.md +++ b/README.md @@ -13,4 +13,4 @@ A guide to simulation and post-processing of the turbulent flow over a periodic A tutorial for running a minimal channel simulation in `Nek5000`, and doing a postprocessing using proper orthogonal decomposition (POD). ### `/pod/` -A python package for POD using Nek5000 data. +A python package for POD and DMD using Nek5000 data. diff --git a/pod/channel/MODULES/DMDmodule.py b/pod/channel/MODULES/DMDmodule.py new file mode 100644 index 0000000..f3eeee6 --- /dev/null +++ b/pod/channel/MODULES/DMDmodule.py @@ -0,0 +1,67 @@ +""" +The module computes DMD modes via +singular value decomposition +""" + +import numpy as np +import numpy.linalg as la + +def hermitian(arr): + return arr.conj().T + +def DMD(Usnp,mvect,nsnap,r,ifsym,deltaT): + """ + Module to compute DMD modes + Args: + - Usnp = snapshots matrix Usnp + - mvect = array with square root of mass weights (nGLLe*ncomponents) + - nsnap = number of snapshots + - ifsym = to mirror data wrt symmetry x axis + - r = reduced rank number + - deltaT = timestep + Returns: + - Phi = DMD Modes + - Lambdat = DMD Eigenvalues + - a1 = DMD Coefficients + + """ + + # Snapshot pairs + X1 = Usnp[:,0:-1] + X2 = Usnp[:,1:] + + # SVD of X1: + # - U = matrix whose columns are the spatial modes in POD expansion + # - Sigma = diagonal matrix of singular values for Usnp + # - Vt = conjugate transpose of matrix with temporal coefficients + U,Sigma,Vt = la.svd(X1,full_matrices=False) + V = hermitian(Vt) + + # Rank r Reduction + U = U[:,:r] + Sigma = np.diag(Sigma[:r]) + V = V[:,:r] + + # Computation of Atilde + Atilde = hermitian(U) @ X2 @ V @ la.inv(Sigma) + + # Computation of Eigenvalues and Eigenmodes + Lambdat, W = la.eig(Atilde) + Phi = X2 @ V @ la.inv(Sigma) @ W + omega = np.log(Lambdat)/deltaT ### continuous time eigenvalues + + # Normalize the DMD Modes + Phi = Phi / np.linalg.norm(Phi, axis=0) + + # DMD Coefficients + a1 = la.lstsq(Phi, Usnp[:,0], rcond=None)[0] + + # Reordering eigenvalues and modes based on mode energies + + mode_energy = np.abs(a1)**2 + idx = np.flip(np.argsort(mode_energy)) + Lambdat = Lambdat[idx] + omega = omega[idx] + Phi = Phi[:, idx] + + return Phi, Lambdat, a1, omega \ No newline at end of file diff --git a/pod/channel/MODULES/PODmodule.py b/pod/channel/MODULES/PODmodule.py index 27c5b69..57a2638 100644 --- a/pod/channel/MODULES/PODmodule.py +++ b/pod/channel/MODULES/PODmodule.py @@ -15,7 +15,7 @@ def POD(Usnp,mvect,nsnap,ifsym): Module to compute POD modes Args: - Usnp = snapshots matrix Usnp - - mvect = array with square of mass weights (nGLLe*ncomponents) + - mvect = array with square root of mass weights (nGLLe*ncomponents) - nsnap = number of snapshots - ifsym = to mirror data wrt symmetry x axis Returns: diff --git a/pod/channel/MODULES/dbMaker.py b/pod/channel/MODULES/dbMaker.py index ab65627..ae9c90f 100644 --- a/pod/channel/MODULES/dbMaker.py +++ b/pod/channel/MODULES/dbMaker.py @@ -8,20 +8,20 @@ -def dbMan_sym(info,iff,infom,iffm,infos,iffs): +def dbMan_sym(info,iff,if3D,infom,iffm,infos,iffs): """ To manage the databases creation. """ - db = dbCreator(info,iff) - db_m = dbCreator(infom,iffm) + db = dbCreator(info,iff,if3D) + db_m = dbCreator(infom,iffm,if3D) # Consistency check mass matrix/snapshots mesh data_ms = db_m['data'][0] if (db['data'][0].nel != data_ms.nel): print('N.elements of mass matrix element != N.elements of velocity fields !!!') exit(0) # Symmetric fields database - db_s = dbCreator(infos,iffs) + db_s = dbCreator(infos,iffs,if3D) if (db['data'][0].nel != db_s['data'][0].nel): print('N.elements of mirrored field != N.elements of original fields !!!') exit(0) @@ -31,13 +31,13 @@ def dbMan_sym(info,iff,infom,iffm,infos,iffs): -def dbMan(info,iff,infom,iffm): +def dbMan(info,iff,if3D,infom,iffm): """ To manage the databases creation. """ - db = dbCreator(info,iff) - db_m = dbCreator(infom,iffm) + db = dbCreator(info,iff,if3D) + db_m = dbCreator(infom,iffm,if3D) # Consistency check mass matrix/snapshots mesh data_ms = db_m['data'][0] if (db['data'][0].nel != data_ms.nel): @@ -51,12 +51,12 @@ def dbMan(info,iff,infom,iffm): -def dbCreator(info,iff): +def dbCreator(info,iff,if3D): """ The database is created for the solution field or the mass matrix. """ if (iff=='field'): - db = dbCreator_v(info) + db = dbCreator_v(info, if3D) elif (iff=='mass'): db = dbCreator_m(info) @@ -66,7 +66,7 @@ def dbCreator(info,iff): -def dbCreator_v(info): +def dbCreator_v(info, if3D): """ The database is created for the solution field (all the snapshots). """ @@ -74,13 +74,35 @@ def dbCreator_v(info): caseName_= info['caseName'] start_ = info['startID'] end_ = info['endID'] - + pre_=caseName_+'0.f' nSnap=0 ll =[] + + if if3D: + vars = ['ux', 'uy', 'uz', 'pressure', 'temperature'] + else: + vars = ['ux', 'uy', 'pressure', 'temperature'] + + if (info['qoiName'] == 'temperature'): + vars.remove('temperature') + elif (info['qoiName'] == 'pressure'): + vars.remove('pressure') + elif (info['qoiName'] == 'uvel'): + vars.remove('uvel') + elif (info['qoiName'] == 'vvel'): + vars.remove('vvel') + elif (info['qoiName'] == 'wvel'): + vars.remove('wvel') + elif (info['qoiName'] == 'velocity'): + vars.remove('uvel') + vars.remove('vvel') + if if3D: + vars.remove('wvel') + for id_ in range(start_,end_+1): nSnap+=1 - data = neksuite.readnek(path_+pre_+str(id_).zfill(5),'float32') + data = neksuite.readnek(path_+pre_+str(id_).zfill(5),'float32', skip_vars=vars) ll.append(data) time_=datetime.now() diff --git a/pod/channel/MODULES/freqAnalysis.py b/pod/channel/MODULES/freqAnalysis.py new file mode 100644 index 0000000..af4cb21 --- /dev/null +++ b/pod/channel/MODULES/freqAnalysis.py @@ -0,0 +1,79 @@ +import numpy as np +import matplotlib.pyplot as plt + +plt.rcParams.update({ + 'font.size': 10, + 'font.family': 'serif', + 'font.serif': 'Times New Roman', + 'mathtext.fontset': 'cm', + 'axes.labelsize': 12, + 'axes.titlesize': 12, + 'axes.grid': True, + 'axes.edgecolor': 'black', + 'axes.linewidth': 0.8, + 'grid.linestyle': '--', + 'grid.alpha': 0.7, + 'grid.color': 'gray', + 'xtick.labelsize': 10, + 'ytick.labelsize': 10, + 'xtick.direction': 'in', + 'ytick.direction': 'in', + 'xtick.major.size': 5, + 'ytick.major.size': 5, + 'lines.linewidth': 1.5, + 'lines.markersize': 5, + 'legend.fontsize': 10, + 'legend.frameon': False, + 'legend.loc': 'best', + 'savefig.dpi': 300, + 'savefig.format': 'png', +}) + + +def freqPOD(A2,maxModes,nplt,timestep,info): + + #A2 -- POD Coefficients (m,n) where m is the number of snapshots and n is the mode + #Mode 0 is the time average. + + fileName = info['outputPath']+'PODFreq.txt' + F = open(fileName,'w') + F.write("# Snapshots POD result \n") + F.write("# Dominant frequencies calculated using FFT \n") + F.write("# According to our convention the mode 0 is the mean value \n") + F.write("# ------------------------------------------------------------------\n") + F.write("# Mode\t F1\t F2\t F3\t F4\t F5\n") + F.write("# ------------------------------------------------------------------\n") + + # Define frequencies for FFT + frequencies = np.fft.rfftfreq(nplt, d=timestep) # Generates only positive frequencies + print(frequencies) + + for i in range(1,maxModes): + # Step 1: Compute FFT + fft_result = np.fft.fft(A2[:, i] - np.mean(A2[:, i])) + power = np.abs(fft_result[:nplt // 2 + 1]) # Compute power spectrum (positive frequencies only) + + # Step 2: Filter the power spectrum + mean_power = np.mean(power) + std_power = np.std(power) + threshold = mean_power + 2 * std_power + filtered_power = np.where(power > threshold, power, 0) + + # Step 3: Identify dominant frequencies + idx = np.argsort(filtered_power) # Sort indices based on filtered power + f = frequencies[np.flip(idx)[0:5]] # Top 5 frequencies + F.write("%g\t%g\t%g\t%g\t%g\t%g \n" % \ + (i, f[0],f[1],f[2],f[3],f[4])) + + plt.figure(figsize=(4,3)) + plt.plot(frequencies, power, label='FFT Power', alpha=0.6) + plt.plot(frequencies, filtered_power, label='Filtered Power', linestyle="solid") + plt.xlabel('Frequency (Hz)') + plt.ylabel('Power') + plt.grid(True) + plt.legend() + plt.tight_layout() + plt.savefig(info['outputPath']+f"Mode{i}.png") + plt.close() + + F.close() diff --git a/pod/channel/MODULES/outpWriter.py b/pod/channel/MODULES/outpWriter.py index b2d7f00..1b9877d 100644 --- a/pod/channel/MODULES/outpWriter.py +++ b/pod/channel/MODULES/outpWriter.py @@ -23,7 +23,8 @@ def modes(nmod,db,info,L,if3D): - if3D = to identify 3D or 2D case """ - + print(np.shape(L)) + for m in range(nmod+1): data = db @@ -80,17 +81,24 @@ def modes(nmod,db,info,L,if3D): neksuite.writenek(info['outputPath']+'PODmod'+info['caseName']+'0.f'+str(m).zfill(5),data_o) print('POD mode number %d has been saved' % (m)) - elif (info['module']=='DMD'): + print('Writing: '+info['outputPath']+'PODmod'+info['caseName']+'.nek5000') + elif (info['module']=='DMD'): # Writing data in nek/visit format + data_o.time = m + neksuite.writenek(info['outputPath']+'DMDmod'+info['caseName']+'0.f'+str(m).zfill(5),data_o) print('DMD mode number %d has been saved' % (m)) - - print('Writing: '+info['outputPath']+'PODmod'+info['caseName']+'.nek5000') - with open(info['outputPath']+'PODmod'+info['caseName']+'.nek5000', "w") as f: - f.write('filetemplate: PODmod' + info['caseName']+'%01d.f%05d\n') - f.write('firsttimestep: 0\n') - f.write('numtimesteps: %d\n' % (nmod+1)) + print('Writing: '+info['outputPath']+'DMDmod'+info['caseName']+'.nek5000') + if (info['module']=='POD'): + modname = 'PODmod' + elif (info['module']=='DMD'): + modname = 'DMDmod' + with open(info['outputPath']+modname+info['caseName']+'.nek5000', "w") as f: + f.write('filetemplate: '+ modname + info['caseName']+'%01d.f%05d\n') + f.write('firsttimestep: 0\n') + f.write('numtimesteps: %d\n' % (nmod+1)) + return diff --git a/pod/channel/MODULES/plotter.py b/pod/channel/MODULES/plotter.py index dac5b20..f6da4bb 100644 --- a/pod/channel/MODULES/plotter.py +++ b/pod/channel/MODULES/plotter.py @@ -25,7 +25,7 @@ -def pplot(A2,Lam2,nplt,iff): +def pplot(A2,Lam2,nplt,iff,info): """ Module to plot eigenvalues spectrum and cumulative sum Args: @@ -53,28 +53,30 @@ def pplot(A2,Lam2,nplt,iff): plt.show() elif (iff=='eigen'): - # Eigenvalues - print('singular values:',Lam2[0:nplt]) - plt.plot(np.arange(nplt),np.log10((Lam2[0:nplt])/sum(Lam2[0:nplt])),'^--r',label=r'$\lambda_k/\sum_{i=0}^m{\lambda_i}$') + # Eigenvalues without mean + print('singular values:',Lam2[1:nplt]) + plt.plot(np.arange(1,nplt),(Lam2[1:nplt])*100/sum(Lam2[1:nplt]),'^--r',label=r'$\lambda_k/\sum_{i=0}^m{\lambda_i}$') # plt.plot(np.cumsum(Lam2[0:nplt])/sum(Lam2[0:nplt]),'o-b',label='$\sum_{i=1}^k\lambda_i/\sum_{i=1}^m{\lambda_i}$') - plt.xlabel(r'$k$') - plt.ylabel('log10 energy') - plt.legend(loc='best') + plt.xlabel(r'Modes ($k$)') + plt.ylabel('Energy %') + #plt.legend(loc='best') plt.grid() - plt.show() + plt.savefig(info['outputPath']+"eigen.png") + plt.close() # here we only show the cummulative sum without the mean plt.plot(np.arange(nplt-1)+1,np.cumsum(Lam2[1:nplt])/sum(Lam2[1:nplt]),'o-b',label='$\sum_{i=1}^k\lambda_i/\sum_{i=1}^m{\lambda_i}$') plt.xlabel(r'$k$') plt.legend(loc='best') plt.grid() - plt.show() + plt.savefig(info['outputPath']+"sum.png") + plt.close() return -def pplotc(Lambdat): +def pplotc(Lambdat, info): """ Module to plot eigenvalues spectrum Args: @@ -86,7 +88,7 @@ def pplotc(Lambdat): thet=np.linspace(0,2*np.pi,100) plt.plot(np.cos(thet),np.sin(thet),'-k') for i in range(len(Lambdat)): - if abs(Lambdat[i])<1: + if abs(Lambdat[i])<=1: plt.plot(Lambdat[i].real,Lambdat[i].imag,'ob',mfc='b',label=r'$|\lambda|<1$') elif abs(Lambdat[i])>1: plt.plot(Lambdat[i].real,Lambdat[i].imag,'or',mfc='r',label=r'$|\lambda|>1$') @@ -95,11 +97,32 @@ def pplotc(Lambdat): plt.ylabel(r'$Im(\tilde{\Lambda})$') plt.grid(alpha=0.4) plt.title(r'Eigenvalues $\tilde{\mathbf{\Lambda}}$') - plt.show() + plt.savefig(info['outputPath']+"eigen.png") + plt.close() return +def saveOmega(omega, info): + """ + Module to save the continuous eigenvalues spectrum + Args: + - omega = time-continuous eigenvalues + - info = case info + """ + fileName = info['outputPath']+'DMDomega.txt' #'./OUTPUT/eigns.txt' + F = open(fileName,'w') + F.write("# Snapshots DMD result \n") + F.write("# Time Continuous Eigenvalues of linear operator A: Ax=b \n") + F.write("# According to our convention the mode 0 is the mean value \n") + F.write("# ------------------------------------------------------------------\n") + F.write("# ------------------------------------------------------------------\n") + F.write("# Real{Eigenvalues}\t Im{Eigenvalues}\t \n") + F.write("# ------------------------------------------------------------------\n") + for i in range(0,np.size(omega)): + F.write("%g\t%g\t \n" % \ + (omega[i].real, omega[i].imag/2/np.pi)) + F.close() def savetxt(Lam2,info): """ @@ -135,11 +158,11 @@ def savetxt(Lam2,info): F.write("# According to our convention the mode 0 is the mean value \n") F.write("# ------------------------------------------------------------------\n") F.write("# ------------------------------------------------------------------\n") - F.write("# Real{Eigenvalues}\t Im{Eigenvalues}\t \n") + F.write("# Real{Eigenvalues}\t Im{Eigenvalues}\t Mag{Eigenvalues}\t Arg{Eigenvalues}\t \n") F.write("# ------------------------------------------------------------------\n") for i in range(0,np.size(Lam2)): - F.write("%g\t%g\t \n" % \ - (Lam2[i].real, Lam2[i].imag )) + F.write("%g\t%g\t%g\t%g\t \n" % \ + (Lam2[i].real, Lam2[i].imag, abs(Lam2[i]), np.arctan2(Lam2[i].imag, Lam2[i].real))) F.close() return diff --git a/pod/channel/MODULES/reader.py b/pod/channel/MODULES/reader.py index 886fc76..b8f7a12 100644 --- a/pod/channel/MODULES/reader.py +++ b/pod/channel/MODULES/reader.py @@ -34,6 +34,7 @@ def read_input(filename): ifPickRead = eval(params['ifPickRead']) r = int(params['rpar']) timeprdc = int(params['tprdc']) + timestep = float(params['deltaT']) # All snapshots are used. @@ -49,7 +50,8 @@ def read_input(filename): 'startID':start_, 'endID':end_, 'variable':variable, - 'qoiName':qoiName} + 'qoiName':qoiName, + 'deltaT':timestep} # - mass matrix ( X-VELOCITY COMPONENT IS USED, SAVED HERE bm1 MATRIX) massName = 'bm1'+caseName info_m = {'dataPath':path, diff --git a/pod/channel/README.md b/pod/channel/README.md index deaee16..294acae 100644 --- a/pod/channel/README.md +++ b/pod/channel/README.md @@ -1,7 +1,7 @@ # Data-Driven Methods for Nek5000 Simulations Python code for relevant data-driven techniques with I/O designed for Nek5000. So far, the available methods are: - Proper Orthogonal Decomposition (POD) - - Dynamic Mode Decomposition (DMD), removed for this short tutorial + - Dynamic Mode Decomposition (DMD) ## General: - Nek5000 V19 diff --git a/pod/channel/ddmMain.py b/pod/channel/ddmMain.py index 3c3d6da..211aad3 100644 --- a/pod/channel/ddmMain.py +++ b/pod/channel/ddmMain.py @@ -12,14 +12,16 @@ """ import sys + # Local modules -sys.path.append('/home/pschlatt/NEK/git/TurbCourse/pod/channel/MODULES/') +sys.path.append('C:/Users/gaurav/Desktop/pod/MODULES') from snapMaker import snpAssembler,snpAssembler_sym from reader import read_input from dbMaker import dbMan,dbMan_sym from PODmodule import POD -#from DMDmodule import DMD -from plotter import pplot,pplotc,savetxt +from DMDmodule import DMD +from plotter import pplot,pplotc,savetxt,saveOmega +from freqAnalysis import freqPOD from outpWriter import modes,snaprcn,prdct from pickManager import pickReader,pickReader_sym @@ -38,7 +40,7 @@ qoiName,nsnap,nplt,r,timeprdc, \ outMod,outSnp,maxMode, \ if3D,ifsym,ifPickSave,ifPickRead, \ -info,info_m,info_s = read_input('/home/pschlatt/NEK/git/TurbCourse/pod/channel/input.txt') +info,info_m,info_s = read_input('C:/Users/gaurav/Desktop/pod/input.txt') @@ -49,7 +51,7 @@ if ifPickRead: # reading pickle Usnp,db,db_s,mvect = pickReader_sym(info,info_s) else: # building data - db,db_m,db_s,data_ms = dbMan_sym(info,'field',info_m,'mass',info_s,'field') # Database generation + db,db_m,db_s,data_ms = dbMan_sym(info,'field',if3D,info_m,'mass',info_s,'field') # Database generation Usnp,mvect = snpAssembler_sym(db,db_s,data_ms,info,nsnap,ifsym,if3D,ifPickSave) # Snapshots matrix assembly else: @@ -57,10 +59,10 @@ if ifPickRead: # reading pickle Usnp,db,mvect = pickReader(info) else: # building data - db,db_m,data_ms = dbMan(info,'field',info_m,'mass') # Database generation + db,db_m,data_ms = dbMan(info,'field',if3D,info_m,'mass') # Database generation Usnp,mvect = snpAssembler(db,data_ms,info,nsnap,ifsym,if3D,ifPickSave) # Snapshots matrix assembly - +print(db['data'][0].var) # ------- CALCULATION @@ -69,20 +71,21 @@ L,Lam2,R,A2 = POD(Usnp,mvect,nsnap,ifsym) elif (info['module']=='DMD'): # DMD - Phi,Lambdat,a1 = DMD(Usnp,mvect,nsnap,r,ifsym) + Phi,Lambdat,a1,omega = DMD(Usnp,mvect,nsnap,r,ifsym, info['deltaT']) # ------- PLOTS if (info['module']=='POD'): - pplot(A2,Lam2,nplt,'coeff') - pplot(A2,Lam2,nplt,'eigen') + # pplot(A2,Lam2,nplt,'coeff') + pplot(A2,Lam2,nplt,'eigen',info) savetxt(Lam2,info) + freqPOD(A2,10,nplt,0.1,info) elif (info['module']=='DMD'): - pplotc(Lambdat) + pplotc(Lambdat, info) savetxt(Lambdat,info) - + saveOmega(omega,info) # ------- OUTPUT @@ -91,6 +94,6 @@ snaprcn(outSnp,db['data'][0],info,L,A2,maxMode,if3D) elif (info['module']=='DMD'): modes(outMod,db['data'][0],info,Phi.real,if3D) - prdct(timeprdc,db['data'][0],info,Phi,a1,Lambdat,if3D) + #prdct(timeprdc,db['data'][0],info,Phi,a1,Lambdat,if3D) diff --git a/pod/channel/input.txt b/pod/channel/input.txt index 8b8e548..4dc0bb0 100644 --- a/pod/channel/input.txt +++ b/pod/channel/input.txt @@ -27,39 +27,41 @@ The reader module takes as input each variable one space after the equal. 12) To read compressed pickle file 13) Others parameters - (FOR DMD) Reduced rank parameter. Truncation factor for DMD construction/prediction -- Number of the initial snapshot +- Number of the initial snapshots +- TIme interval between snapshots #0) - module = POD + module = DMD #1) - path_in = /home/pschlatt/NEK/run/ - path_out = /home/pschlatt/NEK/run/PODmodes/ + path_in = path/to/dns/data + path_out = path/to/dmd/output #2) - caseName = channel + caseName = cyl #3) - nsnap = 1000 + nsnap = 200 #4) - if3D = True + if3D = False #5) - variable = vector - qoiName = velocity + variable = scalar + qoiName = temperature #6) ifsym = False path_m = HERE caseName_m = HERE #7) - nplt = 50 + nplt = 20 #8) - outMod = 50 + outMod = 10 #9) - outSnp = 200 + outSnp = 2 #10) - maxMode = 2 + maxMode = 10 #11) ifPickSave = False #12) ifPickRead = False #13) - rpar = 5000 - tprdc = 999 + rpar = 21 + tprdc = 1 + deltaT = 0.1 \ No newline at end of file