rbmatlab  1.16.09
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');
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
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');
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');
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...
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
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
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;
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;
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',[]);
547  '-t',num2str(ti),'.png']);
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');
575
576  case 13
577  disp('POD-Greedy basis generation, caution takes +- 1 hour!');
578
579  params = [];
580  params.coarse_factor = 4;
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
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);
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]);
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]);
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
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
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
736
737  case 17
738  disp('Effectivity computation for 200 random points');
739
740  % plot: test effectivity over diffusivity parameter
741
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
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