Commit b6163f4f authored by Felix Kong's avatar Felix Kong

vortex DMD: use class only, rm old "script" code

parent fa635adb
...@@ -37,6 +37,7 @@ classdef dmd < handle ...@@ -37,6 +37,7 @@ classdef dmd < handle
if nargin == 2 if nargin == 2
obj.r = varargin{1}; obj.r = varargin{1};
end end
tic;
[obj.U0,obj.S0,obj.V0] = svd(obj.X,'econ'); [obj.U0,obj.S0,obj.V0] = svd(obj.X,'econ');
% truncate % truncate
...@@ -51,6 +52,12 @@ classdef dmd < handle ...@@ -51,6 +52,12 @@ classdef dmd < handle
% Phi projects from weight-space to original data space (e.g. % Phi projects from weight-space to original data space (e.g.
% image space) % image space)
obj.Phi = obj.X2*obj.V*inv(obj.S)*obj.W; obj.Phi = obj.X2*obj.V*inv(obj.S)*obj.W;
% tell user how long it took
regressionTime = toc;
if(obj.verbose)
fprintf('DMD regression() took %f seconds.',regressionTime);
end
end end
function X2hat = reconstruct(obj) function X2hat = reconstruct(obj)
...@@ -72,6 +79,10 @@ classdef dmd < handle ...@@ -72,6 +79,10 @@ classdef dmd < handle
end end
function xhat = predict(obj,initial_condition,num_timesteps) function xhat = predict(obj,initial_condition,num_timesteps)
% predicts num_timesteps timesteps from an initial condition
% like reconstruct, but more general.
% can be used to predict more timesteps than were in the data
tic; tic;
if(isempty(obj.U0) || isempty(obj.S0) || isempty(obj.V0)) if(isempty(obj.U0) || isempty(obj.S0) || isempty(obj.V0))
error('Need to run regression() before running predict().'); error('Need to run regression() before running predict().');
...@@ -88,7 +99,10 @@ classdef dmd < handle ...@@ -88,7 +99,10 @@ classdef dmd < handle
assert(norm(imag(xhat)) <= 1e-6); assert(norm(imag(xhat)) <= 1e-6);
xhat = real(xhat); xhat = real(xhat);
toc; predictTime = toc;
if(obj.verbose)
fprintf('Predict() took %f seconds.\n',predictTime);
end
end end
end end
end end
......
clear all, close all, clc
load data/VortexData.mat % loads variables VORTALL, nx, ny
X = VORTALL(:,1:end-1);
X2 = VORTALL(:,2:end);
[U0,S0,V0] = svd(X,'econ');
%% Compute DMD (Phi are eigenvectors)
r = 21; % truncate at 21 modes
U = U0(:,1:r);
S = S0(1:r,1:r);
V = V0(:,1:r);
Atilde = U'*X2*V*inv(S);
[W,Lambda] = eig(Atilde);
Phi = X2*V*inv(S)*W;
eigs = diag(Lambda);
b0 = Phi\X(:,1);
%% prepare objects for animating the raw data
figure;
subplot(3,1,1);
xIC = reshape(VORTALL(:,1),nx,ny); % vortex field first timestep
im = imagesc(xIC); % plot vorticity field for first timestep
load CCcool.mat
colormap(CC); % use custom colormap
% clean up axes
set(gca,'XTick',[1 50 100 150 200 250 300 350 400 449],'XTickLabel',{'-1','0','1','2','3','4','5','6','7','8'})
set(gca,'YTick',[1 50 100 150 199],'YTickLabel',{'2','1','0','-1','-2'});
set(gcf,'Position',[100 100 300 130])
axis equal
hold on
% plot cylinder
theta = (1:100)/100'*2*pi;
x = 49+25*sin(theta);
y = 99+25*cos(theta);
fill(x,y,[.3 .3 .3]) % place cylinder
plot(x,y,'k','LineWidth',1.2) % cylinder boundary
title('Data')
colorbar;
% add contour lines (positive = solid, negative = dotted);
[~,cneg] = contour(xIC,[-5.5:.5:-.5 -.25 -.125],':k','LineWidth',1.2);
[~,cpos] = contour(xIC,[.125 .25 .5:.5:5.5],'-k','LineWidth',1.2);
%% prepare objects to animate the DMD reconstruction
subplot(3,1,2);
imDmd = imagesc(xIC);
load CCcool.mat
colormap(CC); % use custom colormap
% clean up axes
set(gca,'XTick',[1 50 100 150 200 250 300 350 400 449],'XTickLabel',{'-1','0','1','2','3','4','5','6','7','8'})
set(gca,'YTick',[1 50 100 150 199],'YTickLabel',{'2','1','0','-1','-2'});
set(gcf,'Position',[100 74 604 713])
axis equal
hold on
% plot cylinder
theta = (1:100)/100'*2*pi;
x = 49+25*sin(theta);
y = 99+25*cos(theta);
fill(x,y,[.3 .3 .3]) % place cylinder
plot(x,y,'k','LineWidth',1.2) % cylinder boundary
% add contour lines (positive = solid, negative = dotted);
[~,cnegDmd] = contour(xIC,[-5.5:.5:-.5 -.25 -.125],':k','LineWidth',1.2);
[~,cposDmd] = contour(xIC,[.125 .25 .5:.5:5.5],'-k','LineWidth',1.2);
titlename = sprintf('DMD reconstruction (%d dimensional, sorta)', r);
title(titlename)
colorbar;
%% Reconstruction error
subplot(3,1,3);
imErr = surf(zeros(size(xIC)));
set(gca,'XTick',[1 50 100 150 200 250 300 350 400 449],'XTickLabel',{'-1','0','1','2','3','4','5','6','7','8'})
set(gca,'YTick',[1 50 100 150 199],'YTickLabel',{'2','1','0','-1','-2'});
set(gcf,'Position',[100 74 604 713])
axis equal
hold on
imErr.EdgeAlpha = 0;
zlim([-0.01 0.01]);
% plot cylinder
theta = (1:100)/100'*2*pi;
x = 49+25*sin(theta);
y = 99+25*cos(theta);
fill(x,y,[.3 .3 .3]) % place cylinder
plot(x,y,'k','LineWidth',1.2) % cylinder boundary
title('Reconstruction error')
colorbar;
%% Run the animation
for tstep = 1:size(VORTALL,2) % every column is a timestep
% update ground truth
xTrue = reshape(VORTALL(:,tstep),nx,ny);
im.CData = xTrue;
cneg.ZData = xTrue;
cpos.ZData = xTrue;
% update dmd reconstruction
xkvec = Phi*Lambda^(tstep-1)*b0;
assert(all(abs(imag(xkvec)) <= 1e-6));
xk = reshape(real(xkvec),nx,ny);
imDmd.CData = xk;
cnegDmd.ZData = xk;
cposDmd.ZData = xk;
% update error plot
errork = xTrue-xk;
fprintf('Max abs error: %f \n', max(abs(errork),[],'all'));
imErr.ZData = errork;
view([0 90]);
drawnow;
% pause(0.1);
end
clear all, close all, clc clear all, close all, clc
load data/VortexData.mat % loads variables VORTALL, nx, ny load data/VortexData.mat % loads variables VORTALL, nx, ny
vortex_dmd = dmd(VORTALL); % initialize DMD object
X = VORTALL(:,1:end-1); verbose = true; % print how long things take
X2 = VORTALL(:,2:end); vortex_dmd = dmd(VORTALL,verbose);
[U0,S0,V0] = svd(X,'econ');
%% Compute DMD (Phi are eigenvectors) % perform DMD regression on data
r = 21; % truncate at 21 modes truncation_rank = 21; % rank to truncate reduced dynamics to
U = U0(:,1:r); vortex_dmd.regression(truncation_rank);
S = S0(1:r,1:r);
V = V0(:,1:r);
Atilde = U'*X2*V*inv(S);
[W,Lambda] = eig(Atilde);
Phi = X2*V*inv(S)*W;
eigs = diag(Lambda);
b0 = Phi\X(:,1); % reconstruct VORTALL(:,2:end) using DMD approximation
X2hat = vortex_dmd.reconstruct(); % not used in the script
vortex_dmd.regression(r); % use predict() to reconstruct the data a different way
X2hat = vortex_dmd.reconstruct(); % values may differ very slightly from X2hat
xhat = vortex_dmd.predict(X(:,1),size(VORTALL,2)); initial_condition = VORTALL(:,1);
num_timesteps = size(VORTALL,2); % could use predict() for more timesteps too
xhat = vortex_dmd.predict(initial_condition,num_timesteps);
%% No calculations after this line.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% everything after this line is plotting
% there are no calculations after this line
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% prepare objects for animating the raw data %% prepare objects for animating the raw data
...@@ -73,7 +76,7 @@ plot(x,y,'k','LineWidth',1.2) % cylinder boundary ...@@ -73,7 +76,7 @@ plot(x,y,'k','LineWidth',1.2) % cylinder boundary
% add contour lines (positive = solid, negative = dotted); % add contour lines (positive = solid, negative = dotted);
[~,cnegDmd] = contour(xIC,[-5.5:.5:-.5 -.25 -.125],':k','LineWidth',1.2); [~,cnegDmd] = contour(xIC,[-5.5:.5:-.5 -.25 -.125],':k','LineWidth',1.2);
[~,cposDmd] = contour(xIC,[.125 .25 .5:.5:5.5],'-k','LineWidth',1.2); [~,cposDmd] = contour(xIC,[.125 .25 .5:.5:5.5],'-k','LineWidth',1.2);
titlename = sprintf('DMD reconstruction (%d dimensional dynamics)', r); titlename = sprintf('DMD reconstruction (%d dimensional dynamics)', truncation_rank);
title(titlename) title(titlename)
colorbar; colorbar;
...@@ -101,7 +104,6 @@ title('Reconstruction error') ...@@ -101,7 +104,6 @@ title('Reconstruction error')
colorbar; colorbar;
%% Run the animation %% Run the animation
tic;
for tstep = 1:size(VORTALL,2) % every column is a timestep for tstep = 1:size(VORTALL,2) % every column is a timestep
% update ground truth % update ground truth
xTrue = reshape(VORTALL(:,tstep),nx,ny); xTrue = reshape(VORTALL(:,tstep),nx,ny);
...@@ -110,12 +112,7 @@ for tstep = 1:size(VORTALL,2) % every column is a timestep ...@@ -110,12 +112,7 @@ for tstep = 1:size(VORTALL,2) % every column is a timestep
cpos.ZData = xTrue; cpos.ZData = xTrue;
% update dmd reconstruction % update dmd reconstruction
% xkvec = Phi*Lambda^(tstep-1)*b0; xk = reshape(xhat(:,tstep),nx,ny);
% assert(all(abs(imag(xkvec)) <= 1e-6));
% xk = reshape(real(xkvec),nx,ny);
% reconstruction(:,tstep) = xkvec;
xk = reshape(real(xhat(:,tstep)),nx,ny);
imDmd.CData = xk; imDmd.CData = xk;
cnegDmd.ZData = xk; cnegDmd.ZData = xk;
...@@ -130,10 +127,5 @@ for tstep = 1:size(VORTALL,2) % every column is a timestep ...@@ -130,10 +127,5 @@ for tstep = 1:size(VORTALL,2) % every column is a timestep
drawnow; drawnow;
% pause(0.1); % pause(0.1);
end end
toc;
disp('End of animation.'); disp('End of animation.');
%%
% xhat = vortex_dmd.predict(X(:,1),size(VORTALL,2));
% disp(norm(reconstruction-xhat))
% disp(norm(reconstruction(:,2:end)-vortex_dmd.reconstruct()))
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment