gdl_Nadapt.m

Go to the documentation of this file.

% dgl_adaptivity
%
% script performing the experiments for the N-adaptivity with the
% extended error estimator

% Bernard Haasdonk 7.11.2008

%step = 1 % generate figures
%step = 2 % perform time measurements and plot time over error/estimator/Nav
%step = 3 % profiling non-adaptive
%step = 4 % profiling adaptive
%step = 5 % simulation with longer time
step = 6 % new figures ignoring error

disp('preparing simulations:');
load demo_lin_evol_detailed_data;
offline_data = rb_offline_prep(detailed_data,params);

%params.N_adaptivity = 'off';
%params.N = 10 % set max!
%params.error_norm = 'l1l2';

params.N_mode = 'adapt';
params.N = 100; % set max!
%params.N = size(detailed_data.RB,2); % set max!
%params.error_norm = 'l1l2';
params.error_norm = 'l2';

% the following illustrates non-rigorous error-estimator for adaptivity!!!!
params.N_adapt_epsilon = 1e-4;
%params.N_adapt_epsilon = 2e-5;
%params.N_adapt_epsilon = 1e-8;
params.N_adapt_look_ahead_nt = 20;
params.N_adapt_slope = 10;
%  params.N_adapt_epsilon = 1e-3; => everything over 50% N
%params = set_mu([1,1,0],params); % convection
params = set_mu([1,0,1e-9],params); % convection 
%params = set_mu([1,0,0],params); % convection => super "prediction"
%params = set_mu([1,1,5e-8],params); % diffusion
%params = set_mu([0,1,5e-8],params); % diffusion
reduced_data = rb_online_prep(offline_data, params);

switch step 
  case 1
%1. N adaptivity plot
disp('detailed_simulation:');
%params = set_mu([1,0,1e-9],params); % convection 
params = set_mu([1,0,0],params); % convection 
params.N_adapt_epsilon = 5e-5;
params.N_adapt_look_ahead_nt = 10;
params.N_adapt_slope = 5;
U = detailed_simulation(detailed_data.grid,params);
disp('reduced_simulation:');
simulation_data = rb_simulation(reduced_data,params);
Uappr = rb_reconstruction(detailed_data,simulation_data); 
err = fv_error(Uappr,U,detailed_data.W,params);
%  params.title = 'reduced simulation result';  
%  plot_element_data_sequence(detailed_data.grid,Uappr,params);
if isfield(simulation_data,'N'); 
  Ns = simulation_data.N;
else
  Ns = params.N * ones(1,params.nt+1);
end;
figure; % subplot(2,1,2); 
plot([0:params.nt],Ns);
title('model dimension N');
set(gca,'XLim',[0,params.nt]);
set(gca,'YLim',[0,100]);
xlabel('time index k');
set(gcf,'Position',[472   404   564   238]);
saveas(gcf,'Nadaptive_error','epsc');
%subplot(2,1,1), 
figure, plot([0:params.nt],[simulation_data.Delta;...
              0:(params.N_adapt_epsilon)/params.nt:params.N_adapt_epsilon;err]');
legend({'estimator','target','error'},'location','west');
title('error estimator evolution');
set(gca,'XLim',[0,params.nt]);
xlabel('time index k');
set(gcf,'Position',[472   404   564   238]);
disp('please save figure as eps (automatic saving does not get aspectratio)')
%saveas(gcf,'Nadaptive_Ns','epsc');

%plot_element_data_sequence(detailed_data.grid,Uappr,params);
%plot_element_data_sequence(detailed_data.grid,Uappr-U,params);

params.N_mode = 'fixed';
params.N = 40;
reduced_data = rb_online_prep(offline_data, params);
disp('reduced_simulation:');
simulation_data = rb_simulation(reduced_data,params);
Uappr = rb_reconstruction(detailed_data,simulation_data); 
err = fv_error(Uappr,U,detailed_data.W,params);
%  params.title = 'reduced simulation result';  
%  plot_element_data_sequence(detailed_data.grid,Uappr,params);
if isfield(simulation_data,'N'); 
  Ns = simulation_data.N;
else
  Ns = params.N * ones(1,params.nt+1);
end;
figure, subplot(2,1,2); plot([0:params.nt],Ns);
title('model dimension N');
set(gca,'XLim',[0,params.nt]);
set(gca,'YLim',[0,100]);
xlabel('time index k');
subplot(2,1,1), plot([0:params.nt],[simulation_data.Delta;...
              0:(params.N_adapt_epsilon)/params.nt:params.N_adapt_epsilon;err]');
legend({'estimator','target','error'},'location','northwest');
title('error estimator evolution');
set(gca,'XLim',[0,params.nt]);
xlabel('time index k');
disp('please save figure as eps (automatic saving does not get aspectratio)')


case 2
% time measurement:

params.N = 100;

nrep = 2;
t_exact = zeros(1,nrep);
disp('detailed_simulation:');
for r = 1:nrep
  tic, U = detailed_simulation(detailed_data.grid,params);
  t_exact(r) = toc;
end;
disp('                      time      av(N)     Delta       err ');
disp(['Exact simulation: ', ...
      num2str(mean(t_exact)), '+-',num2str(std(t_exact)),...
      '   ', num2str(detailed_data.grid.nelements),' 0 0 ']);

disp('reduced_simulations, fixed Ns:');
Ns = [100,85,70,55,40,25,10];
t_Nfix = zeros(length(Ns),nrep);
Delta_Nfix  = zeros(size(Ns));
% prevent matrix inversion for fair comparison
params.skip_matrix_inversion = 1; 
params.N_mode = 'fixed';

%dummy simulation
dummy_reduced_data = rb_online_prep(offline_data, params); 
dummy = rb_simulation(dummy_reduced_data,params);

for r = 1:nrep
for n = 1:length(Ns) 
  simulation_data = [];
  params.N = Ns(n);
  reduced_data = rb_online_prep(offline_data, params); 
  tic;
  simulation_data = rb_simulation(reduced_data,params);
  t_Nfix(n,r)= toc;
  Uappr = rb_reconstruction(detailed_data,simulation_data); 
  err = fv_error(Uappr,U,detailed_data.W,params);  
  err_Nfix(n) = max(err);
  Delta_Nfix(n) = simulation_data.Delta(end);
%  disp(['N = ',num2str(Ns(n)),'  ',num2str(t_Nfix(n,r)),'    ',...
%       num2str(Ns(n)),'    ',...
%       num2str(Delta_Nfix(n)),'    ',num2str(err_Nfix(n))]);
end;
end;
disp('overall:')
for n = 1:length(Ns) 
  disp(['N = ',num2str(Ns(n)),'  ',num2str(mean(t_Nfix(n,:))),...
        '+-',num2str(std(t_Nfix(n,:))),'   ',...
        num2str(Ns(n)),'    ',...
        num2str(Delta_Nfix(n)),'    ',num2str(err_Nfix(n))]);
end;

disp('reduced_simulations, adaptive Ns:');
epsilons = [5e-6,1e-5,5e-5,1e-4,5e-4,1e-3];
Delta_Nadaptive  = zeros(size(epsilons));
t_Nadaptive = zeros(length(epsilons),nrep);
Nav_Nadaptive = zeros(size(epsilons));
params.N = 100; %size(detailed_data.RB,2);
reduced_data = rb_online_prep(offline_data, params); 
params.N_mode = 'adapt';

for r = 1:nrep
for e = 1:length(epsilons) 
  simulation_data =[];
  params.N_adapt_epsilon = epsilons(e);
  tic;
  simulation_data = rb_simulation(reduced_data,params);
  t_Nadaptive(e,r)= toc;
  Uappr = rb_reconstruction(detailed_data,simulation_data); 
  err = fv_error(Uappr,U,detailed_data.W,params);  
  err_Nadaptive(e) = max(err);
  Delta_Nadaptive(e) = simulation_data.Delta(end);
  Nav_Nadaptive(e) = mean(simulation_data.N);
%  disp(['epsilon = ',num2str(epsilons(e)),...
%       '  ',num2str(t_Nadaptive(e,r)),'    ',num2str(Nav_Nadaptive(e)),...
%       '    ',...
%       num2str(Delta_Nadaptive(e)),'    ',num2str(err_Nadaptive(e))]);
end;
end;
disp('overall:')
for e = 1:length(epsilons) 
  disp(['epsilon = ',num2str(epsilons(e)),...
        '  ',num2str(mean(t_Nadaptive(e,:))),'+-',...
        num2str(std(t_Nadaptive(e,:))),'    ',num2str(Nav_Nadaptive(e)),...
        '    ',...
        num2str(Delta_Nadaptive(e)),'    ',num2str(err_Nadaptive(e))]);
end;

  figure
  plot(mean(t_Nadaptive,2),Delta_Nadaptive','b-');
  hold on;
  plot(mean(t_Nfix,2),Delta_Nfix','g-');
  hold on;
%  plot(mean(t_Nswitch,2),Delta_Nswitch','r-');
%  hold on;
  xlabel('runtime')
  ylabel('error estimator')
  set(gca,'Yscale','log')
  title('error estimator over runtime')
%  legend('N-adaptive','N-fixed','N-switch')
  legend('N-adaptive','N-fixed')
  
  figure
  plot(mean(t_Nadaptive,2),err_Nadaptive','b-');
  hold on;
  plot(mean(t_Nfix,2),err_Nfix','g-');
%  hold on;
%  plot(mean(t_Nswitch,2),err_Nswitch','r-');
  xlabel('runtime')
  ylabel('error')
  set(gca,'Yscale','log')
  title('error over runtime')
%  legend('N-adaptive','N-fixed','N-switch')
 legend('N-adaptive','N-fixed')

  figure
  plot(Nav_Nadaptive,err_Nadaptive','b-');
  hold on;
  plot(Ns,err_Nfix','g-');
%  hold on;
%  plot(Nav_Nswitch,err_Nswitch','r-');
  xlabel('average N')
  ylabel('error')
  set(gca,'Yscale','log')
  title('error over average N')
%  legend('N-adaptive','N-fixed','N-switch')
  legend('N-adaptive','N-fixed')

 case 3   % profiling non-adaptivity
  % => 56.80 fuer bicgstab, 4.17 fuer res_norm_sqr
  
  profile clear;
  profile on
  disp('without adaptivity:');
  params.N_mode = 'fixed';
  params.skip_matrix_inversion = 1; 
  tic
  for i = 1:100
    disp(i)
    simulation_data = rb_simulation(reduced_data,params);
  end;
  toc
  profview
  profsave(profile('info'),'profile_results_non_adaptive')


 case 4  % profiling non-adaptivity
  % => 39.19 sec fuer bicgstab, 3.57 fuer res_norm_sqr
  profile clear;
  profile on
  disp('with adaptivity:');
  params.N_mode = 'adapt';
  tic
  for i = 1:100
    disp(i)
    simulation_data = rb_simulation(reduced_data,params);
  end;
  toc
  profview
  profsave(profile('info'),'profile_results_adaptive')
 
 case 5  % test: adaptive and fixed N should give identical result!

params.N = 100;
disp('detailed_simulation:');
U = detailed_simulation(detailed_data.grid,params);
% prevent matrix inversion for fair comparison
params.N_mode = 'fixed';
params.skip_matrix_inversion = 1; 
%dummy simulation
simulation_data_fix = [];
reduced_data = rb_online_prep(offline_data, params); 
simulation_data_fix = rb_simulation(reduced_data,params);
Uappr_fix = rb_reconstruction(detailed_data,simulation_data_fix); 
err_fix = fv_error(Uappr_fix,U,detailed_data.W,params);  
disp(['error fixed N=100: ',num2str(err_fix(end))]);

disp('reduced_simulations, adaptive Ns:');
params.N = 100; %size(detailed_data.RB,2);
params.N_mode = 'adapt';
reduced_data = rb_online_prep(offline_data, params); 
simulation_data_adapt =[];
params.N_adapt_epsilon = 1e-6;
simulation_data_adapt = rb_simulation(reduced_data,params);
Uappr_adapt = rb_reconstruction(detailed_data,simulation_data_adapt); 
err_adapt = fv_error(Uappr_adapt,U,detailed_data.W,params);  
disp(['error adapt N=100: ',num2str(err_adapt(end))]);

 case 6 % new figures: no error....

linestyles = {'-.','-','-',':'};  
linecolors = {[0,0,1],[0,0.5,0],[1,0,0],[0,0,0]};  
linewidths = [1,1,2,2];

params = set_mu([1,0,0],params); % convection 
params.N_adapt_epsilon = 5e-5;
params.N_adapt_look_ahead_nt = 20;
params.N_adapt_slope = 10;
disp('reduced_simulation:');
simulation_data = rb_simulation(reduced_data,params);

params.N_mode = 'fixed';
params.N = 20;
reduced_data = rb_online_prep(offline_data, params);
disp('reduced_simulation:');
simulation_data_N1 = rb_simulation(reduced_data,params);
params.N = 60;
reduced_data = rb_online_prep(offline_data, params);
disp('reduced_simulation:');
simulation_data_N2 = rb_simulation(reduced_data,params);
%Uappr = rb_reconstruction(detailed_data,simulation_data); 
if isfield(simulation_data,'N'); 
  Ns = simulation_data.N;
else
  Ns = params.N * ones(1,params.nt+1);
end;
figure; % subplot(2,1,2); 
p = plot([0:params.nt],[ones(size(Ns))*20;ones(size(Ns))*60;Ns]');
for i=1:length(p)
  set(p(i),'Color',linecolors{i},...
           'LineWidth',linewidths(i),...
           'LineStyle',linestyles{i})
end;
%set(p(3),'LineWidth',2);
%set(p(1),'LineStyle','-.');
title('model dimension N');
legend('N-fixed N=20','N-fixed N=60','N-adaptive')
set(gca,'XLim',[0,params.nt]);
set(gca,'YLim',[0,100]);
xlabel('time index k');
saveas(gcf,'Nadaptive_error_large','epsc');
set(gcf,'Position',[472   404   564   238]);
%subplot(2,1,1), 
figure, p = plot([0:params.nt],[simulation_data_N1.Delta; ...
                           simulation_data_N2.Delta;...
                           simulation_data.Delta;...
                    0:(params.N_adapt_epsilon)/params.nt:params.N_adapt_epsilon
                   ]');
%keyboard;
for i=1:length(p)
  set(p(i),'Color',linecolors{i},...
           'LineWidth',linewidths(i),...
           'LineStyle',linestyles{i})
end;
%keyboard;
%set(p(1),'LineStyle','-.');
%set(p(4),'LineStyle',':');
%set(p(3),'LineWidth',2);
title('error estimator evolution');
l = legend({'N-fixed, N=20',    'N-fixed, N=60', ...
        'N-adaptive', ...
        'N-adaptive-target'},'location','northwest');
set(gca,'XLim',[0,params.nt]);
xlabel('time index k');
saveas(gcf,'Nadaptive_Ns_large','epsc');
set(gcf,'Position',[472   404   564   238]);
po = get(l,'Position');
set(l,'Position',po);
disp('please save figures as eps (automatic saving does not get aspectratio)')

 otherwise
  error('step unknown');
end;
%| \docupdate