Contents
function bt_mor_DAE2(problem,lvl,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. (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. % % 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 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 % % 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 = 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.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); 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.B=mat.mat_v.B{lvl}; eqn.C=mat.mat_v.C{lvl}; eqn.st=mat.mat_mg.nv(lvl); st=eqn.st; eqn.haveE=1; n=size(eqn.E_,1); 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 ...
eqn.type='N'; eqn.G=eqn.B; if strcmp(problem,'NSE') && (re>200) eqn.V=-mat.mat_v.Feed_0{lvl}; eqn.U = eqn.B; eqn.haveUV=1; 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); opts.shifts.p=mess_para(eqn,opts,oper); tic; outB = mess_lradi(eqn,opts,oper); toc; if not(istest) figure(); semilogy(outB.res); title('AXM^T + MXA^T = -BB^T'); xlabel('number of iterations'); ylabel('normalized residual norm'); pause(1); end disp('size outB.Z:'); disp(size(outB.Z));
ADI Shifts: -41.6133 -31.2785 -27.3893 -11.4121 -9.4014 ADI step: 1 normalized residual: 1.035178e-01 relative change in Z: 1.000000e+00 ADI step: 2 normalized residual: 1.749181e-02 relative change in Z: 4.304568e-01 ADI step: 3 normalized residual: 4.942211e-03 relative change in Z: 1.789458e-01 ADI step: 4 normalized residual: 3.043153e-03 relative change in Z: 8.604064e-02 ADI step: 5 normalized residual: 2.575881e-03 relative change in Z: 1.379911e-02 updating shifts -56.7048 -9.7891 -281.5743 -24.2170 -137.0822 -15.5985 ADI step: 6 normalized residual: 1.007967e-03 relative change in Z: 2.030694e-02 ADI step: 7 normalized residual: 8.692139e-04 relative change in Z: 5.589770e-03 ADI step: 8 normalized residual: 1.323972e-05 relative change in Z: 1.279104e-02 ADI step: 9 normalized residual: 7.696915e-06 relative change in Z: 1.389208e-03 ADI step: 10 normalized residual: 5.786755e-07 relative change in Z: 1.458978e-03 ADI step: 11 normalized residual: 4.733094e-07 relative change in Z: 1.388050e-04 updating shifts -61.4575 -10.5444 -349.4876 -156.9795 -18.5766 -256.1629 ADI step: 12 normalized residual: 2.184857e-07 relative change in Z: 2.032810e-04 ADI step: 13 normalized residual: 1.923623e-07 relative change in Z: 6.172848e-05 ADI step: 14 normalized residual: 1.462355e-09 relative change in Z: 1.583981e-04 Elapsed time is 0.019494 seconds. size outB.Z: 180 70

eqn.type = 'T'; eqn.G=eqn.C'; if strcmp(problem,'NSE') && (re>200) if (re == 500) % the dataset stores the feedback of the adjoint system transposed % for Reynolds 500 eqn.U=-mat.mat_v.Feed_1{lvl}'; else eqn.U=-mat.mat_v.Feed_1{lvl}; end eqn.V = eqn.C'; eqn.haveUV=1; 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); opts.shifts.p=mess_para(eqn,opts,oper); tic; outC = mess_lradi(eqn, opts, oper); toc; if not(istest) figure(); semilogy(outC.res); title('A^TXM + M^TXA = -C^TC'); xlabel('number of iterations'); ylabel('normalized residual norm'); pause(1); end disp('size outC.Z:'); disp(size(outC.Z));
ADI Shifts: -75.5197 -43.1860 -40.1578 -21.5501 -9.2164 ADI step: 1 normalized residual: 1.208573e-01 relative change in Z: 1.000000e+00 ADI step: 2 normalized residual: 2.367322e-02 relative change in Z: 5.518605e-01 ADI step: 3 normalized residual: 5.199387e-03 relative change in Z: 2.434720e-01 ADI step: 4 normalized residual: 1.537880e-03 relative change in Z: 1.189586e-01 ADI step: 5 normalized residual: 1.254685e-03 relative change in Z: 2.993634e-02 updating shifts -53.6871 -9.7921 -289.6480 -127.0587 -19.1556 -38.2682 ADI step: 6 normalized residual: 5.739858e-04 relative change in Z: 1.498369e-02 ADI step: 7 normalized residual: 5.005320e-04 relative change in Z: 5.048697e-03 ADI step: 8 normalized residual: 5.718206e-06 relative change in Z: 1.158868e-02 ADI step: 9 normalized residual: 7.946930e-07 relative change in Z: 1.384702e-03 ADI step: 10 normalized residual: 6.225670e-07 relative change in Z: 3.830421e-04 ADI step: 11 normalized residual: 3.990819e-07 relative change in Z: 2.276257e-04 updating shifts -53.7786 -377.0595 -9.7915 -134.6499 -19.0765 -251.9305 ADI step: 12 normalized residual: 2.165395e-07 relative change in Z: 1.967278e-04 ADI step: 13 normalized residual: 1.398423e-09 relative change in Z: 1.993795e-04 Elapsed time is 0.016688 seconds. size outC.Z: 180 65

Compute reduced system matrices
Perform Square Root Method (SRM)
% BT tolerance and maximum order for the ROM 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:st,1:st)*TR); ROM.B = TL'*eqn.B(1:st,:); ROM.C = eqn.C(:,1:st)*TR; toc;
Elapsed time is 0.008581 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:
switch lower(problem) case 'stokes' eqn.B=eqn.Borig; eqn.C=eqn.Corig; case 'nse' n = size(eqn.A_,1); eqn.B(st+1:n,:) = zeros(n-st,size(eqn.B,2)); eqn.C(:,st+1:n) = zeros(size(eqn.C,1),n-st); end 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); err = out.err; 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.372866 seconds.


if istest if max(err)>=opts.srm.tol, error('MESS:TEST:accuracy','unexpectedly inaccurate result'); end else figure; semilogy(hsv); title('Computed Hankel singular values'); xlabel('index'); ylabel('magnitude'); end

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.138973 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*TL'*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.138561 seconds.


