burgers_fv.m

Go to the documentation of this file.

% small script demonstrating the burgers equation with explicit fv
% discretization and RB model reduction

% Bernard Haasdonk 14.8.2007

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%% Select here, what is to be performed with the burgers data %%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
step = 1 % single detailed simulation with given data and plot. Run
         % this with varying parameters mu until sure that scheme
         % is stable. Modify dt or the data-functions accordingly,
         % until a nice parameter-domain with uniformly stable
         % detailed scheme is obtained.
%step = 2 % generate colateral reduced basis of L_E operator
%step = 3 % compare single detailed simulation with and without
         % interpolated space operator
%step = 4 % generate dummy reduced basis from single trajectory and check, if
         % ei_interpolation with projection on this space maintains
         % result. A simple reduced simulation can also be
         % performed. All results should be visually identical
%step = 5 % generate reduced basis
%step = 6 % time measurement of reduced simulation and 
         % use reduced basis in rb_demo_gui
%step = 7 % generate error-landscape over varying N and M
          % can take several hours!!!
         
% output-filenames in rbmatlabresult      
CRBfname = 'burgers_fv_CRB_interpol.mat';
detailedfname = 'burgers_fv_detailed_interpol.mat';
         
% grid parameters
params = [];
% rectangular or triangular
params.gridtype = 'triagrid';
params.grid_initfile = '2dsquaretriang';

%params.gridtype = 'rectgrid';
%params.xnumintervals = 20;
%params.ynumintervals = 20;
%params.xrange = [0,1];
%params.yrange = [0,1];

% set all to neuman-boundary by specifying "rectangles", all
% boundary within is set to boundary type
params.bnd_rect_corner1 = [-1; ...
                            -1];
params.bnd_rect_corner2 = [ 2; ...
                            2];

% -1 means dirichlet, -2 means neumann
params.bnd_rect_index =   [-2]; % first is dirichlet, second neuman

% time discretization
%params.T = 1.0;
%params.nt = 200; % guess and adjust later

params.T = 0.4;
params.nt = 80; % guess and adjust later
%params.nt = 100; % guess and adjust later
%params.T = 0.1;
%params.nt = 50; % guess and adjust later
%params.400 => non-bounded u

% output settings
params.verbose = 10;
params.axis_equal = 1;
params.clim = [-1,1];

% data functions
%params.name_init_values = 'blobs'; 
%params.radius = 0.1;
%params.gamma = 100;
%params.beta = 1; % init data weight between 0 and 1

% name of function in RBmatlab/datafunc/init_values
params.name_init_values = 'gradient_box'; 
% parameters for data functions
params.c_init_up = 1.0;
%params.c_init_lo = 1.0;
params.c_init_lo = 0.0;
%params.c_init_lo = -1.0;

% name of function in RBmatlab/datafunc/dirichlet_values
%params.name_dirichlet_values = 'homogeneous';
%params.c_dir = 0;

% name of function in RBmatlab/datafunc/neuman_values
params.name_neuman_values = 'homogeneous';
params.c_neu = 0;

%params.name_neuman_values = 'homogeneous';
%params.c_neu = 0;

% convective flux
params.name_flux = 'burgers_parabola';
% for simplifying computation in case of linear flux:
params.flux_linear = 0;
%params.flux_quad_degree = 2;

% parameter of chosen conv-flux
params.vmax = 0.30;
%params.vrot_angle = 0;
params.vrot_angle = -pi/4;
%params.vrot_angle = -pi/10;
params.flux_pdeg = 2; % burgers exponent

% rb-formulation
params.mu_names = {'c_init_lo','vrot_angle'}; 
params.mu_ranges = {[-1,1],[-pi/4, 0]};
params.rb_problem_type = 'nonlin_evol';
%params.problem_type = 'lin_evol';

% finite volume settings
params.name_diffusive_num_flux = 'none';
params.name_convective_num_flux = 'lax-friedrichs';
params.lxf_lambda = 1.66
%params.lxf_lambda = 1.25;
%params.lxf_lambda = fv_search_max_lxf_lambda([],params)
%params.name_convective_num_flux = 'enquist-osher';

% projection of analytical initial data to discrete function
params.init_values_algorithm = 'fv_init_values';
params.implicit_operators_algorithm = 'fv_operators_implicit';
% get mass matrix for inner product computation from DOF-vectors
params.inner_product_matrix_algorithm = 'fv_inner_product_matrix';
params.data_const_in_time = 1;

% settings for empirical interpolation
params.L_E_local_name = 'fv_conv_explicit_space';
par.range = params.mu_ranges;
params.ei_numintervals = [4,4]; % 8x8 parameter grid
%params.Mmax = 300;
params.Mmax = 150;
params.ei_detailed_savepath = 'burgers_fv_5x5_detailed';
params.ei_operator_savepath = 'burgers_fv_5x5_ei_operator'; 
params.ei_time_indices = 1:params.nt; 
params.CRB_generation_mode = 'param-time-space-grid';
params.ei_target_error = 'interpol';

% settings for reduced basis generation
params.detailed_simulation_algorithm = 'detailed_simulation';
%params.operators_algorithm = 'fv_operators_implicit_explicit';
%params.init_values_algorithm = 'fv_init_values';
%params.inner_product_matrix_algorithm = 'fv_inner_product_matrix';
params.error_algorithm = 'fv_error';
%params.lxf_lambda = 1.0194e+003;
%params.data_const_in_time = 1;
params.error_norm = 'l2';
params.RB_extension_algorithm = 'RB_extension_PCA_fixspace';
%                   'RB_extension_max_error_snapshot'}; 
params.RB_stop_timeout = 2*60*60; % 2 hours             
params.RB_stop_epsilon = 1e-5; 
params.RB_error_indicator = 'error'; % true error
%params.RB_error_indicator = 'estimator'; % Delta from rb_simulation
params.RB_stop_Nmax = 100; 
params.RB_generation_mode = 'uniform_fixed';
params.RB_numintervals = params.ei_numintervals;
params.RB_detailed_train_savepath = params.ei_detailed_savepath;

switch step
 case 1 % single detailed simulation and plot
  disp('performing single detailed simulation')
  grid = construct_grid(params);
  U = detailed_simulation(grid, params);
  plot_element_data_sequence(grid,U,params);
 case 2 % empirical interpolation of space operator
  disp('constructing colateral reduced basis')
  params.RB_generation_mode = 'none';
  detailed_data = rb_detailed_prep(params);
  save(fullfile(rbmatlabresult,CRBfname),...
       'detailed_data','params');
  detailed_data.RB = zeros(detailed_data.grid.nelements,1);
  rb_plot_detailed_data(detailed_data,params);
  close(gcf-2);
 case 3 % compare detailed simulation with and without interpolation
  disp(['comparison between detailed simulation with and without', ...
        ' interpolation']);
  tmp = load(fullfile(rbmatlabresult,CRBfname));
  detailed_data = tmp.detailed_data;
  U = detailed_simulation(detailed_data.grid, params);
  params.title = 'detailed simulation without interpolation';
  plot_element_data_sequence(detailed_data.grid,U,params);
  params.M = size(detailed_data.QM{1},2);
%  params.M = 400;
  Ui =  detailed_ei_nonlin_evol_simulation(detailed_data,params);
  params.title = 'detailed simulation with empirical interpolation';
  plot_element_data_sequence(detailed_data.grid,Ui,params);
  params.title = 'error';
  plot_element_data_sequence(detailed_data.grid,abs(Ui-U),params);
  disp(['maximum l-infty error:',num2str(max(max(abs(Ui-U))))]);
 case 4 % construct dummy reduced basis by single trajectory and simulate
  load(fullfile(rbmatlabresult,CRBfname));
  disp('detailed interpolated simulation for basis construction:')
  params.M = size(detailed_data.QM{1},2);
  Ui =  detailed_ei_nonlin_evol_simulation(detailed_data,params);
  A = feval(params.inner_product_matrix_algorithm,detailed_data.grid,params);
  UON = orthonormalize(Ui,A);
  detailed_data.RB = UON;
  disp('reduced simulation:')
  offline_data = rb_offline_prep(detailed_data,params);
  params.N = size(detailed_data.RB,2);
  params.M = size(detailed_data.QM{1},2);
  reduced_data = rb_online_prep(offline_data, params);
  simulation_data = rb_simulation(reduced_data,params);
  Uappr = rb_reconstruction(detailed_data,simulation_data);

  disp('detailed interpolated and rb-projected simulation:')
  params.M = size(detailed_data.QM{1},2);
  params.N = size(detailed_data.RB,2);
  Uirb =  detailed_ei_rb_proj_nonlin_evol_simulation(detailed_data,params);

  params.title = 'reduced simulation result';   
  plot_element_data_sequence(detailed_data.grid,Uappr,params);
  params.title = 'detailed ei_interpol simulation result';      
  plot_element_data_sequence(detailed_data.grid,Ui,params);
  params.title = 'detailed ei_interpol and rb_projected simulation result';
  plot_element_data_sequence(detailed_data.grid,Uirb,params);
  
 case 5 % reduced basis
  disp('constructing reduced basis')
  tmp = load(fullfile(rbmatlabresult,CRBfname));
  detailed_data = tmp.detailed_data;
  detailed_data = rb_basis_generation(detailed_data, ...
                                      params);
  save(fullfile(rbmatlabresult,detailedfname),...
       'detailed_data','params');
  plot(detailed_data.RB_info.max_err_sequence);
  set(gca,'Yscale','log');
  title('RB-generation error convergence');
 case 6
  load(fullfile(rbmatlabresult,detailedfname));
  disp('reduced simulation:')
  offline_data = rb_offline_prep(detailed_data,params);
  params.N = size(detailed_data.RB,2);
  params.M = size(detailed_data.QM,2);
  reduced_data = rb_online_prep(offline_data, params);
  tic;
  simulation_data = rb_simulation(reduced_data,params); 
  t = toc;
  disp(['time for online phase: t = ',num2str(t)]);
  
  disp('full simulation:')
  tic;
  U = detailed_simulation(detailed_data.grid, params);
  t = toc;
  disp(['time for detailed simulation: t = ',num2str(t)]);
 
  demo_rb_gui(detailed_data,[],params);

 case 7 % training-error landscape
  disp('warning: takes a few hours!');
  load(fullfile(rbmatlabresult,detailedfname));
  
  % runs = {'train_linRB','test_linRB'};
  %  dirnames = {testdir};
  %  testdir = 'chemnitz_test_100';
  
  %  runs = {'train','test'};
  %  dirnames = {params.RB_detailed_train_savepath, testdir};
  
  runs = {'train'};
  dirnames = {params.RB_detailed_train_savepath};

  %runs = {'test'};
  %dirnames = {testdir};

  % generate test data if not existent
  if ismember('test',runs)
    % generate test data if not existent 
    rand('state',123456);
    M = rand_uniform(100,params.mu_ranges);
    save_detailed_simulations(M,testdir,detailed_data.grid,params)
  end;
   
  % load demo_nonlin_evol_detailed_data_LE_on_RB;
  offline_data = rb_offline_prep(detailed_data,params);
%  Nsamples = 3; 
  Nsamples = 20; 
    Nmax = size(detailed_data.RB,2);
  %Nmax = 150;
  Ns = round((0:(Nsamples-1)) * ...
             (Nmax-1)/(Nsamples-1)+1);
  
  %Msamples = 3; 
  Msamples = 25; 
  Mmax = size(detailed_data.BM,2);
  %Mmax = 300; % size(detailed_data.BM,2);
  Ms = round((0:(Msamples-1)) * ...
             (Mmax-1)/(Msamples-1)+1);
  
  for ru = 1:length(runs);
    disp(['starting run ',runs{ru}]);
    % storage for maximum l2-errors (i.e. linfty([0,T],l2)-norm)
    errs = zeros(Msamples,Nsamples);
    inds = zeros(Msamples,Nsamples);
    mu_inds = zeros(Msamples,Nsamples);
    % assuming, that training data is in the following directory
    settings = load(fullfile(rbmatlabtemp,...
                             dirnames{ru},'settings.mat'));
    
    for mu_ind = 1:size(settings.M,2);
      disp(['processing parameter vector ',num2str(mu_ind),...
            '/',num2str(size(settings.M,2))]);
      params = set_mu(settings.M(:,mu_ind),params);
      sim_data = load_detailed_simulation(mu_ind,settings.savepath,params);
      U_H = sim_data.U;
      for Nind = 1:Nsamples
        N = Ns(Nind);
        for Mind = 1:Msamples
          fprintf('.');
          M = Ms(Mind);
          
          % nonlinear simulation
          params.M = M;
          params.N = N;
          reduced_data = rb_online_prep(offline_data,params);
          try
            simulation_data = rb_simulation(reduced_data,params);
            U = rb_reconstruction(detailed_data,simulation_data);
            l2_errors = fv_l2_error(U_H,U,detailed_data.W);
            [err,ind] = max(l2_errors);
            if err(1)>errs(Mind,Nind)
              errs(Mind,Nind) = err(1);
              inds(Mind,Nind) = ind(1);
              mu_inds(Mind,Nind) = mu_ind;
            end;
          catch
            errs(Mind,Nind) = NaN;
            inds(Mind,Nind) = NaN;
            mu_inds(Mind,Nind) = mu_ind;            
          end
        end;
      end;
      fprintf('.');
      disp('saving temporary workspace in burgers_fv_step7_*_tmp.mat');
      save(fullfile(rbmatlabtemp,...
                    ['burgers_fv_step7_result_',...
                     runs{ru},'_tmp.mat']));
    end; 
    % define stability-region as error being smaller than
    % sqrt(diffmax * area), e.g. diffmax = 4, area = 2e-7
    
    stab_limit = 1e-2;
    C = ones(size(errs));
    i = find(errs<stab_limit);
    C(i) = 2;
    
    if isempty(find(C==1)); % i.e. all stable
      f1 = figure, surf(Ms, Ns,errs');
    else % some not stable   
      f1 = figure, surf(Ms, Ns,errs',C');
      shading interp;
    end;
    % figure, pcolor(Ms, Ns,C);
    title([runs{ru},' L-infty([0,T],L2) error']);
    set(gca,'Zscale','log');
    xlabel('M');
    ylabel('N');
    zlabel('error');
    
    f2 = figure, surf(Ms, Ns,inds');
    title([runs{ru},' time index of maximum l2-error']);
    %set(gca,'Zscale','log');
    xlabel('M');
    ylabel('N');
    zlabel('time index');
    
    save(fullfile(rbmatlabresult,...
                  ['burgers_fv_step7_result_',...
                   runs{ru},'.mat']));
  end;

 otherwise
  error('step-number is unknown!');
end;


% TO BE ADJUSTED TO NEW SYNTAX
%| \docupdate