Skip to content

Commit

Permalink
add matlab driver
Browse files Browse the repository at this point in the history
  • Loading branch information
@Ping-Hsuan
Ping-Hsuan committed Jul 7, 2024
1 parent cfc8d71 commit e7d9011
Show file tree
Hide file tree
Showing 4 changed files with 330 additions and 0 deletions.
96 changes: 96 additions & 0 deletions bin/matlab_driver/BDFk_EXTk.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,96 @@
classdef BDFk_EXTk
properties
u
alphas
betas
dt
T_final
nsteps
iostep
ito
ext
hfac
rhs
end
methods
% Constructor
function obj = BDFk_EXTk(nsteps, dt, iostep, nb)
obj.alphas = zeros(3,3); % initialize vectors to hold EXTk coefficients
obj.betas = zeros(4,3); % initialize vectors to hold BDFk coefficients
obj.ext = zeros(nb,3); % initialize vectors to hold extrapolation vectors
obj.u = zeros(nb+1, 3); % initialize vectors to hold ROM coefficients
obj.hfac = [];
obj.dt = dt; % Simultaion time step size
obj.nsteps = nsteps; % Simulation total time steps
obj.iostep = iostep; % Every #steps to store ROM quantities
obj.T_final = obj.nsteps*obj.dt; % Simulation final time
end
function obj = setup(obj)
% Setup BDFk/EXTk coefficients

% Setup EXT1 coefficients
obj.alphas(1,1) = 1.0;

% Setup EXT2 coefficients
obj.alphas(1,2) = 2.0;
obj.alphas(2,2) = -1.0;

% Setup EXT3 coefficients
obj.alphas(1,3) = 3.0;
obj.alphas(2,3) = -3.0;
obj.alphas(3,3) = 1.0;

% BDF1 coefficients
obj.betas(1,1) = 1.0;
obj.betas(2,1) = -1.0;

% BDF2 coefficients
obj.betas(1,2) = 1.5;
obj.betas(2,2) = -2.0;
obj.betas(3,2) = 0.5;

% BDF3 coefficients
obj.betas(1,3) = 11.0/6;
obj.betas(2,3) = -3.0;
obj.betas(3,3) = 1.5;
obj.betas(4,3) = -1.0/3;
end

function obj = setrhs(obj, rom, u)

obj.ext(:,3) = obj.ext(:,2);
obj.ext(:,2) = obj.ext(:,1);
obj.ext(:,1) = rom.setrhs(u, "semi-implicit");

obj.rhs = obj.ext*obj.alphas(:,obj.ito);
% Later move rom.bu into rom.setrhs (need inverse though)
obj.rhs = obj.rhs - rom.bu*(u(2:end,:)*obj.betas(2:end,obj.ito))/obj.dt;
end

function [next] = advance(obj, rom)
if isempty(obj.hfac)
h=rom.bu*obj.betas(1,obj.ito)/obj.dt+rom.mu*rom.au;
obj.hfac=chol(h);
end
next = [1,(obj.hfac\(obj.hfac'\obj.rhs))'];
if obj.ito <= 2
obj.hfac = [];
end
end

end

methods (Static)
function a = shift(a,b,n)
for i=n:-1:2
a(:,i)=a(:,i-1);
end
a(:,1)=b;
end

function rom = collect_statistics(rom, u_new)
rom.ua = rom.ua + u_new';
rom.u2a = rom.u2a + u_new'*u_new;
end
end
end
71 changes: 71 additions & 0 deletions bin/matlab_driver/grom.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,71 @@
%#####################################
%
%# GROM class inherits from NekROM
%# v0.0.0
%
%# Ping-Hsuan Tsai
%# 2024-07-04
%
%#####################################

classdef grom < nekrom
% Class for GROM (Galerkin Based Reduced Order Model)
properties
ua
u2a
mu
Re
end

methods
% Constructor
function obj = grom(path)
% Call the constructor of the superclass (NekROM)
obj@nekrom(path);
end

function classname = str(obj)
classname = upper(class(obj));
end

function obj = get_N_dim_ops(obj, nb)
% Get N dimensional operators and vectors
% : nb: number of modes
% : returns obj: NekROM object with N dimensional operators and vectors
obj = get_N_dim_ops@nekrom(obj, nb);
end


function obj = initialize_vars(obj, Re)
% Initialize other properties as needed
% : returns obj: NeROM object with initialized variables
obj.ua = zeros(obj.nb+1, 1); % initialize ROM averaged velocity coefficients
obj.u2a = zeros(obj.nb+1, obj.nb+1); % initialize ROM averaged velocity squared coefficients
obj.Re = Re;
obj.mu = 1./obj.Re;
end

function [rhs] = setrhs(obj, u, method)
% Set right hand side of the G-ROM
% : u: ROM velocity coefficients
% : returns rhs: right hand side of the G-ROM
if method == "semi-implicit"
rhs = -reshape(obj.cu*u(:,1),obj.nb,obj.nb+1)*u(:,1); % nonlinear term
rhs = rhs-obj.mu*obj.au0(2:end,1); % viscous term of the zeroth mode
end
end

function dump_rom_statistics(rom, nsteps)
ua = rom.ua/(nsteps);
u2a = rom.u2a/(nsteps);

fileID = fopen("./ua_nsteps"+nsteps+"N_"+rom.nb,'w');
fprintf(fileID,"%24.15e\n",ua);
fclose(fileID);

fileID = fopen("./u2a_nsteps"+nsteps+"N_"+rom.nb,'w');
fprintf(fileID,"%24.15e\n",u2a);
fclose(fileID);
end
end
end
73 changes: 73 additions & 0 deletions bin/matlab_driver/main.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,73 @@
%#####################################
%
%# Projection Based ROM driver in Matlab for NekROM
%# v0.0.0
%
%# Currently only support Galerkin ROM
%
%# Ping-Hsuan Tsai
%# 2024-07-04
%
%#####################################

clear all; close all;

% Parameters users need to specify
Re = 204; % Reynolds number
nb = 20; % number of modes
dt = 0.001; % Simultaion time step size
nsteps = 20000; % Simulation total time steps
iostep = 10; % Every #steps to store ROM quantities

rom = grom('../ops/'); % Create an instance of the grom class and load NekROM operators

rom = rom.get_N_dim_ops(nb); % Extract nb-dimensional operators and vectors
rom = rom.initialize_vars(Re); % Initialize variables

timestepper = BDFk_EXTk(nsteps, dt, iostep, nb);
timestepper = timestepper.setup();

ucoef = zeros(rom.nb+1, (timestepper.nsteps/timestepper.iostep));
% Solving ROM using BDFk/EXTk time stepping scheme
timestepper.u(:, 1) = rom.u0;
for istep=1:timestepper.nsteps
timestepper.ito=min(istep, 3);

timestepper = timestepper.setrhs(rom, timestepper.u); % Compute the RHS in BDFk/EXTk scheme
[u_new] = timestepper.advance(rom); % Solve for solution at the next time step
timestepper.u = timestepper.shift(timestepper.u, u_new, 3);

rom = timestepper.collect_statistics(rom, u_new); % Compute mean coefficient and mean squared coefficient

if (mod(istep, timestepper.iostep) == 0)
ucoef(:, istep/timestepper.iostep) = timestepper.u(:, 1);
end

end

% Plot the first five modes behavior in time

% Setup time stamp for ROM
% TODO: Make it cleaner
t_rom = linspace(timestepper.dt, timestepper.T_final, timestepper.nsteps/timestepper.iostep);
for i=2:min(rom.nb+1, 6)
figure(1)
plot(t_rom, ucoef(i, :), 'r'); hold on
legend(rom.str(), 'FontSize', 14);
xlabel('$t$', 'Interpreter', 'latex', 'FontSize', 14);
ylabel(['$u_', num2str(i), '$'], 'Interpreter', 'latex', 'FontSize', 14);
title(['Mode ', num2str(i), ' behavior of ', num2str(rom.nb), '-dimensional ', rom.str()], 'Interpreter', 'latex', 'FontSize', 14);
saveas(gcf, sprintf('u%d.png', i))
close(1)
end

dump_rom_statistics(rom, timestepper.nsteps);

figure(1)
plot(rom.uas(1:rom.nb+1), 'k-o'); hold on
plot(rom.ua/timestepper.nsteps, 'r-x'); hold off
legend(rom.str(), 'FOM', 'FontSize', 14);
xlabel('Mode $i$', 'Interpreter', 'latex', 'FontSize', 14);
ylabel('$u_i$', 'Interpreter', 'latex', 'FontSize', 14);
title(['Averaged coefficients', ' of ', num2str(rom.nb), '-dimensional ', rom.str()], 'Interpreter', 'latex', 'FontSize', 14);
saveas(gcf, sprintf('ua_compare_N%d.png', rom.nb))
90 changes: 90 additions & 0 deletions bin/matlab_driver/nekrom.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,90 @@
%#####################################
%
%# NekROM class
%# v0.0.0
%
%# Ping-Hsuan Tsai
%# 2024-07-04
%
%#####################################

classdef nekrom
% Class for NekROM
properties
mb % Total number of modes generated by NekROM
ns % Number of snapshots
aufull % Stiffness matrix of size (mb+1 x mb+1)
bufull % Mass matrix of size (mb+1 x mb+1)
cufull % Advection tensor of size (mb x mb+1 x mb+1)
u0full % Initial condition of size (mb+1)
ukfull % Snapshot projection matrix of size (mb+1 x ns)
uas % Averaged velocity coefficients of the snapshots (mb+1)
nb % Number of modes
au0 % Stiffness matrix of size (nb+1 x nb+1)
bu0 % Mass matrix of size (nb+1 x nb+1)
cu % Advection tensor of size (nb*(nb+1) x nb+1)
u0 % Initial condition of size (nb+1)
au % Stiffness matrix of size (nb x nb)
bu % Mass matrix of size (nb x nb)
uk % Snapshot projection matrix of size (nb+1 x ns)
ukmin % Minimum value of uk for each row
ukmax % Maximum value of uk for each row
end

methods
% Constructor
function obj = nekrom(path)
if nargin == 0
disp('Path is required');
return;
end
obj = obj.load_nekrom_ops(path);
end

function obj = load_nekrom_ops(obj, path)
% Load NekROM operators and vectors
% : path: path to the NekROM operators and vectors
% : returns obj: NekROM object with operators and vectors loaded

fprintf('Loading NekROM operators and vectors under path: %s\n', path);
% Load total number of modes generated by NekROM
obj.mb = dlmread(fullfile(path + "nb"));
tt = dlmread(path + "au");
obj.aufull = reshape(dlmread(path + "au"), obj.mb+1, obj.mb+1); % Load stiffness matrix
obj.bufull = reshape(dlmread(path + "bu"), obj.mb+1, obj.mb+1); % Load mass matrix
obj.cufull = reshape(dlmread(path + "cu"), obj.mb, obj.mb+1, obj.mb+1); % Load advection tensor
obj.u0full = dlmread(path + "u0"); % Load initial condition
obj.ukfull = reshape(dlmread(path + "uk"), obj.mb+1, []);
obj.ns = dlmread(path + "ns"); % load number of snapshots
obj.uas = dlmread(path + "uas"); % load averaged velocity coefficients of the snapshots
if size(obj.ukfull, 2) ~= obj.ns
error('Number of columns in obj.ukfull is not equal to ns.');
end
fprintf("Done loading NekROM operators ... \n");
end

function obj = get_N_dim_ops(obj, nb)
% Get N dimensional operators and vectors
% : nb: number of modes
% : returns obj: NekROM object with N dimensional operators and vectors

fprintf('Get N dimensional operators and vectors\n');

obj.nb = nb;
index = 1:obj.nb+1; % index including zeroth mode
index1 = 1:obj.nb; % index for the first dimension of the advection operator

obj.au0 = obj.aufull(index, index); % get stiffness matrix of size (nb+1 x nb+1)
obj.bu0 = obj.bufull(index, index); % get mass matrix of size (nb+1 x nb+1)
obj.au = obj.au0(2:end, 2:end); % get stiffness matrix of size (nb x nb)
obj.bu = obj.bu0(2:end, 2:end); % get mass matrix of size (nb x nb)
obj.cu = obj.cufull(index1, index, index); % get advection tensor of size (nb x (nb+1) x nb+1)
obj.cu = reshape(obj.cu, obj.nb*(obj.nb+1), obj.nb+1); % reshape advection tensor to size (nb*(nb+1) x nb+1)
obj.u0 = obj.u0full(index); % get initial condition of size (nb+1)
obj.uk = obj.ukfull(index, 1:obj.ns); % get snapshot projection matrix of size (nb+1 x ns)
obj.ukmin = min(obj.uk, [], 2); % get minimum value of uk for each row
obj.ukmax = max(obj.uk, [], 2); % get maximum value of uk for each row
end

end
end

0 comments on commit e7d9011

Please sign in to comment.