nonlin_symmetry.m

Go to the documentation of this file.

function [ddescr, rdescr, dmodel, rmodel] = nonlin_symmetry(steps, combined, M_by_N_ratio, noof_ei_extensions, random, explicit, num_cpus, params)
% function [descr, rdescr, dmodel, rmodel] = nonlin_symmetry(steps, combined, M_by_N_ratio, noof_ei_extensions, random, explicit, num_cpus, params)
% small script demonstrating the burgers equation with explicit fv
% discretization and RB model reduction
% example is meant to demonstrate automatic symmetry detection by
% the algorithm.
% These results are presented at HYP2008

% Bernard Haasdonk 3.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
%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 demo rb gui
%step = 7 % generate error-landscape over varying N and M.
          % This can take several hours!!!
%step = 8 % generate estimator-landscape over varying N and M and
%           Mstrich. This  can take several hours!!!
%step = 9 % generate different Figures for presentation and time measurement
%step = 10 % start demo_rb_gui
%step = 11 % comparison of error estimator and error

if nargin == 0
  steps = [0,5];
end

%steps = [1,2,3,4,5,7,8];

if nargin <= 1
  combined           = true;
  M_by_N_ratio       = 0;
  noof_ei_extensions = 1;
  random             = false;
  explicit           = true;
  num_cpus           = 4;
end
if nargin <= 7
  params = [];
end

rmodel = [];

for si=1:length(steps)

step = steps(si);

% output-filenames in rbmatlabresult
detailedfname = 'nonlin_symmetry_detailed.mat';
testdir = 'nonlin_symmetry_test_10';
results_fname_base = 'nonlin_symmetry_step7_result';

plot_params.axis_equal = 1;
plot_params.clim = [0,1];
plot_params.plot = @fv_plot;

params.explicit = explicit;
params.model_type = 'nonlin_symmetry';
params.step7_outputfile = fullfile(rbmatlabresult, [params.model_type, '_step7_output']);
params.step8_outputfile = fullfile(rbmatlabresult, [params.model_type, '_step8_output']);

extrapar.Mstrich = 100;

% grid parameters
switch step
  case 0
    if ~isempty(getenv('MACHINEFILE'))
      machinefile = getenv('MACHINEFILE');
      FakeMPI.MPI_Init(machinefile, fullfile(rbmatlabhome, 'fakeMPI', 'run_parCompErrs.sh'));
      FakeMPI.MPI_Run;
    end

    ddescr = nonlin_symmetry_descr(params);
    rdescr.rb_problem_type = ddescr.rb_problem_type;

    if ~combined
      rdescr.Mmax_small = 0;
    else
      rdescr.Mmax_small = 20;
    end

    rdescr.extension_M_by_N_r       = M_by_N_ratio; % 2.5;
    rdescr.refinement_theta         = 0.15;
    rdescr.noof_ei_extensions       = noof_ei_extensions;
    rdescr.train_num_intervals      = 25;
    rdescr.val_size                 = 4;
    rdescr.stop_Mmax                = 600;
    rdescr.stop_Nmax                = 300;
    rdescr.stop_timeout             = 12*60*60;
    rdescr.stop_epsilon             = 1e-9;
    rdescr.epsilon_M_by_N_r         = 1e-1;
    rdescr.stop_max_val_train_ratio = 1.15;

    if random
      rdescr.train_sample_mode        = 'random';
      rdescr.train_num_intervals      = 400;
      rdescr.stop_max_val_train_ratio = inf;
    end

    [rmodel, dmodel] = gen_models(ddescr, rdescr);

    params.mu_ranges = rmodel.mu_ranges;
    params.mu_names  = rmodel.mu_names;

    dmodel.num_cpus   = num_cpus;

    model_data = gen_model_data(dmodel);
 case 1 % single detailed simulation and plot
  step1_detailed_simulation;
 case 2 % empirical interpolation of space operator
  step2_empirical_interpolation;
 case 3 % compare detailed simulation with and without interpolation
  step3_detailed_ei_simulation;
 case 4 % construct dummy reduced basis by single trajectory and simulate
  step4_dummy_reduced_basis;
 case 5 % reduced basis
   detailed_data = gen_detailed_data(rmodel, model_data);
   save(fullfile(rbmatlabresult, detailedfname), 'detailed_data', 'rmodel');
   FakeMPI.MPI_Finalize;
 case 6
  step6_demo_rb_gui;
 case 7
  load(fullfile(rbmatlabresult,detailedfname));
  model.verbose = 1;

  model.Mstrich   = 0;
  % maximum number of reduced basis functions
  model.Nmax      = 135;
  % maximum number of collateral reduced basis functions
  model.Mmax      = size(detailed_data.BM{1},2);
  % fixed ratio between M and N in couples plots.
  model.M_by_N_ratio = 1.5;
  % number of reduced basis sizes to be tested.
  Nsize           = 20;
  % number of collateral reduced basis sizes to be tested.
  Msize           = 25;
  % sample of values for the coupling variable 'c' for which errors are
  % computed.
  cmax        = min(1,model.Mmax/(model.M_by_N_ratio*model.Nmax));
  csample     = cmax/20:cmax/20:cmax;
  % axis description for plots with coupling variable 'c'
  cdescr      = '(N,M) = c * (Nmax, Mmax)';
  % size of mu vector test set.
  mu_set_size = 20;

  step7_error_landscape;
 case 7.5
  load(params.step65_output);
  for i = 1:length(output)
    plot_error_landscape(output{i});
  end
 case 8
  disp('warning: takes a few hours!');
  load(fullfile(rbmatlabresult,detailedfname));

  % maximum number of reduced basis vectors
  model.Nmax = size( detailed_data.RB, 2 );
  % maximum number of collateral reduced basis vectors used for
  model.Mmax = size( detailed_data.BM{1}, 2 ) - extrapar.Mstrich;
  % fixed ratio between M and N in couples plots.
  model.M_by_N_ratio = 1.5;

  % number of reduced basis sizes to be tested.
  Nsize           = 20;
  % number of collateral reduced basis sizes to be tested.
  Msize           = 25;
  % Mstrich values to be tested
  Mstrich_samples = 1:2:extrapar.Mstrich;
  % sample of values for the coupling variable 'c' for which
  % estimators are computed.
  cmax            = min(1,model.Mmax/(model.M_by_N_ratio*model.Nmax));
  csample         = cmax/20:cmax/20:cmax;
  % size of mu vector test set
  mu_set_size     = 20;
  % axis description for plots with coupling variable 'c'
  cdescr          = '(N,M)                                             = c * (Nmax, Mmax)';

  step8_estimator_landscape;
 case 8.5
  load(params.step75_output);
  for i = 1:length(output)
    plot_error_landscape(output{i});
  end
 case 9 % generate figures and runtime table
  load(fullfile(rbmatlabresult,detailedfname));

  % anfangsdaten, Enddaten p=1, Enddaten p=2

  model.flux_pdeg = 1;
  U1 = detailed_simulation(detailed_data.grid,model);
  model.flux_pdeg = 2;
  U2 = detailed_simulation(detailed_data.grid,model);

  model.no_lines = 1;
  figure;
  plot_element_data(detailed_data.grid,U1(:,1),model);
  title('Initial Data');
  set(gcf,'Position',[420   663   565   285]);
  saveas(gcf,'smooth_init_data.eps','epsc');

  figure;
  plot_element_data(detailed_data.grid,U1(:,end),model);
  title('End Data p=1');
  set(gcf,'Position',[420   663   565   285]);
  saveas(gcf,'smooth_end_data_p1.eps','epsc');

  figure;
  plot_element_data(detailed_data.grid,U2(:,end),model);
  title('End Data p=2');
  set(gcf,'Position',[420   663   565   285]);
  saveas(gcf,'smooth_end_data_p2.eps','epsc');

  model.no_lines = 1;
  model.M = length(detailed_data.TM{1});
  rb_plot_interpolation_points(detailed_data,model);
  set(gcf,'Position',[420   663   565   285]);
  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]);
  box on;
  keyboard;
  saveas(gcf,'smooth_interpolation_points.eps','epsc');

  % nicely arranged error landscape figures:

%  load(fullfile(rbmatlabresult,...
%               [results_fname_base,'_train.mat']));
  load(fullfile(rbmatlabresult,...
                [results_fname_base,'_test.mat']));

  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);
  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]);

  keyboard;

  % generate runtime table
  nreps = 10;

  Ns = [10,20,30,40,50,60, 70, 80, 90, 100];
  Ms = [15,30,45,60,75,90,105,120, 135,150];
  t_detailed = zeros(nreps,1);
  for rep = 1:nreps
    model.flux_pdeg = rand(1)+1;
%    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
      model.flux_pdeg = rand(1)+1;
      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);

  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]);

 case 10 % start demo_rb_gui
  load(fullfile(rbmatlabresult,detailedfname));
  % Mmax = 150; Nmax = 100
  model.N       = 100;
  model.M       = 150;
  model.Mstrich = 10;
  model.enable_error_estimator = 1;
%  model.RB_error_indicator = 'estimator';
  demo_rb_gui(model,detailed_data,[],plot_params);

 case 11 % compare error and error estimator
  load(fullfile(rbmatlabresult,detailedfname));
  model.N       = 135;
  model.M       = 290;
%  model.M       = 150;
%  model.Mstrich = 10; % set below in loop
  model.enable_error_estimator = 1;
 
  model_data = gen_model_data(model);
%  detailed_data = gen_detailed_data(model,model_data);
  reduced_data = gen_reduced_data(model, detailed_data);

  model = model.set_mu(model,[1.5]);
%  model = model.set_mu(model,2);
%  model = model.set_mu(model,1);
  sim_data = detailed_simulation(model,model_data);
%  Mstrich = [1,10,150];
  Mstrich = [1,5,10];
  for Mi = 1:length(Mstrich);
    model.Mstrich = Mstrich(Mi);
    rb_sim_data = rb_simulation(model, reduced_data);
    rb_sim_data = rb_reconstruction(model, detailed_data, rb_sim_data);

    De{Mi} = rb_sim_data.Delta;
  end;
  err_data = sim_data;
  err_data.U = sim_data.U - rb_sim_data.U; 
  err = sqrt(max(0,sum(err_data.U' * detailed_data.W * err_data.U)))';
  max_t_err = err;
  
  % somehow replace by matrix version!!!
  for i = 2:length(max_t_err)
    if max_t_err(i) < max_t_err(i-1)
      max_t_err(i) = max_t_err(i-1);
    end;
  end;
  figure;
  plot([De{1}(:),De{2}(:),De{3}(:),max_t_err(:),err(:)]);
  Ylim = get(gca,'Ylim');
  set(gca,'Ylim',[0,max(Ylim)]);
  title('rb state error and estimator');
  legend({'M''= 1, estimator \Delta(t)',...
          'M''= 5, estimator \Delta(t)',...
          'M''= 10, estimator \Delta(t)',...
          '||e||_{L\infty([0,t],L2)}',...
          '||e(t)||_{L2}'},'Location','NorthWest');
 
 otherwise
  error('step-number is unknown!');
end

end

%| \docupdate