I am trying to compare two different snapshot proper orthogonal decomposition methods (eigenvalue problem and svd). I am able to obtain the energy fraction graph for the eigenvalue problem, but the svd method (in yellow) is returning much different data points (see attached figure). I believe I am not coding the E_frame correctly or the E_fraction correctly. The energy fraction graphs should be similar to each other. Can anyone help with this? I attached supporting functions, however, the files of the snapshots are too large to upload as a zip file.
clear
clc
close all
cd('C:\Users\blake\OneDrive\600 and 20')
% savePath='';
autoselect=0;
%% Manual file selection
if autoselect==0
[fileName,PathName,FilterIndex] = uigetfile({'*.tiff';'*.tif'},'MultiSelect','on');
nameIdx=(find(PathName=='\',2,'last'));
expCond=PathName(nameIdx(1)+1:nameIdx(2)-1);
end
cd(PathName)
%% Initialize variables and constants
%constants
Nrel=length(fileName);
Nt = 30;
temp = imread(fileName{1},1);
[Ny,Nx,~] = size(temp);
%Ny represent the y direction in cartesian notation but also represents rows in pictures
%Nx represent the x directoin in cartensian notation but also represents
%columns in pictures
%image crop to remove injector and some
image_p = processimage(temp);
image_p2 = edge(image_p);
JetCenterY = findjet(image_p,30);
imLength = Nx-50; %reduce right most dark edge, chamber's shadow
Nx = imLength - JetCenterY + 1;
%colormap for POD modes plotting
load('cmap_MOD')
% Allocate storage for data
E_fraction = zeros(Nrel+1,Nt+1); %% pcolor doesn't print last col/row, so add 1 on each direction
modes = zeros(Nx*Ny,Nt, Nrel);
images = zeros(Ny,Nx,Nt,Nrel);
% Detect start function
SOIframes = SOI_Determination(fileName); %start of injection
EOIframes = SOIframes + Nt; %end of injection
%% Load Data
for j = 1:Nrel
for i = 1:Nt
image = imread(fileName{j},SOIframes(j)+i-1);
frames = processimage(image);
images(:,:,i,j) = frames(:,JetCenterY:imLength);
end
end
q_mean = mean(images,4);
Q = images - q_mean; % see message below
%negative values imply that the mean does not have that feature, i.e. it is unique to cycle. because
%black is denoted by 0 while white has a value greater than 0 (255 in 8 bit
%notation). so when you subtracted a non-existant feature (background
%color > 0 ) from the spray (value of 0) the result will b a negative value
%
%% Compute POD - Snapshot via eig()
% Compute Space-only POD
[phiSO,lambdaSO,aSO] = SpaceOnlyPOD(Q);
% Mean and total energies
e_mean = sum(sum(sum(abs(q_mean).^2,1),2),3)./(Ny*Nx*Nt);
e_total = sum(lambdaSO) + e_mean;
% Make Space-only POD eigenvalue figure
figure(1);
hold on;
plot((1:Nrel),[e_mean;mean(lambdaSO(1:end-1,:),2)]./e_total,'ko');
set(gca,'Yscale','log');
%title([Case,': ','$e^{\prime} / e_{total}$ = ',num2str(100*sum(lambda)/e_total),'\%'],'Interpreter','Latex');
hold off
box on;
%% Compute POD - Snapshot via svd()
for k= 1:Nt
S=reshape(images(:,:,k,:),[Ny*Nx,Nrel]);
% compute correlation matrix
C=(S'*S)/(Nrel);
%% Solve the eigenvalue problem
[beta, lmd, ~]=svd(C);
% where beta = U (columns are POD modes = normalized unit vectors),
% lmd = E (diagonal matrix with singular values, [sigma]= energy captured in each mode)
% eigenvalues are [lambda] = sigma^2/M
%% Compute the normalized basis function
% basis function represent the extracted flow pattern (i.e. coherent
% structures). Not necesarily reflect a real flow structure. A real
% flow patter only appears prominent in a POD mode when the flow
% structure is nearly the same and in the same location in each flow
% field (i.e. snapshot)
phix = S*beta; % projection of scalar quantity (i.e. intensity) onto the eigenvectors
% Normalize basis
phi = normc(phix);
% double check if this is necessary**
%% Compute coefficients
coeff=S'*phi;
coeff=coeff';
coeffMatrix(:,:,k)=coeff;
coeff_sq=coeff.^2;
%E_frame=sum(coeff_sq(2:end,:))./coeff_sq(1,:);
E_frame=sqrt(sum(coeff_sq(2:end,:)))./sqrt(coeff_sq(1,:));
E_fraction(1:Nrel,k)=E_frame'*100;
% Rearrange basis vectors into single matrix
for n = 1:Nrel
modes(:,k,n)=phix(:,n);
end
end
modes = reshape(modes,[Ny, Nx, Nt, Nrel]);
% Plotting Energy Fractions Together
figure(6);
hold on;
plot((1:Nrel),[e_mean;mean(lambdaSO(1:end-1,:),2)]./e_total,'-o');
plot(E_fraction(:,1),'-o');
set(gca,'Yscale','log');
[t,s] = title('600 & 20 Energy Fraction Graph','Eig() = Orange, SVD() = Yellow');
xlabel('POD Mode #')
ylabel('Energy Fraction')
%title([Case,': ','$e^{\prime} / e_{total}$ = ',num2str(100*sum(lambda)/e_total),'\%'],'Interpreter','Latex');
hold off
box on;
