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 = {
'cone_weight',
'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 = 1;
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);
114 % set mu in model and base_model!!!
115 ds_model = ds_model.set_mu(ds_model,mu);
117 ds_sim_data = detailed_simulation(ds_model,ds_model_data);
121 params.plot_title='detailed implicit lin_ds';
122 lin_evol_from_ds_sim_data.U = ds_sim_data.X;
123 lin_evol_from_ds_sim_data.y = ds_sim_data.Y;
124 plot_sim_data(lin_evol_model,lin_evol_model_data,lin_evol_from_ds_sim_data,params);
126 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
127 %step = 3; % conversion to DS primal model, DS detailed simulation
128 % and comparison with lin_evol
129 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
131 case 3 % initialization of DS model and lin_ds detailed simulation
132 disp('initialization of lin_ds model from lin_evol model and plot');
133 %
verbose(1,'initialization of lin_ds model from lin_evol model and ...
136 lin_evol_model = advection_fv_output_model(params);
137 %lin_evol_model.debug = 1;
138 lin_evol_model_data = gen_model_data(lin_evol_model);
140 % set some additional fields in model which will be copied into
141 % ds model depending on whether constants are known or not:
142 estimate_constants = 0;
143 if estimate_constants
144 %%% if constant estimaton has to be performed:
145 lin_evol_model.estimate_lin_ds_nmu = 4;
146 lin_evol_model.estimate_lin_ds_nX = 10;
147 lin_evol_model.estimate_lin_ds_nt = 4;
150 % the following values are rough bounds, so could be
151 % non-rigorous, then larger search must be performed.
152 lin_evol_model.state_bound_constant_C1 = 1;
153 lin_evol_model.output_bound_constant_C2 = 1.2916;
156 ds_model = lin_ds_from_lin_evol_model(lin_evol_model);
157 % ds_model.theta = 0.0;
158 %ds_model.u_function_ptr = @(params) 0.1; % define new
159 ds_model_data = gen_model_data(ds_model);
161 % set mu in model and base_model!!!
162 ds_model = ds_model.set_mu(ds_model,mu);
164 ds_sim_data = detailed_simulation(ds_model,ds_model_data);
165 lin_evol_model = lin_evol_model.set_mu(lin_evol_model,mu);
166 lin_evol_sim_data = detailed_simulation(lin_evol_model,...
167 lin_evol_model_data);
168 err = lin_evol_sim_data.U - ds_sim_data.X;
169 disp(['max error in DOFs of lin_evol and ds simulation: ',...
170 num2str(max(abs(err(:))))]);
172 params.plot_title='detailed lin_evol';
174 plot_sim_data(lin_evol_model,lin_evol_model_data,lin_evol_sim_data,params);
175 params.plot_title='detailed lin_ds';
176 lin_evol_from_ds_sim_data.U = ds_sim_data.X;
177 lin_evol_from_ds_sim_data.y = ds_sim_data.Y;
178 plot_sim_data(lin_evol_model,lin_evol_model_data,lin_evol_from_ds_sim_data,params);
179 params.plot_title='error explicit lin_ds and lin_evol';
180 error_sim_data = lin_evol_sim_data;
181 error_sim_data.U = error_sim_data.U-ds_sim_data.X;
182 plot_sim_data(lin_evol_model,lin_evol_model_data,error_sim_data,params);
185 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
186 %step = 4 % reduced Basis generation by single trajectories and
187 %comparison of detailed and reduced simulation
188 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
190 case 4 % generate some trajectories
192 ds_model = fast_model(params);
193 % lin_evol_model = advection_fv_output_model(params);
194 %lin_evol_model.debug = 1;
195 % lin_evol_model_data = gen_model_data(lin_evol_model);
197 % set some additional fields in model which will be copied into
198 % ds model depending on whether constants are known or not:
199 % estimate_constants = 0;
200 % if estimate_constants
201 % %%% if constant estimaton has to be performed:
202 % lin_evol_model.estimate_lin_ds_nmu = 4;
203 % lin_evol_model.estimate_lin_ds_nX = 10;
204 % lin_evol_model.estimate_lin_ds_nt = 4;
207 % % the following values are rough bounds, so could be
208 % % non-rigorous, then larger search must be performed.
209 % lin_evol_model.state_bound_constant_C1 = 1;
210 % lin_evol_model.output_bound_constant_C2 = 1.2916;
213 % ds_model = lin_ds_from_lin_evol_model(lin_evol_model);
214 %ds_model.u_function_ptr = @(params) 0.1; % define new
215 ds_model_data = gen_model_data(ds_model);
217 % mu = [0.5,0.5,0.5];
219 mus = {[0,0,0],[1,1,1],[1,1,0.5],[0,0.5,1]};
220 for m = 1:length(mus)
221 disp(['computing trajectory number ',num2str(m)]);
223 ds_model = ds_model.set_mu(ds_model,mu);
226 %disp(
'set theta to 1!!!');
228 ds_sim_data = detailed_simulation(ds_model,ds_model_data);
230 disp([
'time for full simulation: ',num2str(t_det)]);
231 save(fullfile(rbmatlabresult,[
'trajectory',num2str(m)]),...
232 'ds_model',
'ds_model_data',
'ds_sim_data');
235 case 4.5 % compute reduced basis.
236 % make sure that the previous trajectories are computed!!
238 disp(
'PCA of snapshots takes a while...');
239 load(fullfile(rbmatlabresult,
'trajectory1'),...
240 'ds_model',
'ds_model_data',
'ds_sim_data');
241 % smaller epsilon does result in more than 1 vector
for mu=0,0,0
243 RB = PCA_fixspace(ds_sim_data.X,[],ds_model_data.G,[],[],epsilon);
244 %RB = orthonormalize_qr(ds_sim_data.X,ds_model_data.G);
245 % select correct orthogonal
set
246 K = RB
' * ds_model_data.G * RB;
247 i = find(abs(diag(K)-1)>1e-2);
249 disp(['found
',num2str(size(RB,2)),' basis vectors after traj. 1
']);
255 disp(['reduced to
',num2str(size(RB,2)),' basis vectors.
']);
258 load(fullfile(rbmatlabresult,'trajectory2
'),...
259 'ds_model
','ds_model_data
','ds_sim_data
');
260 RB2 = PCA_fixspace(ds_sim_data.X,RB,ds_model_data.G,[],[],epsilon);
261 %RB = orthonormalize_qr(ds_sim_data.X,ds_model_data.G);
262 % select correct orthogonal set
263 K = RB2' * ds_model_data.G * RB2;
264 i = find(abs(diag(K)-1)>1e-2);
269 disp(['found ',num2str(size(RB2,2)),' basis vectors after traj. 2']);
271 RB2 = RB2(:,1:(i-1));
272 disp(['reduced to ',num2str(size(RB2,2)),' basis vectors after traj. 2']);
274 % RB = orthonormalize(RB);
277 W = ds_model_data.G * V;
278 rbfname = fullfile(rbmatlabresult,'advection_fv_output_basis');
281 ds_model.RB_basis_filename = rbfname;
282 ds_model.RB_generation_mode = 'file_load';
284 detailed_data = gen_detailed_data(ds_model,ds_model_data);
287 save(fullfile(rbmatlabresult,'advection_fv_output_detailed_data_2traj'),...
288 'model','detailed_data')
290 case 4.75 % compare reduced and detailed simulation
292 load(fullfile(rbmatlabresult,'advection_fv_output_detailed_data_2traj'),...
293 'model','detailed_data')
294 load(fullfile(rbmatlabresult,'trajectory2'),...
299 %perform single reduced simulation for given parameter
300 reduced_data = gen_reduced_data(ds_model,detailed_data);
301 ds_model.N = reduced_data.N;
303 % ds_model.nt = ds_model.nt/4;
304 % disp('set nt to quarter!!');
305 rb_sim_data = rb_simulation(ds_model,reduced_data);
307 rb_sim_data = rb_reconstruction(ds_model,detailed_data,rb_sim_data)
308 disp(['time for reduced simulation: ',num2str(t_red)]);
310 % plot of reduced trajectory
311 lin_evol_model = ds_model.base_model;
312 lin_evol_model_data = gen_model_data(lin_evol_model);
315 params.show_colorbar = 1;
316 params.axis_equal = 1;
317 params.plot = lin_evol_model.plot;
318 % plot single snapshot
319 %lin_evol_model.plot(rb_sim_data.X(:,1),lin_evol_model_data.grid, ...
323 params.title = 'reduced solution';
327 params.title = 'error';
329 lin_evol_model_data.grid,params)
331 % plot of reduced output
332 figure,plot(rb_sim_data.time,[rb_sim_data.Y;ds_sim_data.Y]);
333 legend('reduced output','detailed output');
337 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
338 case 5 % generate plots of detailed data for paper:
339 % step = 5 % figures of detailed simuations based on step 4
340 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
342 disp('generate plots of detailed_data for morepas poster');
344 load(fullfile(rbmatlabresult,'advection_fv_output_detailed_data_2traj'),...
345 'model','detailed_data')
347 % model = fast_model(params);
348 % rbfname = 'advection_fv_output_basis';
349 % load(fullfile(rbmatlabresult,rbfname),'ds_model','ds_model_data');
351 % make sure, that these trajectories exist.
352 trajectories = [3,4];
354 for tr = 1:length(trajectories);
355 trname = ['trajectory',num2str(trajectories(tr))];
356 % plots of trajectory
357 load(fullfile(rbmatlabresult,trname));
358 lin_evol_model = ds_model.base_model;
359 lin_evol_model_data = gen_model_data(lin_evol_model);
362 params.show_colorbar = 0;
363 params.axis_equal = 1;
364 params.plot = lin_evol_model.plot;
365 lin_evol_model.plot(ds_sim_data.X(:,1),lin_evol_model_data.grid, ...
369 saveas(gcf,fullfile(rbmatlabresult,[trname,'_t1.eps']),'epsc');
372 lin_evol_model.plot(ds_sim_data.X(:,64),lin_evol_model_data.grid, ...
376 saveas(gcf,fullfile(rbmatlabresult,[trname,'_tmiddle.eps']),'epsc');
379 lin_evol_model.plot(ds_sim_data.X(:,end),lin_evol_model_data.grid, ...
383 saveas(gcf,fullfile(rbmatlabresult,[trname,'_tend.eps']),'epsc');
386 p = ds_model.
plot_sim_data(ds_model, ds_model_data, ds_sim_data,params);
387 % Ylim = get(gca,'Ylim');
388 % Ylim = [Ylim(1),Ylim(2)*1.1];
389 set(gca,'Ylim',[0,0.09]);
390 set(p(2),'Color',[0,0.8,0]);
391 saveas(gcf,fullfile(rbmatlabresult,[trname,'_output.eps']),'epsc');
396 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
397 % step = 6 % figures of reduced simuations based on step 4
398 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
399 case 6 % generate plots of reduced model for poster/paper
401 disp('generate plots of reduced simulations output estimation');
403 load(fullfile(rbmatlabresult,'advection_fv_output_detailed_data_2traj'),...
404 'model','detailed_data')
408 % load(fullfile(rbmatlabresult,'trajectory1'));
409 load(fullfile(rbmatlabresult,'trajectory2'));
411 % ds_model.error_estimation = 1;
412 % ds_model.RB_generation_mode = 'file_load'
413 % ds_model.RB_basis_filename = fullfile(...
414 % rbmatlabresult,'advection_fv_output_basis');
415 % detailed_data = gen_detailed_data(ds_model,ds_model_data);
416 reduced_data = gen_reduced_data(ds_model,detailed_data);
420 Ns = Ns(find(reduced_data.N>=Ns));
422 for nind = 1:length(Ns)
423 % perform reduced simulation
424 ds_model.N = Ns(nind);
425 rb_sim_data = rb_simulation(ds_model,reduced_data);
426 rb_sim_data = rb_reconstruction(ds_model,detailed_data,rb_sim_data)
429 params.show_colorbar = 1;
430 params.axis_equal = 1;
432 % plot of reduced output
434 p = plot(rb_sim_data.time,...
436 rb_sim_data.Y+rb_sim_data.DeltaY; ...
437 rb_sim_data.Y-rb_sim_data.DeltaY; ...
439 colors = {[1,0,0],[1,0,0],[1,0,0],[0,1,0]};
440 linestyles = {
'-',
'-.',
'-.',
'-'};
443 set(p(i),'Color',colors{i},
'Linestyle',linestyles{i},...
444 'Linewidth',widths(i));
446 %
set marker on rb_sim_data:
447 Xdata =
get(p(1),
'XData');
448 Ydata =
get(p(1),
'YData');
449 Zdata =
get(p(1),
'ZData');
450 set(p(1),
'XData',Xdata(1:4:end));
451 set(p(1),
'YData',Ydata(1:4:end));
452 set(p(1),
'ZData',Zdata(1:4:end));
453 set(p(1),
'Marker',
'o');
455 legend([
'y\^(t): reduced output'],...
456 'y\^(t)+\Delta_y(t)',...
457 'y\^(t)-\Delta_y(t)',...
458 [
'y(t): detailed output'],
'Location',
'NorthWest');
459 title([
'k = ',num2str(ds_model.N)],
'Fontsize',15);
461 fullfile(rbmatlabresult,...
462 [
'output_estimation',num2str(ds_model.N),
'.eps']),
'epsc');
466 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
467 % step = 7 % time measurements
468 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
470 disp(
'time measurements for morepas poster');
474 t_detailed = zeros(1,nreps);
476 % initialize detailed model and simulate
478 load(fullfile(rbmatlabresult,
'advection_fv_output_detailed_data_2traj'),...
479 'model',
'detailed_data')
482 % lin_evol_model = advection_fv_output_model(params);
483 %lin_evol_model.debug = 1;
484 %lin_evol_model_data = gen_model_data(lin_evol_model);
486 % set some additional fields in model which will be copied into
487 % ds model depending on whether constants are known or not:
488 % estimate_constants = 0;
489 % if estimate_constants
490 % %%% if constant estimaton has to be performed:
491 % lin_evol_model.estimate_lin_ds_nmu = 4;
492 % lin_evol_model.estimate_lin_ds_nX = 10;
493 % lin_evol_model.estimate_lin_ds_nt = 4;
496 % % the following values are rough bounds, so could be
497 % % non-rigorous, then larger search must be performed.
498 % lin_evol_model.state_bound_constant_C1 = 1;
499 % lin_evol_model.output_bound_constant_C2 = 1.2916;
502 % ds_model = lin_ds_from_lin_evol_model(lin_evol_model);
503 % % set accelerated data functions from below
504 % ds_model.A_function_ptr = @(model,model_data) ...
505 % eval_affine_decomp(@A_components,...
506 % @A_coefficients,...
509 % ds_model.B_function_ptr = @(model,model_data) ...
510 % eval_affine_decomp(@B_components,...
511 % @B_coefficients,...
513 % %ds_model.u_function_ptr = @(params) 0.1; % define new
514 ds_model_data = gen_model_data(ds_model);
516 disp('skipping detailed simulations...')
520 disp(['repetition = ',num2str(r)])
521 mu = [1,1,0.5]-rand(1,3)*0.01;
522 ds_model = ds_model.set_mu(ds_model,mu);
524 ds_sim_data = detailed_simulation(ds_model,ds_model_data);
525 t_detailed(1,r) = toc;
530 disp('--------------------------------------------------')
531 disp('detailed simulation')
533 disp(['mean t_detailed=',num2str(mean(t_detailed))]);
534 disp('--------------------------------------------------')
536 % reduced simulation with and without error estimation
538 clear('ds_sim_data');
539 ds_model.RB_generation_mode = 'file_load';
540 ds_model.RB_basis_filename = ...
541 fullfile(rbmatlabresult,'advection_fv_output_basis');
542 detailed_data = gen_detailed_data(ds_model,ds_model_data);
543 reduced_data = gen_reduced_data(ds_model,detailed_data);
545 Ns = [150:-10:10]; %,30,20,10];
546 Ns = Ns(find(reduced_data.N>=Ns));
548 t_reduced_with_err_est = zeros(length(Ns),nreps);
549 t_reduced_without_err_est = zeros(length(Ns),nreps);
551 % reduced simulation with and without error estimation
554 ds_model.error_estimation = err_est;
555 disp(['err_est= ',num2str(err_est)]);
556 for n_ind = 1:length(Ns)
557 disp(['M = ',num2str(Ns(n_ind))]);
558 ds_model.N = Ns(n_ind);
560 disp(['repetition = ',num2str(r)])
561 mu = [1,1,1]-rand(1,3)*0.01;
562 ds_model = ds_model.set_mu(ds_model,mu);
564 rb_sim_data = rb_simulation(ds_model,reduced_data);
565 if size(rb_sim_data.Xr,1)~=Ns(n_ind)
566 disp('N dimension not set correctly!');
570 t_reduced_with_err_est(n_ind,r) = toc;
572 t_reduced_without_err_est(n_ind,r) = toc;
578 % output measurements
580 disp('--------------------------------------------------')
581 disp('detailed simulation')
583 disp(['mean t_detailed=',num2str(mean(t_detailed))]);
585 disp('--------------------------------------------------')
586 disp('reduced simulation with error estimation')
587 disp(t_reduced_with_err_est);
588 disp(['mean t_reduced_with_err_est=',...
589 num2str(mean(t_reduced_with_err_est,2)')]);
590 disp(['Ns =',num2str(Ns)]);
592 disp('--------------------------------------------------')
593 disp('reduced simulation without error estimation')
594 disp(t_reduced_without_err_est);
595 disp(['mean t_reduced_without_err_est=',...
596 num2str(mean(t_reduced_without_err_est,2)')]);
597 disp(['Ns =',num2str(Ns)]);
598 disp('--------------------------------------------------')
600 save(fullfile(rbmatlabresult,'time_measurements'));
602 % dim_x = 65536, nt = 2048
604 %mean t_detailed=64.3806
605 %--------------------------------------------------
606 %reduced simulation with error estimation
607 % 3.4542 3.5969 3.5659 3.4216 3.6170
608 % 3.3648 3.3816 3.4218 3.4391 3.4552
609 % 3.1643 3.1312 3.1244 3.1433 3.1720
610 % 3.0594 3.0496 3.0628 3.0486 3.0387
611 % 2.9993 2.9953 3.0078 2.9413 2.9603
612 % 2.9555 2.9093 2.9091 2.9510 2.9256%%%
614 %mean t_reduced_with_err_est=3.5311 3.4125 3.147 3.0518 2.9808 2.9301
615 %Ns =60 50 40 30 20 10
616 %--------------------------------------------------
617 %reduced simulation without error estimation
618 % 2.3975 2.4220 2.2924 2.4490 2.3302
619 % 2.3322 2.3198 2.3425 2.3056 2.3512
620 % 2.2227 2.2734 2.2210 2.2582 2.2504
621 % 2.2531 2.2180 2.2429 2.2212 2.2612
622 % 2.2283 2.1979 2.2581 2.2241 2.1957
623 % 2.2274 2.2014 2.2577 2.1856 2.2268
625 %mean t_reduced_without_err_est=2.3782 2.3303 2.2451 2.2393 2.2208 2.2198
626 %Ns =60 50 40 30 20 10
627 %--------------------------------------------------
629 case 8 % compare original and accelerated detailed model
631 disp('compare original and accelerated lin_ds models');
634 % initialize detailed model and simulate
635 lin_evol_model = advection_fv_output_model(params);
636 estimate_constants = 0;
637 if estimate_constants
638 %%% if constant estimaton has to be performed:
639 lin_evol_model.estimate_lin_ds_nmu = 4;
640 lin_evol_model.estimate_lin_ds_nX = 10;
641 lin_evol_model.estimate_lin_ds_nt = 4;
644 % the following values are rough bounds, so could be
645 % non-rigorous, then larger search must be performed.
646 lin_evol_model.state_bound_constant_C1 = 1;
647 lin_evol_model.output_bound_constant_C2 = 1.2916;
650 ds_model = lin_ds_from_lin_evol_model(lin_evol_model);
651 %ds_model.u_function_ptr = @(params) 0.1; % define new
652 ds_model_data = gen_model_data(ds_model);
654 % set accelerated data functions
655 ds_model_fast = lin_ds_from_lin_evol_model(lin_evol_model);
656 ds_model_fast.A_function_ptr = @(model,model_data) ...
657 eval_affine_decomp(@A_components,...
661 ds_model_fast.B_function_ptr = @(model,model_data) ...
662 eval_affine_decomp(@B_components,...
666 % ds_model_fast.C_function_ptr = @(model,model_data) ...
667 % eval_affine_decomp(@C_components,...
668 % @C_coefficients,...
671 % disp('skipping detailed simulations...')
674 mu = [1,1,0.5]-rand(1,3)*0.01;
675 ds_model = ds_model.set_mu(ds_model,mu);
676 ds_model_fast = ds_model_fast.set_mu(ds_model_fast,mu);
678 ds_sim_data = detailed_simulation(ds_model,ds_model_data);
681 ds_sim_data_fast = detailed_simulation(ds_model_fast,ds_model_data);
682 t_detailed_fast = toc;
684 disp(['original t = ',num2str(t_detailed),...
685 ', accelerated t = ',num2str(t_detailed_fast)]);
686 if ~isequal(ds_sim_data,ds_sim_data_fast)
687 disp('simulation results of original and accelerated model differ!')
689 disp('simulation results of original and accelerated model equal!')
694 case 9 % generate basis from POD-Greedy
697 ds_model = fast_model(params);
698 ds_model.cone_range = [0.4,0.6];
699 ds_model.base_model.cone_range = ds_model.cone_range;
701 % shrink range to 'single point'
702 % r = [0.4,0.6]; % => 2 -> 0.8 error nach 10 basisvektoren
703 % r = [0.0,1.0]; % => 4.6 -> 1.9 error nach 10 basisvektoren
704 % r = [0.45,0.55]; % => 2.38 -> 0.51 error nach 10 basisvektoren
705 ds_model.mu_ranges = {[0,1],[0.4,0.6],[0.4,0.6]};
706 % ds_model.rb_basis_generation_ptr = ...
707 % @basisgen_POD_greedy_uniform_fixed;
709 % ds_model.RB_generation_mode =
'uniform_fixed';
710 % ds_model.RB_num_intervals = [2,2,2];
712 ds_model_data = gen_model_data(ds_model);
714 disp(
'generating basis from POD-Greedy');
716 %
for basis generation,
set further fields:
718 ds_model.verbose = 5; % => only dots in greedy loop
719 ds_model.RB_generation_mode =
'greedy_uniform_fixed';
720 ds_model.RB_numintervals = [2,1,1]; % => 0.5 in cone_pos obtained
721 % ds_model.RB_numintervals = [4,4,4];
722 ds_model.RB_stop_epsilon = 1e-5;
723 % ds_model.RB_stop_timeout = 3600*6;
724 % ds_model.RB_stop_Nmax =150;
725 ds_model.RB_stop_timeout = 600;
726 ds_model.RB_stop_Nmax =10;
728 ds_model.RB_error_indicator =
'estimator';
730 % ds_model.get_estimator_from_sim_data = @(sim_data) sim_data.DeltaX(end);
732 detailed_data = gen_detailed_data(ds_model,ds_model_data);
735 dfname = fullfile(rbmatlabresult,
'advection_fv_output_lin_ds_detailed');
736 save(dfname,
'model',
'detailed_data');
738 disp(
'successfully generated and saved detailed_data');
740 params.plot_title =
'reduced basis';
741 params.plot = model.base_model.plot;
742 lin_evol_model_data = gen_model_data(model.base_model);
743 plot_sequence(detailed_data.V,lin_evol_model_data.grid,params)
747 case 10 % reduced simulations and error plots
749 %basis_from_step = 9; % step 9
750 basis_from_step = 4; % step 4
752 if basis_from_step == 9
753 dfname = fullfile(rbmatlabresult,'advection_fv_output_lin_ds_detailed');
755 disp('loaded basis from step 9')
757 load(fullfile(rbmatlabresult,'advection_fv_output_detailed_data_2traj'),...
758 'model','detailed_data')
760 % model = fast_model(params);
761 model_data = gen_model_data(model);
762 % dfname = fullfile(rbmatlabresult,'advection_fv_output_detailed_data_146');
763 % model.RB_basis_filename = ...
764 % fullfile(rbmatlabresult,'advection_fv_output_basis');
765 % model.RB_generation_mode = 'file_load';
766 % detailed_data = gen_detailed_data(model,model_data);
767 % dfname = fullfile(rbmatlabresult,'advection_fv_output_detailed');
769 tmp = load(fullfile(rbmatlabresult,'trajectory2'));
770 ds_sim_data = tmp.ds_sim_data;
771 mu_from_file = get_mu(tmp.ds_model);
772 disp('loaded basis from step 4')
773 model = model.set_mu(model,mu_from_file);
775 % trajectory2 above is made with this
776 % model = model.set_mu(model,[1,1,1]);
777 % model = model.set_mu(model,[1,1,1]);
778 % trajectory3 above is made with this:
779 %model = model.set_mu(model,[1,1,0.5]);
780 %model = model.set_mu(model,[0.5,0.6,0.6]);
781 %model = model.set_mu(model,[0.5,0.4,0.4]);
782 %model = model.set_mu(model,[0.5,0.5,0.5]);
783 model.N = size(detailed_data.V,2);
785 reduced_data = gen_reduced_data(model,detailed_data);
786 rb_sim_data = rb_simulation(model,reduced_data);
789 params.show_colorbar = 1;
790 params.axis_equal = 1;
792 % plot of reduced output
794 p = plot(rb_sim_data.time,...
796 rb_sim_data.Y+rb_sim_data.DeltaY; ...
797 rb_sim_data.Y-rb_sim_data.DeltaY]);
798 colors = {[1,0,0],[1,0,0],[1,0,0],[0,1,0]};
799 linestyles = {
'-',
'-.',
'-.',
'-'};
802 set(p(i),'Color',colors{i},
'Linestyle',linestyles{i},...
803 'Linewidth',widths(i));
805 %
set marker on rb_sim_data:
806 Xdata =
get(p(1),
'XData');
807 Ydata =
get(p(1),
'YData');
808 Zdata =
get(p(1),
'ZData');
809 set(p(1),
'XData',Xdata(1:4:end));
810 set(p(1),
'YData',Ydata(1:4:end));
811 set(p(1),
'ZData',Zdata(1:4:end));
812 set(p(1),
'Marker',
'o');
814 legend([
'y\^(t): reduced output'],...
815 'y\^(t)+\Delta_y(t)',...
816 'y\^(t)-\Delta_y(t)',
'Location',
'NorthWest');
817 title([
'k = ',num2str(model.N)],
'Fontsize',15);
819 % table of error, error estimator and effectivity
821 reduced_data = gen_reduced_data(model,detailed_data);
824 Ns = Ns(find(reduced_data.N>=Ns));
826 err = zeros(1,length(Ns));
827 err_estim = zeros(1,length(Ns));
828 effectivity = zeros(1,length(Ns));
830 for ni = 1:length(Ns)
832 rb_sim_data = rb_simulation(model, reduced_data);
833 err(ni) = abs(ds_sim_data.Y(end) - rb_sim_data.Y(end));
834 err_estim(ni) = rb_sim_data.DeltaY(end);
835 effectivity(ni) = err_estim(ni)/err(ni);
841 disp('finished computing effectivities')
844 case 11 % test of matlab ode solver for resulting sparse system
846 params.coarse_factor = 8; % nake model smaller
847 model = fast_model(params);
848 model.RB_generation_mode = 'from_detailed_data';
850 % detailed simulation:
851 model_data = gen_model_data(model);
852 model_data.RB = zeros(size(model_data.G,1),0);
853 detailed_data = gen_detailed_data(model,model_data);
855 model.decomp_mode = 0; % complete
856 model = model.set_time(model,0)
857 model = model.set_mu(model,[0.5,1,1]);
859 A = @(model,model_data,t) ...
860 model.A_function_ptr(model.set_time(model,t),model_data);
862 B = @(model,model_data,t) ...
863 model.B_function_ptr(model.set_time(model,t),model_data);
865 x0 = model.x0_function_ptr(model,model_data);
867 [t,x] = ode15s(@(t,x) A(model,model_data,t)*x + B(model,model_data,t), ...
870 disp(['ode15 takes ',num2str(time_ode15),' sec']);
872 params.plot_title = 'ode15 solution';
873 params.plot = model.base_model.plot;
874 params.axis_equal = 1;
875 lin_evol_model_data = gen_model_data(model.base_model);
879 sim_data = detailed_simulation(model,model_data);
881 disp(['rbmatlab takes ',num2str(time_rbmatlab),' sec']);
882 params.plot_title = 'rbmatlab solution';
883 params.plot = model.base_model.plot;
884 params.axis_equal = 1;
885 lin_evol_model_data = gen_model_data(model.base_model);
888 % further plots for MCMCS paper
891 load(fullfile(rbmatlabresult,'advection_fv_output_detailed_data_2traj'),...
892 'model','detailed_data')
894 % plot of error estimator over parameter variation
895 % model = fast_model(params);
896 % model_data = gen_model_data(model);
897 % dfname = fullfile(rbmatlabresult,'advection_fv_output_detailed_data_146');
898 % model.RB_basis_filename = ...
899 % fullfile(rbmatlabresult,'advection_fv_output_basis');
900 % model.RB_generation_mode = 'file_load';
901 % detailed_data = gen_detailed_data(model,model_data);
902 % save(fullfile(rbmatlabresult,...
903 % 'advection_fv_output_detailed_2trajectories.mat'),...
904 % 'model','detailed_data')
906 % dfname = fullfile(rbmatlabresult,'advection_fv_output_detailed');
908 tmp = load(fullfile(rbmatlabresult,'trajectory2'));
909 ds_sim_data = tmp.ds_sim_data;
910 disp('loaded basis from step 4')
912 % factors = 0:0.05:1;
913 reduced_data = gen_reduced_data(model,detailed_data);
915 Ns = Ns(find(reduced_data.N>=Ns));
917 err_estim = zeros(length(factors),length(Ns));
918 for ni = 1:length(Ns)
920 for i =1:length(factors)
921 model = model.set_mu(model, factors(i)*[1,1,1]);
922 rb_sim_data = rb_simulation(model, reduced_data);
923 % err(ni) = abs(ds_sim_data.Y(end) - rb_sim_data.Y(end));
924 err_estim(i,ni) = rb_sim_data.DeltaY(end);
925 % effectivity(ni) = err_estim(ni)/err(ni);
929 p = plot(factors,err_estim');
930 colors = {[1,0,0],[0,0.6,0],[0,0,1],[0,0,0]};
931 linestyles = {
'-',
'-.',
'-',
'-.'};
933 set(p(i),
'Linewidth',2,
'Color',colors{i},
'Linestyle',linestyles{i});
936 xlabel(
'parameter \mu_1=\mu_2=\mu_3');
937 ylabel(
'\Delta_y(\mu)');
938 title(
'error estimator from parameter sweep')
939 legend({
'k=1',
'k=8',
'k=16',
'k=32'},
'Location',
'NorthWest');
940 saveas(gcf,fullfile(rbmatlabresult,
'parameter_sweep.eps'),
'epsc');
943 disp(
'step number unknown');
946 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
947 % generation of ds_model from lin_evol_model and replace
949 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
951 function ds_model = fast_model(params);
953 lin_evol_model = advection_fv_output_model(params);
954 %lin_evol_model.debug = 1;
955 %lin_evol_model_data = gen_model_data(lin_evol_model);
957 %
set some additional fields in model which will be copied into
958 % ds model depending on whether constants are known or not:
959 estimate_constants = 0;
960 if estimate_constants
961 %%%
if constant estimaton has to be performed:
962 lin_evol_model.estimate_lin_ds_nmu = 4;
963 lin_evol_model.estimate_lin_ds_nX = 10;
964 lin_evol_model.estimate_lin_ds_nt = 4;
967 % the following values are rough bounds, so could be
968 % non-rigorous, then larger search must be performed.
969 lin_evol_model.state_bound_constant_C1 = 1;
970 lin_evol_model.output_bound_constant_C2 = 1.2916;
973 ds_model = lin_ds_from_lin_evol_model(lin_evol_model);
974 %
set accelerated data functions from below
975 ds_model.A_function_ptr = @(model,model_data) ...
976 eval_affine_decomp(@A_components,...
980 ds_model.B_function_ptr = @(model,model_data) ...
981 eval_affine_decomp(@B_components,...
984 %ds_model.u_function_ptr = @(params) 0.1; % define
new
987 %disp(
'set theta to 1!!!');
993 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
994 % auxialiary coefficient functions
for acceleration of
995 % advection model:
explicit implementation of coefficients
996 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
998 function Acomp = A_components(model,model_data)
1000 Acomp = model_data.A_components;
1002 function Acoeff = A_coefficients(model);
1005 %model.base_model.decomp_mode = 2;
1006 %model.base_model.t = model.t;
1007 %mu = get_mu(model);
1008 %model.base_model = model.set_mu(model.base_model,mu);
1009 %[L_I_coeff, L_E_coeff, b_coeff] = ...
1010 % model.base_model.operators_ptr(model.base_model, );
1011 %% L_E = Id + Delta t * A
1012 %Acoeff = L_E_coeff(2:end)/model.base_model.dt;
1014 %Acoeff = -model.base_model.dt * ...
1015 % model.base_model.velocity_coefficients_ptr(model.base_model);
1017 Acoeff = - [model.vx_weight, model.vy_weight]*(1-model.t);
1021 function Bcomp = B_components(model,model_data)
1023 Bcomp = model_data.B_components;
1025 function Bcoeff = B_coefficients(model);
1028 % hmmm the coefficients are now called twice in A and B
1029 % should somehow be cached later...
1030 %model.base_model.decomp_mode = 2;
1031 %model.base_model.t = model.t;
1032 %mu = get_mu(model);
1033 %model.base_model = model.set_mu(model.base_model,mu);
1034 %[L_I_coeff, L_E_coeff, b_coeff] = ...
1035 % model.base_model.operators_ptr(model.base_model, );
1036 %% L_E = Id + Delta t * A
1037 %Bcoeff2 = b_coeff/model.base_model.dt;
1038 vcoeff = - [model.vx_weight, model.vy_weight]*(1-model.t);
1039 Q_0 = model.base_model.cone_number;
1040 dcoeff = zeros(1,Q_0);
1041 max_pos = model.cone_weight * (Q_0-1)+1;
1044 dcoeff(q) = (1-(max_pos-q))*(1-t) * ((max_pos>=q) && (max_pos < q+1)) ...
1045 + (1+(max_pos-q))*(1-t) * ((max_pos>=q-1) && (max_pos < q));
1047 v = -(vcoeff(:)*(dcoeff(:)'));
1053 %
function Ccomp = C_components(model,model_data)
1054 %Ccomp = model_data.C_components;
1056 %
function Ccoeff = C_coefficients(model);
1057 %model.base_model.decomp_mode = 2;
1058 %model.base_model.t = model.t;
1059 %mu = get_mu(model);
1060 %model.base_model = model.set_mu(model.base_model,mu);
1062 % model.base_model.operators_output(model.base_model);