Contents

function IRKA_mor_Stokes(istest)
% Computes a standard ROM by implicitly running IRKA for the equivalent
% projected system on the hidden manifold.
%
% Inputs:
% istest        flag to determine whether this demo runs as a CI test or
%               interactive demo
%               (optional, defaults to 0, i.e. interactive demo)
%

%
% This file is part of the M-M.E.S.S. project
% (http://www.mpi-magdeburg.mpg.de/projects/mess).
% Copyright © 2009-2022 Jens Saak, Martin Koehler, Peter Benner and others.
% All rights reserved.
% License: BSD 2-Clause License (see COPYING)
%

IRKA tolerance and maximum iteration number

opts.irka.maxiter = 150;
opts.irka.r = 30;
opts.irka.flipeig = 0;
opts.irka.h2_tol = 1e-6;

opts.irka.init = 'logspace';

if nargin < 1, istest = 0; end

if istest
    opts.irka.info = 1;
else
    opts.irka.info = 2;
end

oper = operatormanager('dae_2');

Problem data

nin = 5;
nout = 5;
nx = 10;
ny = 10;
[eqn.E_, eqn.A_, eqn.Borig, eqn.Corig] = ...
    stokes_ind2(nin, nout, nx, ny);
n = size(eqn.E_, 1);
eqn.haveE = 1;
st=trace(eqn.E_); % Stokes is FDM discretized, so so this is
                  % the dimension of the velocity space
eqn.st = st;
eqn.B = eqn.Borig(1:st, :);
eqn.C = eqn.Corig(:, 1:st);
degrees of freedom: 
------------------------------
 total         :   280
 velocity      :   180
 pressure      :   100
 n_finite      :    81
------------------------------
Generating FVM matrices...
 -> Laplacians...
 -> gradient and divergence operator...
 -> B1, C1 and v1 velocity nodes...
 -> B2, C2 and v2 velocity nodes...
Setting up system matrices ...

Compute reduced system matrices

t_mess_tangential_irka = tic;
[ROM.E, ROM.A, ROM.B, ROM.C, ~, ~, ~, ~, W] = ...
    mess_tangential_irka(eqn, opts, oper);
t_elapsed1 = toc(t_mess_tangential_irka);
fprintf(1,'mess_tangential_irka took %6.2f seconds \n', t_elapsed1);
Warning: IRKA step 1 : 1 non-stable reduced eigenvalues detected.
 
IRKA step   1, rel. chg. shifts = 1.680119e+00 , rel. H2-norm chg. ROM = 1.000000e+00
Warning: IRKA step 2 : 1 non-stable reduced eigenvalues detected.
 
IRKA step   2, rel. chg. shifts = 1.560014e+00 , rel. H2-norm chg. ROM = Inf
IRKA step   3, rel. chg. shifts = 1.560738e+00 , rel. H2-norm chg. ROM = Inf
IRKA step   4, rel. chg. shifts = 1.210988e+00 , rel. H2-norm chg. ROM = 2.284273e-06
IRKA step   5, rel. chg. shifts = 8.204189e-01 , rel. H2-norm chg. ROM = 1.319332e-06
IRKA step   6, rel. chg. shifts = 1.210158e+00 , rel. H2-norm chg. ROM = 4.763188e-07
IRKA terminated due to relative change of ROMs criterion.

Computing TFMs of original and reduced order systems and MOR errors
 Step  10 / 100 Step  20 / 100 Step  30 / 100 Step  40 / 100 Step  50 / 100 Step  60 / 100 Step  70 / 100 Step  80 / 100 Step  90 / 100 Step 100 / 100

mess_tangential_irka took   1.39 seconds 
t_eval_ROM = tic;

Evaluate the ROM quality

while the Gramians are computed exploiting the DAE structure, due to the construction of the function handles we can not do so for the transfer function. Therefore we need to extend the matrices B and C and call the 'default' usfs for unstructured computation:

eqn.B = eqn.Borig;
eqn.C = eqn.Corig;
oper = operatormanager('default');

if istest
    opts.sigma.info = 0;
else
    opts.sigma.info = 2;
end

opts.sigma.fmin = -3;
opts.sigma.fmax = 4;

out = mess_sigma_plot(eqn, opts, oper, ROM);

t_elapsed2 = toc(t_eval_ROM);
fprintf(1,'evaluation of ROM matrices took %6.2f seconds \n' , t_elapsed2);
Computing TFMs of original and reduced order systems and MOR errors
 Step  10 / 100 Step  20 / 100 Step  30 / 100 Step  40 / 100 Step  50 / 100 Step  60 / 100 Step  70 / 100 Step  80 / 100 Step  90 / 100 Step 100 / 100

evaluation of ROM matrices took   0.57 seconds 
maerr = max(abs(out.err));
mrerr = max(abs(out.relerr));
if istest && (maerr>1e-6 || mrerr>1e-4)
    error('MESS:TEST:accuracy',['unexpectedly inaccurate result.\n' ...
                        'max. abs err: %e (allowed 1e-6)\n' ...
                        'max rel err: %e (allowed 1e-4)'], maerr, mrerr);
end
problem = 'Stokes';
fprintf(['\nComputing open loop step response of original and ', ...
    'reduced-order systems and time domain MOR errors\n']);
open_step(eqn, ROM.A, ROM.B, ROM.C, problem, istest);
Computing open loop step response of original and reduced-order systems and time domain MOR errors
 Implicit Euler step 500 / 5000 Implicit Euler step 1000 / 5000 Implicit Euler step 1500 / 5000 Implicit Euler step 2000 / 5000 Implicit Euler step 2500 / 5000 Implicit Euler step 3000 / 5000 Implicit Euler step 3500 / 5000 Implicit Euler step 4000 / 5000 Implicit Euler step 4500 / 5000 Implicit Euler step 5000 / 5000

implicit euler took   0.25 seconds 
fprintf('\nComputing ROM based feedback\n');
if exist('care', 'file')
    [~, ~, Kr] = care(ROM.A, ROM.B, ROM.C'*ROM.C, eye(size(ROM.B, 2)));
else
    Y = care_nwt_fac([], ROM.A, ROM.B, ROM.C, 1e-12, 50);
    Kr = (Y * ROM.B)' * Y;
end
K = [Kr * W' * eqn.E_(1:st, 1:st), zeros(size(Kr, 1), n-st)];
Computing ROM based feedback
fprintf(['\nComputing closed loop step response of original and ', ...
    'reduced-order systems and time domain MOR errors\n']);
closed_step(eqn, ROM.A, ROM.B, ROM.C, problem, K, Kr, istest);
Computing closed loop step response of original and reduced-order systems and time domain MOR errors
 Implicit Euler step 500 / 5000 Implicit Euler step 1000 / 5000 Implicit Euler step 1500 / 5000 Implicit Euler step 2000 / 5000 Implicit Euler step 2500 / 5000 Implicit Euler step 3000 / 5000 Implicit Euler step 3500 / 5000 Implicit Euler step 4000 / 5000 Implicit Euler step 4500 / 5000 Implicit Euler step 5000 / 5000

implicit euler took   0.41 seconds