1 function res = rb_tutorial_buggy(step)
3 % script
for RB Summerschool in Munich, 16.9.2013-19.9.2013
4 % script used
for generating pictures
for course
6 % B. Haasdonk, 7.9.2013
13 case 1 % thermal block
for different configurations and
18 params.mu_range = [0.1;10];
19 params.numintervals_per_block = 5;
21 model_data = gen_model_data(model);
23 plot_params.axis_equal = 1;
24 plot_params.axis_tight = 1;
26 % plot 1: four identical values
28 model = set_mu(model,[1,1,1,1]);
29 sim_data = detailed_simulation(model,model_data);
32 % plot 2: lower and upper values identical
34 model = set_mu(model,[0.1,0.1,10,10]);
35 sim_data = detailed_simulation(model,model_data);
38 % plot 3: four different values
40 model = set_mu(model,[0.3,0.2,0.2,0.7]);
41 sim_data = detailed_simulation(model,model_data);
44 % plot 4: larger division
45 % mu = rand(params.B1*params.B2,1)+0.1;
47 mu = [ 0.8094 0.8547 0.3760 0.7797 0.7551 0.2626, ...
48 0.2190 0.5984 1.0597 0.4404 0.6853 0.3238, ...
49 0.8513 0.3551 0.6060 0.7991 0.9909 1.0593 ...
50 0.6472 0.2386 0.2493 0.3575 0.9407 0.3543 ...
51 0.9143 0.3435 1.0293 0.4500 0.2966 0.3511 ...
52 0.7160 0.5733 0.4517 0.9308 0.6853 0.6497]
';
56 params.mu_range = [0.1;10];
57 model = thermalblock_model_struct(params);
58 model_data = gen_model_data(model);
59 model = set_mu(model,mu);
60 sim_data = detailed_simulation(model,model_data);
61 plot_sim_data(model,model_data,sim_data,plot_params);
63 case 2 % thermal block detailed gui
67 params.numintervals_per_block = 50;
68 params.mu_range = [0.1;10];
69 params.numintervals_per_block = 5;
70 model = thermalblock_model_struct(params);
71 model_data = gen_model_data(model);
73 plot_params.axis_equal = 1;
74 plot_params.axis_tight = 1;
75 plot_params.yscale_uicontrols = 0.7;
76 demo_detailed_gui(model,model_data,[],plot_params);
77 set(gcf,'Position
',[584 245 553 691]);
79 case 3 % equidistant samples for basis generation, offline versus online time
84 params.numintervals_per_block = 20;
85 params.mu_range = [0.1;10];
86 disp('initializing model
')
87 model = thermalblock_model_struct(params);
89 disp('-----------------------------------------------------
')
90 disp('OFFLINE PHASE:
');
91 disp('generating model data: grid, fem inner product matrix, etc.
');
93 model_data = gen_model_data(model);
94 model.RB_generation_mode = 'lagrangian
';
95 model.RB_mu_list = lin_equidist_samples(model,[4,2,2,2]);
96 disp('generating detailed data: A_q,f_q,l_q components, reduced basis.
');
97 detailed_data = gen_detailed_data(model,model_data);
98 disp('generating reduced data: A_Nq,f_Nq,l_Nq components
');
99 reduced_data = gen_reduced_data(model, detailed_data);
101 disp('Offline time:
')
103 save('thermalblock_2x2_data
');
105 disp('-----------------------------------------------------
')
106 disp('ONLINE PHASE:
');
107 % online phase for single simulation:
108 disp('reduced simulation: assembling system and solve
');
109 model = set_mu(model,[1,1,1,1]);
111 sim_data = rb_simulation(model,reduced_data);
114 disp('Online time:
')
117 disp('starting reduced basis gui:
')
119 plot_params.axis_equal = 1;
120 plot_params.axis_tight = 1;
121 plot_params.yscale_uicontrols = 0.7;
122 demo_rb_gui(model,detailed_data,reduced_data,plot_params);
124 case 4 % 1-parameter example: error, estimator, effectivity bound plot.
126 % generate reduced model
127 disp('generating reduced model
')
131 params.numintervals_per_block = 20;
132 params.mu_range = [0.1;10];
133 model = thermalblock_model_struct(params);
134 model_data = gen_model_data(model);
135 model.RB_generation_mode = 'lagrangian
';
137 mu_list = {[0.1,c,c,c],...
142 model.RB_mu_list = mu_list;
143 detailed_data = gen_detailed_data(model,model_data);
144 reduced_data = gen_reduced_data(model,detailed_data);
146 disp('computing parameter sweep
')
147 mu1s = linspace(0.1,4,1000);
148 % rapidly computable error landscape:
149 Deltas = zeros(length(mu1s),1);
150 for n = 1:length(mu1s)
152 error('fix something here in the following line!
');
154 model = set_mu(model,mu);
155 sim_data = rb_simulation(model,reduced_data);
156 Deltas(n) = sim_data.Delta;
159 plot(mu1s,Deltas,'Linewidth
',2);
160 xlabel('mu1
','Fontsize
',15);
161 ylabel('error/estimator
','Fontsize
',15);
162 title('error estimator by parameter sweep
','Fontsize
',15);
163 set(gca,'Yscale
','log
');
164 legend('estimator \Delta_u(\mu)
');
166 % comparison with expensive true error
168 disp('computing
true error
for some parameters
')
169 mu1s = linspace(0.1,4,40);
170 errs = zeros(length(mu1s),1);
171 for n = 1:length(mu1s)
173 mu = [mu1s(n),c,c,c];
174 model = set_mu(model,mu);
175 sim_data = detailed_simulation(model,model_data);
176 rb_sim_data = rb_simulation(model,reduced_data);
177 % linear combination of reduced vector with basis:
178 rb_sim_data= rb_reconstruction(model,detailed_data,rb_sim_data);
180 error('fix something here in the following line!
');
182 errs(n) = fem_h10_norm(err);
186 plot(mu1s,errs,'r*
','Markersize
',10);
187 legend({'estimator \Delta_u(\mu)
','true error |u-u_N|
'},'Fontsize
',15);
189 %set(gcf,'Position
', [379 492 434 328]);
190 %saveas(gcf,'error_parameter_sweep.png
');
191 %saveas(gcf,'error_parameter_sweep.eps
','epsc
');
193 save('rb_tutorial_step4
');
195 case 5 % Effectivity plot
196 load('rb_tutorial_step4
');
198 disp('computing parameter sweep
')
199 mu1s = linspace(0.1,4,40);
200 % rapidly computable error landscape:
201 Deltas = zeros(length(mu1s),1);
202 for n = 1:length(mu1s)
204 mu = [mu1s(n),0.1,0.1,0.1];
205 model = set_mu(model,mu);
206 sim_data = rb_simulation(model,reduced_data);
207 Deltas(n) = sim_data.Delta;
211 error('fix something here in the following line!
');
214 % remove large values, that occured by division by (almost) zero
217 plot(mu1s,etas,'b*
','Markersize
',10);
218 xlabel('\mu1
','Fontsize
',15);
219 ylabel('effectivity
','Fontsize
',15);
220 title('effectivity over parameter
','Fontsize
',15);
221 % set(gca,'Yscale
','log
');
226 effectivities = gammas/alpha;
227 plot(mu1s,effectivities,'r-
','Linewidth
',2);
229 legend({'effectivity \Delta_u/|u-u_N|
',...
230 'upper bound \gamma(\mu)/\alpha(\mu)
'},...
232 'Location
','NorthWest
');
234 % set(gcf,'Position
', [379 492 434 328]);
235 % saveas(gcf,'effectivity_parameter_sweep.png
');
236 % saveas(gcf,'effectivity_parameter_sweep.eps
','epsc
');
238 case 6 % Convergence plot for test error on equidistant samples
243 params.numintervals_per_block = 20;
244 params.mu_range = [0.5;2];
245 model = thermalblock_model_struct(params);
246 model_data = gen_model_data(model);
247 model.RB_generation_mode = 'model_RB_basisgen
';
248 model.RB_basisgen = @lagrangian_orth_basisgen;
249 Ns = 2:8; % linear dependent for higher values...
251 maxDeltas = zeros(length(Ns),1);
252 maxerrs = zeros(length(Ns),1);
254 disp('Precomputing test snapshots
');
256 [U,test_mus] = filecache_function(@precompute_test_snapshots,...
257 params,model,model_data,ntest);
259 disp('Generating basis and determining max error, estimators
');
261 for Ni = 1:length(Ns)
263 train_mu1s = linspace(params.mu_range(1),params.mu_range(2),N);
264 train_mus = [train_mu1s(:), ones(N,params.B1*params.B2-1)*1];
265 model.RB_mu_list = mat2cell(train_mus,ones(N,1),params.B1*params.B2);
266 detailed_data = gen_detailed_data(model,model_data);
267 reduced_data = gen_reduced_data(model,detailed_data);
269 model = set_mu(model,test_mus(i,:));
270 rb_sim_data = rb_simulation(model, reduced_data);
271 rb_sim_data = rb_reconstruction(model,detailed_data,rb_sim_data);
272 err = rb_sim_data.uh;
273 err.dofs = err.dofs - U(:,i);
274 errnorm = fem_h10_norm(err);
275 if maxerrs(Ni)<errnorm
278 if maxDeltas(Ni)<rb_sim_data.Delta
279 maxDeltas(Ni)=rb_sim_data.Delta;
285 % save('rb_tutorial_step6
');
286 % case 6.1 % plot of error convergenc
287 % load('rb_tutorial_step6
');
288 plot(Ns,maxDeltas,'b-
','Linewidth
',2);
290 plot(Ns,maxerrs,'r:
','Linewidth
',2);
291 set(gca,'Yscale
','log
');
292 legend({'estimator \Delta_u(\mu)
','error |u(\mu)-u_N(\mu)|
'},'Fontsize
',15);
293 title('test error decay
for growing sample size
','Fontsize
',15);
294 xlabel('sample/basis Size N
','Fontsize
',15);
295 ylabel('error/estimator
','Fontsize
',15);
297 % set(gcf,'Position
', [379 492 434 328]);
298 % saveas(gcf,'lagrangian_error_convergence.png
');
299 % saveas(gcf,'lagrangian_error_convergence.eps
','epsc
');
302 load('rb_tutorial_step6
');
304 plot_params.axis_equal = 1;
305 plot_params.axis_tight = 1;
306 plot_params.plot = @plot_vertex_data;
307 plot_params.title='Orthonormal Reduced Basis \Phi_N
';
309 plot_sequence(detailed_data.RB,detailed_data.grid,plot_params);
312 plot_params.show_colorbar = 0;
313 plot_params.no_lines = 1;
317 plot_vertex_data(detailed_data.grid,detailed_data.RB(:,n), ...
319 title(['n=
',num2str(n)],'Fontsize
',15);
323 set(gcf,'Position
',[32 49 1569 906]);
324 % saveas(gcf,'orthonormal_basis.png
');
325 % saveas(gcf,'orthonormal_basis.eps
','epsc
');
327 case 8 % Greedy error convergence
332 params.numintervals_per_block = 20;
333 params.mu_range = [0.5;2];
334 model = thermalblock_model_struct(params);
335 model_data = gen_model_data(model);
336 model.RB_stop_epsilon = 1e-6;
337 model.RB_stop_Nmax = 50;
338 model.RB_train_size = 1000;
340 disp('Generating basis and determining max error, estimators
');
342 detailed_data = gen_detailed_data(model,model_data);
343 reduced_data = gen_reduced_data(model,detailed_data);
345 disp('Precomputing test snapshots
');
347 [U,test_mus] = filecache_function(@precompute_test_snapshots,...
348 params,model,model_data,ntest);
350 Ns = 1:size(detailed_data.RB,2);
351 maxDeltas = zeros(length(Ns),1);
352 maxerrs = zeros(length(Ns),1);
353 for Ni = 1:length(Ns)
356 reduced_data_subset = model.reduced_data_subset(model,reduced_data);
358 model = set_mu(model,test_mus(i,:));
359 rb_sim_data = rb_simulation(model, reduced_data_subset);
360 rb_sim_data = rb_reconstruction(model,detailed_data,rb_sim_data);
361 err = rb_sim_data.uh;
362 err.dofs = err.dofs - U(:,i);
363 errnorm = fem_h10_norm(err);
364 if maxerrs(Ni)<errnorm
367 if maxDeltas(Ni)<rb_sim_data.Delta
368 maxDeltas(Ni)=rb_sim_data.Delta;
373 save('rb_tutorial_step8
');
375 case 8.1 % plot of greedy errors
376 load('rb_tutorial_step8
');
378 plot(Ns,maxDeltas,'b-
','Linewidth
',2);
380 plot(Ns,maxerrs,'r:
','Linewidth
',2);
381 set(gca,'Yscale
','log
');
382 % if wanted, training error can be plotted
383 % plot(Ns(1:end-2),detailed_data.RB_info.max_err_est_sequence(2:end),...
384 % 'g-.
','Linewidth
',2);
385 set(gca,'Yscale
','log
');
386 legend({'test estimator \Delta_u(\mu)
','test error |u(\mu)-u_N(\mu)|
'},'Fontsize
',15);
387 title('Error decay
for Greedy basis
','Fontsize
',15);
388 xlabel('sample/basis Size N
','Fontsize
',15);
389 ylabel('error/estimator
','Fontsize
',15);
391 % set(gcf,'Position
', [379 492 434 328]);
392 % saveas(gcf,'greedy_test_error.png
');
393 % saveas(gcf,'greedy_test_error.eps
','epsc
');
395 % possible further cases:
396 % plot of sample points: not spectacular...
397 % load('rb_tutorial_step8
');
398 % mus = detailed_data.RB_info.max_mu_sequence;
399 % mumat = cell2mat(mus)
400 % plot(mumat(1,:),mumat(2,:),'*
');
402 case 9 % experiment with basis from 4
404 case 10 % convergence curves, training error for different block numbers
407 disp('step number unknown
');
411 function mu_list = lin_equidist_samples(model,numintervals)
412 % generate cell array of parameter vectors
413 % numintervals a vector of number of equal sized intervals per
414 % parameter direction
416 p.range = model.mu_ranges;
417 p.numintervals = numintervals;
419 mu_list = mat2cell(c.vertex,ones(size(c.vertex,1),1),...
422 function [U,test_mus] = precompute_test_snapshots(params,model,model_data,ntest)
423 % function precomputing test snapshots for step 6
424 test_mu1s = rand_uniform(ntest,{params.mu_range});
425 test_mus = [test_mu1s(:), ones(ntest,params.B1*params.B2-1)*1];
426 U = zeros(model_data.df_info.ndofs,0);
429 model = set_mu(model,test_mus(i,:));
430 sim_data = detailed_simulation(model, model_data);
431 U = [U,sim_data.uh.dofs(:)];
435 function RB = lagrangian_orth_basisgen(model, ...
438 for m = 1:length(model.RB_mu_list)
439 model = set_mu(model,model.RB_mu_list{m});
440 sim_data = detailed_simulation(model,detailed_data);
442 Utot = model.get_dofs_from_sim_data(sim_data);
444 U = model.get_dofs_from_sim_data(sim_data);
448 % Gram-Schmidt orthonormalization
449 KN = Utot'*model.get_inner_product_matrix(detailed_data)*Utot;
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.
Customizable implementation of an abstract greedy algorithm.