247 resolvent analysis

.gitignore

@@ -40,6 +40,8 @@ pyLOM/DMD/wrapper.c
 pyLOM/DMD/wrapper.html
 pyLOM/SPOD/wrapper.c
 pyLOM/SPOD/wrapper.html
+pyLOM/RES/wrapper.c
+pyLOM/RES/wrapper.html
 pyLOM/vmmath/*.c
 pyLOM/vmmath/*.html
 pyLOM/inp_out/*.c

Examples/RES/example_RES_jet.py

@@ -0,0 +1,53 @@
+from __future__ import print_function, division
+
+import mpi4py
+mpi4py.rc.recv_mprobe = False
+
+import os, numpy as np
+import pyLOM
+
+
+# Parameters
+DATAFILE = './DATA/jetLES.h5'
+VARIABLES = 'PRESS'
+
+param = 2 * np.pi # Define the normalization parameter for the frequency (in St)
+f = 0.8 # Desired frequency (in St)
+n_modes = 5 # Desired number of modes to save 
+modes = np.arange(1,n_modes+1,dtype=np.int32)
+
+
+# Load the mesh
+m = pyLOM.Mesh.load(DATAFILE)
+pyLOM.pprint(0,'mesh loaded', flush=True)
+
+
+# Load the dataset
+d = pyLOM.Dataset.load(DATAFILE,ptable=m.partition_table)
+X = d[VARIABLES]
+t = d.get_variable('time')
+dt = t[1] - t[0]
+pyLOM.pprint(0,'dataset loaded', flush=True)
+
+
+# Compute the DMD of the case
+muReal, muImag, Phi, bJov = pyLOM.DMD.run(X, r=2e-1, remove_mean=True)
+delta, omega = pyLOM.DMD.frequency_damping(muReal,muImag,dt)
+freq = omega / param
+
+
+# Compute the Resolvent Analysis
+U, S, V = pyLOM.RES.run(Phi, delta, freq, f=f, Q=None)
+
+
+# Extract the desired modes
+U2, V2 = pyLOM.RES.extract_modes(U,V,1,len(d),modes=modes,kind='real') # kind can be 'real', 'imag' or 'abs'
+
+
+# Write the modes to be visualized with paraview
+d.add_field(f'forcing_modes',len(modes),V2)
+d.add_field(f'response_modes',len(modes),U2)
+pyLOM.io.pv_writer(m,d,f'modes_{f}',basedir='./',instants=[0],times=[0.],vars=[f'forcing_modes',f'response_modes'],fmt='vtkh5')
+
+
+pyLOM.cr_info()

Examples/RES/example_RES_plots_jet.py

@@ -0,0 +1,70 @@
+from __future__ import print_function, division
+
+import mpi4py
+mpi4py.rc.recv_mprobe = False
+
+import os, numpy as np
+import matplotlib.pyplot as plt
+import pyLOM
+
+pyLOM.gpu_device(gpu_per_node=4)
+
+
+# Parameters
+DATAFILE = './DATA/jetLES.h5'
+VARIABLE = 'PRESS'
+
+param = 2 * np.pi # Define the normalization parameter for the frequency (in St)
+f_list = [0.2, 0.4, 0.6, 0.8, 1.0] # Desired frequencies (in St)
+n_modes = 5 # Desired number of modes to save 
+modes = np.arange(1,n_modes+1,dtype=np.int32)
+
+
+# Load the mesh
+m = pyLOM.Mesh.load(DATAFILE)
+pyLOM.pprint(0,'mesh loaded', flush=True)
+
+
+# Load the dataset
+d = pyLOM.Dataset.load(DATAFILE,ptable=m.partition_table).to_gpu([VARIABLE])
+X = d[VARIABLE]
+t = d.get_variable('time')
+dt = t[1] - t[0]
+pyLOM.pprint(0,'dataset loaded', flush=True)
+
+
+# Compute the DMD of the case
+muReal, muImag, Phi, bJov = pyLOM.DMD.run(X, r=2e-1, remove_mean=True)
+delta, omega = pyLOM.DMD.frequency_damping(muReal,muImag,dt)
+freq = omega / param
+
+
+S_list = np.empty((0,Phi.shape[1]))
+# Compute the Resolvent Analysis for every frequency
+for f in f_list:
+    U, S, V = pyLOM.RES.run(Phi, delta, freq, f=f, Q=None)
+    S_list = np.vstack([S_list, S])
+
+    # Extract the desired modes
+    U2, V2 = pyLOM.RES.extract_modes(U,V,1,len(d),modes=modes,kind='real') # kind can be 'real', 'imag' or 'abs'
+
+
+    # Write the modes to be visualized with paraview
+    d.add_field(f'forcing_modes',len(modes),V2)
+    d.add_field(f'response_modes',len(modes),U2)
+    pyLOM.io.pv_writer(m,d.to_cpu(['forcing_modes','response_modes']),f'modes_{f}',basedir='./',instants=[0],times=[0.],vars=['forcing_modes','response_modes'],fmt='vtkh5')
+
+
+if pyLOM.utils.is_rank_or_serial(0):
+    # Plot the cumulative energy gains
+    pyLOM.RES.plotEnergy(S_list[0,:])
+    plt.savefig('energy.png', dpi=300)
+
+    # Plot the energy gains vs frequency for the desired modes
+    pyLOM.RES.plotEvW(S_list, f_list, modes)
+    plt.savefig('EvW.png', dpi=300)
+
+
+
+pyLOM.cr_info()
+pyLOM.show_plots()

Makefile

@@ -262,13 +262,14 @@ clean:
 	-@cd pyLOM; rm -rf POD/__pycache__ POD/*.c POD/*.cpp POD/*.html
 	-@cd pyLOM; rm -rf DMD/__pycache__ DMD/*.c DMD/*.cpp DMD/*.html
 	-@cd pyLOM; rm -rf SPOD/__pycache__ SPOD/*.c SPOD/*.cpp SPOD/*.html
+	-@cd pyLOM; rm -rf RES/__pycache__ RES/*.c RES/*.cpp RES/*.html
 	-@cd pyLOM; rm -rf vmmath/__pycache__ vmmath/*.c vmmath/*.cpp vmmath/*.html
 	-@cd pyLOM; rm -rf inp_out/__pycache__ inp_out/*.c inp_out/*.cpp inp_out/*.html
 	-@cd pyLOM; rm -rf NN/__pycache__ NN/architectures/__pycache__ 
 
 cleanall: clean
 	-@rm -rf build
-	-@cd pyLOM; rm vmmath/*.so POD/*.so DMD/*.so SPOD/*.so
+	-@cd pyLOM; rm vmmath/*.so POD/*.so DMD/*.so SPOD/*.so RES/*.so
 
 ifeq ($(USE_MKL),ON)
 uninstall_vector_matrix: uninstall_mkl 

options.cfg

@@ -24,7 +24,7 @@ USE_GESVD        = OFF
 USE_COMPILED     = ON
 # Comma separated list of the modules to be compiled
 # if USE_COMPILED = ON
-MODULES_COMPILED = MATH.MATHS,MATH.AVERAGING,MATH.QR,MATH.SVD,MATH.FFT,MATH.GEOMETRIC,MATH.TRUNCATION,MATH.STATS,MATH.REGRESSION,ROM.POD,ROM.DMD,ROM.SPOD
+MODULES_COMPILED = MATH.MATHS,MATH.AVERAGING,MATH.QR,MATH.SVD,MATH.FFT,MATH.GEOMETRIC,MATH.TRUNCATION,MATH.STATS,MATH.REGRESSION,MATH.LINEAR,MATH.DATAPROCESSING,ROM.POD,ROM.DMD,ROM.SPOD,ROM.RES
 
 
 ## Optimization, host and CPU type
@@ -39,11 +39,10 @@ TUNE = skylake
 PYTHON = python3
 PIP    = pip3
 
-
 ## Versions of the libraries
 #
 ONEAPI_VERS   = 2024.2.0.634
-OPENBLAS_VERS = 0.3.17
+OPENBLAS_VERS = 0.3.34
 LAPACK_VERS   = 3.9.0
 KISSFFT_VERS  = 131.1.0
 FFTW_VERS     = 3.3.8

pyLOM/DMD/__init__.py

@@ -7,7 +7,7 @@
 # Last rev: 30/09/2021
 
 # Functions coming from DMD
-from .wrapper import run, frequency_damping, reconstruction_jovanovic
+from .wrapper import run, frequency_damping, reconstruction_jovanovic, run_new
 from .utils   import extract_modes, save, load
 from .plots   import plotMode, ritzSpectrum, amplitudeFrequency, dampingFrequency, plotResidual, plotSnapshot
 

pyLOM/DMD/wrapper.py

@@ -10,8 +10,7 @@
 import numpy as np
 
 from ..utils.gpu import cp
-from ..vmmath    import vecmat, matmul, temporal_mean, subtract_mean, tsqr_svd, transpose, eigen, cholesky, diag, polar, vandermonde, conj, inv, flip, matmulp, vandermondeTime
-from ..POD       import truncate
+from ..vmmath    import matmul, temporal_mean, subtract_mean, transpose, eigen, cholesky, diag, polar, vandermonde, conj, inv, flip, vandermondeTime, linear_operator, concatenate, separate
 from ..utils     import cr_nvtx as cr, cr_start, cr_stop
 
 
@@ -60,32 +59,16 @@ def run(X, r, remove_mean = True):
 		- b:        Amplitude of the DMD modes
 		- X_DMD:    Reconstructed flow
 	'''
-	# Remove temporal mean or not, depending on the user choice
-	if remove_mean:
-		cr_start('DMD.temporal_mean',0)
-		#Compute temporal mean
-		X_mean = temporal_mean(X)
-		#Subtract temporal mean
-		Y = subtract_mean(X, X_mean)
-		cr_stop('DMD.temporal_mean',0)
+	# Prepare matrices and calculate the linear operator
+	if (type(X) is list):
+		Y, Z = concatenate(X, remove_mean=remove_mean)
+		U, S, VT, Atilde = linear_operator(Y, Z, r)
+
 	else:
-		Y = X.copy()
-
-	# Compute SVD
-	cr_start('DMD.SVD',0)
-	U, S, VT = tsqr_svd(Y[:, :-1])
-	cr_stop('DMD.SVD',0)
-	# Truncate according to residual
-	cr_start('DMD.truncate', 0)
-	U, S, VT = truncate(U, S, VT, r)
-	cr_stop('DMD.truncate', 0)
-
-	# Project A (Jacobian of the snapshots) into POD basis
-	cr_start('DMD.linear_mapping',0)
-	aux1   = matmulp(transpose(U), Y[:, 1:])
-	aux2   = transpose(vecmat(1./S, VT))
-	Atilde = matmul(aux1, aux2)
-	cr_stop('DMD.linear_mapping',0)
+		Y, Z = separate(X, remove_mean=remove_mean)
+		U, S, VT, Atilde = linear_operator(Y, Z, r)
+
+	del U
 
 	# Eigendecomposition of Atilde: Eigenvectors given as complex matrix
 	# NOTE: there is no implementation of eig in cupy yet
@@ -96,12 +79,12 @@ def run(X, r, remove_mean = True):
 	w      = cp.asarray(w)      if type(Atilde) is cp.ndarray else w
 
 	# Mode computation
-	Phi =  matmul(matmul(matmul(Y[:, 1:], transpose(VT)), diag(1/S)), w)/(muReal + muImag*1J)
+	Phi =  matmul(matmul(matmul(Z, transpose(VT)), diag(1/S)), w)/(muReal + muImag*1J)
 	cr_stop('DMD.modes',0)
 
 	# Amplitudes according to: Jovanovic et. al. 2014 DOI: 10.1063
 	cr_start('DMD.amplitudes',0)
-	Vand = vandermonde(muReal, muImag, muReal.shape[0], Y.shape[1]-1)
+	Vand = vandermonde(muReal, muImag, muReal.shape[0], Y.shape[1])
 	P    = matmul(transpose(conj(w)), w)*conj(matmul(Vand, transpose(conj(Vand))))
 	Pl   = cholesky(P)
 	G    = matmul(diag(S), VT)