Contents

function bt_mor_DAE2(problem, level, re, istest)
% Computes a standard ROM by implicitly solving the generalized Lyapunov
% equations for the equivalent projected system on the hidden manifold.
%
% Inputs:
% problem       either 'Stokes' or 'NSE' to choose the Stokes demo or the
%               linearized Navier-Stokes-Equation.
%               (optional, defaults to 'Stokes')
%
% level         discretization level 1 through 5
%               (optional, only used in 'NSE' case, default: 1)
%
% re            Reynolds number 300, 400, or 500
%               (optional, only used in 'NSE' case, default: 500)
%
% istest        flag to determine whether this demo runs as a CI test or
%               interactive demo
%               (optional, defaults to 0, i.e. interactive demo)
%
% Note that the 'NSE' option requires additional data available in a
% separate 270MB archive and at least the 5th discretization level needs a
% considerable amount of main memory installed in your machine.
%
% See:
% P. Benner, J. Saak, M. M. Uddin, Balancing based model reduction for
% structured index-2 unstable descriptor systems with application to flow
% control, Numerical Algebra, Control and Optimization 6 (1) (2016) 1–20.
% https://doi.org/10.3934/naco.2016.6.1

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

ADI tolerance and maximum iteration number

opts.adi.maxiter      = 350;
opts.adi.res_tol      = sqrt(eps);
opts.adi.rel_diff_tol = 1e-16;
opts.adi.info         = 1;
opts.shifts.info      = 1;
opts.norm             = 'fro';

[oper, opts] = operatormanager(opts, 'dae_2');

Problem data

if nargin < 1
    problem = 'stokes';
end
if nargin < 2
    level   = 1;
end
if nargin < 3
    re      = 500;
end
if nargin < 4
    istest  = false;
end

problem = lower(problem);

switch problem
    case 'stokes'
        nin  = 5;
        nout = 5;
        nx   = 10;
        ny   = 10;
        [eqn.E_, eqn.A_, eqn.Borig, eqn.Corig] = ...
            stokes_ind2(nin, nout, nx, ny, opts);
        n = size(eqn.E_, 1);
        eqn.haveE = true;
        n_ode = trace(eqn.E_); % Stokes is FDM discretized, so so this is
        % the dimension of the velocity space
        eqn.manifold_dim = n_ode;
        eqn.B = eqn.Borig(1:n_ode, :);
        eqn.C = eqn.Corig(:, 1:n_ode);
    case 'nse'
        [eqn, K0primal, K0dual] = mess_get_NSE(re, level);
        n_ode = eqn.manifold_dim;
        n = size(eqn.E_, 1);
    otherwise
        mess_err(opts, ...
                 'illegal_input', ...
                 'input ''problem'' must be either ''NSE'' or ''Stokes''');
end
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 ...
eqn.type       = 'N';
% Activate stabilizing (Bernoulli) feedback
if strcmp(problem, 'nse')
    eqn.V      = -K0primal';
    eqn.U      = eqn.B;
    eqn.haveUV = true;
end

opts.shifts.num_desired = 6;
opts.shifts.num_Ritz    = 40;
opts.shifts.num_hRitz   = 40;
opts.shifts.method      = 'projection';

opts.shifts.b0 = ones(size(eqn.A_, 1), 1);

t_mess_lradi = tic;
outB = mess_lradi(eqn, opts, oper);
t_elapsed1 = toc(t_mess_lradi);
mess_fprintf(opts, ...
             'mess_lradi took %6.2f seconds \n', ...
             t_elapsed1);

if not(istest)
    figure('Name', 'Residual history (controllability)');
    semilogy(outB.res, 'LineWidth', 3);
    title('A X E^T + E X A^T = -BB^T');
    xlabel('number of iterations');
    ylabel('normalized residual norm');
    pause(1);
end
[mZ, nZ] = size(outB.Z);
mess_fprintf(opts, 'size outB.Z: %d x %d\n\n', mZ, nZ);
ADI Shifts:
-4.161327e+01  
-3.127846e+01  
-2.738932e+01  
-1.141206e+01  
-9.401353e+00  
ADI step:    1 normalized residual: 1.980092e-01 relative change in Z: 1.000000e+00  
ADI step:    2 normalized residual: 3.345840e-02 relative change in Z: 4.304568e-01  
ADI step:    3 normalized residual: 9.453481e-03 relative change in Z: 1.789458e-01  
ADI step:    4 normalized residual: 5.820955e-03 relative change in Z: 8.604064e-02  
ADI step:    5 normalized residual: 4.927155e-03 relative change in Z: 1.379911e-02  
updating shifts
p:
-5.670482e+01  
-9.789089e+00  
-2.815743e+02  
-2.421695e+01  
-1.370822e+02  
-1.559854e+01  
ADI step:    6 normalized residual: 1.928044e-03 relative change in Z: 2.030694e-02  
ADI step:    7 normalized residual: 1.662636e-03 relative change in Z: 5.589770e-03  
ADI step:    8 normalized residual: 2.532499e-05 relative change in Z: 1.279104e-02  
ADI step:    9 normalized residual: 1.472269e-05 relative change in Z: 1.389208e-03  
ADI step:   10 normalized residual: 1.106893e-06 relative change in Z: 1.458978e-03  
ADI step:   11 normalized residual: 9.053482e-07 relative change in Z: 1.388050e-04  
updating shifts
p:
-6.145754e+01  
-1.054438e+01  
-3.494876e+02  
-1.569795e+02  
-1.857658e+01  
-2.561629e+02  
ADI step:   12 normalized residual: 4.179204e-07 relative change in Z: 2.032810e-04  
ADI step:   13 normalized residual: 3.679514e-07 relative change in Z: 6.172848e-05  
ADI step:   14 normalized residual: 2.797198e-09 relative change in Z: 1.583981e-04  
mess_lradi took   0.08 seconds   
size outB.Z: 180 x 70  
  
eqn.type = 'T';
% Activate stabilizing (Bernoulli) feedback (for the dual system)
if strcmp(problem, 'nse')
    eqn.U      = -K0dual';
    eqn.V      = eqn.C';
    eqn.haveUV = true;
end

opts.shifts.num_desired = 6;
opts.shifts.num_Ritz    = 40;
opts.shifts.num_hRitz   = 40;
opts.shifts.method      = 'projection';

opts.shifts.b0 = ones(size(eqn.A_, 1), 1);

t_mess_lradi = tic;
outC = mess_lradi(eqn, opts, oper);
t_elapsed2 = toc(t_mess_lradi);
mess_fprintf(opts, ...
             'mess_lradi took %6.2f seconds \n', ...
             t_elapsed2);

if not(istest)
    figure('Name', 'Residual history (observability)');
    semilogy(outC.res, 'LineWidth', 3);
    title('A^T X E + E^T X A = -C^T C');
    xlabel('number of iterations');
    ylabel('normalized residual norm');
    pause(1);
end
[mZ, nZ] = size(outC.Z);
mess_fprintf(opts, 'size outC.Z: %d x %d\n\n', mZ, nZ);
ADI Shifts:
-7.551972e+01  
-4.318595e+01  
-4.015777e+01  
-2.155011e+01  
-9.216396e+00  
ADI step:    1 normalized residual: 2.795911e-01 relative change in Z: 1.000000e+00  
ADI step:    2 normalized residual: 5.476559e-02 relative change in Z: 5.518605e-01  
ADI step:    3 normalized residual: 1.202826e-02 relative change in Z: 2.434720e-01  
ADI step:    4 normalized residual: 3.557729e-03 relative change in Z: 1.189586e-01  
ADI step:    5 normalized residual: 2.902586e-03 relative change in Z: 2.993634e-02  
updating shifts
p:
-5.368712e+01  
-9.792074e+00  
-2.896480e+02  
-1.270587e+02  
-1.915558e+01  
-3.826816e+01  
ADI step:    6 normalized residual: 1.327858e-03 relative change in Z: 1.498369e-02  
ADI step:    7 normalized residual: 1.157930e-03 relative change in Z: 5.048697e-03  
ADI step:    8 normalized residual: 1.322849e-05 relative change in Z: 1.158868e-02  
ADI step:    9 normalized residual: 1.838442e-06 relative change in Z: 1.384702e-03  
ADI step:   10 normalized residual: 1.440246e-06 relative change in Z: 3.830421e-04  
ADI step:   11 normalized residual: 9.232355e-07 relative change in Z: 2.276257e-04  
updating shifts
p:
-5.377860e+01  
-3.770595e+02  
-9.791499e+00  
-1.346499e+02  
-1.907646e+01  
-2.519305e+02  
ADI step:   12 normalized residual: 5.009422e-07 relative change in Z: 1.967278e-04  
ADI step:   13 normalized residual: 3.235110e-09 relative change in Z: 1.993795e-04  
mess_lradi took   0.08 seconds   
size outC.Z: 180 x 65  
  

Compute reduced system matrices

Perform Square Root Method (SRM)

% BT tolerance and maximum order for the ROM
t_SRM             = tic;
opts.srm.tol      = 1e-5;
opts.srm.max_ord  = 250;

% SRM verbosity
if istest
    opts.srm.info = 1;
else
    opts.srm.info = 2;
end

% The actual SRM
[TL, TR, hsv] = mess_square_root_method(eqn, opts, oper, outB.Z, outC.Z);
reduced system order: 19  (max possible/allowed: 65/250)  
  
ROM.A = TL' * (eqn.A_(1:n_ode, 1:n_ode) * TR);
ROM.B = TL' * eqn.B(1:n_ode, :);
ROM.C = eqn.C(:, 1:n_ode) * TR;

t_elapsed3 = toc(t_SRM);
mess_fprintf(opts, ...
             'computation of reduced system matrices took %6.2f seconds \n', ...
             t_elapsed3);
computation of reduced system matrices took   0.03 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:

switch lower(problem)
    case 'stokes'
        eqn.B = eqn.Borig;
        eqn.C = eqn.Corig;
    case 'nse'
        n = size(eqn.A_, 1);
        eqn.B(n_ode + 1:n, :) = zeros(n - n_ode, size(eqn.B, 2));
        eqn.C(:, n_ode + 1:n) = zeros(size(eqn.C, 1), n - n_ode);
end
[oper, opts] = operatormanager(opts, 'default');

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

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

opts.tf_plot.type = 'sigma';

out = mess_tf_plot(eqn, opts, oper, ROM);
err = out.err;

t_elapsed4 = toc(t_eval_ROM);
mess_fprintf(opts, ...
             'evaluation of rom quality took %6.2f seconds \n', ...
             t_elapsed4);
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 quality took   0.44 seconds   
if istest
    if max(err) >= opts.srm.tol
        mess_err(opts, 'TEST:accuracy', 'unexpectedly inaccurate result');
    end
else
    figure('Name', 'Computed Hankel singular values');
    semilogy(hsv, 'LineWidth', 3);
    title('Computed Hankel singular values');
    xlabel('index');
    ylabel('magnitude');
end
mess_fprintf(opts, ...
             ['\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.14 seconds 
mess_fprintf(opts, '\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 * TL' * eqn.E_(1:n_ode, 1:n_ode), zeros(size(Kr, 1), n - n_ode)];
Computing ROM based feedback
mess_fprintf(opts, ...
             ['\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.16 seconds