advection_fv_output_morepas.m

Go to the documentation of this file.

% small script implementing a simple advection example
% for producing matrices to be used in the RB-DS framework
% Discretization with FV Functions
%
% mu_names  = {'cone_weight','vx_weight','vy_weight'};
% mu_ranges = [0,1]^3

% Bernard Haasdonk 25.8.2009

% cartesian grid

step = 1; % initialization of model and plot of model data
%step = 2; % detailed simulation of lin_evol
%step = 3; % conversion to DS primal model, DS detailed simulation 
          % and comparison with lin_evol
%step = 4 % reduced Basis generation by single trajectory and
          %comparison of detailed and reduced simulation
% step = 5 % figures of detailed simuations based on step
          % 4
%step = 6 % figures of reduced simuations based on step
          % 
%step = 7 % time measurements
          
switch step
  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %step = 1; % initialization of model and plot of model data
 case 1
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  
  model = advection_fv_output_model([]);
  model_data = gen_model_data(model,[]);
  plot(model_data.grid);
  %    inspect(model_data.grid)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %step = 2; % detailed simulation of lin_evol
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 case 2
  model = advection_fv_output_model([]);
  model_data = gen_model_data(model,[]);
  %    model.debug = 1;
  %    disp('model debugging turned on.')
  model.cone_weight = 0.5;
  model.vx_weight = 1;
  model.vy_weight = 1;
  sim_data = detailed_simulation(model,model_data,[]);
  plot_sim_data(model,model_data,sim_data);

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%step = 3; % conversion to DS primal model, DS detailed simulation 
        % and comparison with lin_evol
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 case 3 % initialization of DS model and lin_ds detailed simulation

        lin_evol_model = advection_fv_output_model([]);
  %lin_evol_model.debug = 1;
  lin_evol_model_data = gen_model_data(lin_evol_model,[]);
  
  % set some additional fields in model which will be copied into
  % ds model depending on whether constants are known or not:
  estimate_constants = 0;
  if estimate_constants
    %%% if constant estimaton has to be performed:
    lin_evol_model.estimate_lin_ds_nmu = 4;
    lin_evol_model.estimate_lin_ds_nX = 10;
    lin_evol_model.estimate_lin_ds_nt = 4;
  else
    %%% otherwise:
    % the following values are rough bounds, so could be 
    % non-rigorous, then larger search must be performed.
    lin_evol_model.state_bound_constant_C1 = 1;
    lin_evol_model.output_bound_constant_C2 = 1.2916;
  end;
  
  ds_model = lin_ds_from_lin_evol_model(lin_evol_model);
  %ds_model.u_function_ptr = @(params) 0.1; % define new
  ds_model_data = gen_model_data(ds_model,[]);
  mu = [0.5,1,1];
  ds_model = ds_model.set_mu(ds_model,mu);
  ds_sim_data = detailed_simulation(ds_model,ds_model_data,[]);
  lin_evol_model = lin_evol_model.set_mu(lin_evol_model,mu);
  lin_evol_sim_data = detailed_simulation(lin_evol_model,...
                                          lin_evol_model_data,[]);
  err = lin_evol_sim_data.U - ds_sim_data.X;
  disp(['max error in DOFs of lin_evol and ds simulation: ',...
       num2str(max(abs(err(:))))]);

  lin_evol_model.plot_title='detailed lin_evol';
  plot_sim_data(lin_evol_model,lin_evol_model_data,lin_evol_sim_data);
  lin_evol_model.plot_title='detailed lin_ds';
  lin_evol_from_ds_sim_data.U = ds_sim_data.X;
  lin_evol_from_ds_sim_data.y = ds_sim_data.Y;
  plot_sim_data(lin_evol_model,lin_evol_model_data,lin_evol_from_ds_sim_data);
  lin_evol_model.plot_title='error';
  error_sim_data = lin_evol_sim_data;
  error_sim_data.U = error_sim_data.U-ds_sim_data.X;
  plot_sim_data(lin_evol_model,lin_evol_model_data,error_sim_data);
  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %step = 4 % reduced Basis generation by single trajectory and
          %comparison of detailed and reduced simulation
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 case 4 % generate and store reduced basis and test reduced simulation   
  lin_evol_model = advection_fv_output_model([]);
  %lin_evol_model.debug = 1;
  lin_evol_model_data = gen_model_data(lin_evol_model,[]);
  
  % set some additional fields in model which will be copied into
  % ds model depending on whether constants are known or not:
  estimate_constants = 0;
  if estimate_constants
    %%% if constant estimaton has to be performed:
    lin_evol_model.estimate_lin_ds_nmu = 4;
    lin_evol_model.estimate_lin_ds_nX = 10;
    lin_evol_model.estimate_lin_ds_nt = 4;
  else
    %%% otherwise:
    % the following values are rough bounds, so could be 
    % non-rigorous, then larger search must be performed.
    lin_evol_model.state_bound_constant_C1 = 1;
    lin_evol_model.output_bound_constant_C2 = 1.2916;
  end;
  
  ds_model = lin_ds_from_lin_evol_model(lin_evol_model);
  %ds_model.u_function_ptr = @(params) 0.1; % define new
  ds_model_data = gen_model_data(ds_model,[]);
  mu = [1,1,0.5];
  %mu = [0,0.5,1];
  ds_model = ds_model.set_mu(ds_model,mu);
  ds_sim_data = detailed_simulation(ds_model,ds_model_data,[]);
  save(fullfile(rbmatlabresult,'trajectory1'),...
       'ds_model','ds_model_data','ds_sim_data');

%  %  mu = [1,1,0.5];
%  mu = [0,0.5,1];
%  ds_model = ds_model.set_mu(ds_model,mu);
%  ds_sim_data = detailed_simulation(ds_model,ds_model_data,[]);
%  save(fullfile(rbmatlabresult,'trajectory2'),...
%       'ds_model','ds_model_data','ds_sim_data');
%    
%  % single trajectory for basis, POD
  
  load(fullfile(rbmatlabresult,'trajectory1'),...
       'ds_model','ds_model_data','ds_sim_data');
  
  rbfname = 'advection_fv_output_basis';        
  RB = PCA_fixspace(ds_sim_data.X,[],ds_model_data.G);
  % select correct orthogonal set
  K = RB' * ds_model_data.G * RB;
  i = find(abs(diag(K)-1)>1e-2);
  i = min(i);
  RB = RB(:,1:(i-1));
  V = RB;
  W = ds_model_data.G * RB;  
  save(fullfile(rbmatlabresult,rbfname),'V','W','ds_model','ds_model_data');
  
  ds_model.RB_filename = 'advection_fv_output_basis';
  detailed_data = gen_detailed_data(ds_model,ds_model_data,[]);
  
  %perform single reduced simulation for given parameter
  reduced_data = gen_reduced_data(ds_model,detailed_data,[]);
  rb_sim_data = rb_simulation(ds_model,reduced_data,[]);
  rb_sim_data = rb_reconstruction(ds_model,detailed_data,rb_sim_data)

  % plot of reduced trajectory
  lin_evol_model = ds_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');
  
  keyboard;
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 case 5 % generate plots of detailed data for poster: 
 % step = 5 % figures of detailed simuations based on step 4
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  

  
  rbfname = 'advection_fv_output_basis';        
  load(fullfile(rbmatlabresult,rbfname),'ds_model','ds_model_data');
    
  % plots of trajectory 1 
  load(fullfile(rbmatlabresult,'trajectory1'));
  lin_evol_model = ds_model.base_model;
  lin_evol_model_data = gen_model_data(lin_evol_model,[]);  
  
  params.no_lines = 1;
  params.show_colorbar = 0;
  params.axis_equal = 1;
  params.plot = lin_evol_model.plot;
  lin_evol_model.plot(ds_sim_data.X(:,1),lin_evol_model_data.grid, ...
                      params);
  set(gca,'Xtick',[])
  set(gca,'Ytick',[])
  saveas(gcf,'trajectory1_t1.eps','epsc');
  close(gcf);
  
  lin_evol_model.plot(ds_sim_data.X(:,end),lin_evol_model_data.grid, ...
                      params);
  set(gca,'Xtick',[])
  set(gca,'Ytick',[])
  saveas(gcf,'trajectory1_tend.eps','epsc');
  close(gcf);

  p = ds_model.plot_sim_data(ds_model, ds_model_data, ds_sim_data);
 % Ylim = get(gca,'Ylim');
 % Ylim = [Ylim(1),Ylim(2)*1.1];
  set(gca,'Ylim',[0,0.09]);
  set(p(2),'Color',[0,0.8,0]);
  saveas(gcf,'trajectory1_output.eps','epsc');
  close(gcf);
  close(gcf);
   
  % plots of trajectory 2 
  load(fullfile(rbmatlabresult,'trajectory2'));
  lin_evol_model = ds_model.base_model;
  lin_evol_model_data = gen_model_data(lin_evol_model,[]);  
  
  params.no_lines = 1;
  params.show_colorbar = 0;
  params.axis_equal = 1;
  lin_evol_model.plot(ds_sim_data.X(:,1),lin_evol_model_data.grid, ...
                      params);
  set(gca,'Xtick',[])
  set(gca,'Ytick',[])
  saveas(gcf,'trajectory2_t1.eps','epsc');
  close(gcf);

  lin_evol_model.plot(ds_sim_data.X(:,end),lin_evol_model_data.grid, ...
                      params);
  set(gca,'Xtick',[])
  set(gca,'Ytick',[])
  saveas(gcf,'trajectory2_tend.eps','epsc');
  close(gcf);
  
  p = ds_model.plot_sim_data(ds_model, ds_model_data, ds_sim_data);
%  Ylim = get(gca,'Ylim');
%  Ylim = [Ylim(1),Ylim(2)*1.1];
  set(gca,'Ylim',[0,0.9]);
  set(p(2),'Color',[0,0.8,0]);
  saveas(gcf,'trajectory2_output.eps','epsc');
  close(gcf);
  close(gcf);

  
  disp('generated 4 pictures of 2 trajectories')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% step = 6 % figures of reduced simuations based on step 4
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 case 6 % generate plots of reduced model for poster 
  
  % as reference
  load(fullfile(rbmatlabresult,'trajectory1'));
  
  ds_model.error_estimation = 1;
  ds_model.RB_filename = 'advection_fv_output_basis';
  detailed_data = gen_detailed_data(ds_model,ds_model_data,[]);
  reduced_data = gen_reduced_data(ds_model,detailed_data,[]);
  
  Ns = [58,50,40,30,20,10];
  
  for nind = 1:length(Ns)
    % perform reduced simulation
    ds_model.N = Ns(nind);
    rb_sim_data = rb_simulation(ds_model,reduced_data,[]);
    rb_sim_data = rb_reconstruction(ds_model,detailed_data,rb_sim_data)
    
    params.no_lines = 1;
    params.show_colorbar = 1;
    params.axis_equal = 1;
    
    % plot of reduced output
    figure;
    p = plot(rb_sim_data.time,...
                [rb_sim_data.Y;...
                 rb_sim_data.Y+rb_sim_data.DeltaY; ...
                 rb_sim_data.Y-rb_sim_data.DeltaY; ...
                 ds_sim_data.Y]);
    colors = {[1,0,0],[1,0,0],[1,0,0],[0,1,0]};
    linestyles = {'-','-.','-.','-'};
    widths = [2,2,2,2];
    for i = 1:length(p)
      set(p(i),'Color',colors{i},'Linestyle',linestyles{i},...
               'Linewidth',widths(i));
    end;
    % set marker on rb_sim_data:
    Xdata = get(p(1),'XData');
    Ydata = get(p(1),'YData');
    Zdata = get(p(1),'ZData');
    set(p(1),'XData',Xdata(1:4:end));
    set(p(1),'YData',Ydata(1:4:end));
    set(p(1),'ZData',Zdata(1:4:end));
    set(p(1),'Marker','o');
    
    legend(['y\^(t): reduced output'],...
           'y\^(t)+\Delta_y(t)',...
           'y\^(t)-\Delta_y(t)',...
           ['y(t): detailed output'],'Location','NorthWest');
    title(['k = ',num2str(ds_model.N)],'Fontsize',15);
    saveas(gcf,['output_estimation',num2str(ds_model.N),'.eps'],'epsc');
    close(gcf);
  end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% step = 7 % time measurements
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 case 7
  
  clear all;
  Ns = [50,40,30,20,10]; %,30,20,10];
  nreps = 5;
  t_detailed = zeros(1,nreps);
  t_reduced_with_err_est = zeros(length(Ns),nreps);
  t_reduced_without_err_est = zeros(length(Ns),nreps);
  
  % initialize detailed model and simulate
  
  lin_evol_model = advection_fv_output_model([]);
  %lin_evol_model.debug = 1;
  %lin_evol_model_data = gen_model_data(lin_evol_model,[]);
  
  % set some additional fields in model which will be copied into
  % ds model depending on whether constants are known or not:
  estimate_constants = 0;
  if estimate_constants
    %%% if constant estimaton has to be performed:
    lin_evol_model.estimate_lin_ds_nmu = 4;
    lin_evol_model.estimate_lin_ds_nX = 10;
    lin_evol_model.estimate_lin_ds_nt = 4;
  else
    %%% otherwise:
    % the following values are rough bounds, so could be 
    % non-rigorous, then larger search must be performed.
    lin_evol_model.state_bound_constant_C1 = 1;
    lin_evol_model.output_bound_constant_C2 = 1.2916;
  end;
  
  ds_model = lin_ds_from_lin_evol_model(lin_evol_model);
  %ds_model.u_function_ptr = @(params) 0.1; % define new
  ds_model_data = gen_model_data(ds_model,[]);
  
%  disp('skipping detailed simulations...')
%  if 0
    
  for r = 1:nreps
    disp(['repetition = ',num2str(r)])
    mu = [1,1,0.5]-rand(1,3)*0.01;
    ds_model = ds_model.set_mu(ds_model,mu);
    tic,
    ds_sim_data = detailed_simulation(ds_model,ds_model_data,[]);
    t_detailed(1,r) = toc;
  end;

%  end;
  
  disp('--------------------------------------------------')
  disp('detailed simulation')
  disp(t_detailed);
  disp(['mean t_detailed=',num2str(mean(t_detailed))]);
  disp('--------------------------------------------------')
  
  % reduced simulation with and without error estimation  
  
  clear('ds_sim_data');
  ds_model.RB_filename = 'advection_fv_output_basis';
  detailed_data = gen_detailed_data(ds_model,ds_model_data,[]);
  reduced_data = gen_reduced_data(ds_model,detailed_data,[]);
  
  %   reduced simulation with and without error estimation
  
  for err_est = 0:1;
    ds_model.error_estimation = err_est;
    disp(['err_est= ',num2str(err_est)]);
    for n_ind = 1:length(Ns)
      disp(['M = ',num2str(Ns(n_ind))]);
      ds_model.N = Ns(n_ind);
      for r = 1:nreps
        disp(['repetition = ',num2str(r)])
        mu = [1,1,0.5]-rand(1,3)*0.01;
        ds_model = ds_model.set_mu(ds_model,mu);
        tic;
        rb_sim_data = rb_simulation(ds_model,reduced_data,[]);
        if err_est == 1
          t_reduced_with_err_est(n_ind,r) = toc;
        else
          t_reduced_without_err_est(n_ind,r) = toc;
        end;
      end;
    end;
  end;
    
  % output measurements

  disp('--------------------------------------------------')
  disp('detailed simulation')
  disp(t_detailed);
  disp(['mean t_detailed=',num2str(mean(t_detailed))]);

  disp('--------------------------------------------------')
  disp('reduced simulation with error estimation')
  disp(t_reduced_with_err_est);
  disp(['mean t_reduced_with_err_est=',...
        num2str(mean(t_reduced_with_err_est,2)')]);
  disp(['Ns =',num2str(Ns)]);
  
  disp('--------------------------------------------------')
  disp('reduced simulation without error estimation')
  disp(t_reduced_without_err_est);
  disp(['mean t_reduced_without_err_est=',...
        num2str(mean(t_reduced_without_err_est,2)')]);
  disp(['Ns =',num2str(Ns)]);
  disp('--------------------------------------------------')

  save('time_measurements')
 otherwise
  disp('step number unknown');
end;