lin_ds_adjoint_experiments.m

Go to the documentation of this file.

% script performing different experiments for linear dynamical
% system

step = 1; % initialize model, compute some trajectories
%step = 2; % generate reduced basis by POD
%step = 3; % comparison detailed and reduced simulation

params = [];
model = [];
%params.coarse_factor = 2; % fact 2 => ndofs = 16384, nt = 1024, 15s detailed
%params.coarse_factor = 1; % fact 1 => ndofs = 32768, nt = 2048, 60s detailed
params.coarse_factor = 1; 
%verbose(1);

% some filenames for storing data
trajectory_fnbase = 'advection_fv_output_vconst_trajectory';
detailed_data_fname = 'advection_fv_output_vconst_detailed_data';
rb_fname = 'advection_fv_output_vconst_rb';     

switch step
 case 1
  
  model = lin_ds_advection_vconst_fast_model(params);
  
  model_data = gen_model_data(model);

  mus = {[0,0,0],[1,1,1],[1,1,0.5],[0,0.5,1]};

  for m = 1:length(mus)
    disp(['computing trajectory number ',num2str(m)]);
    mu = mus{m};
    model = model.set_mu(model,mu);
    ds_sim_data = detailed_simulation(model,model_data);
    save(fullfile(rbmatlabtemp,[trajectory_fnbase,num2str(m)]),...
         'model','model_data','ds_sim_data');
  end;
  
 case 2
  
  disp('PCA of snapshots takes a while...');
  load(fullfile(rbmatlabtemp,[trajectory_fnbase,'1']),...
       'model','model_data','ds_sim_data');
  % smaller epsilon does result in more than 1 vector for mu=0,0,0
  epsilon = 1e-4;
  RB = PCA_fixspace(ds_sim_data.X,[],model_data.G,[],[],epsilon);
  %RB = orthonormalize_qr(ds_sim_data.X,model_data.G);
  % select correct orthogonal set
  K = RB' * model_data.G * RB;
  i = find(abs(diag(K)-1)>1e-2);
  i = min(i);
  disp(['found ',num2str(size(RB,2)),' basis vectors after traj. 1']);
  if isempty(i)
    i = size(RB,2)+1;
  end;
  %  keyboard;
  RB = RB(:,1:(i-1));
  disp(['reduced to ',num2str(size(RB,2)),' basis vectors.']);    
  %  keyboard;
  
  load(fullfile(rbmatlabtemp,[trajectory_fnbase,'2']),...
       'model','model_data','ds_sim_data');
  RB2 = PCA_fixspace(ds_sim_data.X,RB,model_data.G,[],[],epsilon);
  %RB = orthonormalize_qr(ds_sim_data.X,model_data.G);
  % select correct orthogonal set
  K = RB2' * model_data.G * RB2;
  i = find(abs(diag(K)-1)>1e-2);
  i = min(i);
  if isempty(i)
    i = size(RB2,2)+1;
  end;
  disp(['found ',num2str(size(RB2,2)),' basis vectors after traj. 2']);
  %  keyboard;
  RB2 = RB2(:,1:(i-1));
  disp(['reduced to ',num2str(size(RB2,2)),' basis vectors after traj. 2']);
  RB = [RB, RB2];
  %  RB = orthonormalize(RB);
  
  V = RB;
  W = model_data.G * V;  
  save(fullfile(rbmatlabtemp,rb_fname),'RB');
  
  model.RB_basis_filename = fullfile(rbmatlabtemp,rb_fname);
  model.RB_generation_mode = 'file_load'; 
  
  detailed_data = gen_detailed_data(model,model_data);
  
  model = model;
  save(fullfile(rbmatlabtemp,detailed_data_fname),...
       'model','detailed_data')

 case 3 % compare detailed and reduced simulation
    
  load(fullfile(rbmatlabtemp,detailed_data_fname),...
       'model','detailed_data')
  load(fullfile(rbmatlabtemp,[trajectory_fnbase,'2']),...
       'ds_sim_data');

  model = model;
  
  %perform single reduced simulation for given parameter
  reduced_data = gen_reduced_data(model,detailed_data);
  model.N = reduced_data.N;
  tic,
%  model.nt = model.nt/4;
%  disp('set nt to quarter!!');
  rb_sim_data = rb_simulation(model,reduced_data);
  t_red = toc;
  rb_sim_data = rb_reconstruction(model,detailed_data,rb_sim_data)
  disp(['time for reduced simulation: ',num2str(t_red)]);
  
  % plot of reduced trajectory
  lin_evol_model = model.base_model;
  lin_evol_model_data = gen_model_data(lin_evol_model);  
  
  params.no_lines = 1;
  params.show_colorbar = 1;
  params.axis_equal = 1;
  params.plot = lin_evol_model.plot;
  % plot single snapshot
  %lin_evol_model.plot(rb_sim_data.X(:,1),lin_evol_model_data.grid, ...
  %                   params);
  
  % plot sequence:
  params.title = 'reduced solution';
  plot_sequence(rb_sim_data.X,lin_evol_model_data.grid,params)  

  % plot sequence:
  params.title = 'error';
  plot_sequence((rb_sim_data.X-ds_sim_data.X),...
                lin_evol_model_data.grid,params)  

  % plot of reduced output
  figure,plot(rb_sim_data.time,[rb_sim_data.Y;ds_sim_data.Y]); 
  legend('reduced output','detailed output');

 case 4 % einfache Basis durch balanciertes Abschneiden erzeugen
     % A = ...
     % B = ...
     % k = 10;
            
     %[V,W] = balanced_truncation(A,B,C,D,k)
     
 case 5 % test der Basis anahand von balancierten System    
     
 otherwise 
  error('step unknown');
end;