-
Notifications
You must be signed in to change notification settings - Fork 0
/
svd_data.m
58 lines (43 loc) · 1.25 KB
/
svd_data.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
clc
clear
close all
% load data
velocity = load('velocity.dat');
pressure = load('pressure.dat');
%%
% transpose the data to get each snapshot as a column vector
velocity = velocity';
pressure = pressure';
%%
% determine the mean of the quanitity under consideration
mean_vel = mean(velocity,2);
mean_pres = mean(pressure,2);
% determine the dynamic part of the quantity
dynamic_vel = velocity - repmat(mean_vel, 1, 500);
dynamic_pres = pressure - repmat(mean_pres, 1, 500);
% set the values below the tolerance as zero to avoid precision errors
dynamic_vel(abs(dynamic_vel) < 1e-15) = 0;
dynamic_pres(abs(dynamic_pres) < 1e-15) = 0;
% Perform SVD on the data
[Uv,Sv,Vv] = svd(dynamic_vel,"econ");
[Up,Sp,Vp] = svd(dynamic_pres,"econ");
% Write the first n modes to the files
Uvr = Uv(:,1:10);
Upr = Up(:,1:10);
writematrix(Uvr, 'v_pod_modes.txt')
writematrix(Upr, 'p_pod_modes.txt')
% Plot the energy content of the modes
eig_v = diag(Sv);
eig_p = diag(Sp);
energy_v = zeros(1,10);
energy_p = zeros(1,10);
n = zeros(1,11);
for i = 1:11
n(i) = i-1;
energy_p(i) = sum(eig_p(1:i-1))/sum(eig_p);
energy_v(i) = sum(eig_v(1:i-1))/sum(eig_v);
end
figure(1)
plot(n, energy_v)
figure(2)
plot(n, energy_p)