chemnitz_gdl_experiments.m

Go to the documentation of this file.

% script collecting the commands for running the chemnitz-experiments
% i.e. for the paper 
%
%B. Haasdonk, M. Ohlberger, M., G. Rozza: A Reduced Basis Method for Evolution 
%Schemes with Parameter-Dependent Explicit Operators.  Preprint 09/07 - N, 
%FB 10 , University of Muenster, Preprint Nr. 17/2007, 
%Institute of Mathematics, University of Freiburg, 2007, Accepted by ETNA. 

% Bernard Haasdonk 3.9.2007

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%% Select here, what is to be performed with the 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 collateral reduced basis of L_E operator
step = 3 % compare single detailed simulation with and without
         % interpolated space operator vary the parameters to ensure that
         % the empirical interpolation is stable
%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 of reduced basis in rb_demo_gui
%step = 7 % generate train and/or test-error-landscape by doing 
         % reduced simulations for these 
                 
% load linear gdl example, but turn it into nonlinear in the following
load demo_lin_evol_params

%CRBfname = 'chemnitz_CRB_new_interpol';
%detailed_data_fname = 'chemnitz_detailed_data_new_interpol';
%CRBfname = 'chemnitz_CRB_new';
%detailed_data_fname = 'chemnitz_detailed_data_new';
CRBfname = 'chemnitz_CRB_new_interpol';
detailed_data_fname = 'chemnitz_detailed_data_new_interpol';

% finite volume settings
%params.operators_algorithm = 'fv_operators_implicit_explicit';
%params.init_values_algorithm = 'fv_init_values';
params.implicit_operators_algorithm = 'fv_operators_implicit';
params.inner_product_matrix_algorithm = 'fv_inner_product_matrix';
params.data_const_in_time = 1;
params.L_E_local_name = 'fv_conv_explicit_space';
params.dt = params.T/params.nt;

% change to nonlin problem type
params.mu_names = {'beta','c_init'}; 
%params.mu_ranges = {[-1,1],[0, 0.4]};
params.mu_ranges = {[0,1],[0, 1]};
params.rb_problem_type = 'nonlin_evol';
params.k = 0;
params.beta = 0;
params.c_init = 1;

% settings for empirical interpolation
par.range = params.mu_ranges;
params.ei_numintervals = [4,4]; % 5x5 parameter grid
%params.Mmax = 300;
%params.Mmax = 500;
%params.Mmax = 150;
params.Mmax = 100;
params.ei_detailed_savepath = 'chemnitz_5x5_detailed';
params.ei_operator_savepath = 'chemnitz_5x5_ei_operator_new'; 
params.ei_time_indices = 1:params.nt; 
% spaeter einbauen: params.ei_stop_epsilon = 1e-12;
params.CRB_generation_mode = 'param-time-space-grid';
%params.ei_target_error = 'approx';
%params.ei_target_error = 'interpol';
%disp('taking interpolation error as target in empirical interpolation!!');

% settings for reduced basis generation
params.detailed_simulation_algorithm = 'detailed_simulation';
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 = 3*60*60; % 3 hours             
params.RB_stop_epsilon = 1e-9; 
params.RB_error_indicator = 'error'; % true error
%params.RB_error_indicator = 'estimator'; % Delta from rb_simulation
params.RB_stop_Nmax = 150; 
params.RB_generation_mode = 'uniform_fixed';
params.RB_numintervals = params.ei_numintervals;
params.RB_detailed_train_savepath = params.ei_detailed_savepath;

% caching of velocity field for acceleration of GDL-example
params.filecache_velocity_matrixfile_extract = 2;

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,'chemnitz_CRB.mat'));
  tmp = load(fullfile(rbmatlabresult,CRBfname));
  detailed_data = tmp.detailed_data;
  params = tmp.params;
  params.k = 0;
  params.beta = 0.4;
  params.c_init = 1;
  U = detailed_simulation(detailed_data.grid, params);
  params.title = 'detailed simulation without interpolation';
  plot_element_data_sequence(detailed_data.grid,U,params);
  M = 100;
  %M = round(size(detailed_data.QM,2)/2)
  %  M = round(size(detailed_data.QM,2)/2)
  params.M = M;
%  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 = round(size(detailed_data.QM{1},2)/2);
  params.M = size(detailed_data.QM{1},2);
%  params.nt = params.nt /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;
%  params = tmp.params;
  detailed_data = rb_basis_generation(detailed_data, ...
                                      params);
  save(fullfile(rbmatlabresult,detailed_data_fname),...
       'detailed_data','params');
 case 6
  load(fullfile(rbmatlabresult,detailed_data_fname));
  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);
  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 % train-error landscape and test-error landscape
  disp('warning: takes a few hours!');
  load(fullfile(rbmatlabresult,detailed_data_fname));

  % replace RB with dummy basis from linear case!!
  %disp('replacing RB with linear basis!! remove this afterwards!!');
  %tmp = load('demo_lin_evol_detailed_data');
  %detailed_data.RB = tmp.detailed_data.RB;  
  
 % 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 = 11; 
  %  Nmax = size(detailed_data.RB,2);
  Nmax = 50;
  Ns = round((0:(Nsamples-1)) * ...
             (Nmax-1)/(Nsamples-1)+1);
  
  %Msamples = 3; 
  Msamples = 11; 
  %  Mmax = size(detailed_data.BM,2);
  Mmax = 100; % 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);
          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;
        end;
      end;
      fprintf('.');
    end; 
    % define stability-region as error being smaller than
    % sqrt(diffmax * area), e.g. diffmax = 4, area = 2e-7
    
    stab_limit = 1e-3;
    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,...
                  ['chemnitz_gdl_experiment_step7_result_',...
                   runs{ru},'.mat']));
  end;
  
 otherwise
  error('step-number is unknown!');
end;

% TO BE ADJUSTED TO NEW SYNTAX
%| \docupdate