weak_strong_POD_Greedy_comp.m

Go to the documentation of this file.

function res = weak_strong_POD_Greedy_comp(step)
% script comparing weak and strong POD-Greedy on lin_evol default
% example
%
% expected result:
% * strong POD-Greedy: longer runtime, better quality than
%      weak-POD-Greedy for same training set
% * weak-POD-Greedy on larger train set: 
%      faster than coarser strong POD-Greedy
%      better tes quality than coarser strong POD-Greedy

% B. Haasdonk 1.10.2012

model = convdiff_model;

if nargin<1
  step = 1;
end;

strong_PODGreedy_fn = 'strong_PODGreedy_fn.mat';
weak_coarse_PODGreedy_fn = 'weak_coarse_PODGreedy_fn.mat';
weak_fine_PODGreedy_fn = 'weak_fine_PODGreedy_fn.mat';
weak_finer_PODGreedy_fn = 'weak_finer_PODGreedy_fn.mat';

coarse_Nmax = 85;
fine_Nmax = 100;
finer_Nmax = 100;

switch step
 case 1 % basis generation by strong POD-Greedy

  disp('--------------------------------------')
  disp('performing strong POD-Greedy procedure');
  model.RB_numintervals = [2,2,2];
  % compute 150 basis vectors
  model.RB_stop_Nmax = coarse_Nmax;
  model.RB_stop_epsilon = 0;
  model.RB_error_indicator = 'projection-error';
  model.RB_detailed_train_savepath = ...
      fullfile('PODgreedy_exp_train');
  % delete if existing
  if exist(fullfile(rbmatlabtemp,model.RB_detailed_train_savepath))
    disp('deleting old training data for correct timing measurement');
    delete(fullfile(rbmatlabtemp,model.RB_detailed_train_savepath,'*.*'));
    rmdir(fullfile(rbmatlabtemp,model.RB_detailed_train_savepath));
  end;
  
  model_data = gen_model_data(model);
  tic;
  detailed_data = filecache_function(@gen_detailed_data,model,model_data);
%  detailed_data.RB_info.computation_time
  t = toc;
  disp(['basis_generation time: ',num2str(t),' sec.']);
  detailed_data.RB_info.basis_generation_time = t;
  save(strong_PODGreedy_fn,'model','detailed_data');
  
 case 2 % basis generation by weak coarse and weak fine POD-Greedy
  
  % using beta = const, gamma = const
  
  disp('-------------------------------------------')
  disp('performing weak coarse POD-Greedy procedure');
  % compute 125 trajectories
  model.RB_numintervals = [2,2,2];
  % compute 150 basis vectors
  model.RB_stop_Nmax = fine_Nmax;
  model.RB_stop_epsilon = 0;
  model.RB_error_indicator = 'estimator';
  model.error_norm = 'l2l2';
  model.inf_sup_beta_lower_bound = @(model) 1; % note, this is wrong!
  disp('warning: inf-sup-constant set to 1, unrealistic! ');
  
  model_data = gen_model_data(model);
  tic;
  detailed_data = filecache_function(@gen_detailed_data,model,model_data);
%  detailed_data.RB_info.computation_time
  t = toc;
  disp(['basis_generation time: ',num2str(t),' sec.']);
  detailed_data.RB_info.basis_generation_time = t;
  save(weak_coarse_PODGreedy_fn,'model','detailed_data');

  % fine set:
  disp('-------------------------------------------')
  disp('performing weak fine POD-Greedy procedure');
  detailed_data = [];
  model.RB_numintervals = [4,4,4];
  tic;
  detailed_data = gen_detailed_data(model,model_data);
%  detailed_data.RB_info.computation_time
  t = toc;
  disp(['basis_generation time: ',num2str(t),' sec.']);
  detailed_data.RB_info.basis_generation_time = t;
  save(weak_fine_PODGreedy_fn,'model','detailed_data');

 case 2.1
  
  % using beta = const, gamma = const
  
  disp('-------------------------------------------')
  disp('performing weak finer POD-Greedy procedure');
  % compute 9^3 trajectories
  model.RB_numintervals = [6,6,6];
  model.RB_stop_Nmax = finer_Nmax;
  model.RB_stop_epsilon = 0;
  model.RB_error_indicator = 'estimator';
  model.error_norm = 'l2l2';
  model.inf_sup_beta_lower_bound = @(model) 1; % note, this is wrong!
  disp('warning: inf-sup-constant set to 1, unrealistic! ');
  
  model_data = gen_model_data(model);
  tic;
  detailed_data = gen_detailed_data(model,model_data);
%  detailed_data.RB_info.computation_time
  t = toc;
  disp(['basis_generation time: ',num2str(t),' sec.']);
  detailed_data.RB_info.basis_generation_time = t;
  save(weak_finer_PODGreedy_fn,'model','detailed_data');

 case 3 % compare offline runtimes
  
  disp('--------------------------------------')
  disp('comparing runtimes');
  tmp = load(strong_PODGreedy_fn);
  tstrong = tmp.detailed_data.RB_info.basis_generation_time
  tmp = load(weak_coarse_PODGreedy_fn);
  tweakcoarse = tmp.detailed_data.RB_info.basis_generation_time
  tmp = load(weak_fine_PODGreedy_fn);
  tweakfine = tmp.detailed_data.RB_info.basis_generation_time    
  tmp = load(weak_finer_PODGreedy_fn);
  tweakfiner = tmp.detailed_data.RB_info.basis_generation_time    
  
 case 4 % determine relative RB-test-error decay: 50 random points
  
  % 50 random trajectories
  Ntest = 50;

  RandStream.setDefaultStream(RandStream('mt19937ar','seed',100000));
  Mtest = rand_uniform(Ntest,model.mu_ranges);
  test_savepath = 'PODGreedy_test';
  
  fns = {strong_PODGreedy_fn,weak_coarse_PODGreedy_fn,...
         weak_fine_PODGreedy_fn, weak_finer_PODGreedy_fn};
  Nmax = [coarse_Nmax,coarse_Nmax,fine_Nmax,finer_Nmax];
  
  model_data = gen_model_data(model);

  test_errs = cell(1,length(fns));
  Nss = cell(1,length(fns));
  for s = 1:length(fns)
    Nss{s} = 5:5:Nmax(s);
    Nss{s} = [1,Nss{s}];
    tmp = load(fns{s});
    Nrb = size(tmp.detailed_data.RB,2);
    if Nss{s}(end)~=Nrb
      Nss{s} = [Nss{s}, Nrb];
    end;
    test_errs{s} = ... 
        filecache_function(@compute_test_errs,fns{s},Nss{s},Mtest,test_savepath);    
  end;
  save('PODGreedy_test_step4');
  
 case 4.5
  
  linestyles = {'b','r','r','r'};
  load('PODGreedy_test_step4');

  % plot results
%  for s = 1:length(fns);
  for s = 3;
    plot(Nss{s},test_errs{s},linestyles{s});
    hold on;
  end;
%  legend({'strong POD-Greedy, coarse','weak POD-Greedy, coarse',...
%         'weak POD-Greedy, fine', 'weak POD-Greedy, finer'});
  legend({'weak POD-Greedy'});
  axis tight;
  
  set(gca,'Yscale','log');
  disp('please adjust plot');
  title('maximum relative test rb-error','Fontsize',15);
  xlabel('dimension N','Fontsize',15);
  ylabel('rb-error','Fontsize',15);
  c = get(gca,'children');
  set(c(:),'Linewidth',2)
  set(gca,'Fontsize',15)
  save('PODGreedy_test')
  saveas(gcf,'PODGreedy_relative_test_rb_error.png')
  keyboard;

  figure;
  load(fns{3});
  plot_detailed_data(model,detailed_data);

  plot_params.no_lines = 1;
  plot_params.show_colorbar = 0;
  plot_params.axis_equal = 1;
  for i=97:100
    model.plot(detailed_data.grid,detailed_data.RB(:,i),plot_params); 
    axis off
    saveas(gcf,['basis_func',num2str(i),'.png']);
    pause;
  end;
  
  
 case 5 % determine relative projection-test-error decay: 50 random points
  
  Ntest = 50;
  model_data = gen_model_data(model);
  
  RandStream.setDefaultStream(RandStream('mt19937ar','seed',100000));
  Mtest = rand_uniform(Ntest,model.mu_ranges);
  test_savepath = 'PODGreedy_test';
  
%  fns = {strong_PODGreedy_fn,weak_coarse_PODGreedy_fn,...
%        weak_fine_PODGreedy_fn, weak_finer_PODGreedy_fn};
%  Nmax = [coarse_Nmax,coarse_Nmax,fine_Nmax,finer_Nmax];
  fns = {weak_fine_PODGreedy_fn};
  Nmax = [fine_Nmax];
  linestyles = {'r'};
  
  test_errs = cell(1,length(fns));
  Nss = cell(1,length(fns));
  for s = 1:length(fns)
    Nss{s} = 5:5:Nmax(s);
    Nss{s} = [1,Nss{s}];
    tmp = load(fns{s});
    Nrb = size(tmp.detailed_data.RB,2);
    if Nss{s}(end)~=Nrb
      Nss{s} = [Nss{s}, Nrb];
    end;
    test_errs{s} = ... 
        filecache_function(@compute_test_projection_errs,...
                           model,model_data,fns{s},Nss{s},Mtest);    
  end;
  save('PODGreedy_test_step5');
  
 case 5.5
  
  load('PODGreedy_test_step5');
  % plot results
  for s = 1:length(fns);
    plot(Nss{s},test_errs{s},linestyles{s});
    hold on;
  end;
  %  legend({'strong POD-Greedy, coarse','weak POD-Greedy, coarse',...
  %       'weak POD-Greedy, fine', 'weak POD-Greedy, finer'});  
  legend({'weak POD-Greedy'});  
  set(gca,'Yscale','log');
  c = get(gca,'children');
  set(c(:),'Linewidth',2)
  set(gca,'Fontsize',15)
  title('maximum relative test projection error ','Fontsize',15);
  xlabel('dimension N','Fontsize',15);
  ylabel('projection error','Fontsize',15);
  disp('please adjust plot');
  save('PODGreedy_projection_test')
  saveas(gcf,'PODGreedy_relative_test_projection_error.png')
  keyboard;
    
  %case 6 % determine train-error decay: only reasonable if identical
  %       % train set
  
  % Ns = 5:5:150;
  % Ns = [1,Ns];
  
 otherwise 
  error('step unknown')
end;

function errs = compute_test_errs(fn,Ns,Mtest,test_savepath)
tmp = load(fn);
model = tmp.model;
model.relative_error = 1;
model.error_norm = 'l2l2';
model.inf_sup_beta_lower_bound = @(model) 1; % note, this is wrong!
detailed_data = tmp.detailed_data;
reduced_data = gen_reduced_data(model,detailed_data);
errs = zeros(1,length(Ns));
for Ni = 1:length(Ns)
  tmp_model = model;
  tmp_model.N = Ns(Ni);
  tmp_reduced_data = model.reduced_data_subset(tmp_model,reduced_data);
  test_errs = rb_test_error(model,detailed_data,...
                            tmp_reduced_data,Mtest,test_savepath);
  errs(Ni) = max(test_errs);
end;  


function errs = compute_test_projection_errs(model,model_data,fn,Ns,Mtest);
tmp = load(fn);
detailed_data = tmp.detailed_data;

%model.relative_error = 1;
%model.error_norm = 'l2l2';
%model.inf_sup_beta_lower_bound = @(model) 1; % note, this is wrong!
%detailed_data = tmp.detailed_data;
%reduced_data = gen_reduced_data(model,detailed_data);
errs = zeros(1,length(Ns));
for i = 1:length(Mtest);
  model = set_mu(model,Mtest(:,i));
  sim_data = filecache_function(@detailed_simulation,model, ...
                                model_data);
  test_errs = zeros(length(Ns),length(Mtest));
  for Ni = 1:length(Ns)
    proj_sol = detailed_data.RB(:,1:Ns(Ni))*...
        (detailed_data.RB(:,1:Ns(Ni))' * detailed_data.W) * sim_data.U;
    er = fv_l2l2_error(sim_data.U,proj_sol, ...
                       detailed_data.W,model);
    test_errs(Ni,i) = er(end);
  end;  
end;  
errs = max(test_errs,[],2);