From 348f0db7adf8fca35a39222c9735495abd77005b Mon Sep 17 00:00:00 2001 From: Weibo Date: Tue, 12 Nov 2024 16:52:32 +0800 Subject: [PATCH] POD analysis for ramp --- source/pod_ramp.py | 217 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 217 insertions(+) create mode 100644 source/pod_ramp.py diff --git a/source/pod_ramp.py b/source/pod_ramp.py new file mode 100644 index 0000000..5dbbb88 --- /dev/null +++ b/source/pod_ramp.py @@ -0,0 +1,217 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Wed Nov 6 16:00:05 2024 + pod for ramp cases +@author: weibo +""" + +# %% Load necessary module +import matplotlib +import matplotlib.pyplot as plt +import numpy as np +import pandas as pd +import pod as pod +import plt2pandas as p2p +from timer import timer +from scipy.interpolate import griddata +import os +import sys + +tsize = 13 +nsize = 10 +matplotlib.rcParams["xtick.direction"] = "out" +matplotlib.rcParams["ytick.direction"] = "out" +matplotlib.rc("font", size=tsize) +plt.close("All") +plt.rc("text", usetex=True) +font = { + "family": "Times New Roman", # 'color' : 'k', + "weight": "normal", + "size": "large", +} +cm2in = 1 / 2.54 + + +# %% prep data +path = "/media/weibo/VID2/AAS/ramp_st14_2nd/" +p2p.create_folder(path) +pathP = path + "probes/" +pathF = path + "Figures/" +pathM = path + "MeanFlow/" +pathS = path + "SpanAve/" +pathT = path + "TimeAve/" +pathI = path + "Instant/" +pathD = path + "Domain/" +pathPOD = path + "POD/" +pathSS = path + 'TP_2D_Z_001/' +timepoints = np.arange(600, 900 + 0.25, 0.25) +dirs = sorted(os.listdir(pathSS)) +if np.size(dirs) != np.size(timepoints): + sys.exit("The NO of snapshots are not equal to the NO of timespoints!!!") +DataFrame = pd.read_hdf(pathSS + dirs[0]) +ind0 = (DataFrame['y'] == 0.0) +DataFrame['walldist'][ind0] = 0.0 + +grouped = DataFrame.groupby(['x', 'y']) +DataFrame = grouped.mean().reset_index() +NewFrame = DataFrame.query("x>=-160.0 & x<=60.0 & walldist>=0.0 & y<=36.0") + +ind = NewFrame.index.values +xval = DataFrame['x'][ind] +yval = DataFrame['y'][ind] +x, y = np.meshgrid(np.unique(xval), np.unique(yval)) +x1 = -160.0 +x2 = 60.0 +y1 = 0.0 +y2 = 36.0 + +var0 = 'u' +var1 = 'v' +var2 = 'p' +col = [var0, var1, var2] +fa = 1 +FirstFrame = DataFrame[col].values +Snapshots = FirstFrame[ind].ravel(order='F') + +m, = np.shape(Snapshots) +n = np.size(timepoints) +o = np.size(col) +if (m % o != 0): + sys.exit("Dimensions of snapshots are wrong!!!") +m = int(m/o) +varset = {var0: [0, m], + var1: [m, 2*m], + var2: [2*m, 3*m] + } +# %% load data +with timer("Load Data"): + for i in range(np.size(dirs)-1): + TempFrame = pd.read_hdf(pathSS + dirs[i+1]) + grouped = TempFrame.groupby(['x', 'y']) + TempFrame = grouped.mean().reset_index() + if np.shape(TempFrame)[0] != np.shape(DataFrame)[0]: + sys.exit('The input snapshots does not match!!!') + NextFrame = TempFrame[col].values + Snapshots = np.vstack((Snapshots, NextFrame[ind].ravel(order='F'))) + DataFrame += TempFrame +Snapshots = Snapshots.T + +AveFlow = DataFrame/np.size(dirs) +meanflow = AveFlow.query("x>=-160.0 & x<=60.0 & walldist>=0.0 & y<=36.0") + +# %% POD compute +if np.size(dirs) != np.size(timepoints): + sys.exit("The NO of snapshots are not equal to the NO of timespoints!!!") + +with timer("POD computing"): + eigval, eigvec, phi, coeff = \ + pod.pod(Snapshots, fluc=True, method='svd') + +meanflow.to_hdf(pathPOD + 'Meanflow.h5', 'w', format='fixed') +np.save(pathPOD + 'eigval', eigval) +np.save(pathPOD + 'eigvec', eigvec) +np.save(pathPOD + 'phi1', phi[:, :600]) +np.save(pathPOD + 'phi2', phi[:, 600:]) +np.save(pathPOD + 'coeff', coeff) + +# %% optional: load POD results +eigval = np.load(pathPOD + 'eigval.npy') +eigvec = np.load(pathPOD + 'eigvec.npy') +coeff = np.load(pathPOD + 'coeff.npy') +# phi1 = np.load(pathPOD + 'phi1.npy') +# phi2 = np.load(pathPOD + 'phi2.npy') +phi = np.hstack((np.load(pathPOD + 'phi1.npy'), np.load(pathPOD + 'phi2.npy'))) +# %%############################################################################ +""" + Plot eigenvalue spectrum +""" +EFrac, ECumu, N_modes = pod.pod_eigspectrum(97, eigval) +np.savetxt(pathPOD+'EnergyFraction97.dat', EFrac, fmt='%1.7e', delimiter='\t') +N_modes = -N_modes + 1201 + 2 +print("the NO of modes is ", N_modes) +var = var0 +matplotlib.rc('font', size=tsize) +fig1, ax1 = plt.subplots(figsize=(10*cm2in, 8*cm2in)) +xaxis = np.arange(0, N_modes) +ax1.scatter( + xaxis[1:], + EFrac[1:N_modes], + c='black', + marker='o', + s=10.0, +) # fraction energy of every eigval mode +# ax1.legend('E_i') +ax1.set_ylim(bottom=0) +ax1.set_xlabel('Mode', fontsize=tsize) +ax1.set_ylabel(r'$E_i$', fontsize=tsize) +ax1.grid(visible=True, which='both', linestyle=':') +ax1.tick_params(labelsize=nsize) +ax2 = ax1.twinx() # cumulation energy of first several modes +# ax2.fill_between(xaxis, ECumu[:N_modes], color='grey', alpha=0.5) +# ESum = np.zeros(N_modes) +ESum = ECumu[:N_modes] +ax2.plot(xaxis, ESum, color='grey', label=r'$ES_i$') +ax2.set_ylim([60, 100]) +ax2.set_ylabel(r'$ES_i$', fontsize=tsize) +ax2.tick_params(labelsize=nsize) +plt.tight_layout(pad=0.5, w_pad=0.2, h_pad=1) +plt.savefig(pathPOD+str(N_modes)+'_PODEigSpectrum80.svg', bbox_inches='tight') +plt.show() + +# %% specific mode in space +ind = 20 +var = var0 +fa = 1.0 # 1.7*1.7*1.4 # +x, y = np.meshgrid(np.unique(xval), np.unique(yval)) +newflow = phi[:, ind]*coeff[ind, 0] +modeflow = newflow[varset[var][0]:varset[var][1]] +print("The limit value: ", np.min(modeflow)*fa, np.max(modeflow)*fa) +u = griddata((xval, yval), modeflow, (x, y))*fa +corner = (x < 0.0) & (y < 0.0) +u[corner] = np.nan +matplotlib.rc('font', size=tsize) +fig, ax = plt.subplots(figsize=(14*cm2in, 5*cm2in)) +c1 = -0.02 # -0.01 # -0.006 +c2 = -c1 # 0.063 + +lev1 = np.linspace(c1, c2, 11) +lev2 = np.linspace(c1, c2, 6) +cbar = ax.contourf(x, y, u, cmap='coolwarm', levels=lev1, extend='both') +# ax.contour(x, y, u, levels=lev2, colors='k', linewidths=0.8, extend='both') +# colors=('#66ccff', '#e6e6e6', '#ff4d4d'), extend='both') +ax.set_xlim(-100, 60) +ax.set_ylim(0, 25) +ax.set_xticks(np.linspace(-100, 60, 5)) +ax.set_yticks(np.linspace(0, 25, 6)) +ax.tick_params(labelsize=nsize) +# cbar.cmap.set_under('#053061') +# cbar.cmap.set_over('#67001f') +ax.set_xlabel(r'$x/l_r$', fontsize=tsize) +ax.set_ylabel(r'$y/l_r$', fontsize=tsize) +# add colorbar +rg2 = np.linspace(c1, c2, 3) +cbaxes = fig.add_axes([0.2, 0.76, 0.30, 0.07]) # x, y, width, height +cbar1 = plt.colorbar(cbar, cax=cbaxes, orientation="horizontal", + ticks=rg2, extendrect='False') +cbar1.formatter.set_powerlimits((-2, 2)) +cbar1.ax.xaxis.offsetText.set_fontsize(nsize) +cbar1.update_ticks() +cbar1.set_label(r'$\varphi_{}$'.format(var), rotation=0, + x=-0.18, labelpad=-20, fontsize=tsize) +cbaxes.tick_params(labelsize=nsize) +# Add shock wave +shock = pd.read_csv(pathM + "ShockLineFit.dat", skipinitialspace=True) +ax.plot(shock.x, shock.y, 'gray', linewidth=1.0) +# Add sonic line +sonic = pd.read_csv(pathM+'SonicLine.dat', skipinitialspace=True) +ax.plot(sonic.x, sonic.y, 'g--', linewidth=1.0) +# Add boundary layer +# boundary = np.loadtxt(pathM+'BoundaryEdge.dat', skiprows=1) +# ax.plot(boundary[:, 0], boundary[:, 1], 'k', linewidth=1.0) +# Add dividing line(separation line) +# dividing = np.loadtxt(pathM+'BubbleLine.dat', skiprows=1) +# ax.plot(dividing[:, 0], dividing[:, 1], 'k--', linewidth=1.0) + +plt.savefig(pathPOD+var+'_PODMode'+str(ind)+'.svg', bbox_inches='tight') +plt.show() \ No newline at end of file