2 % small script implementing a simple advection example
3 %
for producing matrices to be used in the RB-DS framework
4 % Discretization with FV Functions.
6 % mu_names = {
'vx_weight',
'vy_weight'};
11 % step = 1; % initialization of model and plot of model data
12 % step = 2; % detailed simulation of lin_evol
13 % step = 2.5; % conversion to DS primal model, lin_ds detailed simulation
14 % step = 3; % conversion to DS primal model, DS detailed simulation
15 % and comparison with lin_evol
16 % step = 4 % compute several trajectories.
17 % step = 4.5 % compute reduced basis based on trajecories.
18 % step = 4.75 % comparison of detailed and reduced simulation
19 % step = 5 % figures of detailed simulations based on step 4
20 % step = 6 % figures of reduced simuations based on step 4
21 % step = 7 % time measurements based on step 4
22 % step = 8; % compare original and accelerated ds_model
23 % step = 9; % generate greedy basis on subset of parameter domain
24 % step = 10; % reduced simulations and error plots
25 % step = 11; % test of matlab ode solver
for resulting sparse system
26 % step = 12; % further plots
for MCMDS paper
28 % Bernard Haasdonk 25.8.2009
38 %params.coarse_factor = 2; % fact 2 => ndofs = 16384, nt = 1024, 15s detailed
39 %params.coarse_factor = 1; % fact 1 => ndofs = 32768, nt = 2048, 60s detailed
40 params.coarse_factor = 2;
45 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
46 %step = 1; % initialization of model and plot of model data
48 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
50 %
verbose(1,
'initializaton of lin_evol_model and plot model_data');
51 disp(
'initializaton of lin_evol_model and plot model_data');
52 model = advection_fv_output_model(params);
53 model_data = gen_model_data(model);
54 params.axis_equal = 1;
55 params.axis_tight = 1;
56 plot(model_data.grid,params);
57 % inspect(model_data.grid)
58 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
59 %step = 2; % detailed simulation of lin_evol
60 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
62 %
verbose(1,
'detailed simulation of lin_evol_model and plot');
63 disp(
'detailed simulation of lin_evol_model and plot');
64 model = advection_fv_output_model(params);
65 model_data = gen_model_data(model);
67 % disp(
'model debugging turned on.')
68 model.cone_weight = 0.0;
69 model.vx_weight = 0.0;
70 model.vy_weight = 0.0;
73 model.cone_weight = 1.0;
74 % model.vx_weight = 1;
75 % model.vy_weight = 0;
76 sim_data = detailed_simulation(model,model_data);
80 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
81 %step = 2.5; % conversion to DS primal model, DS detailed simulation
82 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
84 disp('initialization of lin_ds model from lin_evol model');
85 %
verbose(1,'initialization of lin_ds model from lin_evol model and ...
88 lin_evol_model = advection_fv_output_model(params);
89 %lin_evol_model.debug = 1;
90 lin_evol_model_data = gen_model_data(lin_evol_model);
92 % set some additional fields in model which will be copied into
93 % ds model depending on whether constants are known or not:
94 estimate_constants = 0;
96 %%% if constant estimaton has to be performed:
97 lin_evol_model.estimate_lin_ds_nmu = 4;
98 lin_evol_model.estimate_lin_ds_nX = 10;
99 lin_evol_model.estimate_lin_ds_nt = 4;
102 % the following values are rough bounds, so could be
103 % non-rigorous, then larger search must be performed.
104 lin_evol_model.state_bound_constant_C1 = 1;
105 lin_evol_model.output_bound_constant_C2 = 1.2916;
108 ds_model = lin_ds_from_lin_evol_model(lin_evol_model);
109 %ds_model.u_function_ptr = @(params) 0.1; % define new
110 % ds_model.theta = 1;
112 ds_model_data = gen_model_data(ds_model);
115 % set mu in model and base_model!!!
116 ds_model = ds_model.set_mu(ds_model,mu);
118 ds_sim_data = detailed_simulation(ds_model,ds_model_data);
122 params.plot_title='detailed implicit lin_ds';
123 lin_evol_from_ds_sim_data.U = ds_sim_data.X;
124 lin_evol_from_ds_sim_data.y = ds_sim_data.Y;
125 plot_sim_data(lin_evol_model,lin_evol_model_data,lin_evol_from_ds_sim_data,params);
127 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
128 %step = 3; % conversion to DS primal model, DS detailed simulation
129 % and comparison with lin_evol
130 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
132 case 3 % initialization of DS model and lin_ds detailed simulation
133 disp('initialization of lin_ds model from lin_evol model and plot');
134 %
verbose(1,'initialization of lin_ds model from lin_evol model and ...
137 lin_evol_model = advection_fv_output_model(params);
138 %lin_evol_model.debug = 1;
139 lin_evol_model_data = gen_model_data(lin_evol_model);
141 % set some additional fields in model which will be copied into
142 % ds model depending on whether constants are known or not:
143 estimate_constants = 0;
144 if estimate_constants
145 %%% if constant estimaton has to be performed:
146 lin_evol_model.estimate_lin_ds_nmu = 4;
147 lin_evol_model.estimate_lin_ds_nX = 10;
148 lin_evol_model.estimate_lin_ds_nt = 4;
151 % the following values are rough bounds, so could be
152 % non-rigorous, then larger search must be performed.
153 lin_evol_model.state_bound_constant_C1 = 1;
154 lin_evol_model.output_bound_constant_C2 = 1.2916;
157 ds_model = lin_ds_from_lin_evol_model(lin_evol_model);
158 % ds_model.theta = 0.0;
159 %ds_model.u_function_ptr = @(params) 0.1; % define new
160 ds_model_data = gen_model_data(ds_model);
163 % set mu in model and base_model!!!
164 ds_model = ds_model.set_mu(ds_model,mu);
166 ds_sim_data = detailed_simulation(ds_model,ds_model_data);
167 lin_evol_model = lin_evol_model.set_mu(lin_evol_model,mu);
168 lin_evol_sim_data = detailed_simulation(lin_evol_model,...
169 lin_evol_model_data);
170 err = lin_evol_sim_data.U - ds_sim_data.X;
171 disp(['max error in DOFs of lin_evol and ds simulation: ',...
172 num2str(max(abs(err(:))))]);
174 params.plot_title='detailed lin_evol';
176 plot_sim_data(lin_evol_model,lin_evol_model_data,lin_evol_sim_data,params);
177 params.plot_title='detailed lin_ds';
178 lin_evol_from_ds_sim_data.U = ds_sim_data.X;
179 lin_evol_from_ds_sim_data.y = ds_sim_data.Y;
180 plot_sim_data(lin_evol_model,lin_evol_model_data,lin_evol_from_ds_sim_data,params);
181 params.plot_title='error explicit lin_ds and lin_evol';
182 error_sim_data = lin_evol_sim_data;
183 error_sim_data.U = error_sim_data.U-ds_sim_data.X;
184 plot_sim_data(lin_evol_model,lin_evol_model_data,error_sim_data,params);
186 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
187 %step = 4 % reduced Basis generation by single trajectories and
188 %comparison of detailed and reduced simulation
189 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
191 case 4 % generate some trajectories
193 ds_model = fast_model(params);
194 ds_model.enable_error_estimator = 0;
195 % lin_evol_model = advection_fv_output_model(params);
196 %lin_evol_model.debug = 1;
197 % lin_evol_model_data = gen_model_data(lin_evol_model);
199 % set some additional fields in model which will be copied into
200 % ds model depending on whether constants are known or not:
201 % estimate_constants = 0;
202 % if estimate_constants
203 % %%% if constant estimaton has to be performed:
204 % lin_evol_model.estimate_lin_ds_nmu = 4;
205 % lin_evol_model.estimate_lin_ds_nX = 10;
206 % lin_evol_model.estimate_lin_ds_nt = 4;
209 % % the following values are rough bounds, so could be
210 % % non-rigorous, then larger search must be performed.
211 % lin_evol_model.state_bound_constant_C1 = 1;
212 % lin_evol_model.output_bound_constant_C2 = 1.2916;
215 % ds_model = lin_ds_from_lin_evol_model(lin_evol_model);
216 %ds_model.u_function_ptr = @(params) 0.1; % define new
217 ds_model_data = gen_model_data(ds_model);
219 % mu = [0.5,0.5,0.5];
221 %mus = {[0,0,0],[1,1,1],[1,1,0.5],[0,0.5,1]};
222 mus = {[0,0],[1,1],[1,0.5],[0.5,1]};
223 for m = 1:length(mus)
224 disp(['computing trajectory number ',num2str(m)]);
226 ds_model = ds_model.set_mu(ds_model,mu);
229 %disp(
'set theta to 1!!!');
231 ds_sim_data = detailed_simulation(ds_model,ds_model_data);
233 disp([
'time for full simulation: ',num2str(t_det)]);
234 save(fullfile(rbmatlabresult,[
'trajectory',num2str(m)]),...
235 'ds_model',
'ds_model_data',
'ds_sim_data');
238 case 4.5 % compute reduced basis.
239 % make sure that the previous trajectories are computed!!
241 disp(
'PCA of snapshots takes a while...');
242 load(fullfile(rbmatlabresult,
'trajectory1'),...
243 'ds_model',
'ds_model_data',
'ds_sim_data');
244 % smaller epsilon does result in more than 1 vector
for mu=0,0,0
246 RB = PCA_fixspace(ds_sim_data.X,[],ds_model_data.G,[],[],epsilon);
247 %RB = orthonormalize_qr(ds_sim_data.X,ds_model_data.G);
248 % select correct orthogonal set
249 K = RB
' * ds_model_data.G * RB;
250 i = find(abs(diag(K)-1)>1e-2);
252 disp(['found
',num2str(size(RB,2)),' basis vectors after traj. 1
']);
258 disp(['reduced to
',num2str(size(RB,2)),' basis vectors.
']);
261 load(fullfile(rbmatlabresult,'trajectory2
'),...
262 'ds_model
','ds_model_data
','ds_sim_data
');
263 RB2 = PCA_fixspace(ds_sim_data.X,RB,ds_model_data.G,[],[],epsilon);
264 %RB = orthonormalize_qr(ds_sim_data.X,ds_model_data.G);
265 % select correct orthogonal set
266 K = RB2' * ds_model_data.G * RB2;
267 i = find(abs(diag(K)-1)>1e-2);
272 disp(['found ',num2str(size(RB2,2)),' basis vectors after traj. 2']);
274 RB2 = RB2(:,1:(i-1));
275 disp(['reduced to ',num2str(size(RB2,2)),' basis vectors after traj. 2']);
277 % RB = orthonormalize(RB);
280 W = ds_model_data.G * V;
281 rbfname = fullfile(rbmatlabresult,'advection_fv_output_basis');
284 ds_model.RB_basis_filename = rbfname;
285 ds_model.RB_generation_mode = 'file_load';
287 detailed_data = gen_detailed_data(ds_model,ds_model_data);
290 save(fullfile(rbmatlabresult,'advection_fv_output_detailed_data_2traj'),...
291 'model','detailed_data')
293 case 4.75 % compare reduced and detailed simulation
295 load(fullfile(rbmatlabresult,'advection_fv_output_detailed_data_2traj'),...
296 'model','detailed_data')
297 load(fullfile(rbmatlabresult,'trajectory2'),...
302 %perform single reduced simulation for given parameter
303 reduced_data = gen_reduced_data(ds_model,detailed_data);
304 ds_model.N = reduced_data.N;
306 % ds_model.nt = ds_model.nt/4;
307 % disp('set nt to quarter!!');
308 rb_sim_data = rb_simulation(ds_model,reduced_data);
310 rb_sim_data = rb_reconstruction(ds_model,detailed_data,rb_sim_data)
311 disp(['time for reduced simulation: ',num2str(t_red)]);
313 % plot of reduced trajectory
314 lin_evol_model = ds_model.base_model;
315 lin_evol_model_data = gen_model_data(lin_evol_model);
318 params.show_colorbar = 1;
319 params.axis_equal = 1;
320 params.plot = lin_evol_model.plot;
321 % plot single snapshot
322 %lin_evol_model.plot(rb_sim_data.X(:,1),lin_evol_model_data.grid, ...
326 params.title = 'reduced solution';
330 params.title = 'error';
332 lin_evol_model_data.grid,params)
334 % plot of reduced output
335 figure,plot(rb_sim_data.time,[rb_sim_data.Y;ds_sim_data.Y]);
336 legend('reduced output','detailed output');
340 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
341 case 5 % generate plots of detailed data for paper:
342 % step = 5 % figures of detailed simuations based on step 4
343 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
345 disp('generate plots of detailed_data for morepas poster');
347 load(fullfile(rbmatlabresult,'advection_fv_output_detailed_data_2traj'),...
348 'model','detailed_data')
350 % model = fast_model(params);
351 % rbfname = 'advection_fv_output_basis';
352 % load(fullfile(rbmatlabresult,rbfname),'ds_model','ds_model_data');
354 % make sure, that these trajectories exist.
355 trajectories = [3,4];
357 for tr = 1:length(trajectories);
358 trname = ['trajectory',num2str(trajectories(tr))];
359 % plots of trajectory
360 load(fullfile(rbmatlabresult,trname));
361 lin_evol_model = ds_model.base_model;
362 lin_evol_model_data = gen_model_data(lin_evol_model);
365 params.show_colorbar = 0;
366 params.axis_equal = 1;
367 params.plot = lin_evol_model.plot;
368 lin_evol_model.plot(lin_evol_model_data.grid,ds_sim_data.X(:,1), ...
372 saveas(gcf,fullfile(rbmatlabresult,[trname,'_t1.eps']),'epsc');
375 lin_evol_model.plot(lin_evol_model_data.grid,ds_sim_data.X(:,64), ...
379 saveas(gcf,fullfile(rbmatlabresult,[trname,'_tmiddle.eps']),'epsc');
382 lin_evol_model.plot(lin_evol_model_data.grid,ds_sim_data.X(:,end), ...
386 saveas(gcf,fullfile(rbmatlabresult,[trname,'_tend.eps']),'epsc');
389 p = ds_model.
plot_sim_data(ds_model, ds_model_data, ds_sim_data,params);
390 % Ylim = get(gca,'Ylim');
391 % Ylim = [Ylim(1),Ylim(2)*1.1];
392 set(gca,'Ylim',[0,0.09]);
393 set(p(2),'Color',[0,0.8,0]);
394 saveas(gcf,fullfile(rbmatlabresult,[trname,'_output.eps']),'epsc');
399 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
400 % step = 6 % figures of reduced simuations based on step 4
401 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
402 case 6 % generate plots of reduced model for poster/paper
404 disp('generate plots of reduced simulations output estimation');
406 load(fullfile(rbmatlabresult,'advection_fv_output_detailed_data_2traj'),...
407 'model','detailed_data')
411 % load(fullfile(rbmatlabresult,'trajectory1'));
412 load(fullfile(rbmatlabresult,'trajectory2'));
413 ds_model.enable_error_estimator = 1;
414 % ds_model.error_estimation = 1;
415 % ds_model.RB_generation_mode = 'file_load'
416 % ds_model.RB_basis_filename = fullfile(...
417 % rbmatlabresult,'advection_fv_output_basis');
418 % detailed_data = gen_detailed_data(ds_model,ds_model_data);
419 reduced_data = gen_reduced_data(ds_model,detailed_data);
423 Ns = Ns(find(reduced_data.N>=Ns));
425 for nind = 1:length(Ns)
426 % perform reduced simulation
427 ds_model.N = Ns(nind);
428 rb_sim_data = rb_simulation(ds_model,reduced_data);
429 rb_sim_data = rb_reconstruction(ds_model,detailed_data,rb_sim_data)
432 params.show_colorbar = 1;
433 params.axis_equal = 1;
435 % plot of reduced output
437 p = plot(rb_sim_data.time,...
439 rb_sim_data.Y+rb_sim_data.DeltaY; ...
440 rb_sim_data.Y-rb_sim_data.DeltaY; ...
442 colors = {[1,0,0],[1,0,0],[1,0,0],[0,1,0]};
443 linestyles = {
'-',
'-.',
'-.',
'-'};
446 set(p(i),'Color',colors{i},
'Linestyle',linestyles{i},...
447 'Linewidth',widths(i));
449 % set marker on rb_sim_data:
450 Xdata =
get(p(1),
'XData');
451 Ydata =
get(p(1),
'YData');
452 Zdata =
get(p(1),
'ZData');
453 set(p(1),
'XData',Xdata(1:4:end));
454 set(p(1),
'YData',Ydata(1:4:end));
455 set(p(1),
'ZData',Zdata(1:4:end));
456 set(p(1),
'Marker',
'o');
458 legend([
'y\^(t): reduced output'],...
459 'y\^(t)+\Delta_y(t)',...
460 'y\^(t)-\Delta_y(t)',...
461 [
'y(t): detailed output'],
'Location',
'NorthWest');
462 title([
'k = ',num2str(ds_model.N)],
'Fontsize',15);
464 fullfile(rbmatlabresult,...
465 [
'output_estimation',num2str(ds_model.N),
'.eps']),
'epsc');
469 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
470 % step = 7 % time measurements
471 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
473 disp(
'time measurements for morepas poster');
477 t_detailed = zeros(1,nreps);
479 % initialize detailed model and simulate
481 load(fullfile(rbmatlabresult,
'advection_fv_output_detailed_data_2traj'),...
482 'model',
'detailed_data')
485 % lin_evol_model = advection_fv_output_model(params);
486 %lin_evol_model.debug = 1;
487 %lin_evol_model_data = gen_model_data(lin_evol_model);
489 % set some additional fields in model which will be copied into
490 % ds model depending on whether constants are known or not:
491 % estimate_constants = 0;
492 % if estimate_constants
493 % %%% if constant estimaton has to be performed:
494 % lin_evol_model.estimate_lin_ds_nmu = 4;
495 % lin_evol_model.estimate_lin_ds_nX = 10;
496 % lin_evol_model.estimate_lin_ds_nt = 4;
499 % % the following values are rough bounds, so could be
500 % % non-rigorous, then larger search must be performed.
501 % lin_evol_model.state_bound_constant_C1 = 1;
502 % lin_evol_model.output_bound_constant_C2 = 1.2916;
505 % ds_model = lin_ds_from_lin_evol_model(lin_evol_model);
506 % % set accelerated data functions from below
507 % ds_model.A_function_ptr = @(model,model_data) ...
508 % eval_affine_decomp(@A_components,...
509 % @A_coefficients,...
512 % ds_model.B_function_ptr = @(model,model_data) ...
513 % eval_affine_decomp(@B_components,...
514 % @B_coefficients,...
516 % %ds_model.u_function_ptr = @(params) 0.1; % define new
517 ds_model_data = gen_model_data(ds_model);
519 disp('skipping detailed simulations...')
523 disp(['repetition = ',num2str(r)])
524 %mu = [1,1,0.5]-rand(1,3)*0.01;
525 mu = [1,0.5]-rand(1,2)*0.01;
526 ds_model = ds_model.set_mu(ds_model,mu);
528 ds_sim_data = detailed_simulation(ds_model,ds_model_data);
529 t_detailed(1,r) = toc;
534 disp('--------------------------------------------------')
535 disp('detailed simulation')
537 disp(['mean t_detailed=',num2str(mean(t_detailed))]);
538 disp('--------------------------------------------------')
540 % reduced simulation with and without error estimation
542 clear('ds_sim_data');
543 ds_model.RB_generation_mode = 'file_load';
544 ds_model.RB_basis_filename = ...
545 fullfile(rbmatlabresult,'advection_fv_output_basis');
546 detailed_data = gen_detailed_data(ds_model,ds_model_data);
547 reduced_data = gen_reduced_data(ds_model,detailed_data);
549 Ns = [150:-10:10]; %,30,20,10];
550 Ns = Ns(find(reduced_data.N>=Ns));
552 t_reduced_with_err_est = zeros(length(Ns),nreps);
553 t_reduced_without_err_est = zeros(length(Ns),nreps);
555 % reduced simulation with and without error estimation
558 ds_model.error_estimation = err_est;
559 disp(['err_est= ',num2str(err_est)]);
560 for n_ind = 1:length(Ns)
561 disp(['M = ',num2str(Ns(n_ind))]);
562 ds_model.N = Ns(n_ind);
564 disp(['repetition = ',num2str(r)])
565 %mu = [1,1,1]-rand(1,3)*0.01;
566 mu = [1,1]-rand(1,2)*0.01;
567 ds_model = ds_model.set_mu(ds_model,mu);
569 rb_sim_data = rb_simulation(ds_model,reduced_data);
570 if size(rb_sim_data.Xr,1)~=Ns(n_ind)
571 disp('N dimension not set correctly!');
575 t_reduced_with_err_est(n_ind,r) = toc;
577 t_reduced_without_err_est(n_ind,r) = toc;
583 % output measurements
585 disp('--------------------------------------------------')
586 disp('detailed simulation')
588 disp(['mean t_detailed=',num2str(mean(t_detailed))]);
590 disp('--------------------------------------------------')
591 disp('reduced simulation with error estimation')
592 disp(t_reduced_with_err_est);
593 disp(['mean t_reduced_with_err_est=',...
594 num2str(mean(t_reduced_with_err_est,2)')]);
595 disp(['Ns =',num2str(Ns)]);
597 disp('--------------------------------------------------')
598 disp('reduced simulation without error estimation')
599 disp(t_reduced_without_err_est);
600 disp(['mean t_reduced_without_err_est=',...
601 num2str(mean(t_reduced_without_err_est,2)')]);
602 disp(['Ns =',num2str(Ns)]);
603 disp('--------------------------------------------------')
605 save(fullfile(rbmatlabresult,'time_measurements'));
607 % dim_x = 65536, nt = 2048
609 %mean t_detailed=64.3806
610 %--------------------------------------------------
611 %reduced simulation with error estimation
612 % 3.4542 3.5969 3.5659 3.4216 3.6170
613 % 3.3648 3.3816 3.4218 3.4391 3.4552
614 % 3.1643 3.1312 3.1244 3.1433 3.1720
615 % 3.0594 3.0496 3.0628 3.0486 3.0387
616 % 2.9993 2.9953 3.0078 2.9413 2.9603
617 % 2.9555 2.9093 2.9091 2.9510 2.9256%%%
619 %mean t_reduced_with_err_est=3.5311 3.4125 3.147 3.0518 2.9808 2.9301
620 %Ns =60 50 40 30 20 10
621 %--------------------------------------------------
622 %reduced simulation without error estimation
623 % 2.3975 2.4220 2.2924 2.4490 2.3302
624 % 2.3322 2.3198 2.3425 2.3056 2.3512
625 % 2.2227 2.2734 2.2210 2.2582 2.2504
626 % 2.2531 2.2180 2.2429 2.2212 2.2612
627 % 2.2283 2.1979 2.2581 2.2241 2.1957
628 % 2.2274 2.2014 2.2577 2.1856 2.2268
630 %mean t_reduced_without_err_est=2.3782 2.3303 2.2451 2.2393 2.2208 2.2198
631 %Ns =60 50 40 30 20 10
632 %--------------------------------------------------
634 case 8 % compare original and accelerated detailed model
636 disp('compare original and accelerated lin_ds models');
639 % initialize detailed model and simulate
640 lin_evol_model = advection_fv_output_model(params);
641 estimate_constants = 0;
642 if estimate_constants
643 %%% if constant estimaton has to be performed:
644 lin_evol_model.estimate_lin_ds_nmu = 4;
645 lin_evol_model.estimate_lin_ds_nX = 10;
646 lin_evol_model.estimate_lin_ds_nt = 4;
649 % the following values are rough bounds, so could be
650 % non-rigorous, then larger search must be performed.
651 lin_evol_model.state_bound_constant_C1 = 1;
652 lin_evol_model.output_bound_constant_C2 = 1.2916;
655 ds_model = lin_ds_from_lin_evol_model(lin_evol_model);
656 ds_model.cone_weight = lin_evol_model.cone_weight;
657 %ds_model.u_function_ptr = @(params) 0.1; % define new
658 ds_model_data = gen_model_data(ds_model);
660 % set accelerated data functions
661 ds_model_fast = lin_ds_from_lin_evol_model(lin_evol_model);
662 ds_model_fast.cone_weight = lin_evol_model.cone_weight;
663 ds_model_fast.A_function_ptr = @(model,model_data) ...
664 eval_affine_decomp(@A_components,...
668 ds_model_fast.B_function_ptr = @(model,model_data) ...
669 eval_affine_decomp(@B_components,...
673 % ds_model_fast.C_function_ptr = @(model,model_data) ...
674 % eval_affine_decomp(@C_components,...
675 % @C_coefficients,...
678 % disp('skipping detailed simulations...')
681 %mu = [1,1,0.5]-rand(1,3)*0.01;
682 mu = [1,0.5]-rand(1,2)*0.01;
683 ds_model = ds_model.set_mu(ds_model,mu);
684 ds_model_fast = ds_model_fast.set_mu(ds_model_fast,mu);
686 ds_sim_data = detailed_simulation(ds_model,ds_model_data);
689 ds_sim_data_fast = detailed_simulation(ds_model_fast,ds_model_data);
690 t_detailed_fast = toc;
692 disp(['original t = ',num2str(t_detailed),...
693 ', accelerated t = ',num2str(t_detailed_fast)]);
694 if ~isequal(ds_sim_data,ds_sim_data_fast)
695 disp('simulation results of original and accelerated model differ!')
697 disp('simulation results of original and accelerated model equal!')
702 case 9 % generate basis from POD-
Greedy
705 ds_model = fast_model(params);
706 ds_model.cone_range = [0.4,0.6];
707 ds_model.base_model.cone_range = ds_model.cone_range;
709 % shrink range to 'single point'
710 % r = [0.4,0.6]; % => 2 -> 0.8 error nach 10 basisvektoren
711 % r = [0.0,1.0]; % => 4.6 -> 1.9 error nach 10 basisvektoren
712 % r = [0.45,0.55]; % => 2.38 -> 0.51 error nach 10 basisvektoren
713 %ds_model.mu_ranges = {[0,1],[0.4,0.6],[0.4,0.6]};
714 ds_model.mu_ranges = {[0.4,0.6],[0.4,0.6]};
715 % ds_model.rb_basis_generation_ptr = ...
716 % @basisgen_POD_greedy_uniform_fixed;
718 % ds_model.RB_generation_mode =
'uniform_fixed';
719 % ds_model.RB_num_intervals = [2,2,2];
721 ds_model_data = gen_model_data(ds_model);
723 disp(
'generating basis from POD-Greedy');
725 %
for basis generation, set further fields:
727 ds_model.verbose = 5; % => only dots in greedy loop
728 ds_model.RB_generation_mode =
'greedy_uniform_fixed';
729 ds_model.RB_numintervals = [2,1,1]; % => 0.5 in cone_pos obtained
730 ds_model.RB_numintervals = [1,1];
731 % ds_model.RB_numintervals = [4,4,4];
732 ds_model.RB_stop_epsilon = 1e-5;
733 % ds_model.RB_stop_timeout = 3600*6;
734 % ds_model.RB_stop_Nmax =150;
735 ds_model.RB_stop_timeout = 600;
736 ds_model.RB_stop_Nmax =10;
738 ds_model.RB_error_indicator =
'estimator';
739 ds_model.enable_error_estimator = 1;
741 % ds_model.get_estimator_from_sim_data = @(sim_data) sim_data.DeltaX(end);
743 detailed_data = gen_detailed_data(ds_model,ds_model_data);
746 dfname = fullfile(rbmatlabresult,
'advection_fv_output_lin_ds_detailed');
747 save(dfname,
'model',
'detailed_data');
749 disp(
'successfully generated and saved detailed_data');
751 params.plot_title =
'reduced basis';
752 params.plot = model.base_model.plot;
753 lin_evol_model_data = gen_model_data(model.base_model);
754 plot_sequence(detailed_data.V,lin_evol_model_data.grid,params)
758 case 10 % reduced simulations and error plots
760 %basis_from_step = 9; % step 9
761 basis_from_step = 4; % step 4
763 if basis_from_step == 9
764 dfname = fullfile(rbmatlabresult,'advection_fv_output_lin_ds_detailed');
766 disp('loaded basis from step 9')
768 load(fullfile(rbmatlabresult,'advection_fv_output_detailed_data_2traj'),...
769 'model','detailed_data')
771 % model = fast_model(params);
772 model.enable_error_estimator = 1;
773 model_data = gen_model_data(model);
774 % dfname = fullfile(rbmatlabresult,'advection_fv_output_detailed_data_146');
775 % model.RB_basis_filename = ...
776 % fullfile(rbmatlabresult,'advection_fv_output_basis');
777 % model.RB_generation_mode = 'file_load';
778 % detailed_data = gen_detailed_data(model,model_data);
779 % dfname = fullfile(rbmatlabresult,'advection_fv_output_detailed');
781 tmp = load(fullfile(rbmatlabresult,'trajectory2'));
782 ds_sim_data = tmp.ds_sim_data;
783 mu_from_file = get_mu(tmp.ds_model);
784 disp('loaded basis from step 4')
785 model = model.set_mu(model,mu_from_file);
787 % trajectory2 above is made with this
788 % model = model.set_mu(model,[1,1,1]);
789 % model = model.set_mu(model,[1,1,1]);
790 % trajectory3 above is made with this:
791 %model = model.set_mu(model,[1,1,0.5]);
792 %model = model.set_mu(model,[0.5,0.6,0.6]);
793 %model = model.set_mu(model,[0.5,0.4,0.4]);
794 %model = model.set_mu(model,[0.5,0.5,0.5]);
795 model.N = size(detailed_data.V,2);
797 reduced_data = gen_reduced_data(model,detailed_data);
798 rb_sim_data = rb_simulation(model,reduced_data);
801 params.show_colorbar = 1;
802 params.axis_equal = 1;
804 % plot of reduced output
806 p = plot(rb_sim_data.time,...
808 rb_sim_data.Y+rb_sim_data.DeltaY; ...
809 rb_sim_data.Y-rb_sim_data.DeltaY]);
810 colors = {[1,0,0],[1,0,0],[1,0,0],[0,1,0]};
811 linestyles = {
'-',
'-.',
'-.',
'-'};
814 set(p(i),'Color',colors{i},
'Linestyle',linestyles{i},...
815 'Linewidth',widths(i));
817 % set marker on rb_sim_data:
818 Xdata =
get(p(1),
'XData');
819 Ydata =
get(p(1),
'YData');
820 Zdata =
get(p(1),
'ZData');
821 set(p(1),
'XData',Xdata(1:4:end));
822 set(p(1),
'YData',Ydata(1:4:end));
823 set(p(1),
'ZData',Zdata(1:4:end));
824 set(p(1),
'Marker',
'o');
826 legend([
'y\^(t): reduced output'],...
827 'y\^(t)+\Delta_y(t)',...
828 'y\^(t)-\Delta_y(t)',
'Location',
'NorthWest');
829 title([
'k = ',num2str(model.N)],
'Fontsize',15);
831 % table of error, error estimator and effectivity
833 reduced_data = gen_reduced_data(model,detailed_data);
836 Ns = Ns(find(reduced_data.N>=Ns));
838 err = zeros(1,length(Ns));
839 err_estim = zeros(1,length(Ns));
840 effectivity = zeros(1,length(Ns));
842 for ni = 1:length(Ns)
844 rb_sim_data = rb_simulation(model, reduced_data);
845 err(ni) = abs(ds_sim_data.Y(end) - rb_sim_data.Y(end));
846 err_estim(ni) = rb_sim_data.DeltaY(end);
847 effectivity(ni) = err_estim(ni)/err(ni);
853 disp('finished computing effectivities')
856 case 11 % test of matlab ode solver for resulting sparse system
858 params.coarse_factor = 8; % nake model smaller
859 model = fast_model(params);
860 model.RB_generation_mode = 'from_detailed_data';
862 % detailed simulation:
863 model_data = gen_model_data(model);
864 model_data.RB = zeros(size(model_data.G,1),0);
865 detailed_data = gen_detailed_data(model,model_data);
867 model.decomp_mode = 0; % complete
868 model = model.set_time(model,0)
869 %model = model.set_mu(model,[0.5,1,1]);
870 model = model.set_mu(model,[1,1]);
872 A = @(model,model_data,t) ...
873 model.A_function_ptr(model.set_time(model,t),model_data);
875 B = @(model,model_data,t) ...
876 model.B_function_ptr(model.set_time(model,t),model_data);
878 x0 = model.x0_function_ptr(model,model_data);
880 [t,x] = ode15s(@(t,x) A(model,model_data,t)*x + B(model,model_data,t), ...
883 disp(['ode15 takes ',num2str(time_ode15),' sec']);
885 params.plot_title = 'ode15 solution';
886 params.plot = model.base_model.plot;
887 params.axis_equal = 1;
888 lin_evol_model_data = gen_model_data(model.base_model);
892 sim_data = detailed_simulation(model,model_data);
894 disp(['rbmatlab takes ',num2str(time_rbmatlab),' sec']);
895 params.plot_title = 'rbmatlab solution';
896 params.plot = model.base_model.plot;
897 params.axis_equal = 1;
898 lin_evol_model_data = gen_model_data(model.base_model);
901 % further plots for MCMCS paper
904 load(fullfile(rbmatlabresult,'advection_fv_output_detailed_data_2traj'),...
905 'model','detailed_data')
907 % plot of error estimator over parameter variation
908 % model = fast_model(params);
909 % model_data = gen_model_data(model);
910 % dfname = fullfile(rbmatlabresult,'advection_fv_output_detailed_data_146');
911 % model.RB_basis_filename = ...
912 % fullfile(rbmatlabresult,'advection_fv_output_basis');
913 % model.RB_generation_mode = 'file_load';
914 % detailed_data = gen_detailed_data(model,model_data);
915 % save(fullfile(rbmatlabresult,...
916 % 'advection_fv_output_detailed_2trajectories.mat'),...
917 % 'model','detailed_data')
919 % dfname = fullfile(rbmatlabresult,'advection_fv_output_detailed');
921 tmp = load(fullfile(rbmatlabresult,'trajectory2'));
922 ds_sim_data = tmp.ds_sim_data;
923 disp('loaded basis from step 4')
924 model.enable_error_estimator = 1;
925 % factors = 0:0.05:1;
926 reduced_data = gen_reduced_data(model,detailed_data);
928 Ns = Ns(find(reduced_data.N>=Ns));
930 err_estim = zeros(length(factors),length(Ns));
931 for ni = 1:length(Ns)
933 for i =1:length(factors)
934 %model = model.set_mu(model, factors(i)*[1,1,1]);
935 model = model.set_mu(model, factors(i)*[1,1]);
936 rb_sim_data = rb_simulation(model, reduced_data);
937 % err(ni) = abs(ds_sim_data.Y(end) - rb_sim_data.Y(end));
938 err_estim(i,ni) = rb_sim_data.DeltaY(end);
939 % effectivity(ni) = err_estim(ni)/err(ni);
943 p = plot(factors,err_estim');
944 colors = {[1,0,0],[0,0.6,0],[0,0,1],[0,0,0]};
945 linestyles = {
'-',
'-.',
'-',
'-.'};
946 for i = 1:min(length(Ns),4);
947 set(p(i),
'Linewidth',2,
'Color',colors{i},
'Linestyle',linestyles{i});
950 xlabel(
'parameter \mu_1=\mu_2=\mu_3');
951 ylabel(
'\Delta_y(\mu)');
952 title(
'error estimator from parameter sweep')
953 legend({
'k=1',
'k=8',
'k=16',
'k=32'},
'Location',
'NorthWest');
954 saveas(gcf,fullfile(rbmatlabresult,
'parameter_sweep.eps'),
'epsc');
957 disp(
'step number unknown');
960 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
961 % generation of ds_model from lin_evol_model and replace
963 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
965 function ds_model = fast_model(params);
967 lin_evol_model = advection_fv_output_model(params);
968 %lin_evol_model.debug = 1;
969 %lin_evol_model_data = gen_model_data(lin_evol_model);
971 % set some additional fields in model which will be copied into
972 % ds model depending on whether constants are known or not:
973 estimate_constants = 0;
974 if estimate_constants
975 %%%
if constant estimaton has to be performed:
976 lin_evol_model.estimate_lin_ds_nmu = 4;
977 lin_evol_model.estimate_lin_ds_nX = 10;
978 lin_evol_model.estimate_lin_ds_nt = 4;
981 % the following values are rough bounds, so could be
982 % non-rigorous, then larger search must be performed.
983 lin_evol_model.state_bound_constant_C1 = 1;
984 lin_evol_model.output_bound_constant_C2 = 1.2916;
987 ds_model = lin_ds_from_lin_evol_model(lin_evol_model);
988 ds_model.cone_weight = lin_evol_model.cone_weight;
989 % set accelerated data functions from below
990 ds_model.A_function_ptr = @(model,model_data) ...
991 eval_affine_decomp(@A_components,...
995 ds_model.B_function_ptr = @(model,model_data) ...
996 eval_affine_decomp(@B_components,...
999 %ds_model.u_function_ptr = @(params) 0.1; % define
new
1001 %ds_model.theta = 1;
1002 %disp(
'set theta to 1!!!');
1008 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1009 % auxialiary coefficient functions
for acceleration of
1010 % advection model:
explicit implementation of coefficients
1011 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1013 function Acomp = A_components(model,model_data)
1015 Acomp = model_data.A_components;
1017 function Acoeff = A_coefficients(model);
1020 %model.base_model.decomp_mode = 2;
1021 %model.base_model.t = model.t;
1022 %mu = get_mu(model);
1023 %model.base_model = model.set_mu(model.base_model,mu);
1024 %[L_I_coeff, L_E_coeff, b_coeff] = ...
1025 % model.base_model.operators_ptr(model.base_model, );
1026 %% L_E = Id + Delta t * A
1027 %Acoeff = L_E_coeff(2:end)/model.base_model.dt;
1029 %Acoeff = -model.base_model.dt * ...
1030 % model.base_model.velocity_coefficients_ptr(model.base_model);
1032 Acoeff = - [model.vx_weight, model.vy_weight]*(1-model.t);
1036 function Bcomp = B_components(model,model_data)
1038 Bcomp = model_data.B_components;
1040 function Bcoeff = B_coefficients(model);
1043 % hmmm the coefficients are now called twice in A and B
1044 % should somehow be cached later...
1045 %model.base_model.decomp_mode = 2;
1046 %model.base_model.t = model.t;
1047 %mu = get_mu(model);
1048 %model.base_model = model.set_mu(model.base_model,mu);
1049 %[L_I_coeff, L_E_coeff, b_coeff] = ...
1050 % model.base_model.operators_ptr(model.base_model, );
1051 %% L_E = Id + Delta t * A
1052 %Bcoeff2 = b_coeff/model.base_model.dt;
1053 vcoeff = - [model.vx_weight, model.vy_weight]*(1-model.t);
1054 Q_0 = model.base_model.cone_number;
1055 dcoeff = zeros(1,Q_0);
1056 max_pos = model.cone_weight * (Q_0-1)+1;
1059 dcoeff(q) = (1-(max_pos-q))*(1-t) * ((max_pos>=q) && (max_pos < q+1)) ...
1060 + (1+(max_pos-q))*(1-t) * ((max_pos>=q-1) && (max_pos < q));
1062 v = -(vcoeff(:)*(dcoeff(:)
'));
1068 %function Ccomp = C_components(model,model_data)
1069 %Ccomp = model_data.C_components;
1071 %function Ccoeff = C_coefficients(model);
1072 %model.base_model.decomp_mode = 2;
1073 %model.base_model.t = model.t;
1074 %mu = get_mu(model);
1075 %model.base_model = model.set_mu(model.base_model,mu);
1077 % model.base_model.operators_output(model.base_model);
function r = verbose(level, message, messageId)
This function displays messages depending on a message-id and/or a level. Aditionally you can set/res...
function [ RBext , dummy ] = RB_extension_PCA_fixspace(model, detailed_data)
function computing a RB basis extension for given parameters by the POD-Greedy algorithm.
function advection_fv_output(step)
small script implementing a simple advection example for producing matrices to be used in the RB-DS f...
function p = plot_sequence(varargin)
plotting a sequence of data slices on polygonal 2d grid (constructed from params if empty) and provid...
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.