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 program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, see <http://www.gnu.org/licenses/>.
%
% Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others
%               2009-2019
%

% 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 = 1;
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 = full(sum(diag(eqn.E_)));
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

tic;
[ROM.E, ROM.A, ROM.B, ROM.C, S, b, c, V, W] = ...
    mess_tangential_irka(eqn, opts, oper);
toc;
Warning: IRKA step 1 : 2 non-stable reduced eigenvalues detected.
 
IRKA step   1, rel. chg. shifts = 3.156664e+02 , rel. H2-norm chg. ROM = 1.000000e+00
Warning: IRKA step 2 : 1 non-stable reduced eigenvalues detected.
 
IRKA step   2, rel. chg. shifts = 4.188718e+02 , rel. H2-norm chg. ROM = Inf
Warning: IRKA step 3 : 1 non-stable reduced eigenvalues detected.
 
IRKA step   3, rel. chg. shifts = 5.524551e+02 , rel. H2-norm chg. ROM = Inf
IRKA step   4, rel. chg. shifts = 1.620365e+02 , rel. H2-norm chg. ROM = Inf
Warning: IRKA step 5 : 1 non-stable reduced eigenvalues detected.
 
IRKA step   5, rel. chg. shifts = 2.057972e+02 , rel. H2-norm chg. ROM = Inf
IRKA step   6, rel. chg. shifts = 2.710761e+02 , rel. H2-norm chg. ROM = 1.106528e-06
IRKA step   7, rel. chg. shifts = 3.164476e+03 , rel. H2-norm chg. ROM = 1.168920e-06
IRKA step   8, rel. chg. shifts = 3.481692e+03 , rel. H2-norm chg. ROM = 7.188294e-06
IRKA step   9, rel. chg. shifts = 3.179503e+02 , rel. H2-norm chg. ROM = 7.175401e-06
IRKA step  10, rel. chg. shifts = 5.023936e+02 , rel. H2-norm chg. ROM = 9.867234e-07
Elapsed time is 0.399755 seconds.
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. Therfore 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);

toc;
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

Elapsed time is 0.404786 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

Elapsed time is 0.127836 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

Elapsed time is 0.146299 seconds.