Contents

function LQR_DAE2(problem,lvl,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. (required)
%
% lvl           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 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-2020
%

Set operations

oper = operatormanager('dae_2');

Problem data

if nargin<1, problem='stokes'; end
if nargin<2, lvl=1; end
if nargin<3, re=500; end
if nargin<4, istest=0; 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);
    eqn.haveE=1;
    st=full(sum(diag(eqn.E_)));
    eqn.st=st;
    eqn.B=eqn.B(1:st,:);
    eqn.C=eqn.C(:,1:st);
  case 'nse'
    try
      load(sprintf('%s/../models/NSE/mat_nse_re_%d',...
                    fileparts(mfilename('fullpath')),re),'mat');
    catch
      error(['The files mat_nse_re_300.mat, mat_nse_re_400.mat and ', ...
          'mat_nse_re_500.mat are available for dowload in a ', ...
          'separate archive (270MB each). Please fetch them from the ', ...
          'MESS download page and unpack them into the ', ...
          'DEMOS/models/NSE folder.']);
    end
    eqn.A_=mat.mat_v.fullA{lvl};
    eqn.E_=mat.mat_v.E{lvl};
    eqn.haveE=1;
    eqn.B=mat.mat_v.B{lvl};
    if re>200
      opts.nm.K0=mat.mat_v.Feed_0{lvl}';
      opts.radi.K0 = opts.nm.K0;
    end
    eqn.C=mat.mat_v.C{lvl};
    eqn.st=mat.mat_mg.nv(lvl);
    otherwise
      error('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=0;
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=1;
%opts.nm.projection.meth='care_nwt_fac';
opts.nm.res=struct('maxiter',10,'tol',1e-6,'orth',0);
opts.nm.linesearch = 1;
opts.nm.inexact = 'superlinear';
opts.nm.tau = 0.1;
opts.nm.accumulateRes = 1;

use low-rank Newton-Kleinman-ADI

tic;
outnm = mess_lrnm(eqn, opts, oper);
toc;

if not(istest)
    figure(1);
    disp(outnm.res);
    semilogy(outnm.res);
    title('0= C^TC + A^TXM + M^TXA -M^TXBB^TXM');
    xlabel('number of newton iterations');
    ylabel('normalized residual norm');
    pause(1);
end
disp('size outnm.Z:');
disp(size(outnm.Z));
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.459721e-03 relative change in Z: 6.921311e-02
		normalized outer residual: 1.459679e-03
ADI step:    7 normalized residual: 9.263121e-05 relative change in Z: 5.492722e-02
		normalized outer residual: 9.253355e-05

outer tolerance: 3.863989e-04

NM step:    2  normalized residual: 9.253355e-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.832583e-05 relative change in Z: 5.135531e-03
		normalized outer residual: 2.832582e-05
ADI step:    7 normalized residual: 2.435343e-05 relative change in Z: 4.159017e-03
		normalized outer residual: 2.435343e-05
ADI step:    8 normalized residual: 2.563350e-07 relative change in Z: 1.651384e-03
		normalized outer residual: 2.563342e-07

outer tolerance: 3.304769e-06

NM step:    3  normalized residual: 2.563342e-07
               relative change in K: 1.145646e-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.406516e-05 relative change in Z: 1.521109e-02
		normalized outer residual: 8.406342e-05
ADI step:    7 normalized residual: 9.211191e-06 relative change in Z: 1.267083e-02
		normalized outer residual: 9.211191e-06
ADI step:    8 normalized residual: 3.101855e-06 relative change in Z: 1.302219e-03
		normalized outer residual: 3.101855e-06
ADI step:    9 normalized residual: 2.053261e-08 relative change in Z: 1.875062e-03
		normalized outer residual: 2.053261e-08
ADI step:   10 normalized residual: 6.170187e-09 relative change in Z: 8.933795e-05
		normalized outer residual: 6.170187e-09
ADI step:   11 normalized residual: 1.864742e-09 relative change in Z: 6.988242e-05
		normalized outer residual: 1.864742e-09
ADI step:   12 normalized residual: 1.354434e-09 relative change in Z: 2.042627e-05
		normalized outer residual: 1.354434e-09
ADI step:   13 normalized residual: 9.506925e-11 relative change in Z: 4.068725e-05
		normalized outer residual: 9.506925e-11

outer tolerance: 3.943603e-09

NM step:    4  normalized residual: 9.506925e-11
               relative change in K: 1.797500e-07
               number of ADI steps: 13 

Elapsed time is 0.115253 seconds.
    0.0035
    0.0001
    0.0000
    0.0000

size outnm.Z:
   180    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 = 1; % Turned on for numerical stability reasons.
opts.radi.get_ZZt = 0;
opts.radi.maxiter = opts.adi.maxiter;
opts.radi.res_tol = opts.nm.res_tol;
opts.radi.rel_diff_tol = 0;
opts.radi.info = 1;


tic;
outradi = mess_lrradi(eqn, opts, oper);
toc;

if not(istest)
    figure();
    semilogy(outradi.res);
    title('0= C^TC + A^TXM + M^TXA -M^TXBB^TXM');
    xlabel('number of iterations');
    ylabel('normalized residual norm');
end
RADI step:    1 normalized residual: 7.035353e-02
RADI step:    2 normalized residual: 1.197302e-02
RADI step:    3 normalized residual: 3.148166e-03
RADI step:    4 normalized residual: 2.370135e-03
RADI step:    5 normalized residual: 6.534372e-05
RADI step:    6 normalized residual: 5.589811e-05
RADI step:    7 normalized residual: 2.820919e-06
RADI step:    8 normalized residual: 2.015254e-07
RADI step:    9 normalized residual: 1.213884e-07
RADI step:   10 normalized residual: 7.433566e-09
RADI step:   11 normalized residual: 4.192365e-09
RADI step:   12 normalized residual: 7.152029e-10
RADI step:   13 normalized residual: 5.434915e-10
RADI step:   14 normalized residual: 4.511120e-10
RADI step:   15 normalized residual: 2.662895e-10
RADI step:   16 normalized residual: 2.161544e-10
RADI step:   17 normalized residual: 1.216147e-10
RADI step:   18 normalized residual: 6.488593e-11
Elapsed time is 0.083617 seconds.

compare

if istest
    if min(outnm.res)>=opts.nm.res_tol, error('MESS:TEST:accuracy','unexpectedly inaccurate result'); end
    if min(outradi.res)>=opts.radi.res_tol, error('MESS: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-');
    title('0= C^TC + A^TXM + M^TXA -M^TXBB^TXM');
    xlabel('number of solves with A+p*M');
    ylabel('normalized residual norm');
    legend('LR-NM','RADI');
end