Contents

function LQR_DAE2(problem, level, re, istest)
% Computes a stabilizing feedback by implicitly solving the generalized
% Riccati equation 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.

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

Set operations

opts = struct();
[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

switch lower(problem)
    case 'stokes'
        nin = 5;
        nout = 5;
        nx = 10;
        ny = 10;
        [eqn.E_, eqn.A_, eqn.B, eqn.C] = ...
            stokes_ind2(nin, nout, nx, ny, opts);
        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.B(1:n_ode, :);
        eqn.C = eqn.C(:, 1:n_ode);
    case 'nse'
        [eqn, K0, ~] = mess_get_NSE(re, level);
        opts.nm.K0 = K0;
        opts.radi.K0 = K0;
    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 ...

First we run the Newton-ADI Method

opts.norm = 2;

% ADI tolerances and maximum iteration number
opts.adi.maxiter       = 300;
opts.adi.res_tol       = 1e-12;
opts.adi.rel_diff_tol  = 1e-16;
opts.adi.info          = 1;
opts.adi.LDL_T         = false;
eqn.type = 'T';
n = size(eqn.A_, 1);
opts.shifts.num_desired = 5; % *nout;
opts.shifts.num_Ritz    = 50;
opts.shifts.num_hRitz   = 25;
opts.shifts.method      = 'projection';
opts.shifts.b0          = ones(n, 1);

Newton tolerances and maximum iteration number

opts.nm.maxiter          = 20;
opts.nm.res_tol          = 1e-10;
opts.nm.rel_diff_tol     = 1e-16;
opts.nm.info             = 1;
opts.nm.projection.freq  = 0;
opts.nm.projection.ortho = true;
% in case you want to e.g. specify the factored Newton solver for
% the projected equations uncomment the following
% opts.nm.projection.meth = 'care_nwt_fac';
opts.nm.res              = struct('maxiter', 10, ...
                                  'tol', 1e-6, ...
                                  'orth', 0);
opts.nm.linesearch       = true;
opts.nm.inexact          = 'superlinear';
opts.nm.tau              = 0.1;
opts.nm.accumulateRes    = true;

use low-rank Newton-Kleinman-ADI

t_mess_lrnm = tic;
outnm = mess_lrnm(eqn, opts, oper);
t_elapsed1 = toc(t_mess_lrnm);
mess_fprintf(opts, 'mess_lrnm took %6.2f seconds \n\n', t_elapsed1);
if not(istest)
    figure(1);
    semilogy(outnm.res, 'LineWidth', 3);
    title('0 = C^T C + A^T X E + E^T X A - E^T X BB^T X E');
    xlabel('number of newton iterations');
    ylabel('normalized residual norm');
    pause(1);
end
[mZ, nZ] = size(outnm.Z);
mess_fprintf(opts, 'size outnm.Z: %d x %d\n\n', mZ, nZ);
ADI step:    1 normalized residual: 4.749604e-01 relative change in Z: 1.000000e+00  
		 normalized outer residual: 4.748502e-01  
ADI step:    2 normalized residual: 2.313480e-01 relative change in Z: 6.115557e-01  
		 normalized outer residual: 2.309720e-01  
ADI step:    3 normalized residual: 1.160895e-01 relative change in Z: 4.086903e-01  
		 normalized outer residual: 1.154641e-01  
ADI step:    4 normalized residual: 3.274932e-02 relative change in Z: 3.384842e-01  
		 normalized outer residual: 3.186747e-02  
ADI step:    5 normalized residual: 4.408328e-03 relative change in Z: 2.031583e-01  
		 normalized outer residual: 3.477590e-03  
  
 outer tolerance: 1.188528e-02  
  
 NM step:    1 normalized residual: 	3.477590e-03  
               relative change in K: 	1.000000e+00  
               number of ADI steps: 	5   
  
ADI step:    1 normalized residual: 4.751369e-01 relative change in Z: 1.000000e+00  
		 normalized outer residual: 4.746418e-01  
ADI step:    2 normalized residual: 2.310230e-01 relative change in Z: 6.109172e-01  
		 normalized outer residual: 2.308408e-01  
ADI step:    3 normalized residual: 1.156490e-01 relative change in Z: 4.073025e-01  
		 normalized outer residual: 1.155924e-01  
ADI step:    4 normalized residual: 3.274385e-02 relative change in Z: 3.354002e-01  
		 normalized outer residual: 3.273882e-02  
ADI step:    5 normalized residual: 4.365684e-03 relative change in Z: 2.003457e-01  
		 normalized outer residual: 4.365677e-03  
ADI step:    6 normalized residual: 1.459719e-03 relative change in Z: 6.921313e-02  
		 normalized outer residual: 1.459678e-03  
ADI step:    7 normalized residual: 9.263129e-05 relative change in Z: 5.492719e-02  
		 normalized outer residual: 9.253362e-05  
  
 outer tolerance: 3.863989e-04  
  
 NM step:    2 normalized residual: 	9.253362e-05  
               relative change in K: 	1.322216e-02  
               number of ADI steps: 	7   
  
ADI step:    1 normalized residual: 5.908896e-01 relative change in Z: 1.000000e+00  
		 normalized outer residual: 5.902469e-01  
ADI step:    2 normalized residual: 2.801793e-01 relative change in Z: 6.987959e-01  
		 normalized outer residual: 2.799246e-01  
ADI step:    3 normalized residual: 7.784030e-02 relative change in Z: 5.237500e-01  
		 normalized outer residual: 7.780849e-02  
ADI step:    4 normalized residual: 6.236314e-04 relative change in Z: 3.257720e-01  
		 normalized outer residual: 6.234901e-04  
ADI step:    5 normalized residual: 5.863780e-05 relative change in Z: 3.252883e-02  
		 normalized outer residual: 5.863776e-05  
ADI step:    6 normalized residual: 2.832581e-05 relative change in Z: 5.135530e-03  
		 normalized outer residual: 2.832580e-05  
ADI step:    7 normalized residual: 2.435341e-05 relative change in Z: 4.159018e-03  
		 normalized outer residual: 2.435341e-05  
ADI step:    8 normalized residual: 2.563225e-07 relative change in Z: 1.651385e-03  
		 normalized outer residual: 2.563217e-07  
  
 outer tolerance: 3.304772e-06  
  
 NM step:    3 normalized residual: 	2.563217e-07  
               relative change in K: 	1.145647e-04  
               number of ADI steps: 	8   
  
ADI step:    1 normalized residual: 5.952306e-01 relative change in Z: 1.000000e+00  
		 normalized outer residual: 5.945829e-01  
ADI step:    2 normalized residual: 3.403191e-01 relative change in Z: 6.609782e-01  
		 normalized outer residual: 3.399864e-01  
ADI step:    3 normalized residual: 1.078645e-01 relative change in Z: 5.650374e-01  
		 normalized outer residual: 1.078086e-01  
ADI step:    4 normalized residual: 2.082073e-02 relative change in Z: 3.428347e-01  
		 normalized outer residual: 2.081777e-02  
ADI step:    5 normalized residual: 1.759708e-04 relative change in Z: 1.767522e-01  
		 normalized outer residual: 1.759642e-04  
ADI step:    6 normalized residual: 8.406515e-05 relative change in Z: 1.521109e-02  
		 normalized outer residual: 8.406341e-05  
ADI step:    7 normalized residual: 9.211192e-06 relative change in Z: 1.267083e-02  
		 normalized outer residual: 9.211192e-06  
ADI step:    8 normalized residual: 3.101854e-06 relative change in Z: 1.302219e-03  
		 normalized outer residual: 3.101854e-06  
ADI step:    9 normalized residual: 2.053260e-08 relative change in Z: 1.875062e-03  
		 normalized outer residual: 2.053260e-08  
ADI step:   10 normalized residual: 6.170173e-09 relative change in Z: 8.933790e-05  
		 normalized outer residual: 6.170173e-09  
ADI step:   11 normalized residual: 1.864738e-09 relative change in Z: 6.988233e-05  
		 normalized outer residual: 1.864738e-09  
ADI step:   12 normalized residual: 1.354432e-09 relative change in Z: 2.042626e-05  
		 normalized outer residual: 1.354432e-09  
ADI step:   13 normalized residual: 9.506907e-11 relative change in Z: 4.068722e-05  
		 normalized outer residual: 9.506907e-11  
  
 outer tolerance: 3.943411e-09  
  
 NM step:    4 normalized residual: 	9.506907e-11  
               relative change in K: 	1.797487e-07  
               number of ADI steps: 	13   
  
mess_lrnm took   0.37 seconds   
  
size outnm.Z: 180 x 54  
  

Lets try the RADI method and compare

opts.norm               = 2;
% RADI-MESS settings
opts.shifts.history     = opts.shifts.num_desired * size(eqn.C, 1);
opts.shifts.num_desired = opts.shifts.num_desired;

% choose either of the three shift methods, here
opts.shifts.method = 'gen-ham-opti';
%     opts.shifts.method = 'heur';
%     opts.shifts.method = 'projection';

opts.shifts.naive_update_mode = false; % .. Suggest false
% (smart update is faster;
%  convergence is the same).
opts.shifts.info              = 0;
opts.radi.compute_sol_fac     = true; % Turned on for numerical stability reasons.
opts.radi.get_ZZt             = false;
opts.radi.maxiter             = opts.adi.maxiter;
opts.radi.res_tol             = opts.nm.res_tol;
opts.radi.rel_diff_tol        = 0;
opts.radi.info                = 1;

t_mess_lrradi = tic;
outradi = mess_lrradi(eqn, opts, oper);
t_elapsed2 = toc(t_mess_lrradi);
mess_fprintf(opts, 'mess_lrradi took %6.2f seconds \n', t_elapsed2);

if not(istest)
    figure();
    semilogy(outradi.res, 'LineWidth', 3);
    title('0 = C^TC + A^T X E + E^T X A - E^T X BB^T X E');
    xlabel('number of iterations');
    ylabel('normalized residual norm');
end
RADI step:    1 normalized residual: 1.190239e-01  
RADI step:    2 normalized residual: 2.025592e-02  
RADI step:    3 normalized residual: 5.326059e-03  
RADI step:    4 normalized residual: 4.009789e-03  
RADI step:    5 normalized residual: 1.105483e-04  
RADI step:    6 normalized residual: 9.456828e-05  
RADI step:    7 normalized residual: 4.772424e-06  
RADI step:    8 normalized residual: 3.409402e-07  
RADI step:    9 normalized residual: 2.053645e-07  
RADI step:   10 normalized residual: 1.257609e-08  
RADI step:   11 normalized residual: 7.092632e-09  
RADI step:   12 normalized residual: 1.209978e-09  
RADI step:   13 normalized residual: 5.526603e-11  
mess_lrradi took   0.14 seconds   

compare

if istest
    if min(outnm.res) >= opts.nm.res_tol
        mess_err(opts, 'TEST:accuracy', ...
                 'unexpectedly inaccurate result');
    end
    if min(outradi.res) >= opts.radi.res_tol
        mess_err(opts, 'TEST:accuracy', ...
                 'unexpectedly inaccurate result');
    end
else
    figure();
    ls_nm = [outnm.adi.niter];
    ls_radi = 1:outradi.niter;

    semilogy(cumsum(ls_nm), outnm.res, 'k--', ...
             ls_radi, outradi.res, 'b-', ...
             'LineWidth', 3);
    title('0 = C^T C + A^T X E + E^T X A - E^T X BB^T X E');
    xlabel('number of solves with A+p*M');
    ylabel('normalized residual norm');
    legend('LR-NM', 'RADI');
end