Skip reading of variables other than QOI from Snapshots and Frequency Analysis for POD

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.

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

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:

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,21 +66,43 @@ def dbCreator(info,iff):
 
 
 
-def dbCreator_v(info):
+def dbCreator_v(info, if3D):
     """
     The database is created for the solution field (all the snapshots).
     """
     path_    = info['dataPath']
     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()

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()

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