POD_matlab程序

% % Create matrix will all fluctuating velocity components for each snapshot in a column
% %
% % Do POD analysis
% n 速度场的总的节点数
% m 时间文件的个数，生成矩阵的维数
% Uall 将所有文件的速度U和速度Ｖ重构得到的POD分析的矩阵
% R 为二阶相关函数构成的矩阵
% eValue 为矩阵R特征值
% eVec 为矩阵R特征向量
% menergy为各特征值所对应的能量百分比
% phi为得到的基函数
clear;clc;

kx=92;
ky=99;
kz=22;
lx1=11;
lx2=82;
lz1=4;
lz2=19;
lx=lx2-lx1+1;
lz=lz2-lz1+1;
n=lx*ky;
dr=0.6872;
x=-35.5357:dr:-35.5357+dr*（kx-1）;
y=-31.3456:dr:-31.3456+dr*（ky-1）;
z=0.6096:dr:0.6096+dr*（kz-1）;
dt=0.002;
t0=fix（lx*0.75）;
lt=200;
% file1=t0+1;
% file2=t0+1+lt;
file1=1;
file2=t0+lt+t0;
lf=file2-file1+1;

m=importdata（‘mean velocity_flipxy.dat‘）;
m1=m.data（: 4）;
m2=reshape（m1 kx ky kz）;
m3=m2（lx1:lx2 : :）;

U=zeros（ n lf）;
for iz=4
    mu = reshape（m3（: : iz） [] 1）;
    for i=file1:file2
        fname  = sprintf（‘d%04d.dat‘ i）;
        a      = importdata（fname）;
        temp   = a.data（:4）;
        temp   = reshape（temp kx ky kz）;
        iu     = -temp（: : iz）;
        iu     = flipud（iu）;
        iu     = fliplr（iu）;
        iu     = iu（lx1:lx2 :）;
%         iu     = iu‘;
        iu     = reshape（iu [] 1）;
        
        U（: i-file1+1） = iu-mu;
        disp（fname）;
    end
    
    R=U‘*U;  % solve: eV is eigenvectors D is eigenvalues in diagonal matrix
    [eVD]=eig（R）; % sort eigenvalues in ascending order-I is sorted index vector
    disp（‘eigenvalues have been done.‘）;
    [LI]=sort（diag（D））;
    for i=1:length（D）
        eValue（length（D）+1-i）=L（i）;   % Eigenvalues sorted in descending order
        eVec（:length（D）+1-i）=eV（:I（i））;  % Eigenvectors sorted in the same order
    end;
    for  k=1:length（eValue）
        if eValue（k）<=0
            eValue（k）=0;              % last eigenvalue should be zero
        end
    end
    menergy=eValue/sum（eValue）;       % relative energy associated with mode m
    
    %计算能量百分比
    for j=1:lf
        cumulant（j）=sum（menergy（1:j））;
    end
    fname=sprintf（‘pod_energy_z%02d.dat‘ iz）;
    fenergy=fopen（fname ‘w‘）;
    fprintf（fenergy‘VARIABLES = “mode“ “Percent“ “Cumulant Percent “\r\n‘）;
    for i=1:length（eValue）
        fprintf（fenergy‘%d\t%g\t%g\r\n‘ i menergy（i） cumulant（i））;
    end
    fclose（fenergy）;
    disp（‘energy percent have been done.‘）;
    
    %计算归一化的基函数
    for i=1:lf
        tmp=U*eVec（:i）;
        phi（:i）=tmp/norm（tmp）;
    end;
    disp（‘normalizations have been done.‘）;
    
    % 输出基函数的等值线图
    fname=sprintf（‘pod_phi_z%02d.dat‘ iz）;
    fphi=fopen（fname‘w‘）;
    fprintf（fphi‘VARIABLES = “x“ “y“ “phi u“\r\n‘）;
    for k=1:50
        fprintf（fphi‘\r\nZone T=“mode %d“i=%d j=%d\r\n‘ k lx ky）;
        for iy=1:ky
            for ix=1:lx
                fprintf（fp