Commit fa635adb authored by Felix Kong's avatar Felix Kong

tested DMD working on vortex data

parent e294c32f
classdef dmd
classdef dmd < handle
%DMD Summary of this class goes here
% Detailed explanation goes here
......@@ -26,27 +26,26 @@ classdef dmd
obj.X = data(:,1:end-1);
obj.X2 = data(:,2:end);
obj.r = size(obj.X,2); % no truncation by default
if nargin == 1
if nargin == 2
obj.verbose = varargin{1};
else
obj.verbose = false;
end
obj.regression(); % it's fast, just run it on construction
end
function obj = regression(obj,varargin) % argument is truncation rank r
if nargin == 1
function regression(obj,varargin) % argument is truncation rank r
if nargin == 2
obj.r = varargin{1};
end
[obj.U0,obj.S0,obj.V0] = svd(obj.X,'econ');
% truncate
obj.U = obj.U0(:,1:obj_r);
obj.U = obj.U0(:,1:obj.r);
obj.S = obj.S0(1:obj.r,1:obj.r);
obj.V = obj.V0(:,obj.r);
obj.V = obj.V0(:,1:obj.r);
% construct weight-space dynamics
obj.Atilde = obj.U'*obj.X2*obj.V/(obj.S);
obj.Atilde = obj.U'*obj.X2*obj.V*inv(obj.S);
[obj.W,obj.Lambda] = eig(obj.Atilde);
% Phi projects from weight-space to original data space (e.g.
......@@ -61,19 +60,35 @@ classdef dmd
% useful for getting a feel for how accurate the low-rank
% approximation is.
if(isempty(obj.U0) || isempty(obj.S0) || isempty(obj.V0))
error('Need to run regression() before running reconstruct().');
end
% UNTESTED, NOT SURE IF THIS IS MATHEMATICALLY CORRECT
X2hat = obj.Phi*obj.Lambda*obj.Phi\obj.X;
X2hatComplex = obj.Phi*obj.Lambda*(obj.Phi\obj.X);
assert(norm(imag(X2hatComplex)) <= 1e-6);
X2hat = real(X2hatComplex);
end
function xhat = predict(obj,initial_condition,num_timesteps)
assert(size(initial_condition) == size(obj.X,1));
tic;
if(isempty(obj.U0) || isempty(obj.S0) || isempty(obj.V0))
error('Need to run regression() before running predict().');
end
assert(all(size(initial_condition) == [size(obj.X,1),1]));
xhat = nan(length(initial_condition),num_timesteps);
xhat(:,1) = obj.Phi\initial_condition;
xhat(:,1) = initial_condition;
b0 = obj.Phi\initial_condition;
for k = 2:num_timesteps
xhat(:,k) = obj.Phi*obj.Lambda*obj.Phi\xhat(:,k-1);
xhat(:,k) = obj.Phi*obj.Lambda^(k-1)*b0;
end
assert(norm(imag(xhat)) <= 1e-6);
xhat = real(xhat);
toc;
end
end
end
......
clear all, close all, clc
load data/VortexData.mat % loads variables VORTALL, nx, ny
vortex_dmd = dmd(VORTALL);
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);
vortex_dmd.regression(r);
X2hat = vortex_dmd.reconstruct();
xhat = vortex_dmd.predict(X(:,1),size(VORTALL,2));
%% 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 data/colormap.mat % loads variable 'CC'
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);
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 dynamics)', 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
tic;
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);
% reconstruction(:,tstep) = xkvec;
xk = reshape(real(xhat(:,tstep)),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
toc;
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