rbmatlab  1.16.09
rb_tutorial_standalone.m
1 function res = rb_tutorial_standalone(step)
2 %function res = rb_tutorial_standalone(step)
3 %
4 % step : index of experiment to be performed
5 % 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 100
6 %
7 % This program allows to reproduce the thermalblock experiments of
8 % the RB-Tutorial use in the book chapter.
9 %
10 % B. Haasdonk: Reduced Basis Methods for Parametrized PDEs --
11 % A Tutorial Introduction for Stationary and Instationary Problems.
12 % Chapter in P. Benner, A. Cohen, M. Ohlberger and K. Willcox (eds.):
13 % "Model Reduction and Approximation: Theory and Algorithms", SIAM, Philadelphia, 2016
14 %
15 % A Preprint version of this chapter is available at:
16 %
17 % http://www.simtech.uni-stuttgart.de/publikationen/prints.php?ID=938
18 %
19 % This is a modified version of the original rb_tutorial.m in
20 % RBmatlab, which is available for download at www.morepas.org.
21 % The modifications are such that this experiments can be run
22 % standalone, i.e. without the RBmatlab classes. But only the stationary
23 % experiments have been transferred, i.e. the instationary experiments can
24 % be found in rb_tutorial.m
25 % For making this file independent of RBmatlab,
26 % the complete FEM-discretization has been precomputed in
27 % a file data_rb_tutorial_standalone.mat, which is used here.
28 % (This data can be (re-)generated in step 100 if RBmatlab is installed.)
29
30 % B. Haasdonk 10.4.2015
31
32 if nargin < 1
33  help rb_tutorial_standalone
34  disp('Now running step 1:')
35  step = 1
36 end;
37
38 switch step
39
40  case 1
41  disp('thermal block for different configurations and corresponding plots');
42  params = [];
43  params.B1 = 2;
44  params.B2 = 2;
45  params.mu_range = [0.1;10];
46  params.numintervals_per_block = 5;
47  model = standalone_thermalblock_model(params);
48  model_data = gen_model_data(model);
49  plot_params = [];
50  plot_params.axis_equal = 1;
51  plot_params.axis_tight = 1;
52
53  % plot 1: four identical values
54  figure;
55  model = set_mu(model,[1,1,1,1]);
56  sim_data = detailed_simulation(model,model_data);
57  plot_sim_data(model,model_data,sim_data,plot_params);
58
59  % plot 2: lower and upper values identical
60  figure;
61  model = set_mu(model,[0.1,0.1,10,10]);
62  sim_data = detailed_simulation(model,model_data);
63  plot_sim_data(model,model_data,sim_data,plot_params);
64
65  % plot 3: four different values
66  figure;
67  model = set_mu(model,[0.3,0.2,0.2,0.7]);
68  sim_data = detailed_simulation(model,model_data);
69  plot_sim_data(model,model_data,sim_data,plot_params);
70
71  % plot 4: larger division
72  % mu = rand(params.B1*params.B2,1)+0.1;
73  figure;
74  mu = [ 0.8094 0.8547 0.3760 0.7797 0.7551 0.2626, ...
75  0.2190 0.5984 1.0597 0.4404 0.6853 0.3238, ...
76  0.8513 0.3551 0.6060 0.7991 0.9909 1.0593 ...
77  0.6472 0.2386 0.2493 0.3575 0.9407 0.3543 ...
78  0.9143 0.3435 1.0293 0.4500 0.2966 0.3511 ...
79  0.7160 0.5733 0.4517 0.9308 0.6853 0.6497]';
80  params = [];
81  params.B1 = 6;
82  params.B2 = 6;
83  params.mu_range = [0.1;10];
84  model = standalone_thermalblock_model(params);
85  model_data = gen_model_data(model);
86  model = set_mu(model,mu);
87  sim_data = detailed_simulation(model,model_data);
88  plot_sim_data(model,model_data,sim_data,plot_params);
89
90  case 2
91  disp('equidistant samples for basis generation, offline versus online time');
92
93  params = [];
94  params.B1 = 2;
95  params.B2 = 2;
96  params.numintervals_per_block = 20;
97  params.mu_range = [0.1;10];
98  disp('initializing model')
99  model = standalone_thermalblock_model(params);
100
101  disp('-----------------------------------------------------')
102  disp('OFFLINE PHASE:');
103  disp('generating model data: grid, fem inner product matrix, etc.');
104  tic
105  model_data = gen_model_data(model);
106  model.RB_generation_mode = 'lagrangian';
107  model.RB_mu_list = lin_equidist_samples(model,[4,2,2,2]);
108  disp('generating detailed data: A_q,f_q,l_q + error estimator components, reduced basis.');
109  detailed_data = gen_detailed_data(model,model_data);
110  disp('generating reduced data: A_Nq,f_Nq,l_Nq + error estimator components');
111  reduced_data = gen_reduced_data(model, detailed_data);
112  t_offline = toc;
113  disp('Offline time: ')
114  t_offline
115 % save('thermalblock_2x2_data');
116
117  disp('-----------------------------------------------------')
118  disp('ONLINE PHASE:');
119  % online phase for single simulation:
120  disp('reduced simulation: assembling system and solve');
121  model = set_mu(model,[1,1,1,1]);
122  tic
123  rb_sim_data = rb_simulation(model,reduced_data);
124  t_online = toc;
125
126  disp('Online time: ')
127  t_online
128
129  case 3
130  disp('1-parameter example: error, estimator, effectivity bound plot.');
131
132  % generate reduced model
133  disp('generating reduced model')
134  params = [];
135  params.B1 = 2;
136  params.B2 = 2;
137  params.numintervals_per_block = 20;
138  params.mu_range = [0.1;10];
139  model = standalone_thermalblock_model(params);
140  model_data = gen_model_data(model);
141  model.RB_generation_mode = 'lagrangian';
142  c = 0.1;
143  mu_list = [[0.1,c,c,c];...
144  [0.5,c,c,c];...
145  [0.9,c,c,c];...
146  [1.3,c,c,c];...
147  [1.7,c,c,c]];
148  model.RB_mu_list = mu_list;
149  detailed_data = gen_detailed_data(model,model_data);
150  reduced_data = gen_reduced_data(model,detailed_data);
151
152  disp('computing parameter sweep')
153  mu1s = linspace(0.1,4,1000);
154  % rapidly computable error landscape:
155  Deltas = zeros(length(mu1s),1);
156  for n = 1:length(mu1s)
157  fprintf('.');
158  mu = [mu1s(n),c,c,c];
159  model = set_mu(model,mu);
160  sim_data = rb_simulation(model,reduced_data);
161  Deltas(n) = sim_data.Delta;
162  end;
163  fprintf('\n');
164  plot(mu1s,Deltas,'Linewidth',2);
165  xlabel('mu1','Fontsize',15);
166  ylabel('error/estimator','Fontsize',15);
167  title('error estimator by parameter sweep','Fontsize',15);
168  set(gca,'Yscale','log');
169  legend('estimator \Delta_u(\mu)');
170
171  % comparison with expensive true error
172
173  disp('computing true error for some parameters')
174  mu1s = linspace(0.1,4,40);
175  errs = zeros(length(mu1s),1);
176  for n = 1:length(mu1s)
177  fprintf('.');
178  mu = [mu1s(n),c,c,c];
179  model = set_mu(model,mu);
180  sim_data = detailed_simulation(model,model_data);
181  rb_sim_data = rb_simulation(model,reduced_data);
182  % linear combination of reduced vector with basis:
183  rb_sim_data= rb_reconstruction(model,detailed_data,rb_sim_data);
184  err = sim_data.uh - rb_sim_data.uh;
185  errs(n) = sqrt(err'*detailed_data.K*err);
186  end;
187  fprintf('\n');
188  hold on;
189  plot(mu1s,errs,'r*','Markersize',10);
190  legend({'estimator \Delta_u(\mu)','true error |u-u_N|'},'Fontsize',15);
191
192  %set(gcf,'Position', [379 492 434 328]);
193  %saveas(gcf,'error_parameter_sweep.png');
194  %saveas(gcf,'error_parameter_sweep.eps','epsc');
195
196  save('rb_tutorial_standalone_step3');
197
198  case 4
199  disp('Effectivity plot');
201
202  disp('computing parameter sweep')
203  mu1s = linspace(0.1,4,40);
204  % rapidly computable error landscape:
205  Deltas = zeros(length(mu1s),1);
206  for n = 1:length(mu1s)
207  fprintf('.');
208  mu = [mu1s(n),0.1,0.1,0.1];
209  model = set_mu(model,mu);
210  sim_data = rb_simulation(model,reduced_data);
211  Deltas(n) = sim_data.Delta;
212  end;
213  fprintf('\n');
214
215  etas = Deltas.*(errs.^(-1));
216  i = find(errs<1e-7);
217  etas(i)=NaN;
218  plot(mu1s,etas,'b*','Markersize',10);
219  xlabel('\mu1','Fontsize',15);
220  ylabel('effectivity','Fontsize',15);
221  title('effectivity over parameter','Fontsize',15);
222 % set(gca,'Yscale','log');
223
224  hold on;
225  gammas = mu1s;
226  alpha = 0.1;
227  effectivities = gammas/alpha;
228  plot(mu1s,effectivities,'r-','Linewidth',2);
229
230  legend({'effectivity \Delta_u/|u-u_N|',...
231  'upper bound \gamma(\mu)/\alpha(\mu)'},...
232  'Fontsize',15,...
233  'Location','NorthWest');
234
235 % set(gcf,'Position', [379 492 434 328]);
236 % saveas(gcf,'effectivity_parameter_sweep.png');
237 % saveas(gcf,'effectivity_parameter_sweep.eps','epsc');
238
239  case 5
240  disp('Convergence plot for test error on equidistant samples');
241
242  params = [];
243  params.B1 = 3;
244  params.B2 = 3;
245  params.numintervals_per_block = 20;
246  params.mu_range = [0.5;2];
247  model = standalone_thermalblock_model(params);
248  model_data = gen_model_data(model);
249 % model.RB_generation_mode = 'greedy_rand_uniform_with_orthogonalization';
250  model.RB_generation_mode = 'lagrangian_with_orthonormalization';
251  Ns = 2:8; % linear dependent for higher values...
252
253  maxDeltas = zeros(length(Ns),1);
254  maxerrs = zeros(length(Ns),1);
255
256  disp('Precomputing test snapshots');
257  ntest = 100;
258  [U,test_mus] = precompute_test_snapshots(params,model,model_data,ntest);
259
260  disp('Generating basis and determining max error, estimators');
261  tic;
262  for Ni = 1:length(Ns)
263  N = Ns(Ni)
264  train_mu1s = linspace(params.mu_range(1),params.mu_range(2),N);
265  train_mus = [train_mu1s(:), ones(N,params.B1*params.B2-1)*1];
266  model.RB_mu_list = train_mus;
267  %mat2cell(train_mus,ones(N,1),params.B1*params.B2);
268
269  detailed_data = gen_detailed_data(model,model_data);
270  reduced_data = gen_reduced_data(model,detailed_data);
271  for i = 1:ntest
272  model = set_mu(model,test_mus(i,:));
273  rb_sim_data = rb_simulation(model, reduced_data);
274  rb_sim_data = rb_reconstruction(model,detailed_data,rb_sim_data);
275  err = rb_sim_data.uh - U(:,i);
276  errnorm = sqrt(err' * detailed_data.K * err);
277  if maxerrs(Ni)<errnorm
278  maxerrs(Ni)=errnorm;
279  end;
280  if maxDeltas(Ni)<rb_sim_data.Delta
281  maxDeltas(Ni)=rb_sim_data.Delta;
282  end;
283  end; % i
284  end; % Ni
285
286  time_step_5 = toc
287  save('rb_tutorial_standalone_step5');
288  % case 6.1 % plot of error convergenc
290  plot(Ns,maxDeltas,'b-','Linewidth',2);
291  hold on
292  plot(Ns,maxerrs,'r:','Linewidth',2);
293  set(gca,'Yscale','log');
294  legend({'estimator \Delta_u(\mu)','error |u(\mu)-u_N(\mu)|'},'Fontsize',15);
295  title('test error decay for growing sample size','Fontsize',15);
296  xlabel('sample/basis Size N','Fontsize',15);
297  ylabel('error/estimator','Fontsize',15);
298
299  % set(gcf,'Position', [379 492 434 328]);
300  % saveas(gcf,'lagrangian_error_convergence.png');
301  % saveas(gcf,'lagrangian_error_convergence.eps','epsc');
302
303  case 6
304  disp('plot basis');
305
307  plot_params = [];
308  plot_params.axis_equal = 1;
309  plot_params.axis_tight = 1;
310  plot_params.plot = @plot_vertex_data;
311  plot_params.title='Orthonormal Reduced Basis \Phi_N';
312
313  % or in single plot
314  plot_params.show_colorbar = 0;
315  plot_params.no_lines = 1;
316  figure;
317  for n = 1:8
318  subplot(2,4,n);
319  phi_ext = zeros(length(detailed_data.grid.X),1);
320  % expand vector to full size including dirichlet values
321  phi_ext(detailed_data.grid.non_dirichlet_indices) = detailed_data.RB(:,n);
322  plot_vertex_data(detailed_data.grid,phi_ext, ...
323  plot_params);
324  title(['n=',num2str(n)],'Fontsize',10);
325  set(gca,'Xtick',[]);
326  set(gca,'Ytick',[]);
327  end;
328  %set(gcf,'Position',[ 32 254 1554 701 ]);
329  %saveasshown(gcf,'orthonormal_basis.png');
330  %saveasshown(gcf,'orthonormal_basis.eps','epsc');
331
332 case 7
333  disp('Greedy error convergence');
334
335  params = [];
336  params.B1 = 2;
337  params.B2 = 2;
338  params.numintervals_per_block = 20;
339  params.mu_range = [0.5;2];
340  model = standalone_thermalblock_model(params);
341  model_data = gen_model_data(model);
342  model.RB_generation_mode = 'greedy_rand_uniform_with_orthogonalization';
343  model.RB_stop_epsilon = 1e-6;
344  model.RB_stop_Nmax = 50;
345  model.RB_train_size = 1000;
346
347  disp('Generating basis and determining max error, estimators');
348  tic;
349  detailed_data = gen_detailed_data(model,model_data);
350  reduced_data = gen_reduced_data(model,detailed_data);
351
352  disp('Precomputing test snapshots');
353  ntest = 100;
354 % [U,test_mus] = filecache_function(@precompute_test_snapshots,...
355 % params,model,model_data,ntest);
356  [U,test_mus] = precompute_test_snapshots(params,model,model_data,ntest);
357
358  Ns = 1:size(detailed_data.RB,2);
359  maxDeltas = zeros(length(Ns),1);
360  maxerrs = zeros(length(Ns),1);
361  for Ni = 1:length(Ns)
362  N = Ns(Ni)
363  model.N = N;
364  reduced_data_subset = reduced_data_subset(model,reduced_data);
365  for i = 1:ntest
366  model = set_mu(model,test_mus(i,:));
367  rb_sim_data = rb_simulation(model, reduced_data_subset);
368  rb_sim_data = rb_reconstruction(model,detailed_data,rb_sim_data);
369  err = rb_sim_data.uh - U(:,i);
370  errnorm = sqrt(err'*detailed_data.K*err);
371  if maxerrs(Ni)<errnorm
372  maxerrs(Ni)=errnorm;
373  end;
374  if maxDeltas(Ni)<rb_sim_data.Delta
375  maxDeltas(Ni)=rb_sim_data.Delta;
376  end;
377  end; % i
378  end; % Ni
379  time_step_6 = toc
380  save('rb_tutorial_standalone_step7');
381
382  case 8
383  disp('plot of greedy errors');
385
386  plot(Ns,maxDeltas,'b-','Linewidth',2);
387  hold on
388  plot(Ns,maxerrs,'r:','Linewidth',2);
389  set(gca,'Yscale','log');
390  % if wanted, training error can be plotted
391  % plot(Ns(1:end-2),detailed_data.RB_info.max_err_est_sequence(2:end),...
392 % 'g-.','Linewidth',2);
393  set(gca,'Yscale','log');
394  legend({'test estimator \Delta_u(\mu)','test error |u(\mu)-u_N(\mu)|'},'Fontsize',15);
395  title('Error decay for Greedy basis','Fontsize',15);
396  xlabel('sample/basis Size N','Fontsize',15);
397  ylabel('error/estimator','Fontsize',15);
398
399 % set(gcf,'Position', [379 492 434 328]);
400 % saveas(gcf,'greedy_test_error.png');
401 % saveas(gcf,'greedy_test_error.eps','epsc');
402
403  % possible further cases:
404  % plot of sample points: not spectacular...
406  % mus = detailed_data.RB_info.max_mu_sequence;
407  % mumat = cell2mat(mus)
408  % plot(mumat(1,:),mumat(2,:),'*');
409
410  case 9
411  disp('experiment with basis from step 3');
412
414
415  mus= {[1,0.2,0.2,0.2],[1,0.1,1,0.1]};
416  for i = 1:length(mus)
417  disp('-------------------------------------------')
418  disp(['Experiment for mu = (',num2str(mus{i}),')']);
419  model = set_mu(model,mus{i});
420  sim_data = detailed_simulation(model,model_data);
421  rb_sim_data = rb_simulation(model,reduced_data);
422  rb_sim_data = rb_reconstruction(model,detailed_data,rb_sim_data);
423  disp(['output s(mu) = ', num2str(sim_data.s)]);
424  disp(['output s_N(mu) = ', num2str(rb_sim_data.s)]);
425  disp(['output error s(mu)-s_N(mu) = ', ...
426  num2str(sim_data.s-rb_sim_data.s)]);
427  disp(['output error bound Delta_s(\mu) = ', num2str(rb_sim_data.Delta_s)]);
428  err = sim_data.uh - rb_sim_data.uh;
429  errnorm = sqrt(err'*detailed_data.K*err);
430  disp(['field error |u-u_N| = ',num2str(errnorm)]);
431  disp(['field error bound Delta_u(\mu) = ', num2str(rb_sim_data.Delta)]);
432  end;
433
434  case 10
435  disp('convergence curves, training error for different block numbers');
436
437  B1S = [2,3,4];
438  B2S = [2,3,4];
439
440  max_err_est_sequences = {};
441
442  for bi = 1:length(B1S)
443  params = [];
444  params.B1 = B1S(bi);
445  params.B2 = B2S(bi);
446  disp('--------------------------------------------------');
447  disp(['generating basis for (B1,B2) = (',...
448  num2str(params.B1),',',num2str(params.B2),')']);
449  params.numintervals_per_block = 20;
450  params.mu_range = [0.5;2];
451  model = standalone_thermalblock_model(params);
452  model_data = gen_model_data(model);
453  model.RB_stop_epsilon = 1e-6;
454  model.RB_stop_Nmax = 50;
455  model.RB_train_size = 1000;
456  detailed_data = gen_detailed_data(model,model_data);
457  max_err_est_sequences = [max_err_est_sequences,...
458  {detailed_data.RB_info.max_err_est_sequence}]
459  end;
460
461 % save('rb_tutorial_standalone_step10');
462 % case 10.1 % plot training error curves
464
465  figure;
466  hold on;
467  legstr = {};
468  styles = {'b-','r-','g-','k-'}
469  for i=1:length(max_err_est_sequences)
470  seq = max_err_est_sequences{i};
471  plot(0:length(seq)-1,seq,styles{i});
472  legstr = [legstr,{['(B1,B2)= (',...
473  num2str(B1S(i)),',',num2str(B2S(i)),')']}]
474  end;
475  legend(legstr,'Fontsize',15);
476  title('greedy training error decay','Fontsize',15);
477  xlabel('basis size N','Fontsize',15);
478  ylabel('training error','Fontsize',15);
479  set(gca,'Yscale','log');
480
481  case 100
482  disp('generate data file data_rb_tutorial_standalone.mat (with RBmatlab!)');
483
484  if ~exist('rbmatlabhome','file')
485  error('This step can only be executed with RBmatlab installed')
486  error(['This step is not required, if the data_rb_tutorial_standalone.mat file is', ...
487  'provided.']);
488  end;
489
490  save_varlist = {};
491  params.mu_range = [0.1;10];
492  B1_list = [2,6,2,3,4];
493  B2_list = [2,6,2,3,4];
494  numintervals_per_block_list = [5,5,20,20,20];
495  for i = 1:length(B1_list);
496  params.B1 = B1_list(i);
497  params.B2 = B2_list(i);
498  params.numintervals_per_block = numintervals_per_block_list(i);
499  model = thermalblock_model_struct(params);
500  % model_data = lin_stat_gen_model_data(model);
501  model_data = [];
502  md.grid = construct_grid(model);
503  md.df_info = feminfo(model,md.grid);
504  model.decomp_mode = 1;
505  [A_comp,f_comp] = model.operators(model,md);
506  % only Neumann in our case:
507  f = f_comp{3};
508  % only Diffusion components:
509  A_comp = A_comp(1:(model.B1*model.B2));
510  l_comp = model.operators_output(model,md);
511  % cut down matrices to inner dofs:
512  dirichlet_indices = md.df_info.dirichlet_gids;
513  dirichlet_flag = zeros(md.grid.nvertices,1);
514  dirichlet_flag(dirichlet_indices)=1;
515  non_dirichlet_indices = find(dirichlet_flag==0);
516  for i = 1:length(A_comp);
517  A_comp_i = A_comp{i};
518  A_comp_i = A_comp_i(non_dirichlet_indices,non_dirichlet_indices);
519  A_comp{i} = A_comp_i;
520  end;
521  f = f(non_dirichlet_indices);
522  l = l_comp{1};
523  l = l(non_dirichlet_indices);
524  % export only structures, not classes:
525  model_data.grid.X = md.grid.X;
526  model_data.grid.Y = md.grid.Y;
527  model_data.grid.VI = md.grid.VI;
528  model_data.grid.dirichlet_indices = md.df_info.dirichlet_gids;
529  model_data.grid.non_dirichlet_indices = non_dirichlet_indices;
530  model_data.A_comp = A_comp;
531  model_data.f = f;
532  model_data.l = l;
533  model_data.K = md.df_info.regularized_h10_inner_product_matrix(...
534  non_dirichlet_indices,non_dirichlet_indices);
535
536  model_data_varname = ['model_data_',num2str(params.B1),'_',num2str(params.B2),'_',...
537  num2str(params.numintervals_per_block)];
538
539  eval([model_data_varname, ' = model_data;']);
540  save_varlist = [save_varlist, {model_data_varname}];
541  end;
542  save('data_rb_tutorial_standalone.mat',save_varlist{:});
543
544  otherwise
545  error('step unknown!');
546 end;
547
548 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
549 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
550 %%%%%%%%%%%%%%%%%%%%%%%%% auxiliary functions %%%%%%%%%%%%%%%%%%%%%%%%%%%%
551 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
552 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
553
554 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
555 %%%%%%%%%%%%%%%%%%%%%%%%% model generation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
556 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
557
558 function model = standalone_thermalblock_model(params)
559 % Thermal Block example
560 %
561 % i.e.
562 % - div ( a(x) grad u(x)) = q(x) on Omega
563 % u(x)) = g_D(x) on Gamma_D
564 % a(x) (grad u(x)) n) = g_N(x) on Gamma_N
565 % s = l(u) linear output functional
566 %
567 % where Omega = [0,1]^2 is divided into B1 * B2 blocks
568 % QA := B1*B2. The heat conductivities are given by mu_i:
569 %
570 % -------------------------
571 % | ... | ..... | mu_QA |
572 % -------------------------
573 % | .... | ... | ..... |
574 % -------------------------
575 % | mu_1 | ... | mu_B1 |
576 % -------------------------
577 %
578 % Gamma_D = upper edge
579 % Gamma_N = boundary(Omega)\Gamma_D
580 %
581 % a(x) = mu_i(x) if x\in B_i
582 % q(x) = 0
583 % g_D = 0 on top
584 % g_N = 1 on lower edge, 0 otherwise
585 % l(u) = average over lower edge
586 %
587 % The discretization information is fully cut out here and loaded
588 % from disc in gen_model_data
589
590 model = [];
591 if nargin == 0
592  params.B1=3;
593  params.B2=3;
594 end
595 model.B1 = params.B1;
596 model.B2 = params.B2;
597 model.number_of_blocks = params.B1*params.B2;
598
599 if ~isfield(params,'numintervals_per_block')
600  params.numintervals_per_block = 5;
601 end;
602 model.numintervals_per_block = params.numintervals_per_block;
603
604 mu_names = {};
605 mu_ranges = {};
606 if ~isfield(params,'mu_range')
607  params.mu_range = [0.1,10];
608 end;
609 mu_range = params.mu_range;
610 for p = 1:model.number_of_blocks
611  mu_names = [mu_names,{['mu',num2str(p)]}];
612  mu_ranges = [mu_ranges,{mu_range}];
613 end;
614 model.mu_names = mu_names;
615 model.mu_ranges = mu_ranges;
616
617 %default values 1 everywhere
618 model.mus = ones(model.number_of_blocks,1);
619
620 % basis generation settings
621 %model.RB_train_rand_seed = 100;
622 model.RB_train_size = 1000;
623 model.RB_stop_epsilon = 1e-3;
624 model.RB_stop_Nmax = 100;
625 model.RB_generation_mode = 'greedy_rand_uniform_with_orthogonalization';
626
627 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
628 %%%%%%%%%%%%%%%%%%%%%%%%% high level functions %%%%%%%%%%%%%%%%%%%%%%%%%%%
629 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
630
631
632 function model = set_mu(model,mu)
633 if length(mu)~=model.number_of_blocks
634  error('length of mu does not fit to number of blocks!');
635 end;
636 model.mus = [mu(:)];
637
638
639 function model_data = gen_model_data(model)
640 % simply load precomputed matrices and grid information
641 model_data_varname = ['model_data_',num2str(model.B1),'_',num2str(model.B2),'_',...
642  num2str(model.numintervals_per_block)];
644 eval(['model_data =', model_data_varname,';']);
645
646
647 function detailed_data = gen_detailed_data(model,model_data)
648 detailed_data = model_data;
649 detailed_data = rb_basis_generation(model,detailed_data);
650
651
652 function detailed_data = rb_basis_generation(model,detailed_data)
653 switch model.RB_generation_mode
654
655  case 'lagrangian'
656  Utot = [];
657  for m = 1:size(model.RB_mu_list,1)
658  model = set_mu(model,model.RB_mu_list(m,:));
659  sim_data = detailed_simulation(model,detailed_data);
660  if isempty(Utot)
661  Utot = sim_data.uh;
662  else
663  Utot = [Utot sim_data.uh];
664  end;
665  end;
666  detailed_data.RB = Utot;
667
668  case 'lagrangian_with_orthonormalization'
669
670  Utot = [];
671  for m = 1:size(model.RB_mu_list,1)
672  model = set_mu(model,model.RB_mu_list(m,:));
673  sim_data = detailed_simulation(model,detailed_data);
674  if isempty(Utot)
675  Utot = sim_data.uh;
676  else
677  Utot = [Utot, sim_data.uh];
678  end;
679  end;
680  % Gram-Schmidt orthonormalization
681  KN = Utot'* detailed_data.K *Utot;
682  KN = 0.5*(KN + KN');
683
684  % in case of singular KN produce less RB vectors:
685  % Ltransp = chol(KN);
686  % RB = Utot / Ltransp;
687
688  %[Ltransp,p] = chol(KN);
689  %if p==0
690  % r = size(Utot,2);
691  %else
692  % r = p-1;
693  %end;
694  %RB = Utot(:,1:r) / Ltransp;
695  %detailed_data.RB = RB;
696
697  RB = orthonormalize_qr(Utot,detailed_data.K,1e-12);
698  detailed_data.RB = RB;
699
700  case 'greedy_rand_uniform_with_orthogonalization'
701
702  sim_data = detailed_simulation(model, detailed_data);
703  n = sqrt(sim_data.uh'*detailed_data.K*sim_data.uh);
704  detailed_data.RB = sim_data.uh / n;
705  reduced_data = gen_reduced_data(model, detailed_data);
706
707  mus = rand_uniform(model.RB_train_size,model.mu_ranges);
708  max_err_est = 1e10;
709  max_err_est_sequence = [];
710  max_mu_sequence = {};
711
712  while( max_err_est> model.RB_stop_epsilon) && ...
713  (size(detailed_data.RB,2)< model.RB_stop_Nmax)
714  max_err_est = 0;
715  mu_max = [];
716  for i = 1:size(mus,2);
717  model = set_mu(model, mus(:,i));
718  rb_sim_data = rb_simulation(model,reduced_data);
719  % disp(['rb sim: Delta = ',num2str(rb_sim_data.Delta)]);
720  if rb_sim_data.Delta > max_err_est
721  mu_max = mus(:,i);
722  max_err_est = rb_sim_data.Delta;
723  end;
724  end;
725  disp(['max error estimator: ',num2str(max_err_est),...
726  ' for mu = ',num2str(mu_max')]);
727  max_err_est_sequence = [max_err_est_sequence,max_err_est];
728  max_mu_sequence = [max_mu_sequence,{mu_max}];
729  model = set_mu(model,mu_max);
730  sim_data = detailed_simulation(model,detailed_data);
731  pr = sim_data.uh - ...
732  detailed_data.RB * (detailed_data.RB' * detailed_data.K * ...
733  sim_data.uh);
734  n = sqrt(max(pr'*detailed_data.K*pr,0));
735  detailed_data.RB = [detailed_data.RB,pr/n];
736  reduced_data = gen_reduced_data(model,detailed_data);
737  disp(['new basis size: ',num2str(size(detailed_data.RB,2))]);
738  end;
739
740  RB_info.max_mu_sequence = max_mu_sequence;
741  RB_info.max_err_est_sequence = max_err_est_sequence;
742  detailed_data.RB_info = RB_info;
743
744  otherwise
745  error('RB_generation_mode unknown!!!');
746
747 end;
748
749
750 function sim_data = detailed_simulation(model,model_data)
751 sim_data =[];
752 A_coeff = model.mus;
753 Q_A = length(A_coeff);
754 A = model_data.A_comp{1} * A_coeff(1);
755 for q=2:Q_A
756  if A_coeff(q)~=0
757  A = A + A_coeff(q)*model_data.A_comp{q};
758  end;
759 end;
760 uh = A\model_data.f;
761 sim_data.uh = uh;
762 sim_data.s = (model_data.l(:)') * sim_data.uh;
763
764
765 function p = plot_sim_data(model,detailed_data,sim_data, ...
766  plot_params)
767 if nargin<4
768  plot_params = [];
769 end;
770
771 % expand vector to full size including dirichlet values
772 uh_ext = zeros(length(detailed_data.grid.X),1);
773 uh_ext(detailed_data.grid.non_dirichlet_indices) = sim_data.uh(:);
774
775 p = plot_vertex_data(detailed_data.grid,uh_ext,plot_params);
776 hold on;
777
778 % plot coarse mesh
779 if ~isfield(plot_params,'plot_blocks')
780  plot_params.plot_blocks = 1;
781 end;
782 if plot_params.plot_blocks
783  X = [0:1/model.B1:1];
784  Y = [0:1/model.B2:1];
785  l1 = line([X;X],...
786  [zeros(1,model.B1+1);...
787  ones(1,model.B1+1)]);
788  set(l1,'color',[0,0,0],'linestyle','-.');
789  %keyboard;
790  l2 = line([zeros(1,model.B2+1);...
791  ones(1,model.B2+1)],...
792  [Y;Y]);
793  set(l2,'color',[0,0,0],'linestyle','-.');
794  p = [p(:);l1(:);l2(:)];
795 end;
796
797
798 function p = plot_vertex_data(grid,data,plot_params);
799 if nargin<3
800  plot_params = [];
801 end;
802 if ~isfield(plot_params,'title');
803  plot_params.title = '';
804 end;
805 if ~isfield(plot_params,'axis_equal');
806  plot_params.axis_equal = 0;
807 end;
808 if ~isfield(plot_params,'no_lines');
809  plot_params.no_lines = 1;
810 end;
811 if ~isfield(plot_params,'show_colorbar');
812  plot_params.show_colorbar = 1;
813 end;
814 if ~isfield(plot_params,'colorbar_location');
815  plot_params.colorbar_location = 'EastOutside';
816 end;
817 % expand vector to full size including dirichlet values
818 XX = grid.X(grid.VI(:));
819 XX = reshape(XX,size(grid.VI)); % nelements*nneigh matrix
820 YY = grid.Y(grid.VI(:));
821 YY = reshape(YY,size(grid.VI)); % nelements*nneigh matrix
822 CC = data(grid.VI);
823 p = patch(XX',YY',0*CC', CC');
824 c = jet(256);
825 colormap(c);
826 if plot_params.axis_equal
827  axis equal;
828  axis tight;
829 end;
830 if plot_params.no_lines
831  set(p,'linestyle','none');
832 end;
833 if plot_params.show_colorbar
834  if isfield(plot_params,'clim')
835  set(gca,'Clim',plot_params.clim)
836  end;
837  colorbar(plot_params.colorbar_location);
838 end;
839 title(plot_params.title);
840
841
842 function reduced_data = gen_reduced_data(model,detailed_data)
843 reduced_data = [];
844 N = size(detailed_data.RB,2);
845 reduced_data.N = N;
846 A_comp = detailed_data.A_comp;
847 f = detailed_data.f;
848 Q_A = length(A_comp);
849 reduced_data.Q_A = Q_A;
850 reduced_data.AN_comp = zeros(N^2,Q_A);
851 for q = 1:Q_A
852  AN_comp_q = detailed_data.RB'*A_comp{q}*detailed_data.RB;
853  reduced_data.AN_comp(:,q) = AN_comp_q(:);
854 end;
855 % f,l non parametric, hence simple projection
856 Q_f = 1;
857 reduced_data.fN = detailed_data.RB' * detailed_data.f;
858 reduced_data.lN = detailed_data.RB' * detailed_data.l;
859
860 % plus error estimation quantities
861 % G = (v_r^q, v_r^q) = v_r^q' * K * v_r^q
862 % with {v_r^q}_q = (v_f^q, v_a^q'n)_{q'n}
863 % G = [Gff, Gfa; Gfa', Gaa];
864 %
865 % (v_f^q,v_f^q') = v_f^q' * K * v_f^q
866 % matrices of coefficient vectors of Riesz-representers:
867 % K * v_f^q = f^q (coefficient vector equations)
868 K = detailed_data.K;
869 K_v_f = detailed_data.f;
870 v_f = K \ K_v_f;
871 Q_r = Q_f + N * Q_A;
872 K_v_r = zeros(length(detailed_data.f),Q_r);
873 K_v_r = [detailed_data.f];
874 for q = 1:Q_A
875  K_v_r = [K_v_r, detailed_data.A_comp{q}*detailed_data.RB];
876 end;
877 v_r = K \ K_v_r;
878 G = v_r' * K_v_r;
879 reduced_data.G = G;
880
881
882 function reduced_data_ss = reduced_data_subset(model,reduced_data)
883 reduced_data_ss = reduced_data;
884
885 if reduced_data.N ~= model.N
886  % extract correct N-sized submatrices from reduced_data
887  if model.N > reduced_data.N
888  error('N too large for current size of reduced basis!');
889  end;
890  N = model.N;
891  dummy = zeros(reduced_data.N);
892  dummy(1:N,1:N)= 1;
893  ind = find(dummy);
894  reduced_data_ss.AN_comp = reduced_data.AN_comp(ind,:);
895  reduced_data_ss.fN = reduced_data.fN(1:N);
896  reduced_data_ss.lN = reduced_data.lN(1:N);
897  % inner product matrix of riesz-representers of residual components:
898  Q_f = 1;
899  Q_a = size(reduced_data_ss.AN_comp,2);
900
901  ind = reshape(1:reduced_data.N*Q_a,reduced_data.N,Q_a)';
902  ind = ind(:,1:N)';
903  ind = [1; ind(:) + 1];
904  reduced_data_ss.G = reduced_data.G(ind,ind);
905
906  reduced_data_ss.N = N;
907
908 end
909
910
911 function rb_sim_data = rb_reconstruction(model,detailed_data,rb_sim_data)
912 rb_sim_data.uh = ...
913  detailed_data.RB(:,1:length(rb_sim_data.uN)) * rb_sim_data.uN;
914
915
916 function rb_sim_data = rb_simulation(model,reduced_data)
917 rb_sim_data =[];
918 A_coeff = model.mus(:);
919 N = reduced_data.N;
920 AN = reshape(reduced_data.AN_comp * A_coeff, N , N);
921 uN = AN\reduced_data.fN;
922 rb_sim_data.uN = uN;
923 rb_sim_data.s = reduced_data.lN(:)' * rb_sim_data.uN;
924 % for elliptic compliant case, X-norm (=H10-norm) error estimator:
925 Q_r = size(reduced_data.G,1);
926 neg_auN_coeff = -uN * A_coeff';
927 f_coeff = 1;
928 res_coeff = [f_coeff; neg_auN_coeff(:)];
929 res_norm_sqr = res_coeff' * reduced_data.G * res_coeff;
930 % prevent possibly negative numerical disturbances:
931 res_norm_sqr = max(res_norm_sqr,0);
932 res_norm = sqrt(res_norm_sqr);
933 alpha_LB = min(model.mus);
934 rb_sim_data.Delta = ...
935  res_norm/alpha_LB;
936 rb_sim_data.Delta_s = ...
937  res_norm_sqr/alpha_LB;
938
939
940 function mu_list = lin_equidist_samples(model,numintervals)
941 mu_list = [];
942 for i = 1:length(model.mus);
943  mu_range = model.mu_ranges{i};
944  n = numintervals(i)+1;
945  mi = linspace(mu_range(1),mu_range(2),n);
946  if i==1
947  mmi = mi;
948  else
949  mmi = ones(size(mu_list,1),1) * mi(:)';
950  end;
951  mmi_vec = mmi(:);
952  mu_list = [repmat(mu_list,n,1), mmi_vec];
953 end;
954
955
956 function [U,test_mus] = precompute_test_snapshots(params,model,model_data,ntest)
957 % function precomputing test snapshots for step 5
958 test_mu1s = rand_uniform(ntest,{params.mu_range});
959 test_mus = [test_mu1s(:), ones(ntest,params.B1*params.B2-1)*1];
960 n = length(model_data.grid.non_dirichlet_indices);
961 U = zeros(n,0);
962 for i = 1:ntest
963  fprintf('.');
964  model = set_mu(model,test_mus(i,:));
965  sim_data = detailed_simulation(model, model_data);
966  U = [U,sim_data.uh(:)];
967 end;
968 fprintf('\n');
969
970 function Xon = orthonormalize_qr(X,A,epsilon)
971 % function Xon = orthonormalize_qr(X[,A,epsilon])
972 %
973 % orthonormalization of vectors in columns of
974 % matrix X to Xon by qr-decomposition of vectors and
975 % the inner product matrix A can be specified additionally, i.e.
976 % <x,x> = x' * A * x. Then the Cholesky factoriziation of A is
977 % performed A = R_M' * R_M
978 %
979 % Then a QR-decomposition of R_M * X is performed which yields the
980 % desired new vectors
981 if nargin < 3
982  epsilon = 1e-7; % => for nonlinear evolution this value is required.
983 end;
984 if isempty(X)
985  Xon = zeros(size(X));
986  return;
987 end;
988 if nargin<2
989  A = 1;
990 end;
991 % check on identity of vectors (can happen, that numerics does not detect
992 % this afterwards due to rounding errors!!)
993 for i=1:(size(X,2)-1)
994  for j=(i+1):size(X,2)
995  if isequal(X(:,i),X(:,j))
996  X(:,j) = 0;
997  end;
998  end;
999 end;
1000 Xon = X;
1001 % incomplete cholesky of inner-product matrix:
1002 if matlab_version > 7.9
1003  opts.shape = 'upper';
1004  opts.type = 'ict';
1005  R_M = ichol(sparse(A), opts);
1006 else
1007  R_M = cholinc(sparse(A),'inf');
1008 end;
1009 % qr decomposition of R_M * X with permutation indices E
1010 [Q, R, E]= qr(R_M * Xon,0 );
1011 % sort such that linear independent first and Q * R = R_M * Xon
1012 Xon = Xon(:,E);
1013 % search nonvanishing diagonal entries of R
1014 ind = find(abs(diag(R))>epsilon);
1015 Xon = Xon(:,ind);
1016 Rind = R(ind,ind);
1017 Xon = Xon(:,ind) / Rind;
1018 % resort vectors such that original order is recovered most:
1019 K = abs(X'* A * Xon);
1020 % search in each row maximum entry
1021 ind = [];
1022 for i = 1:size(Xon,2)
1023  [max_val, max_i] = max(K(i,:));
1024  if length(max_i)>1
1025  error('found multiple matching vectors!');
1026  end;
1027  if max_val > 0
1028  ind = [ind, max_i];
1029  end;
1030  K(:,max_i) = 0; % "delete" from further search
1031 end;
1032 all_ind = ones(1,size(Xon,2));
1033 all_ind(ind) = 0;
1034 rest_ind = find(all_ind);
1035 % resort:
1036 if ~isempty(rest_ind)
1037  Xon = [Xon(:,ind),Xon(:,rest_ind)];
1038 else
1039  Xon = Xon(:,ind);
1040 end;
1041
1042 function M = rand_uniform(N,intervals)
1043 %function M = rand_uniform(N,intervals)
1044 %
1045 % function generating uniformly distributed random data in a hypercube
1046 % N vectors are generated as columns of the matrix M, intervals is a
1047 % cell array indicating the borders of the intervals
1048 %
1049 % example: rand_uniform(100,{[0,1],[100,1000]})
1050  mu_dim = length(intervals);
1051  M = rand(mu_dim,N);
1052  for i=1:mu_dim
1053  mu_range = intervals{i};
1054  M(i,:) = M(i,:)*(mu_range(2)-mu_range(1))+mu_range(1);
1055  end;
1056
1057 function res = matlab_version
1058 %function res = matlab_version
1059 %
1060 % return numerical value of MATLAB version. Currently used to
1061 % switch between ichol and cholinc.
1062 v = ver('MATLAB');
1063 res = str2num(v.Version);
function varargout = filecache_function(funcptr, varargin)
function used for file-caching other function calls.
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
function res = thermalblock(step, params)
thermalblock example
Definition: thermalblock.m:17