3 % small script implementing a simple advection example
4 %
for producing matrices to be used in the RB-DS framework
5 % Discretization with FV Functions. The steps generate the
6 % results for the at-Automatisierungstechnik 2010 Paper
8 % mu_names = {
'cone_weight',
'vx_weight',
'vy_weight'};
13 % step = 1; % initialization of model and plot of model data
14 % step = 2; % detailed simulation of lin_evol
15 % step = 2.5; % conversion to DS primal model, lin_ds detailed simulation
16 % step = 3; % conversion to DS primal model, DS detailed simulation
17 % and comparison with lin_evol
18 % step = 4 % compute several trajectories.
19 % step = 4.5 % compute reduced basis based on trajectories.
20 % step = 4.75 % comparison of detailed and reduced simulation
21 % step = 5 % figures of detailed simulations based on step 4
22 % step = 6 % figures of reduced simulations based on step 4
23 % step = 7 % time measurements based on step 4
24 % step = 8; % compare original and accelerated ds_model
25 % step = 9; % generate greedy basis on subset of parameter domain
26 % step = 10; % reduced simulations and error plots
27 % step = 11; % test of matlab ode solver
for resulting sparse system
28 % step = 12; % further plots
for at-Automatisierungstechnik paper
30 % Bernard Haasdonk 25.8.2009
32 % structured triangular grid
40 %params.coarse_factor = 2; % fact 2 => ndofs = 16384, nt = 1024, 15s detailed
41 %params.coarse_factor = 1; % fact 1 => ndofs = 32768, nt = 2048, 60s detailed
42 params.coarse_factor = 1;
44 trajectory_fnbase =
'advection_fv_output_vconst_trajectory';
45 detailed_data_fname =
'advection_fv_output_vconst_detailed_data';
46 rb_fname =
'advection_fv_output_vconst_rb';
50 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
51 %step = 1; % initialization of model and plot of model data
53 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
55 %
verbose(1,
'initializaton of lin_evol_model and plot model_data');
56 disp(
'initializaton of lin_evol_model and plot model_data');
57 model = advection_fv_output_vconst_model(params);
58 model_data = gen_model_data(model);
59 params.axis_equal = 1;
60 params.axis_tight = 1;
61 plot(model_data.grid,params);
62 % inspect(model_data.grid)
63 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
64 %step = 2; % detailed simulation of lin_evol
65 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
67 %
verbose(1,
'detailed simulation of lin_evol_model and plot');
68 disp(
'detailed simulation of lin_evol_model and plot');
69 % params.coarse_factor = 8;
70 model = advection_fv_output_vconst_model(params);
71 model_data = gen_model_data(model);
73 % disp(
'model debugging turned on.')
74 % model.cone_weight = 0.5;
75 % model.vx_weight = 0.0;
76 % model.vy_weight = 0.0;
77 % model.vx_weight = 1.0;
78 % model.vy_weight = 1.0;
79 % model.cone_amplitude = 100;
80 % model.vx_weight = 1;
81 % model.vy_weight = 0;
82 sim_data = detailed_simulation(model,model_data);
85 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
86 %step = 2.5; % conversion to DS primal model, DS detailed simulation
87 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
89 disp(
'initialization of lin_ds model from lin_evol model');
90 %
verbose(1,
'initialization of lin_ds model from lin_evol model and ...
93 lin_evol_model = advection_fv_output_vconst_model(params);
94 %lin_evol_model.debug = 1;
95 lin_evol_model_data = gen_model_data(lin_evol_model);
97 %
set some additional fields in model which will be copied into
98 % ds model depending on whether constants are known or not:
99 estimate_constants = 0;
100 if estimate_constants
101 %%%
if constant estimaton has to be performed:
102 lin_evol_model.estimate_lin_ds_nmu = 4;
103 lin_evol_model.estimate_lin_ds_nX = 10;
104 lin_evol_model.estimate_lin_ds_nt = 4;
107 % the following values are rough bounds, so could be
108 % non-rigorous, then larger search must be performed.
109 lin_evol_model.state_bound_constant_C1 = 1;
110 lin_evol_model.output_bound_constant_C2 = 1.2916;
113 ds_model = lin_ds_from_lin_evol_model(lin_evol_model);
114 %ds_model.u_function_ptr = @(params) 0.1; % define
new
115 % ds_model.theta = 1;
117 ds_model_data = gen_model_data(ds_model);
119 %
set mu in model and base_model!!!
120 ds_model = ds_model.set_mu(ds_model,mu);
122 ds_sim_data = detailed_simulation(ds_model,ds_model_data);
126 params.plot_title=
'detailed implicit lin_ds';
127 lin_evol_from_ds_sim_data.U = ds_sim_data.X;
128 lin_evol_from_ds_sim_data.y = ds_sim_data.Y;
129 plot_sim_data(lin_evol_model,lin_evol_model_data,lin_evol_from_ds_sim_data,params);
131 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
132 %step = 3; % conversion to DS primal model, DS detailed simulation
133 % and comparison with lin_evol
134 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
136 case 3 % initialization of DS model and lin_ds detailed simulation
137 disp(
'initialization of lin_ds model from lin_evol model and plot');
138 %
verbose(1,
'initialization of lin_ds model from lin_evol model and ...
141 lin_evol_model = advection_fv_output_vconst_model(params);
142 %lin_evol_model.debug = 1;
143 lin_evol_model_data = gen_model_data(lin_evol_model);
145 %
set some additional fields in model which will be copied into
146 % ds model depending on whether constants are known or not:
147 estimate_constants = 0;
148 if estimate_constants
149 %%%
if constant estimaton has to be performed:
150 lin_evol_model.estimate_lin_ds_nmu = 4;
151 lin_evol_model.estimate_lin_ds_nX = 10;
152 lin_evol_model.estimate_lin_ds_nt = 4;
155 % the following values are rough bounds, so could be
156 % non-rigorous, then larger search must be performed.
157 lin_evol_model.state_bound_constant_C1 = 1;
158 lin_evol_model.output_bound_constant_C2 = 1.2916;
161 ds_model = lin_ds_from_lin_evol_model(lin_evol_model);
162 % ds_model.theta = 0.0;
163 %ds_model.u_function_ptr = @(params) 0.1; % define
new
164 ds_model_data = gen_model_data(ds_model);
166 %
set mu in model and base_model!!!
167 ds_model = ds_model.set_mu(ds_model,mu);
169 ds_sim_data = detailed_simulation(ds_model,ds_model_data);
170 lin_evol_model = lin_evol_model.set_mu(lin_evol_model,mu);
171 lin_evol_sim_data = detailed_simulation(lin_evol_model,...
172 lin_evol_model_data);
173 err = lin_evol_sim_data.U - ds_sim_data.X;
174 disp([
'max error in DOFs of lin_evol and ds simulation: ',...
175 num2str(max(abs(err(:))))]);
177 params.plot_title=
'detailed lin_evol';
179 plot_sim_data(lin_evol_model,lin_evol_model_data,lin_evol_sim_data,params);
181 params.plot_title=
'detailed lin_ds';
182 % lin_evol_from_ds_sim_data.U = ds_sim_data.X;
183 % lin_evol_from_ds_sim_data.y = ds_sim_data.Y;
184 %
plot_sim_data(lin_evol_model,lin_evol_model_data,lin_evol_from_ds_sim_data,params);
185 % lin_ds_from_lin_evol_plot_sim_data(ds_model,ds_model_data,ds_sim_data,params);
188 params.plot_title=
'error explicit lin_ds and lin_evol';
189 error_sim_data = lin_evol_sim_data;
190 error_sim_data.U = error_sim_data.U-ds_sim_data.X;
191 plot_sim_data(lin_evol_model,lin_evol_model_data,error_sim_data,params);
194 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
195 %step = 4 % reduced Basis generation by single trajectories and
196 %comparison of detailed and reduced simulation
197 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
199 case 4 % generate some trajectories
201 ds_model = fast_model(params);
202 % lin_evol_model = advection_fv_output_vconst_model(params);
203 %lin_evol_model.debug = 1;
204 % lin_evol_model_data = gen_model_data(lin_evol_model);
206 %
set some additional fields in model which will be copied into
207 % ds model depending on whether constants are known or not:
208 % estimate_constants = 0;
209 %
if estimate_constants
210 % %%%
if constant estimaton has to be performed:
211 % lin_evol_model.estimate_lin_ds_nmu = 4;
212 % lin_evol_model.estimate_lin_ds_nX = 10;
213 % lin_evol_model.estimate_lin_ds_nt = 4;
216 % % the following values are rough bounds, so could be
217 % % non-rigorous, then larger search must be performed.
218 % lin_evol_model.state_bound_constant_C1 = 1;
219 % lin_evol_model.output_bound_constant_C2 = 1.2916;
222 % ds_model = lin_ds_from_lin_evol_model(lin_evol_model);
223 %ds_model.u_function_ptr = @(params) 0.1; % define
new
224 ds_model_data = gen_model_data(ds_model);
226 % mu = [0.5,0.5,0.5];
228 mus = {[0,0,0],[1,1,1],[1,1,0.5],[0,0.5,1]};
229 for m = 1:length(mus)
230 disp(['computing trajectory number ',num2str(m)]);
232 ds_model = ds_model.set_mu(ds_model,mu);
235 %disp(
'set theta to 1!!!');
237 ds_sim_data = detailed_simulation(ds_model,ds_model_data);
239 disp([
'time for full simulation: ',num2str(t_det)]);
240 save(fullfile(rbmatlabresult,[trajectory_fnbase,num2str(m)]),...
241 'ds_model',
'ds_model_data',
'ds_sim_data');
244 case 4.5 % compute reduced basis.
245 % make sure that the previous trajectories are computed!!
247 disp(
'PCA of snapshots takes a while...');
248 load(fullfile(rbmatlabresult,[trajectory_fnbase,
'1']),...
249 'ds_model',
'ds_model_data',
'ds_sim_data');
250 % smaller epsilon does result in more than 1 vector
for mu=0,0,0
252 RB = PCA_fixspace(ds_sim_data.X,[],ds_model_data.G,[],[],epsilon);
253 %RB = orthonormalize_qr(ds_sim_data.X,ds_model_data.G);
254 % select correct orthogonal
set
255 K = RB
' * ds_model_data.G * RB;
256 i = find(abs(diag(K)-1)>1e-2);
258 disp(['found
',num2str(size(RB,2)),' basis vectors after traj. 1
']);
264 disp(['reduced to
',num2str(size(RB,2)),' basis vectors.
']);
267 load(fullfile(rbmatlabresult,[trajectory_fnbase,'2
']),...
268 'ds_model
','ds_model_data
','ds_sim_data
');
269 RB2 = PCA_fixspace(ds_sim_data.X,RB,ds_model_data.G,[],[],epsilon);
270 %RB = orthonormalize_qr(ds_sim_data.X,ds_model_data.G);
271 % select correct orthogonal set
272 K = RB2' * ds_model_data.G * RB2;
273 i = find(abs(diag(K)-1)>1e-2);
278 disp(['found ',num2str(size(RB2,2)),' basis vectors after traj. 2']);
280 RB2 = RB2(:,1:(i-1));
281 disp(['reduced to ',num2str(size(RB2,2)),' basis vectors after traj. 2']);
283 % RB = orthonormalize(RB);
286 W = ds_model_data.G * V;
287 save(fullfile(rbmatlabresult,rb_fname),'RB');
289 ds_model.RB_basis_filename = rb_fname;
290 ds_model.RB_generation_mode = 'file_load';
292 detailed_data = gen_detailed_data(ds_model,ds_model_data);
295 save(fullfile(rbmatlabresult,detailed_data_fname),...
296 'model','detailed_data')
298 case 4.75 % compare reduced and detailed simulation
300 load(fullfile(rbmatlabresult,detailed_data_fname),...
301 'model','detailed_data')
302 load(fullfile(rbmatlabresult,[trajectory_fnbase,'2']),...
307 %perform single reduced simulation for given parameter
308 reduced_data = gen_reduced_data(ds_model,detailed_data);
309 ds_model.N = reduced_data.N;
311 % ds_model.nt = ds_model.nt/4;
312 % disp('set nt to quarter!!');
313 rb_sim_data = rb_simulation(ds_model,reduced_data);
315 rb_sim_data = rb_reconstruction(ds_model,detailed_data,rb_sim_data)
316 disp(['time for reduced simulation: ',num2str(t_red)]);
318 % plot of reduced trajectory
319 lin_evol_model = ds_model.base_model;
320 lin_evol_model_data = gen_model_data(lin_evol_model);
323 params.show_colorbar = 1;
324 params.axis_equal = 1;
325 params.plot = lin_evol_model.plot;
326 % plot single snapshot
327 %lin_evol_model.plot(rb_sim_data.X(:,1),lin_evol_model_data.grid, ...
331 params.title = 'reduced solution';
335 params.title = 'error';
337 lin_evol_model_data.grid,params)
339 % plot of reduced output
340 figure,plot(rb_sim_data.time,[rb_sim_data.Y;ds_sim_data.Y]);
341 legend('reduced output','detailed output');
345 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
346 case 5 % generate plots of detailed data for paper:
347 % step = 5 % figures of detailed simuations based on step 4
348 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
350 load(fullfile(rbmatlabresult,detailed_data_fname),...
351 'model','detailed_data')
353 % model = fast_model(params);
354 % rb_fname = 'advection_fv_output_basis';
355 % load(fullfile(rbmatlabresult,rb_fname),'ds_model','ds_model_data');
357 % make sure, that these trajectories exist.
358 trajectories = [2,3,4];
360 for tr = 1:length(trajectories);
361 trname = [trajectory_fnbase,num2str(trajectories(tr))];
362 % plots of trajectory
363 load(fullfile(rbmatlabresult,trname));
364 lin_evol_model = ds_model.base_model;
365 lin_evol_model_data = gen_model_data(lin_evol_model);
368 params.show_colorbar = 0;
369 params.axis_equal = 1;
370 params.plot = lin_evol_model.plot;
371 lin_evol_model.plot(ds_sim_data.X(:,1),lin_evol_model_data.grid, ...
375 saveas(gcf,fullfile(rbmatlabresult,[trname,'_t1.eps']),'epsc');
378 lin_evol_model.plot(ds_sim_data.X(:,64),lin_evol_model_data.grid, ...
382 saveas(gcf,fullfile(rbmatlabresult,[trname,'_tmiddle.eps']),'epsc');
385 lin_evol_model.plot(ds_sim_data.X(:,end),lin_evol_model_data.grid, ...
389 saveas(gcf,fullfile(rbmatlabresult,[trname,'_tend.eps']),'epsc');
392 p = ds_model.
plot_sim_data(ds_model, ds_model_data, ds_sim_data,params);
393 % Ylim = get(gca,'Ylim');
394 % Ylim = [Ylim(1),Ylim(2)*1.1];
395 set(gca,'Ylim',[0,0.09]);
396 set(p(2),'Color',[0,0.8,0]);
397 saveas(gcf,fullfile(rbmatlabresult,[trname,'_output.eps']),'epsc');
402 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
403 % step = 6 % figures of reduced simuations based on step 4
404 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
405 case 6 % generate plots of reduced model for poster/paper
407 disp('generate plots of reduced simulations output estimation');
409 load(fullfile(rbmatlabresult,detailed_data_fname),...
410 'model','detailed_data')
414 % load(fullfile(rbmatlabresult,'trajectory1'));
415 load(fullfile(rbmatlabresult,[trajectory_fnbase,'2']));
417 % ds_model.error_estimation = 1;
418 % ds_model.RB_generation_mode = 'file_load'
419 % ds_model.RB_basis_filename = fullfile(...
420 % rbmatlabresult,'advection_fv_output_basis');
421 % detailed_data = gen_detailed_data(ds_model,ds_model_data);
422 reduced_data = gen_reduced_data(ds_model,detailed_data);
426 Ns = Ns(find(reduced_data.N>=Ns));
428 for nind = 1:length(Ns)
429 % perform reduced simulation
430 ds_model.N = Ns(nind);
431 rb_sim_data = rb_simulation(ds_model,reduced_data);
432 rb_sim_data = rb_reconstruction(ds_model,detailed_data,rb_sim_data)
435 params.show_colorbar = 1;
436 params.axis_equal = 1;
438 % plot of reduced output
440 p = plot(rb_sim_data.time,...
442 rb_sim_data.Y+rb_sim_data.DeltaY; ...
443 rb_sim_data.Y-rb_sim_data.DeltaY; ...
445 colors = {[1,0,0],[1,0,0],[1,0,0],[0,1,0]};
446 linestyles = {
'-',
'-.',
'-.',
'-'};
449 set(p(i),'Color',colors{i},
'Linestyle',linestyles{i},...
450 'Linewidth',widths(i));
452 %
set marker on rb_sim_data:
453 Xdata =
get(p(1),
'XData');
454 Ydata =
get(p(1),
'YData');
455 Zdata =
get(p(1),
'ZData');
456 set(p(1),
'XData',Xdata(1:4:end));
457 set(p(1),
'YData',Ydata(1:4:end));
458 set(p(1),
'ZData',Zdata(1:4:end));
459 set(p(1),
'Marker',
'o');
461 l = legend([
'y\^(t): reduced output'],...
462 'y\^(t)+\Delta_y(t)',...
463 'y\^(t)-\Delta_y(t)',...
464 [
'y(t): detailed output'],
'Location',
'NorthWest');
465 set(l,
'Fontsize',12);
466 title([
'k = ',num2str(ds_model.N)],
'Fontsize',20);
467 xlabel(
'time t',
'Fontsize',15);
468 ylabel(
'output',
'Fontsize',15);
471 fullfile(rbmatlabresult,...
472 [
'advection_fv_output_vconst_output_estimation',num2str(ds_model.N),
'.eps']),
'epsc');
476 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
477 % step = 7 % time measurements
478 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
483 t_detailed = zeros(1,nreps);
485 % initialize detailed model and simulate
487 load(fullfile(rbmatlabresult,detailed_data_fname),...
488 'model',
'detailed_data')
491 % lin_evol_model = advection_fv_output_vconst_model(params);
492 %lin_evol_model.debug = 1;
493 %lin_evol_model_data = gen_model_data(lin_evol_model);
495 % set some additional fields in model which will be copied into
496 % ds model depending on whether constants are known or not:
497 % estimate_constants = 0;
498 % if estimate_constants
499 % %%% if constant estimaton has to be performed:
500 % lin_evol_model.estimate_lin_ds_nmu = 4;
501 % lin_evol_model.estimate_lin_ds_nX = 10;
502 % lin_evol_model.estimate_lin_ds_nt = 4;
505 % % the following values are rough bounds, so could be
506 % % non-rigorous, then larger search must be performed.
507 % lin_evol_model.state_bound_constant_C1 = 1;
508 % lin_evol_model.output_bound_constant_C2 = 1.2916;
511 % ds_model = lin_ds_from_lin_evol_model(lin_evol_model);
512 % % set accelerated data functions from below
513 % ds_model.A_function_ptr = @(model,model_data) ...
514 % eval_affine_decomp(@A_components,...
515 % @A_coefficients,...
518 % ds_model.B_function_ptr = @(model,model_data) ...
519 % eval_affine_decomp(@B_components,...
520 % @B_coefficients,...
522 % %ds_model.u_function_ptr = @(params) 0.1; % define new
523 ds_model_data = gen_model_data(ds_model);
525 disp('skipping detailed simulations...')
529 disp(['repetition = ',num2str(r)])
530 mu = [1,1,0.5]-rand(1,3)*0.01;
531 ds_model = ds_model.set_mu(ds_model,mu);
533 ds_sim_data = detailed_simulation(ds_model,ds_model_data);
534 t_detailed(1,r) = toc;
539 disp('--------------------------------------------------')
540 disp('detailed simulation')
542 disp(['mean t_detailed=',num2str(mean(t_detailed))]);
543 disp('--------------------------------------------------')
545 % reduced simulation with and without error estimation
547 clear('ds_sim_data');
548 ds_model.RB_generation_mode = 'file_load';
549 ds_model.RB_basis_filename = ...
550 fullfile(rbmatlabresult,rb_fname);
551 detailed_data = gen_detailed_data(ds_model,ds_model_data);
552 reduced_data = gen_reduced_data(ds_model,detailed_data);
554 Ns = [150:-10:10]; %,30,20,10];
555 Ns = Ns(find(reduced_data.N>=Ns));
557 t_reduced_with_err_est = zeros(length(Ns),nreps);
558 t_reduced_without_err_est = zeros(length(Ns),nreps);
560 % reduced simulation with and without error estimation
563 ds_model.error_estimation = err_est;
564 disp(['err_est= ',num2str(err_est)]);
565 for n_ind = 1:length(Ns)
566 disp(['M = ',num2str(Ns(n_ind))]);
567 ds_model.N = Ns(n_ind);
569 disp(['repetition = ',num2str(r)])
570 mu = [1,1,1]-rand(1,3)*0.01;
571 ds_model = ds_model.set_mu(ds_model,mu);
573 rb_sim_data = rb_simulation(ds_model,reduced_data);
574 if size(rb_sim_data.Xr,1)~=Ns(n_ind)
575 disp('N dimension not set correctly!');
579 t_reduced_with_err_est(n_ind,r) = toc;
581 t_reduced_without_err_est(n_ind,r) = toc;
587 % output measurements
589 disp('--------------------------------------------------')
590 disp('detailed simulation')
592 disp(['mean t_detailed=',num2str(mean(t_detailed))]);
594 disp('--------------------------------------------------')
595 disp('reduced simulation with error estimation')
596 disp(t_reduced_with_err_est);
597 disp(['mean t_reduced_with_err_est=',...
598 num2str(mean(t_reduced_with_err_est,2)')]);
599 disp(['Ns =',num2str(Ns)]);
601 disp('--------------------------------------------------')
602 disp('reduced simulation without error estimation')
603 disp(t_reduced_without_err_est);
604 disp(['mean t_reduced_without_err_est=',...
605 num2str(mean(t_reduced_without_err_est,2)')]);
606 disp(['Ns =',num2str(Ns)]);
607 disp('--------------------------------------------------')
609 save(fullfile(rbmatlabresult,'time_measurements'));
611 % dim_x = 65536, nt = 2048
613 %mean t_detailed=64.3806
614 %--------------------------------------------------
615 %reduced simulation with error estimation
616 % 3.4542 3.5969 3.5659 3.4216 3.6170
617 % 3.3648 3.3816 3.4218 3.4391 3.4552
618 % 3.1643 3.1312 3.1244 3.1433 3.1720
619 % 3.0594 3.0496 3.0628 3.0486 3.0387
620 % 2.9993 2.9953 3.0078 2.9413 2.9603
621 % 2.9555 2.9093 2.9091 2.9510 2.9256%%%
623 %mean t_reduced_with_err_est=3.5311 3.4125 3.147 3.0518 2.9808 2.9301
624 %Ns =60 50 40 30 20 10
625 %--------------------------------------------------
626 %reduced simulation without error estimation
627 % 2.3975 2.4220 2.2924 2.4490 2.3302
628 % 2.3322 2.3198 2.3425 2.3056 2.3512
629 % 2.2227 2.2734 2.2210 2.2582 2.2504
630 % 2.2531 2.2180 2.2429 2.2212 2.2612
631 % 2.2283 2.1979 2.2581 2.2241 2.1957
632 % 2.2274 2.2014 2.2577 2.1856 2.2268
634 %mean t_reduced_without_err_est=2.3782 2.3303 2.2451 2.2393 2.2208 2.2198
635 %Ns =60 50 40 30 20 10
636 %--------------------------------------------------
638 case 8 % compare original and accelerated detailed model
640 disp('compare original and accelerated lin_ds models');
643 % initialize detailed model and simulate
644 lin_evol_model = advection_fv_output_vconst_model(params);
645 estimate_constants = 0;
646 if estimate_constants
647 %%% if constant estimaton has to be performed:
648 lin_evol_model.estimate_lin_ds_nmu = 4;
649 lin_evol_model.estimate_lin_ds_nX = 10;
650 lin_evol_model.estimate_lin_ds_nt = 4;
653 % the following values are rough bounds, so could be
654 % non-rigorous, then larger search must be performed.
655 lin_evol_model.state_bound_constant_C1 = 1;
656 lin_evol_model.output_bound_constant_C2 = 1.2916;
659 ds_model = lin_ds_from_lin_evol_model(lin_evol_model);
660 %ds_model.u_function_ptr = @(params) 0.1; % define new
661 ds_model_data = gen_model_data(ds_model);
663 % set accelerated data functions
664 ds_model_fast = lin_ds_from_lin_evol_model(lin_evol_model);
665 ds_model_fast.A_function_ptr = @(model,model_data) ...
666 eval_affine_decomp(@A_components,...
670 ds_model_fast.B_function_ptr = @(model,model_data) ...
671 eval_affine_decomp(@B_components,...
675 % ds_model_fast.C_function_ptr = @(model,model_data) ...
676 % eval_affine_decomp(@C_components,...
677 % @C_coefficients,...
680 % disp('skipping detailed simulations...')
683 mu = [1,1,0.5]-rand(1,3)*0.01;
684 ds_model = ds_model.set_mu(ds_model,mu);
685 ds_model_fast = ds_model_fast.set_mu(ds_model_fast,mu);
687 ds_sim_data = detailed_simulation(ds_model,ds_model_data);
690 ds_sim_data_fast = detailed_simulation(ds_model_fast,ds_model_data);
691 t_detailed_fast = toc;
693 disp(['original t = ',num2str(t_detailed),...
694 ', accelerated t = ',num2str(t_detailed_fast)]);
695 if ~isequal(ds_sim_data,ds_sim_data_fast)
696 disp('simulation results of original and accelerated model differ!')
698 disp('simulation results of original and accelerated model equal!')
703 case 9 % generate basis from POD-Greedy
706 ds_model = fast_model(params);
707 ds_model.cone_range = [0.4,0.6];
708 ds_model.base_model.cone_range = ds_model.cone_range;
710 % shrink range to 'single point'
711 % r = [0.4,0.6]; % => 2 -> 0.8 error nach 10 basisvektoren
712 % r = [0.0,1.0]; % => 4.6 -> 1.9 error nach 10 basisvektoren
713 % r = [0.45,0.55]; % => 2.38 -> 0.51 error nach 10 basisvektoren
714 ds_model.mu_ranges = {[0,1],[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 = [1,1,1]; % => 0.5 in cone_pos obtained
730 % ds_model.RB_numintervals = [4,4,4];
731 ds_model.RB_stop_epsilon = 1e-5;
732 ds_model.RB_stop_timeout = 3600;
733 ds_model.RB_stop_Nmax =150;
734 % ds_model.RB_stop_timeout = 600;
735 % ds_model.RB_stop_Nmax =10;
737 ds_model.RB_error_indicator =
'estimator';
739 % ds_model.get_estimator_from_sim_data = @(sim_data) sim_data.DeltaX(end);
741 detailed_data = gen_detailed_data(ds_model,ds_model_data);
744 dfname = fullfile(rbmatlabresult,[detailed_data_fname,
'_greedy']);
745 save(dfname,
'model',
'detailed_data');
747 disp(
'successfully generated and saved detailed_data');
749 params.plot_title =
'reduced basis';
750 params.plot = model.base_model.plot;
751 lin_evol_model_data = gen_model_data(model.base_model);
752 plot_sequence(detailed_data.V,lin_evol_model_data.grid,params)
756 case 10 % reduced simulations and error plots
758 %basis_from_step = 9; % step 9
759 basis_from_step = 4; % step 4
761 if basis_from_step == 9
762 dfname = fullfile(rbmatlabresult,[detailed_data_fname,'_greedy']);
764 disp('loaded basis from step 9')
766 load(fullfile(rbmatlabresult,detailed_data_fname),...
767 'model','detailed_data')
769 % model = fast_model(params);
770 model_data = gen_model_data(model);
771 % dfname = fullfile(rbmatlabresult,'advection_fv_output_detailed_data_146');
772 % model.RB_basis_filename = ...
773 % fullfile(rbmatlabresult,'advection_fv_output_basis');
774 % model.RB_generation_mode = 'file_load';
775 % detailed_data = gen_detailed_data(model,model_data);
776 % dfname = fullfile(rbmatlabresult,'advection_fv_output_detailed');
778 tmp = load(fullfile(rbmatlabresult,[trajectory_fnbase,'2']));
779 ds_sim_data = tmp.ds_sim_data;
780 mu_from_file = get_mu(tmp.ds_model);
781 disp('loaded basis from step 4')
782 model = model.set_mu(model,mu_from_file);
784 % trajectory2 above is made with this
785 % model = model.set_mu(model,[1,1,1]);
786 % model = model.set_mu(model,[1,1,1]);
787 % trajectory3 above is made with this:
788 %model = model.set_mu(model,[1,1,0.5]);
789 %model = model.set_mu(model,[0.5,0.6,0.6]);
790 %model = model.set_mu(model,[0.5,0.4,0.4]);
791 %model = model.set_mu(model,[0.5,0.5,0.5]);
792 model.N = size(detailed_data.V,2);
794 reduced_data = gen_reduced_data(model,detailed_data);
795 rb_sim_data = rb_simulation(model,reduced_data);
798 params.show_colorbar = 1;
799 params.axis_equal = 1;
801 % plot of reduced output
803 p = plot(rb_sim_data.time,...
805 rb_sim_data.Y+rb_sim_data.DeltaY; ...
806 rb_sim_data.Y-rb_sim_data.DeltaY]);
807 colors = {[1,0,0],[1,0,0],[1,0,0],[0,1,0]};
808 linestyles = {
'-',
'-.',
'-.',
'-'};
811 set(p(i),'Color',colors{i},
'Linestyle',linestyles{i},...
812 'Linewidth',widths(i));
814 %
set marker on rb_sim_data:
815 Xdata =
get(p(1),
'XData');
816 Ydata =
get(p(1),
'YData');
817 Zdata =
get(p(1),
'ZData');
818 set(p(1),
'XData',Xdata(1:4:end));
819 set(p(1),
'YData',Ydata(1:4:end));
820 set(p(1),
'ZData',Zdata(1:4:end));
821 set(p(1),
'Marker',
'o');
823 legend([
'y\^(t): reduced output'],...
824 'y\^(t)+\Delta_y(t)',...
825 'y\^(t)-\Delta_y(t)',
'Location',
'NorthWest');
826 title([
'k = ',num2str(model.N)],
'Fontsize',15);
828 % table of error, error estimator and effectivity
830 reduced_data = gen_reduced_data(model,detailed_data);
833 Ns = Ns(find(reduced_data.N>=Ns));
835 err = zeros(1,length(Ns));
836 err_estim = zeros(1,length(Ns));
837 effectivity = zeros(1,length(Ns));
839 for ni = 1:length(Ns)
841 rb_sim_data = rb_simulation(model, reduced_data);
842 err(ni) = abs(ds_sim_data.Y(end) - rb_sim_data.Y(end));
843 err_estim(ni) = rb_sim_data.DeltaY(end);
844 effectivity(ni) = err_estim(ni)/err(ni);
850 disp('finished computing effectivities')
853 case 11 % test of matlab ode solver for resulting sparse system
855 params.coarse_factor = 8; % nake model smaller
856 model = fast_model(params);
857 model.RB_generation_mode = 'from_detailed_data';
859 % detailed simulation:
860 model_data = gen_model_data(model);
861 model_data.RB = zeros(size(model_data.G,1),0);
862 detailed_data = gen_detailed_data(model,model_data);
864 model.decomp_mode = 0; % complete
865 model = model.set_time(model,0)
866 model = model.set_mu(model,[0.5,1,1]);
868 A = @(model,model_data,t) ...
869 model.A_function_ptr(model.set_time(model,t),model_data);
871 B = @(model,model_data,t) ...
872 model.B_function_ptr(model.set_time(model,t),model_data);
874 x0 = model.x0_function_ptr(model,model_data);
876 [t,x] = ode15s(@(t,x) A(model,model_data,t)*x + B(model,model_data,t), ...
879 disp(['ode15 takes ',num2str(time_ode15),' sec']);
881 params.plot_title = 'ode15 solution';
882 params.plot = model.base_model.plot;
883 params.axis_equal = 1;
884 lin_evol_model_data = gen_model_data(model.base_model);
888 sim_data = detailed_simulation(model,model_data);
890 disp(['rbmatlab takes ',num2str(time_rbmatlab),' sec']);
891 params.plot_title = 'rbmatlab solution';
892 params.plot = model.base_model.plot;
893 params.axis_equal = 1;
894 lin_evol_model_data = gen_model_data(model.base_model);
897 % further plots for MCMCS paper
900 load(fullfile(rbmatlabresult,detailed_data_fname),...
901 'model','detailed_data')
903 % plot of error estimator over parameter variation
904 % model = fast_model(params);
905 % model_data = gen_model_data(model);
906 % dfname = fullfile(rbmatlabresult,'advection_fv_output_detailed_data_146');
907 % model.RB_basis_filename = ...
908 % fullfile(rbmatlabresult,'advection_fv_output_basis');
909 % model.RB_generation_mode = 'file_load';
910 % detailed_data = gen_detailed_data(model,model_data);
911 % save(fullfile(rbmatlabresult,...
912 % 'advection_fv_output_detailed_2trajectories.mat'),...
913 % 'model','detailed_data')
915 % dfname = fullfile(rbmatlabresult,'advection_fv_output_detailed');
917 tmp = load(fullfile(rbmatlabresult,[trajectory_fnbase,'2']));
918 ds_sim_data = tmp.ds_sim_data;
919 disp('loaded basis from step 4')
921 % factors = 0:0.05:1;
922 reduced_data = gen_reduced_data(model,detailed_data);
924 Ns = Ns(find(reduced_data.N>=Ns));
926 err_estim = zeros(length(factors),length(Ns));
927 for ni = 1:length(Ns)
929 for i =1:length(factors)
930 model = model.set_mu(model, factors(i)*[1,1,1]);
931 rb_sim_data = rb_simulation(model, reduced_data);
932 % err(ni) = abs(ds_sim_data.Y(end) - rb_sim_data.Y(end));
933 err_estim(i,ni) = rb_sim_data.DeltaY(end);
934 % effectivity(ni) = err_estim(ni)/err(ni);
938 p = plot(factors,err_estim');
939 colors = {[1,0,0],[0,0.6,0],[0,0,1],[0,0,0]};
940 linestyles = {
'-',
'-.',
'-',
'-.'};
942 set(p(i),
'Linewidth',2,
'Color',colors{i},
'Linestyle',linestyles{i});
945 xlabel(
'parameter \mu_1=\mu_2=\mu_3');
946 ylabel(
'\Delta_y(\mu)');
947 title(
'error estimator from parameter sweep')
948 legend({
'k=1',
'k=8',
'k=16',
'k=32',
'k=64'},
'Location',
'NorthWest');
949 saveas(gcf,fullfile(rbmatlabresult,
'parameter_sweep.eps'),
'epsc');
952 disp(
'step number unknown');
955 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
956 % generation of ds_model from lin_evol_model and replace
958 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
960 function ds_model = fast_model(params);
962 lin_evol_model = advection_fv_output_vconst_model(params);
963 %lin_evol_model.debug = 1;
964 %lin_evol_model_data = gen_model_data(lin_evol_model);
966 %
set some additional fields in model which will be copied into
967 % ds model depending on whether constants are known or not:
968 estimate_constants = 0;
969 if estimate_constants
970 %%%
if constant estimaton has to be performed:
971 lin_evol_model.estimate_lin_ds_nmu = 4;
972 lin_evol_model.estimate_lin_ds_nX = 10;
973 lin_evol_model.estimate_lin_ds_nt = 4;
976 % the following values are rough bounds, so could be
977 % non-rigorous, then larger search must be performed.
978 lin_evol_model.state_bound_constant_C1 = 1;
979 lin_evol_model.output_bound_constant_C2 = 1.2916;
982 ds_model = lin_ds_from_lin_evol_model(lin_evol_model);
983 %
set accelerated data functions from below
984 ds_model.A_function_ptr = @(model,model_data) ...
985 eval_affine_decomp(@A_components,...
989 ds_model.B_function_ptr = @(model,model_data) ...
990 eval_affine_decomp(@B_components,...
993 %ds_model.u_function_ptr = @(params) 0.1; % define
new
996 %disp(
'set theta to 1!!!');
1002 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1003 % auxialiary coefficient functions
for acceleration of
1004 % advection model:
explicit implementation of coefficients
1005 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1007 function Acomp = A_components(model,model_data)
1009 Acomp = model_data.A_components;
1011 function Acoeff = A_coefficients(model);
1013 %model.base_model.decomp_mode = 2;
1014 %model.base_model.t = model.t;
1015 %mu = get_mu(model);
1016 %model.base_model = model.set_mu(model.base_model,mu);
1017 %[L_I_coeff, L_E_coeff, b_coeff] = ...
1018 % model.base_model.operators_ptr(model.base_model, );
1019 %% L_E = Id + Delta t * A
1020 %Acoeff = L_E_coeff(2:end)/model.base_model.dt;
1022 %Acoeff = -model.base_model.dt * ...
1023 % model.base_model.velocity_coefficients_ptr(model.base_model);
1025 Acoeff = - [model.vx_weight, model.vy_weight]*0.5; %*(1-model.t);
1029 function Bcomp = B_components(model,model_data)
1031 Bcomp = model_data.B_components;
1033 function Bcoeff = B_coefficients(model);
1035 % hmmm the coefficients are now called twice in A and B
1036 % should somehow be cached later...
1037 %model.base_model.decomp_mode = 2;
1038 %model.base_model.t = model.t;
1039 %mu = get_mu(model);
1040 %model.base_model = model.set_mu(model.base_model,mu);
1041 %[L_I_coeff, L_E_coeff, b_coeff] = ...
1042 % model.base_model.operators_ptr(model.base_model, );
1043 %% L_E = Id + Delta t * A
1044 %Bcoeff2 = b_coeff/model.base_model.dt;
1045 vcoeff = - [model.vx_weight, model.vy_weight]*0.5;%*(1-model.t);
1046 Q_0 = model.base_model.cone_number;
1047 dcoeff = zeros(1,Q_0);
1048 max_pos = model.base_model.cone_weight * (Q_0-1)+1;
1051 dcoeff(q) = (1-(max_pos-q))*(1-t) * ((max_pos>=q) && (max_pos < q+1)) ...
1052 + (1+(max_pos-q))*(1-t) * ((max_pos>=q-1) && (max_pos < q));
1054 v = -(vcoeff(:)*(dcoeff(:)'));
1055 Bcoeff = v(:)*model.cone_amplitude;
1060 %
function Ccomp = C_components(model,model_data)
1061 %Ccomp = model_data.C_components;
1063 %
function Ccoeff = C_coefficients(model);
1064 %model.base_model.decomp_mode = 2;
1065 %model.base_model.t = model.t;
1066 %mu = get_mu(model);
1067 %model.base_model = model.set_mu(model.base_model,mu);
1069 % model.base_model.operators_output(model.base_model);