Add Matlab driver for NekROM and add integration test to it

.github/workflows/integ_t_matlab.yaml

@@ -0,0 +1,30 @@
+name: MATLAB Driver Test
+
+on: [push]
+
+jobs:
+  test:
+    runs-on: ubuntu-latest
+
+    env:
+      MOR_DIR: /home/runner/work/NekROM/NekROM
+
+    steps:
+    - name: Checkout repository
+      uses: actions/checkout@v2
+
+    - name: Set up MATLAB
+      uses: matlab-actions/setup-matlab@v1
+      with:
+        release: R2023b  # Replace with your MATLAB release
+
+    - name: Install gdown
+      run: pip install gdown
+
+    - name: Install unzip
+      if: runner.os == 'Linux'
+      run: sudo apt-get install -y unzip
+
+    - name: Run MATLAB tests
+      run: |
+        $MOR_DIR/tests/matlab_test/test.sh

bin/matlab_driver/BDFk_EXTk.m

@@ -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

bin/matlab_driver/grom.m

@@ -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

bin/matlab_driver/main.m

@@ -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))

bin/matlab_driver/nekrom.m

@@ -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