riemann_burgers.m

Go to the documentation of this file.

% script demonstrating the burgers equation with explicit fv
% discretization, empirical interpolation and RB model reduction
%
% a step needs to be selected in this script, default is 1.
%
% example is meant to demonstrate automatic space dimension reduction
% of the algorithm.
% The interpolation basis functions are nicely vertical stripes.
% The error convergence in ei however is constant until 100 basis
% functions then drops to zero. Interesting fact, but clearly M not
% variable afterwards!!!
% As a reduced basis, two shock-trajectories are used.

% Bernard Haasdonk 16.6.2008

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%% 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
         % takes still 20 minutes if snapshots are already computed!!
%step = 3 % compare single detailed simulation with and without
         % interpolated space operator
%step = 4 % generate dummy reduced basis from trajectories 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 

%step = 7 % generate error-landscape over varying N (and eventually M)
          % can take several hours!!!
%step = 8 % generate different Figures for presentation and time measurement
%step = 9 % start demo_rb_gui with error estimator        
%step = 10 % start demo_rb_gui without error estimator    

% output-filenames in rbmatlabtemp        
CRBfname = 'riemann_burgers_CRB.mat';
detailedfname = 'riemann_burgers_detailed.mat';
testdir = 'riemann_burgers_test_10';
results_fname_base = 'riemann_burgers_step7_result';
ei_detailed_savepath = 'riemann_burgers_ei_detailed';
ei_operator_savepath = 'riemann_burgers_detailed';

% old settings, now in riemann_burgers_model:

model = riemann_burgers_model([]);

%% grid parameters
%params = [];
%% rectangular or triangular
%%params.gridtype = 'triagrid';
%%params.grid_initfile = '2dsquaretriang';

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

% set upper and lower to neuman noflow
% set left and right to dirichlet

%%epsilon = 1e-6;
%params.bnd_rect_corner1 = [-1, 0-eps, 1-eps; ...
%                          -1, 0-eps, 0-eps];
%params.bnd_rect_corner2 = [ 2, 0+eps, 1+eps; ...
%                           2, 1+eps, 1+eps];

% -1 means dirichlet, -2 means neuman
%params.bnd_rect_index =   [-2,-1,-1]; 

% time discretization
%params.T = 0.5;
%params.nt = 100; % 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 = 9;
%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 = 'waveproduct'; 
%params.name_init_values = 'leftright'; 
%params.name_init_values = 'as_dirichlet'; 

%params.c_init_max = 1.0;
%params.c_init_phase_x = 0.0;
%params.c_init_phase_y = 0.0;
%params.c_init_freq_x = 2*pi;
%params.c_init_freq_y = 2*pi;

%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 = 'leftright';
%params.c_dir_left = -1.0;
%params.c_dir_right = 0.0;
%params.dir_middle = 0.5;

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

% convective flux
%params.name_flux = 'burgers_parabola';
%params.name_flux = 'burgers';
% for simplifying computation in case of linear flux:
%params.flux_linear = 0;
% horizontal right velocity field
%params.flux_vx = 1;
%params.flux_vy = 0;
%params.flux_quad_degree = 2;
%params.flux_pdeg = 2; % burgers exponent

% rb-formulation
%params.mu_names = {'c_init_lo','vrot_angle'}; 
% themiddle parameter makes u_init not decomposable...
%params.mu_names = {'c_dir_left','c_dir_right','dir_middle'}; 
%params.mu_names = {'c_dir_left','c_dir_right','flux_vx'}; 
%params.mu_ranges = {[0,1],[0,1],[0,1]};
%params.mu_ranges = {[-1,1],[-1,1],[-1,1]};
%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 = 0.35355;
%params.lxf_lambda = fv_search_max_lxf_lambda([],params);
%keyboard;
%params.name_convective_num_flux = 'engquist-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,4]; % 5^3 parameter grid
%params.ei_numintervals = [2,2,2]; % 5^3 parameter grid
%params.ei_numintervals = [4,4]; % 5^2 parameter grid
%params.Mmax = 300;
%params.Mmax = 150;
%params.ei_detailed_savepath = ei_detailed_savepath;
%params.ei_operator_savepath = ei_operator_savepath;
%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_generation_mode = 'PCA_trajectories';
%params.RB_mu_list = {[1,0,1],[0,1,-1]}; % two shock-trajectories in
                                        % different direction
%params.RB_numintervals = params.ei_numintervals;
%params.RB_detailed_train_savepath = params.ei_detailed_savepath;

% such that demo_rb_gui looks nicer:
%params.yscale_uicontrols = 0.6;
%params.xscale_gui = 0.5;

model_data = gen_model_data(model);

switch step
 case 1 % single detailed simulation and plot
  disp('performing single detailed simulation')
%  grid = construct_grid(model);
%  shock moving to left:
%  model = set_mu(model,[-0.5,1,-1]); % cdir_left, c_dir_right, vx
%  non-moving shock:
%  model = set_mu(model,[-1,1,-1]); % cdir_left, c_dir_right, vx
%  symmetric rarefaction wave:
%  model = set_mu(model,[1,-1,-1]); % cdir_left, c_dir_right, vx
%  moving shock with velocity 0.5:
  model = set_mu(model,[0,1,-1]); % cdir_left, c_dir_right, vx
  sim_data = detailed_simulation(model, model_data);
  plot_params = [];
  plot_sim_data(model,model_data,sim_data,plot_params);
 case 2 % empirical interpolation of space operator
  tic;
  disp('constructing collateral reduced basis')
  model.RB_generation_mode = 'none';
  detailed_data = gen_detailed_data(model,model_data);
  save(fullfile(rbmatlabtemp,CRBfname),...
       'detailed_data','model');
  detailed_data.RB = zeros(detailed_data.grid.nelements,1);
  plot_params = [];
  plot_params.plot = model.plot;
  plot_detailed_data(model,detailed_data,plot_params);
  close(gcf-2);
  t = toc;
  disp(['computation time for collateral basis: ',num2str(t)]);
 case 3 % compare detailed simulation with and without interpolation
  disp(['comparison between detailed simulation with and without', ...
        ' interpolation']);
  tmp = load(fullfile(rbmatlabtemp,CRBfname));
  detailed_data = tmp.detailed_data;
%  params.flux_pdeg = 2;
  sim_data = detailed_simulation(model,model_data);
  plot_params.title = 'detailed simulation without interpolation';
  plot_sim_data(model,model_data,sim_data,plot_params);
  model.M = size(detailed_data.QM{1},2);
%  model.M = size(detailed_data.QM,2);
%  params.M = 400;
  sim_data_interpol =  nonlin_evol_detailed_ei_simulation(model,detailed_data);
  plot_params.title = 'detailed simulation with empirical interpolation';
  plot_sim_data(model,model_data,sim_data_interpol,plot_params);
  plot_params.title = 'error';
  err = sim_data;
  err.U = err.U-sim_data_interpol.U;
  plot_sim_data(model,model_data,err,plot_params);
  disp(['maximum l-infty error:',num2str(max(max(abs(err.U))))]);
 
 case 4 % construct dummy reduced basis by single trajectory and simulate
  load(fullfile(rbmatlabtemp,CRBfname));
  model = riemann_burgers_model;
  % set mu values new, as may be different in stored file
  model.c_dir_left = 1.0;
  model.c_dir_right = 0.0;
  model.dir_middle = 0.5;
  
  disp('detailed interpolated simulation for basis construction:')
  model.M = size(detailed_data.QM{1},2);
  sim_data_interpol =  ...
      nonlin_evol_detailed_ei_simulation(model,detailed_data);
%  A = feval(params.inner_product_matrix_algorithm, ...
%           detailed_data.grid,params);
  UON = model.orthonormalize(model,model_data,sim_data_interpol.U,1e-4);
%  UON = model.orthonormalize(model,model_data,sim_data_interpol.U,1e-4);
  detailed_data.RB = UON;
  model.N = size(detailed_data.RB,2);
  model.M = size(detailed_data.QM{1},2);
  save(fullfile(rbmatlabtemp,'riemann_burgers_temp'),'model','detailed_data')
%  keyboard;
  disp('reduced simulation:')
  reduced_data = gen_reduced_data(model,detailed_data);
  sim_data = rb_simulation(model,reduced_data);
  sim_data = rb_reconstruction(model,detailed_data,sim_data);

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

  plot_params.title = 'reduced simulation result';      
  plot_sim_data(model,model_data,sim_data,plot_params);
  plot_params.title = 'detailed ei_interpol simulation result'; 
  plot_sim_data(model,model_data,sim_data_interpol,plot_params);
  plot_params.title = 'detailed ei_interpol and rb_projected simulation result';
  plot_sim_data(model,model_data,sim_data_ei_rb_proj,plot_params);
  
 case 5 % reduced basis construction
  disp('constructing reduced basis')
  tmp = load(fullfile(rbmatlabtemp,CRBfname));
  model = riemann_burgers_model;
  detailed_data = tmp.detailed_data;
  detailed_data = rb_basis_generation(model,detailed_data);
  % set these fields such that demo_rb_gui works nicely
  model.N = size(detailed_data.RB,2);
  model.M = size(detailed_data.QM{1},2);
  save(fullfile(rbmatlabtemp,detailedfname),...
       'detailed_data','model');

  if isfield(detailed_data.RB_info,'max_err_sequence')
    plot(detailed_data.RB_info.max_err_sequence);
    set(gca,'Yscale','log');
    title('RB-generation error convergence');
  else
    disp(['skipping plot of error decrease as basis generation does' ...
          ' not provide it']);
  end;
 
 case 6 % time measurement of reduced simulation 

  load(fullfile(rbmatlabtemp,detailedfname));
  disp('reduced simulation:')
  model.N = size(detailed_data.RB,2);
  model.M = size(detailed_data.QM{1},2);
  reduced_data = gen_reduced_data(model,detailed_data);
  tic;
  sim_data = rb_simulation(model,reduced_data); 
  t = toc;
  disp(['time for online phase: t = ',num2str(t)]);
  
  disp('full simulation:')
  tic;
  sim_data = detailed_simulation(model,model_data);
  t = toc;
  disp(['time for detailed simulation: t = ',num2str(t)]);
 
 case 7 % training-error landscape
  disp('to be adjusted to new syntax!!!');
  disp('warning: takes a few hours!');
  load(fullfile(rbmatlabtemp,detailedfname));
  
  % runs = {'train_linRB','test_linRB'};
  %  dirnames = {testdir};
  %  testdir = 'chemnitz_test_100';
  
  %runs = {'train','test'};
  runs = {'test'};
  %dirnames = {params.RB_detailed_train_savepath, testdir};
  dirnames = {testdir};
  
  %  runs = {'train'};
  %  dirnames = {params.RB_detailed_train_savepath};
  %  runs = {'test'};
  %  dirnames = {testdir};

  %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(10,model.mu_ranges);
    save_detailed_simulations(M,testdir,detailed_data.grid,model)
  end;
   
  % load demo_nonlin_evol_detailed_data_LE_on_RB;
  offline_data = rb_offline_prep(detailed_data,model);
%  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{1},2);
  %Mmax = 300; % size(detailed_data.BM{1},2);

  if Msamples > 1
    Ms = round((0:(Msamples-1)) * ...
               (Mmax-1)/(Msamples-1)+1);
  else
    Ms = Mmax;
  end;
  
  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 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))]);
      model = set_mu(settings.M(:,mu_ind),model);
      sim_data = load_detailed_simulation(mu_ind,settings.savepath,model);
      U_H = sim_data.U;
      for Nind = 1:Nsamples
        N = Ns(Nind);
        for Mind = 1:Msamples
          fprintf('.');
          M = Ms(Mind);
          
          % nonlinear simulation
          model.M = M;
          model.N = N;
          reduced_data = rb_online_prep(offline_data,model);
          try
            simulation_data = rb_simulation(reduced_data,model);
            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');
      save(fullfile(rbmatlabtemp,...
                  [results_fname_base,'_',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;

    f1 = figure;
    if Msamples > 1
      error_axis = 'z';
      n_axis = 'y';
      if isempty(find(C==1)); % i.e. all stable
        surf(Ms, Ns,errs');
      else % some not stable   
        surf(Ms, Ns,errs',C');
        shading interp;
      end;
    else
      error_axis = 'y';
      n_axis = 'x';
      plot(Ns,errs);
    end;

%    keyboard;
    
    % figure, pcolor(Ms, Ns,C);
    title([runs{ru},' L-infty([0,T],L2) error']);
    set ( gca, [error_axis, 'scale'], 'log' );
    xlabel('M');
    %ylabel('N');
    eval([n_axis,'label(''N'');']);
    eval([error_axis,'label(''error'');']);
    
    f2 = figure;
    if Msamples > 1
      surf(Ms, Ns,inds');
    else
      plot(Ns, inds);
    end;
    title([runs{ru},' time index of maximum l2-error']);
    %set(gca,'Zscale','log');
    xlabel('M');
    %ylabel('N');
    eval([n_axis,'label(''N'');']);
    %zlabel('time index');
    eval([error_axis,'label(''time index'');']);
    
    save(fullfile(rbmatlabtemp,...
                  [results_fname_base,'_',runs{ru},'.mat']));
  end;

 case 8 % generate figures and runtime table
  
  % detailed_simulation: shock to right, v = 0.9

%  if 0
    
  model.c_dir_left = 0;  
  model.c_dir_right = 0.5;  
  model.flux_vx = 0.9;
  U1 = detailed_simulation(detailed_data.grid,model);

  model.c_dir_left = 0.5;  
  model.c_dir_right = 1;  
  model.flux_vx = -0.5;
  U2 = detailed_simulation(detailed_data.grid,model);

  model.no_lines = 1;
  figure;
  plot_element_data(detailed_data.grid,U1(:,1),model);
  title('Initial Data');
  saveas(gcf,'shock_initial_data1','epsc')
  figure;
  plot_element_data(detailed_data.grid,U2(:,1),model);
  title('Initial Data');
  saveas(gcf,'shock_initial_data2','epsc')
  figure;
  plot_element_data(detailed_data.grid,U1(:,end),model);
  title('End Data');
  saveas(gcf,'shock_end_data1','epsc')
  figure;
  plot_element_data(detailed_data.grid,U2(:,end),model);
  title('End Data');
  saveas(gcf,'shock_end_data2','epsc')
  
  % interpolation points
  
  model.no_lines = 1;
  model.M = length(detailed_data.TM);
  rb_plot_interpolation_points(detailed_data,model);  
  cmap = gray(256);
  cmap = cmap(1:150,:); % take dark shades 
  cmap(1,:) = [1.0,1.0,1.0]; % light yellow
  colormap(cmap);
  caxis([0,model.M]);
  %axis off
  box on
  saveas(gcf,'shock_interpolation_points','epsc')

  offline_data = rb_offline_prep(detailed_data,model);  
  figure, plot(offline_data.grid_local_ext);
  box on
  axis equal
  axis tight
  set(gca,'YLim',[0,1]);
  saveas(gcf,'shock_subgrid','epsc')
  
  % figure of error convergence

  %  load(fullfile(rbmatlabtemp,detailedfname));
  load(fullfile(rbmatlabtemp,...
                      [results_fname_base,'_test.mat']));
  
  f1 = figure;
  error_axis = 'z';
  n_axis = 'y';
  %plot(Ns,errs);  
  surf(Ms, Ns,errs');
%  figure, pcolor(Ms, Ns,C');
  title([runs{ru},' L-infty([0,T],L2) error']);
  set( gca, [error_axis, 'scale'], 'log' );
  xlabel('M');
  ylabel('N');
  eval([n_axis,'label(''N'');']);
  eval([error_axis,'label(''error'');']);
  set(gcf,'Position',[ 120   524   496   370]);
  view(120,30)
  tmp = load('redgreencolmap.mat');
%  cmap = jet(1000);
%  cmap = cmap(end:-1:1,:);
  colormap(tmp.c);
  %  disp('adjust colorbar!');
  set(gcf,'Color',[1.0,1.0,1.0]);
  saveas(gcf,'shock_error_convergence','epsc')

%    f1 = figure, surf(Ms, Ns,errs');
%  title(['test L^\infty([0,T],L^2) error'],'Fontsize',15);
%  title(['test error'],'Fontsize',15); %L^\infty([0,T],L^2) error'],'Fontsize',15);
%  set(gca,'Zscale','log');
%  xlabel('M','Fontsize',15);
%  ylabel('N','Fontsize',15);
%  zlabel('error','Fontsize',15);

  % local grid 
  
  %offline_data = rb_offline_prep(detailed_data,model);
  %plot(offline_data.grid_local_ext);
  %axis off
  %axis equal
  %axis tight
  %set(gcf,'Color',[1,1,1]);

%  end; % if 0
  
  % generate runtime table

  %  load(fullfile(rbmatlabtemp,detailedfname));
  load(fullfile(rbmatlabtemp,...
                      [results_fname_base,'_test.mat']));

  nreps = 10;
  
  Ns = [10,20,30,40,50,60, 70, 80, 90, 100];
  %  Ms = [15,30,45,60,75,90,105,120, 135,150];
  Ms = 100 * ones(size(Ns));
  t_detailed = zeros(nreps,1);
  for rep = 1:nreps
%    model.flux_pdeg = rand(1)+1;
    mu = rand_uniform(1,model.mu_ranges);
    model = set_mu(mu,model);
    %    offline_data = rb_offline_prep(detailed_data,model);
    tic;
    u = detailed_simulation(model);
    t_detailed(rep) = toc;
    disp(['runtime for detailed: ',num2str(t_detailed(rep)),' sec.']);
  end;
  
  t_reduced = zeros(nreps,length(Ns));
  offline_data = rb_offline_prep(detailed_data,model);    
  for rn = 1:length(Ns)
    model.N = Ns(rn);
    model.M = Ms(rn);
    reduced_data = rb_online_prep(offline_data,model);
    for rep = 1:nreps
      mu = rand_uniform(1,model.mu_ranges);
      model = set_mu(mu,model);
      tic;
      simulation_data = rb_simulation(reduced_data,model);
      t_reduced(rep,rn) = toc;
      disp(['N= ',num2str(model.N),', M= ',num2str(model.M)]);
      disp(['runtime for reduced: ',num2str(t_reduced(rep,rn)),' sec.']);
    end;
  end ;
  disp('detailed :');
  disp(' t_mean    t_std ');
  disp([num2str(mean(t_detailed)),'  ',num2str(std(t_detailed))]);
  
  disp('reduced :');
  disp('N     M       t_mean     t_std');
  o = [Ns(:), Ms(:), mean(t_reduced,1)', std(t_reduced,1)'];
  s = sprintf('%2.0f    %2.0f   %2.2f     %2.2f    \n ',o');
  disp(s);
  
  if 0
  % nicely arranged error landscape figures:
  
%  load(fullfile(rbmatlabtemp,...
%               [results_fname_base,'_train.mat']));
  load(fullfile(rbmatlabtemp,...
                [results_fname_base,'_test.mat']));
  
  end; % if 0

 case 9 % start demo_rb_gui with error estimation
  load(fullfile(rbmatlabtemp,detailedfname));
%  model = riemann_burgers_model;
  %  model.show_colorbar = 1;
  model.c_dir_left = -1;
  model.c_dir_right = 0;
  model.flux_vx = -1;
  model.M = 99;
  model.Mstrich = 1;
  model.N = 100;
  model.enable_error_estimator = 1;
  demo_rb_gui(model,detailed_data);

 case 10 % start demo_rb_gui without error estimation
  load(fullfile(rbmatlabtemp,detailedfname));
%  model = riemann_burgers_model;
  %  model.show_colorbar = 1;
  model.M = 100;
  model.Mstrich = 0;
  model.N = 100;
  model.enable_error_estimator = 0;
  demo_rb_gui(model,detailed_data);

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