rbmatlab  1.16.09
 All Classes Namespaces Files Functions Variables Modules Pages
rb_tutorial.m
1 function res = rb_tutorial(step)
2 %function res = rb_tutorial(step)
3 %
4 % step : index of experiment to be performed
5 % 1 - 10: Thermal block, elliptic PDE
6 % 11 - 17: Advection-Diffusion, parabolic PDE
7 %
8 % This program allows to reproduce the thermalblock and advection-diffusion
9 % experiments of the RB-Summerschool and the pictures are used in the
10 % RB-Tutorial book chapter.
11 %
12 % B. Haasdonk: Reduced Basis Methods for Parametrized PDEs --
13 % A Tutorial Introduction for Stationary and Instationary Problems.
14 % Chapter in P. Benner, A. Cohen, M. Ohlberger and K. Willcox (eds.):
15 % "Model Reduction and Approximation: Theory and Algorithms", SIAM, Philadelphia, 2016
16 %
17 % A Preprint version of this chapter is available at:
18 %
19 % http://www.simtech.uni-stuttgart.de/publikationen/prints.php?ID=938
20 
21 % B. Haasdonk, 7.9.2013
22 
23 if nargin<1
24  step = 1
25  help rb_tutorial;
26  disp('now running step 1');
27 end;
28 
29 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
30 %PART 1: Stationary Thermalblock example
31 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
32 
33 switch step
34  case 1
35  disp('thermal block for different configurations and corresponding plots');
36  params = [];
37  params.B1 = 2;
38  params.B2 = 2;
39  params.mu_range = [0.1;10];
40  params.numintervals_per_block = 5;
41  model = thermalblock_model_struct(params);
42  model_data = gen_model_data(model);
43  plot_params = [];
44  plot_params.axis_equal = 1;
45  plot_params.axis_tight = 1;
46 
47  % plot 1: four identical values
48  figure;
49  model = set_mu(model,[1,1,1,1]);
50  sim_data = detailed_simulation(model,model_data);
51  plot_sim_data(model,model_data,sim_data,plot_params);
52 
53  % plot 2: lower and upper values identical
54  figure;
55  model = set_mu(model,[0.1,0.1,10,10]);
56  sim_data = detailed_simulation(model,model_data);
57  plot_sim_data(model,model_data,sim_data,plot_params);
58 
59  % plot 3: four different values
60  figure;
61  model = set_mu(model,[0.3,0.2,0.2,0.7]);
62  sim_data = detailed_simulation(model,model_data);
63  plot_sim_data(model,model_data,sim_data,plot_params);
64 
65  % plot 4: larger division
66  % mu = rand(params.B1*params.B2,1)+0.1;
67  figure;
68  mu = [ 0.8094 0.8547 0.3760 0.7797 0.7551 0.2626, ...
69  0.2190 0.5984 1.0597 0.4404 0.6853 0.3238, ...
70  0.8513 0.3551 0.6060 0.7991 0.9909 1.0593 ...
71  0.6472 0.2386 0.2493 0.3575 0.9407 0.3543 ...
72  0.9143 0.3435 1.0293 0.4500 0.2966 0.3511 ...
73  0.7160 0.5733 0.4517 0.9308 0.6853 0.6497]';
74  params = [];
75  params.B1 = 6;
76  params.B2 = 6;
77  params.mu_range = [0.1;10];
78  model = thermalblock_model_struct(params);
79  model_data = gen_model_data(model);
80  model = set_mu(model,mu);
81  sim_data = detailed_simulation(model,model_data);
82  plot_sim_data(model,model_data,sim_data,plot_params);
83 
84  case 2
85  disp('thermal block detailed gui');
86  params = [];
87  params.B1 = 2;
88  params.B2 = 2;
89  params.numintervals_per_block = 50;
90  params.mu_range = [0.1;10];
91  params.numintervals_per_block = 5;
92  model = thermalblock_model_struct(params);
93  model_data = gen_model_data(model);
94  plot_params = [];
95  plot_params.axis_equal = 1;
96  plot_params.axis_tight = 1;
97  plot_params.yscale_uicontrols = 0.7;
98  demo_detailed_gui(model,model_data,[],plot_params);
99  set(gcf,'Position',[584 245 553 691]);
100 
101  case 3
102  disp('equidistant samples for basis generation, offline versus online time');
103 
104  params = [];
105  params.B1 = 2;
106  params.B2 = 2;
107  params.numintervals_per_block = 20;
108  params.mu_range = [0.1;10];
109  disp('initializing model')
110  model = thermalblock_model_struct(params);
111 
112  disp('-----------------------------------------------------')
113  disp('OFFLINE PHASE:');
114  disp('generating model data: grid, fem inner product matrix, etc.');
115  tic
116  model_data = gen_model_data(model);
117  model.RB_generation_mode = 'lagrangian';
118  model.RB_mu_list = lin_equidist_samples(model,[4,2,2,2]);
119  disp('generating detailed data: A_q,f_q,l_q components, reduced basis.');
120  detailed_data = gen_detailed_data(model,model_data);
121  disp('generating reduced data: A_Nq,f_Nq,l_Nq components');
122  reduced_data = gen_reduced_data(model, detailed_data);
123  t_offline = toc;
124  disp('Offline time: ')
125  t_offline
126  save('thermalblock_2x2_data');
127 
128  disp('-----------------------------------------------------')
129  disp('ONLINE PHASE:');
130  % online phase for single simulation:
131  disp('reduced simulation: assembling system and solve');
132  model = set_mu(model,[1,1,1,1]);
133  tic
134  sim_data = rb_simulation(model,reduced_data);
135  t_online = toc;
136 
137  disp('Online time: ')
138  t_online
139 
140  disp('starting reduced basis gui:')
141  plot_params = [];
142  plot_params.axis_equal = 1;
143  plot_params.axis_tight = 1;
144  plot_params.yscale_uicontrols = 0.7;
145  demo_rb_gui(model,detailed_data,reduced_data,plot_params);
146 
147  case 4
148  disp('1-parameter example: error, estimator, effectivity bound plot.');
149 
150  % generate reduced model
151  disp('generating reduced model')
152  params = [];
153  params.B1 = 2;
154  params.B2 = 2;
155  params.numintervals_per_block = 20;
156  params.mu_range = [0.1;10];
157  model = thermalblock_model_struct(params);
158  model_data = gen_model_data(model);
159  model.RB_generation_mode = 'lagrangian';
160  c = 0.1;
161  mu_list = {[0.1,c,c,c],...
162  [0.5,c,c,c],...
163  [0.9,c,c,c],...
164  [1.3,c,c,c],...
165  [1.7,c,c,c]};
166  model.RB_mu_list = mu_list;
167  detailed_data = gen_detailed_data(model,model_data);
168  reduced_data = gen_reduced_data(model,detailed_data);
169 
170  disp('computing parameter sweep')
171  mu1s = linspace(0.1,4,1000);
172  % rapidly computable error landscape:
173  Deltas = zeros(length(mu1s),1);
174  for n = 1:length(mu1s)
175  fprintf('.');
176  mu = [mu1s(n),c,c,c];
177  model = set_mu(model,mu);
178  sim_data = rb_simulation(model,reduced_data);
179  Deltas(n) = sim_data.Delta;
180  end;
181  fprintf('\n');
182  plot(mu1s,Deltas,'Linewidth',2);
183  xlabel('mu1','Fontsize',15);
184  ylabel('error/estimator','Fontsize',15);
185  title('error estimator by parameter sweep','Fontsize',15);
186  set(gca,'Yscale','log');
187  legend('estimator \Delta_u(\mu)');
188 
189  % comparison with expensive true error
190 
191  disp('computing true error for some parameters')
192  mu1s = linspace(0.1,4,40);
193  errs = zeros(length(mu1s),1);
194  for n = 1:length(mu1s)
195  fprintf('.');
196  mu = [mu1s(n),c,c,c];
197  model = set_mu(model,mu);
198  sim_data = detailed_simulation(model,model_data);
199  rb_sim_data = rb_simulation(model,reduced_data);
200  % linear combination of reduced vector with basis:
201  rb_sim_data= rb_reconstruction(model,detailed_data,rb_sim_data);
202  err = sim_data.uh;
203  err.dofs = err.dofs - rb_sim_data.uh.dofs;
204  errs(n) = fem_h10_norm(err);
205  end;
206  fprintf('\n');
207  hold on;
208  plot(mu1s,errs,'r*','Markersize',10);
209  legend({'estimator \Delta_u(\mu)','true error |u-u_N|'},'Fontsize',15);
210 
211  %set(gcf,'Position', [379 492 434 328]);
212  %saveas(gcf,'error_parameter_sweep.png');
213  %saveas(gcf,'error_parameter_sweep.eps','epsc');
214 
215  save('rb_tutorial_step4');
216 
217  case 5
218  disp('Effectivity plot');
219  load('rb_tutorial_step4');
220 
221  disp('computing parameter sweep')
222  mu1s = linspace(0.1,4,40);
223  % rapidly computable error landscape:
224  Deltas = zeros(length(mu1s),1);
225  for n = 1:length(mu1s)
226  fprintf('.');
227  mu = [mu1s(n),0.1,0.1,0.1];
228  model = set_mu(model,mu);
229  sim_data = rb_simulation(model,reduced_data);
230  Deltas(n) = sim_data.Delta;
231  end;
232  fprintf('\n');
233 
234  etas = Deltas.*(errs.^(-1));
235  i = find(errs<1e-7);
236  %i = find(etas>20);
237  etas(i)=NaN;
238  plot(mu1s,etas,'b*','Markersize',10);
239  xlabel('\mu1','Fontsize',15);
240  ylabel('effectivity','Fontsize',15);
241  title('effectivity over parameter','Fontsize',15);
242 % set(gca,'Yscale','log');
243 
244  hold on;
245  gammas = mu1s;
246  alpha = 0.1;
247  effectivities = gammas/alpha;
248  plot(mu1s,effectivities,'r-','Linewidth',2);
249 
250  legend({'effectivity \Delta_u/|u-u_N|',...
251  'upper bound \gamma(\mu)/\alpha(\mu)'},...
252  'Fontsize',15,...
253  'Location','NorthWest');
254 
255 % set(gcf,'Position', [379 492 434 328]);
256 % saveas(gcf,'effectivity_parameter_sweep.png');
257 % saveas(gcf,'effectivity_parameter_sweep.eps','epsc');
258 
259  case 6
260  disp('Convergence plot for test error on equidistant samples');
261 
262  params = [];
263  params.B1 = 3;
264  params.B2 = 3;
265  params.numintervals_per_block = 20;
266  params.mu_range = [0.5;2];
267  model = thermalblock_model_struct(params);
268  model_data = gen_model_data(model);
269  model.RB_generation_mode = 'model_RB_basisgen';
270  model.RB_basisgen = @lagrangian_orth_basisgen;
271  Ns = 2:8; % linear dependent for higher values...
272 
273  maxDeltas = zeros(length(Ns),1);
274  maxerrs = zeros(length(Ns),1);
275 
276  disp('Precomputing test snapshots');
277  ntest = 100;
278  [U,test_mus] = filecache_function(@precompute_test_snapshots,...
279  params,model,model_data,ntest);
280 
281  disp('Generating basis and determining max error, estimators');
282  tic;
283  for Ni = 1:length(Ns)
284  N = Ns(Ni)
285  train_mu1s = linspace(params.mu_range(1),params.mu_range(2),N);
286  train_mus = [train_mu1s(:), ones(N,params.B1*params.B2-1)*1];
287  model.RB_mu_list = mat2cell(train_mus,ones(N,1),params.B1*params.B2);
288  detailed_data = gen_detailed_data(model,model_data);
289  reduced_data = gen_reduced_data(model,detailed_data);
290  for i = 1:ntest
291  model = set_mu(model,test_mus(i,:));
292  rb_sim_data = rb_simulation(model, reduced_data);
293  rb_sim_data = rb_reconstruction(model,detailed_data,rb_sim_data);
294  err = rb_sim_data.uh;
295  err.dofs = err.dofs - U(:,i);
296  errnorm = fem_h10_norm(err);
297  if maxerrs(Ni)<errnorm
298  maxerrs(Ni)=errnorm;
299  end;
300  if maxDeltas(Ni)<rb_sim_data.Delta
301  maxDeltas(Ni)=rb_sim_data.Delta;
302  end;
303  end; % i
304  end; % Ni
305 
306  time_step_6 = toc
307  save('rb_tutorial_step6');
308  % case 6.1 % plot of error convergenc
309  % load('rb_tutorial_step6');
310  plot(Ns,maxDeltas,'b-','Linewidth',2);
311  hold on
312  plot(Ns,maxerrs,'r:','Linewidth',2);
313  set(gca,'Yscale','log');
314  legend({'estimator \Delta_u(\mu)','error |u(\mu)-u_N(\mu)|'},'Fontsize',15);
315  title('test error decay for growing sample size','Fontsize',15);
316  xlabel('sample/basis Size N','Fontsize',15);
317  ylabel('error/estimator','Fontsize',15);
318 
319  % set(gcf,'Position', [379 492 434 328]);
320  % saveas(gcf,'lagrangian_error_convergence.png');
321  % saveas(gcf,'lagrangian_error_convergence.eps','epsc');
322 
323  case 7
324  disp('plot basis');
325  load('rb_tutorial_step6');
326  plot_params = [];
327  plot_params.axis_equal = 1;
328  plot_params.axis_tight = 1;
329  plot_params.plot = @plot_vertex_data;
330  plot_params.title='Orthonormal Reduced Basis \Phi_N';
331  % plot as sequence
332  plot_sequence(detailed_data.RB,detailed_data.grid,plot_params);
333 
334  % or in single plot
335  plot_params.show_colorbar = 0;
336  plot_params.no_lines = 1;
337  figure;
338  for n = 1:8
339  subplot(2,4,n);
340  plot_vertex_data(detailed_data.grid,detailed_data.RB(:,n), ...
341  plot_params);
342  title(['n=',num2str(n)],'Fontsize',10);
343  set(gca,'Xtick',[]);
344  set(gca,'Ytick',[]);
345  end;
346  set(gcf,'Position',[ 32 254 1554 701 ]);
347  saveasshown(gcf,'orthonormal_basis.png');
348  saveasshown(gcf,'orthonormal_basis.eps','epsc');
349 
350  case 8
351  disp('Greedy error convergence');
352 
353  params = [];
354  params.B1 = 2;
355  params.B2 = 2;
356  params.numintervals_per_block = 20;
357  params.mu_range = [0.5;2];
358  model = thermalblock_model_struct(params);
359  model_data = gen_model_data(model);
360  model.RB_stop_epsilon = 1e-6;
361  model.RB_stop_Nmax = 50;
362  model.RB_train_size = 1000;
363 
364  disp('Generating basis and determining max error, estimators');
365  tic;
366  detailed_data = gen_detailed_data(model,model_data);
367  reduced_data = gen_reduced_data(model,detailed_data);
368 
369  disp('Precomputing test snapshots');
370  ntest = 100;
371  [U,test_mus] = filecache_function(@precompute_test_snapshots,...
372  params,model,model_data,ntest);
373 
374  Ns = 1:size(detailed_data.RB,2);
375  maxDeltas = zeros(length(Ns),1);
376  maxerrs = zeros(length(Ns),1);
377  for Ni = 1:length(Ns)
378  N = Ns(Ni)
379  model.N = N;
380  reduced_data_subset = model.reduced_data_subset(model,reduced_data);
381  for i = 1:ntest
382  model = set_mu(model,test_mus(i,:));
383  rb_sim_data = rb_simulation(model, reduced_data_subset);
384  rb_sim_data = rb_reconstruction(model,detailed_data,rb_sim_data);
385  err = rb_sim_data.uh;
386  err.dofs = err.dofs - U(:,i);
387  errnorm = fem_h10_norm(err);
388  if maxerrs(Ni)<errnorm
389  maxerrs(Ni)=errnorm;
390  end;
391  if maxDeltas(Ni)<rb_sim_data.Delta
392  maxDeltas(Ni)=rb_sim_data.Delta;
393  end;
394  end; % i
395  end; % Ni
396  time_step_6 = toc
397  save('rb_tutorial_step8');
398 
399  case 8.1
400  disp('plot of greedy errors');
401  load('rb_tutorial_step8');
402 
403  plot(Ns,maxDeltas,'b-','Linewidth',2);
404  hold on
405  plot(Ns,maxerrs,'r:','Linewidth',2);
406  set(gca,'Yscale','log');
407  % if wanted, training error can be plotted
408  % plot(Ns(1:end-2),detailed_data.RB_info.max_err_est_sequence(2:end),...
409 % 'g-.','Linewidth',2);
410  set(gca,'Yscale','log');
411  legend({'test estimator \Delta_u(\mu)','test error |u(\mu)-u_N(\mu)|'},'Fontsize',15);
412  title('Error decay for Greedy basis','Fontsize',15);
413  xlabel('sample/basis Size N','Fontsize',15);
414  ylabel('error/estimator','Fontsize',15);
415 
416 % set(gcf,'Position', [379 492 434 328]);
417 % saveas(gcf,'greedy_test_error.png');
418 % saveas(gcf,'greedy_test_error.eps','epsc');
419 
420  % possible further cases:
421  % plot of sample points: not spectacular...
422  % load('rb_tutorial_step8');
423  % mus = detailed_data.RB_info.max_mu_sequence;
424  % mumat = cell2mat(mus)
425  % plot(mumat(1,:),mumat(2,:),'*');
426 
427  case 9
428  disp('experiment with basis from 4');
429 
430  load('rb_tutorial_step4');
431 
432  mus= {[1,0.2,0.2,0.2],[1,0.1,1,0.1]};
433  for i = 1:length(mus)
434  disp('-------------------------------------------')
435  disp(['Experiment for mu = (',num2str(mus{i}),')']);
436  model = set_mu(model,mus{i});
437  sim_data = detailed_simulation(model,model_data);
438  rb_sim_data = rb_simulation(model,reduced_data);
439  rb_sim_data = rb_reconstruction(model,detailed_data,rb_sim_data);
440  disp(['output s(mu) = ', num2str(sim_data.s)]);
441  disp(['output s_N(mu) = ', num2str(rb_sim_data.s)]);
442  disp(['output error s(mu)-s_N(mu) = ', ...
443  num2str(sim_data.s-rb_sim_data.s)]);
444  disp(['output error bound Delta_s(\mu) = ', num2str(rb_sim_data.Delta_s)]);
445  err = sim_data.uh;
446  err.dofs = err.dofs - rb_sim_data.uh.dofs;
447  disp(['field error |u-u_N| = ',num2str(fem_h10_norm(err))]);
448  disp(['field error bound Delta_u(\mu) = ', num2str(rb_sim_data.Delta)]);
449  end;
450 
451  case 10
452  disp('convergence curves, training error for different block numbers');
453 
454  B1S = [2,3,4];
455  B2S = [2,3,4];
456 
457  max_err_est_sequences = {};
458 
459  for bi = 1:length(B1S)
460  params = [];
461  params.B1 = B1S(bi);
462  params.B2 = B2S(bi);
463  disp('--------------------------------------------------');
464  disp(['generating basis for (B1,B2) = (',...
465  num2str(params.B1),',',num2str(params.B2),')']);
466  params.numintervals_per_block = 20;
467  params.mu_range = [0.5;2];
468  model = thermalblock_model_struct(params);
469  model_data = gen_model_data(model);
470  model.RB_stop_epsilon = 1e-6;
471  model.RB_stop_Nmax = 50;
472  model.RB_train_size = 1000;
473  detailed_data = gen_detailed_data(model,model_data);
474  max_err_est_sequences = [max_err_est_sequences,...
475  {detailed_data.RB_info.max_err_est_sequence}]
476  end;
477 
478 % save('rb_tutorial_step10');
479 % case 10.1 % plot training error curves
480 % load('rb_tutorial_step10');
481 
482  figure;
483  hold on;
484  legstr = {};
485  styles = {'b-','r-','g-','k-'}
486  for i=1:length(max_err_est_sequences)
487  seq = max_err_est_sequences{i};
488  plot(0:length(seq)-1,seq,styles{i});
489  legstr = [legstr,{['(B1,B2)= (',...
490  num2str(B1S(i)),',',num2str(B2S(i)),')']}]
491  end;
492  legend(legstr,'Fontsize',15);
493  title('greedy training error decay','Fontsize',15);
494  xlabel('basis size N','Fontsize',15);
495  ylabel('training error','Fontsize',15);
496  set(gca,'Yscale','log');
497 
498 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
499 %PART 2: Instationary Advection-diffusion example
500 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
501 
502  case 11
503  disp('full simulation of advection-diffusion problem.');
504 
505  params = [];
506  params.coarse_factor = 4;
507  model = advection_diffusion_fv_model(params);
508  model_data = gen_model_data(model);
509  plot_params = [];
510  plot_params.axis_equal = 1;
511  plot_params.axis_tight = 1;
512  plot_params.no_lines = 1;
513 
514  model = set_mu(model,[1,0.1]);
515  tic; sim_data = detailed_simulation(model,model_data);
516  disp('time for full simulation (seconds): ')
517  t = toc
518  plot_sim_data(model,model_data,sim_data,plot_params);
519 
520  case 12
521  disp('generate plots of varying mu and t:');
522 
523  params = [];
524  params.coarse_factor = 4;
525  model = advection_diffusion_fv_model(params);
526  model_data = gen_model_data(model);
527  plot_params = [];
528  plot_params.axis_equal = 1;
529  plot_params.axis_tight = 1;
530  plot_params.no_lines = 1;
531  plot_params.show_colorbar = 0;
532  plot_params.clim = [0,1];
533 
534  mus = {[0,0],[1,0],[1,1]};
535  ts = [1,129,257];
536  snapshots = cell(3,3);
537  for mi = 1:length(mus)
538  model = set_mu(model,mus{mi});
539  sim_data = detailed_simulation(model,model_data);
540  for ti = 1:length(ts);
541  snapshots{mi,ti} = sim_data.U(:,ts(ti));
542  plot_element_data(model_data.grid,...
543  sim_data.U(:,ts(ti)),plot_params);
544  set(gca,'Ytick',[]);
545  set(gca,'Xtick',[]);
546  saveas(gcf,['advection_diffusion_snapshot-m',num2str(mi),...
547  '-t',num2str(ti),'.png']);
548  saveas(gcf,['advection_diffusion_snapshot-m',num2str(mi),...
549  '-t',num2str(ti),'.eps'],'epsc');
550  end;
551  end;
552 
553  % single joint plot
554  figure;
555  for mi = 1:length(mus)
556  for ti = 1:length(ts);
557  subplot(3,3,(ti-1)*3 + mi);
558  plot_element_data(model_data.grid,...
559  snapshots{mi,ti} ,plot_params);
560  set(gca,'Ytick',[]);
561  set(gca,'Xtick',[]);
562  if ti == 1
563  mu = mus{mi};
564  title(['\mu = (',num2str(mu(1)),',',num2str(mu(2)),')^T'],...
565  'Fontsize',15);
566  end;
567  end;
568  end;
569  pos = [354 335 1241 613]
570  set(gcf,'Position',pos);
571  disp('please do a manual => File => Export Setup => Export');
572  disp('as files advection_diffusion_snapshots.png and *.eps');
573 % saveas(gcf,'advection_diffusion_snapshots.png');
574 % saveas(gcf,'advection_diffusion_snapshots.eps','epsc');
575 
576  case 13
577  disp('POD-Greedy basis generation, caution takes +- 1 hour!');
578 
579  params = [];
580  params.coarse_factor = 4;
581  model = advection_diffusion_fv_model(params);
582  model_data = gen_model_data(model);
583 
584  % POD-Greedy settings:
585 
586  % choose initial data components as initial basis:
587  model.rb_init_data_basis = @RB_init_data_basis;
588  % or choose empty initial basis:
589  %model.rb_init_data_basis = @RB_init_basis_empty;
590 
591  model.RB_generation_mode = 'greedy_uniform_fixed';
592  model.RB_extension_algorithm = @RB_extension_PCA_fixspace;
593  model.RB_stop_epsilon = 1e-2;
594  model.RB_stop_timeout = 2*60*60; % 2 hours
595  model.RB_stop_Nmax = 100;
596  % decide, whether estimator or true error is error-indicator for greedy
597  % params.RB_error_indicator = 'error'; % true error, i.e. strong POD-greedy
598  model.RB_error_indicator = 'estimator'; % Delta from rb_simulation
599 
600  % choose a 10x10 uniform grid
601  model.RB_numintervals = 9 * ones(size(model.mu_names));
602 
603  detailed_data = gen_detailed_data(model,model_data);
604 
605  save rb_tutorial_step_13.mat;
606 
607  case 14
608  disp('plot first 16 snapshots and point distribution');
609 
610  load rb_tutorial_step_13.mat;
611 
612  plot_params = [];
613  plot_params.axis_equal = 1;
614  plot_params.axis_tight = 1;
615  plot_params.no_lines = 1;
616  plot_params.show_colorbar = 0;
617 
618  plot_detailed_data(model,detailed_data);
619 
620  disp('inspect diagrams and press enter to continue')
621  pause;
622 
623  close(4); % runtime plot not too interesting
624  title('maximal POD-Greedy training estimator','Fontsize',15);
625  xlabel('basis vector number N','Fontsize',15);
626  ylabel('estimator','Fontsize',15);
627  saveas(gcf,'PODgreedy_train_estimator_advection_diffusion.png');
628  saveas(gcf,'PODgreedy_train_estimator_advection_diffusion.eps','epsc');
629  close(3); % training error decay only partially interesting
630  view(2);
631  title('POD-Greedy parameter selection','Fontsize',15);
632  xlabel('mu1','Fontsize',15);
633  ylabel('mu2','Fontsize',15);
634  set(gcf,'Color',[1,1,1]);
635 % export_fig('PODgreedy_points_advection_diffusion.png',gcf);
636  saveas(gcf,'PODgreedy_points_advection_diffusion.png');
637  saveas(gcf,'PODgreedy_points_advection_diffusion.eps','epsc');
638 
639  figure;
640  for n = 1:16
641  subplot(4,4,n);
642  plot_element_data(detailed_data.grid,detailed_data.RB(:,n), ...
643  plot_params);
644  title(['n=',num2str(n)],'Fontsize',15);
645  set(gca,'Xtick',[]);
646  set(gca,'Ytick',[]);
647  end;
648  set(gcf,'Position',[32 49 1569 906]);
649  saveas(gcf,'basis_advection_diffusion.png');
650  saveas(gcf,'basis_advection_diffusion.eps','epsc');
651 
652 % plot_params = [];
653 % plot_params.axis_equal = 1;
654 % plot_params.axis_tight = 1;
655 % plot_params.no_lines = 1;
656 
657  case 15
658  disp('reduced simulation');
659 
660  load rb_tutorial_step_13.mat;
661  reduced_data = gen_reduced_data(model,detailed_data);
662  tic
663 
664  model = set_mu(model,[1,0.1]);
665  rb_sim_data = rb_simulation(model,reduced_data);
666  disp('time for reduced simulation (in seconds): ')
667  t = toc
668  rb_sim_data = rb_reconstruction(model,detailed_data,rb_sim_data);
669 
670  plot_params = [];
671  plot_params.axis_equal = 1;
672  plot_params.axis_tight = 1;
673  plot_params.no_lines = 1;
674  plot_params.show_colorbar = 0;
675 
676  plot_sim_data(model,model_data,rb_sim_data,plot_params);
677  disp('press enter for rb simulation gui')
678  pause;
679 
680  demo_rb_gui(model,detailed_data,[],plot_params,'rb simulation');
681 
682  disp('press enter for rb error simulation gui')
683  pause;
684 
685  plot_params.show_colorbar = 1;
686  demo_rb_error_gui(model,detailed_data,[],plot_params,'rb error simulation');
687 
688  case 16
689  disp('error, estimator by parameter sweep');
690 
691  % sweep along the diagonal of the parameter space, i.e.
692  % some training points, but ┬┤mostly test-points
693 
694  load rb_tutorial_step_13.mat;
695  reduced_data = gen_reduced_data(model,detailed_data);
696 
697  ntest = 101; % higher later, take care that 10x10 training grid is contained
698 
699  % RandStream.setDefaultStream(RandStream('mt19937ar','seed', ...
700  % 12345));
701  % mu_test = rand_uniform(ntest,model.mu_ranges);
702 
703  mu_test = [1;1]*(0:(ntest-1))/(ntest-1)
704 
705  errs_max_t = zeros(ntest,1);
706  errs_final_t = zeros(ntest,1);
707  Deltas_final_t = zeros(ntest,1);
708 
709  W = model.get_inner_product_matrix(detailed_data);
710  for mi = 1:ntest
711  disp(['processing parameter ',num2str(mi),'/',num2str(ntest)]);
712  model = model.set_mu(model,mu_test(:,mi));
713  sim_data = detailed_simulation(model,model_data);
714  rb_sim_data = rb_simulation(model,reduced_data);
715  rb_sim_data = rb_reconstruction(model,detailed_data, ...
716  rb_sim_data);
717  errs = rb_sim_data.U(:,model.save_time_indices+1) -sim_data.U;
718  err_norms = sqrt(sum(errs.* (W * errs)));
719  errs_final_t(mi) = err_norms(end);
720  errm = max(err_norms);
721  errs_max_t(mi) = errm(1);
722  Deltas_final_t(mi) = rb_sim_data.Delta(end);
723  end;
724  p = plot(mu_test(2,:),[errs_max_t(:),Deltas_final_t(:)]);
725  set(p,'Linewidth',2);
726  set(p(2),'Linestyle','-.');
727  title('error and error estimator at final time t=T','Fontsize',15);
728  legend({'error','estimator'},'Fontsize',15);
729  set(gca,'Yscale','log');
730  ylabel('error/estimator','Fontsize',15);
731  xlabel('parameter scaling factor s','Fontsize',15);
732 
733  saveas(gcf,'advection-diffusion-error-parameter-sweep.png');
734  saveas(gcf,'advection-diffusion-error-parameter-sweep.eps','epsc');
735  saveas(gcf,'advection-diffusion-error-parameter-sweep.fig');
736 
737  case 17
738  disp('Effectivity computation for 200 random points');
739 
740  % plot: test effectivity over diffusivity parameter
741 
742  load rb_tutorial_step_13.mat;
743  reduced_data = gen_reduced_data(model,detailed_data);
744 
745  if matlab_version> 7.9
746  s = RandStream('mt19937ar','Seed',12345)
747  RandStream.setGlobalStream(s);
748  else
749  RandStream.setDefaultStream(RandStream('mt19937ar','seed', ...
750  12345));
751  end;
752 
753  ntest = 200; % higher later
754  mu_test = rand_uniform(ntest,model.mu_ranges);
755  errs_max_t = zeros(ntest,1);
756  errs_final_t = zeros(ntest,1);
757  Deltas_final_t = zeros(ntest,1);
758 
759  for mi = 1:ntest
760  disp(['processing parameter ',num2str(mi),'/',num2str(ntest)]);
761  model = model.set_mu(model,mu_test(:,mi));
762  sim_data = detailed_simulation(model,model_data);
763  rb_sim_data = rb_simulation(model,reduced_data);
764  rb_sim_data = rb_reconstruction(model,detailed_data, ...
765  rb_sim_data);
766  errs = rb_sim_data.U(:,model.save_time_indices+1) -sim_data.U;
767  err_norms = sqrt(sum(errs.* (model.get_inner_product_matrix(detailed_data) * errs)));
768  errs_final_t(mi) = err_norms(end);
769  errm = max(err_norms);
770  errs_max_t(mi) = errm(1);
771  Deltas_final_t(mi) = rb_sim_data.Delta(end);
772  end;
773 % keyboard;
774 % [mu_sort,ind ] = sort(mu_test(2,:));
775 % plot(mu_test(2,ind),[Deltas_final_t(ind)./errs_max_t(ind)]);
776 
777  [mu_sort,ind ] = sort(mu_test(1,:));
778  plot(mu_test(1,ind),[Deltas_final_t(ind)./errs_max_t(ind)],'b*');
779  ylim = get(gca,'Ylim');
780  ylim(1) = 0;
781  set(gca,'Ylim',ylim);
782 
783  % plot: maximum test error over Basis size
784  title('effectivities at final time t=T over test set','Fontsize',15);
785  ylabel('effectivity','Fontsize',15);
786  xlabel('velocity weight mu1','Fontsize',15);
787 
788  saveas(gcf,'advection-diffusion-effectivity.png');
789  saveas(gcf,'advection-diffusion-effectivity.eps','epsc');
790  saveas(gcf,'advection-diffusion-effectivity.fig');
791 
792  otherwise
793  disp('step number unknown');
794 end;
795 
796 
797 function mu_list = lin_equidist_samples(model,numintervals)
798 % generate cell array of parameter vectors
799 % numintervals a vector of number of equal sized intervals per
800 % parameter direction
801 p = [];
802 p.range = model.mu_ranges;
803 p.numintervals = numintervals;
804 c = cubegrid(p);
805 mu_list = mat2cell(c.vertex,ones(size(c.vertex,1),1),...
806  size(c.vertex,2));
807 
808 function [U,test_mus] = precompute_test_snapshots(params,model,model_data,ntest)
809 % function precomputing test snapshots for step 6
810 test_mu1s = rand_uniform(ntest,{params.mu_range});
811 test_mus = [test_mu1s(:), ones(ntest,params.B1*params.B2-1)*1];
812 U = zeros(model_data.df_info.ndofs,0);
813 for i = 1:ntest
814  fprintf('.');
815  model = set_mu(model,test_mus(i,:));
816  sim_data = detailed_simulation(model, model_data);
817  U = [U,sim_data.uh.dofs(:)];
818 end;
819 fprintf('\n');
820 
821 function RB = lagrangian_orth_basisgen(model, ...
822  detailed_data)
823 Utot = [];
824 for m = 1:length(model.RB_mu_list)
825  model = set_mu(model,model.RB_mu_list{m});
826  sim_data = detailed_simulation(model,detailed_data);
827  if isempty(Utot)
828  Utot = model.get_dofs_from_sim_data(sim_data);
829  else
830  U = model.get_dofs_from_sim_data(sim_data);
831  Utot = [Utot U];
832  end;
833 end;
834 % Gram-Schmidt orthonormalization
835 %KN = Utot'*model.get_inner_product_matrix(detailed_data)*Utot;
836 %KN = 0.5*(KN + KN');
837 % This does no longer work with recent matlab:
838 %Ltransp = chol(KN);
839 %RB = Utot / Ltransp;
840 RB = orthonormalize_qr(Utot,model.get_inner_product_matrix(detailed_data),1e-12);
function model = thermalblock_model_struct(params)
Thermal Block example similar as described in the book of A.T. patera and G. Rozza (just one paramete...
function p = plot_sim_data(model, model_data, sim_data, plot_params)
function performing the plot of the simulation results as specified in model.
Definition: plot_sim_data.m:17
Customizable implementation of an abstract greedy algorithm.
Definition: DuneRBLeafNode.m:1
function res = thermalblock(step, params)
thermalblock example
Definition: thermalblock.m:17