1 function res = hmm_micro_local(step)
2 %
function performing the multiscale micro local experiments
4 % step 1: simply call Frederikes test
6 % step 2: generate rbmatlab micromodel and run detailed simulation
7 % as nonlin_evol. Should give same result and not too bad
8 % slowing down compared to step 1
12 % step 3: generate approximation: empirical interpolation
14 % very coarse grid: only corners of parameter square
15 % just show feasibility error to detailed one and time measurement
16 % step 3.1: plot ei results
17 % step 4: check that ei approximation is good.
18 % step 5: greedy basis generation
19 % ( step 5.1: trajectory based basis generation: better than step 5)
20 % step 5.2: trajectory based basis generation with svd & gram schmidt:
21 % better than step 5.2
22 %
case 6: compare detailed, detailed_ei, detailed_ri_rb_proj and
23 % rb simulation. Should all be identical => perfect!!
24 %
case 7:
demo_rb_gui with basis from 5, 5.1. or 5.2
25 % => Nmax = 85, but no visual difference until N=46, M = 26
26 %
case 8: compare rb_sim_data with original
27 % time measurement of detailed versus reduced
28 % (
case 9: compare single evolution step with local and global
29 % evaluation of
operator for different M:
for debugging only.)
32 % step 13: as step 3 but large model
34 % very coarse grid: only corners of parameter square
35 % detailed one and time measurement
38 % -Welche Flussfunktion konkret? => u^3
39 % -Dein
test_lokal liefert oszillationen am rechten Rand, ist das
40 % in Ordnung? Insbesondere sehen die Oszillationen bei anderem n_x
41 % wesentlich anders aus => ja, Anfangsphaenomen
42 % -mu = [1,-5];% ul und ur => Oszillationen, nur 1/4 des
43 % Zeitintervalls vernuenftig, ebenso [1,-4]
44 % => ur auf [-3,-1] eingeschränkt
45 % -mu = [4,-1];% Loesung vernuenftig? Oszillatonen links
46 % -Momentane aufloesung: Ort 300, Zeit 7000.
47 % Koennen wir Ortsaufloesung hochschrauben z.B. auf 1000 oder
48 % 10000? Und die Zeitschrittanzahl d.h. Endzeit runterschrauben?
49 % -Echt lokale Auswertung oder Ueberdeckendes Intervall?
50 % -Was sind vernuenftige Parameter und Parameterbereice?
51 % epsilon und alpha sind problematisch weil sie das Gitter
52 % veraendern oder kriegen wir das
"skaliert" hin?
53 % Kriterium: nt und dx fest ueber ganzen parameterbereich
54 % - Was ist
"Kriterium" fuer auftreten von nichtklassischem Schock?
55 % -Was ist ein globales dt, das fuer den gesamten Parameterbereich
56 % funktioniert? Einfache Berechnung anhand der maximalwerte fuer ul und ur?
57 % -Wie Verwendung von HMM_UM_GESCHW?
58 % Bestimmung mittlerer Zustand um und Geschwindigkeit
59 % [um,v]=HMM_UM_GESCHW(u_kontrolle,dx,n_x,ul);
60 % -quantitative Beschreibung was passiert bei kleinem/grossem alpha?
63 % -nice 1D grid
class => done in simple way
64 % -insert
"stencil_size" in model and insert in gen_reduced_data and
65 % reduced_data_subset. stencil_size =1
for one-neigbour schemes,
66 % stencil_size = 2
for two neighbour stept
68 % B. Haasdonk 26.5.2010
75 %params.Time = 2.24*10^(-5);
76 %CRBfname =
'hmm_micro_local_CRB';
77 %detailedfname =
'hmm_micro_local_POD_single_trajectory';
81 params.Time = 2.24*10^(-5)/8400*200; % 100 timesteps
83 params.dt_factor = 0.05;
84 CRBfname =
'hmm_micro_local_CRB_small';
85 detailedfname =
'hmm_micro_local_POD_single_trajectory_small';
89 % params.Time = 2.24*10^(-5); % 100 timesteps
91 % params.dt_factor = 0.05;
92 CRBfname =
'hmm_micro_local_CRB_large';
93 detailedfname =
'hmm_micro_local_POD_single_trajectory_large';
96 model = hmm_micro_local_model(params);
98 model_data = gen_model_data(model);
99 plot_params = merge_model_plot_params(model,[]);
102 case 1 % step 1: simply call Frederikes test
104 [u,dx,fehler]=
test_lokal(2.24*10^(-5),4,-2,4,300);
105 % [u,dx,fehler]=
test_lokal(2.24*10^(-5),4,-5,4,300);
107 disp([
'computation time ',num2str(t),
' sec.']);
109 % [u,dx,fehler]=
test_lokal(2.24*10^(-5),4,-2,4,700);
111 % disp([
'computation time ',num2str(t),
' sec.']);
112 case 2 % call micromodel as nonlin-evol model, same result as step 1
113 % mu = [4,-2];% ul und ur : Simulation OK.
114 % mu = [1,-3];% ul und ur : Simulation OK
115 % mu = [1,-2];% ul und ur : Simulation OK
116 mu = [4,-3];% ul und ur : Simulation OK
117 % mu = [1,-1];% ul und ur : Simulation OK
118 % mu = [4,-5];% ul und ur : Simulation produces nans!
119 model = model.set_mu(model,mu);
121 sim_data = detailed_simulation(model,model_data);
124 disp([
'computation time ',num2str(t),
' sec.']);
126 case {3,13} % generate collateral reduced basis
for operator interpol
127 disp(
'constructing collateral reduced basis')
128 old_generation_mode = model.RB_generation_mode;
129 model.RB_generation_mode = 'none';
131 detailed_data = gen_detailed_data(model, model_data);
133 disp(['computation time for empirical interpolation:',num2str(t),' sec.']);
134 model.RB_generation_mode = old_generation_mode;
135 detailed_data.ei_info{1}.elapsed_time = t;
137 if isfield(model, 'model_type')
138 infix = model.model_type;
142 CRBfname = [model.name, infix, '_CRB_interpol.mat'];
144 save(fullfile(rbmatlabresult,CRBfname),...
145 'detailed_data','model');
146 hmm_micro_local(step+0.1); % plot_detailed_data;
149 load(fullfile(rbmatlabresult,CRBfname),...
150 'detailed_data',
'model');
151 detailed_data.RB = zeros(model.nx,1);
152 plot_detailed_data(model,detailed_data,plot_params);
155 case 4 %check that ei_interpol_detailed_simulation is as good as full
156 disp([
'comparison between detailed simulation with and without', ...
158 tmp = load(fullfile(rbmatlabresult,CRBfname));
159 detailed_data = tmp.detailed_data;
162 model = model.set_mu(model, mu_default);
163 sim_data = detailed_simulation(model, detailed_data);
164 plot_params.title =
'detailed simulation without interpolation';
166 model.M = size(detailed_data.QM{1},2);
168 model = model.set_mu(model, mu_default);
169 ei_sim_data = detailed_ei_simulation(model, detailed_data);
170 plot_params.title =
'detailed simulation with global empirical interpolation';
171 plot_sim_data(model, detailed_data, ei_sim_data, plot_params);
172 plot_params.title =
'error of global ei';
173 diff_data = sim_data;
174 diff_data.U = abs(diff_data.U - ei_sim_data.U);
175 diff_plot_params = plot_params;
176 diff_plot_params.range_lim = [min(diff_data.U(:)),max(diff_data.U(:))];
177 diff_plot_params.clim = [0, max(max(diff_data.U))];
178 plot_sim_data(model, detailed_data, diff_data, diff_plot_params);
179 disp([
'maximum l-infty error:',num2str(max(max(diff_data.U)))]);
181 ei_local_sim_data = nonlin_evol_detailed_local_ei_simulation(model, detailed_data);
182 plot_params.title =
'detailed simulation with global empirical interpolation';
183 plot_sim_data(model, detailed_data, ei_local_sim_data, plot_params);
184 plot_params.title =
'error of local ei';
185 diff_data = sim_data;
186 diff_data.U = abs(diff_data.U - ei_sim_data.U);
187 diff_plot_params = plot_params;
188 diff_plot_params.range_lim = [min(diff_data.U(:)),max(diff_data.U(:))];
189 diff_plot_params.clim = [0, max(max(diff_data.U))];
190 plot_sim_data(model, detailed_data, diff_data, diff_plot_params);
191 disp([
'maximum l-infty error:',num2str(max(max(diff_data.U)))]);
193 case {5,15} % generate reduced basis by greedy
198 % params.M = 296; % full model
202 % online_data = gen_online_data(model,detailed_data);
203 % sim_data = detailed_simulation(model,model_data);
204 % rb_sim_data = rb_simulation(model,online_data);
205 % rb_sim_data = rb_reconstruction(model,detailed_data,rb_sim_data);
206 %
plot_sim_data(model, model_data, rb_sim_data, plot_params);
208 case 5.1 % generate dummy RB by POD of one trajectory
210 tmp = load(fullfile(rbmatlabresult,CRBfname));
211 detailed_data = tmp.detailed_data;
215 model = model.set_mu(model, mu_default);
217 sim_data = detailed_simulation(model, detailed_data);
219 % the foloowing has not so good quality wrt initial data:
220 % RB = model.orthonormalize(model, model_data, sim_data.U);
221 % the following is better, but not well approximating riemann initdata:
222 % [RB,S,V] = svds(sim_data.U,300); % should be exact
223 % combination of inclusion of Riemann data components + svd of remainder
226 detailed_data.RB_info = [];
227 detailed_data.RB_info.M_train = [1 1;-3,-1];
228 RBinit = model.rb_init_data_basis(model, detailed_data);
229 Uorth = sim_data.U - RBinit * (RBinit
' * sim_data.U);
230 [RB,S,V] = svds(Uorth,201); % should be almost exact...
231 % perform extra orthonormalization, as SVD can destroy
232 % orthogonality to inital data
233 RB = RB - RBinit * (RBinit' * RB);
234 [RB,S,V] = svds(RB); % should be almost exact...
235 i = find(diag(S)>=eps)
239 % cut out lower triangle
240 Kerr = K - eye(size(K));
242 imat = (1:N)'*ones(1,N);
244 triamask = imat<=jmat;
245 errvec = sum((Kerr.*triamask).^2);
246 i = find(errvec>1e-10);
250 % if max(max((K-eye(size(K))).^2)>1e-3
252 % error('problem in orthogonality of reduced basis, check!')
254 disp(['generated Nmax = ',num2str(size(RB,2)),' RB vectors']);
256 plot_params.title = 'reduced basis by svd of single trajectory';
257 detailed_data.RB = [RBinit,RB];
258 plot_params.range_lim = 'auto';
259 plot_detailed_data(model,detailed_data,plot_params);
261 % detailedfname = 'hmm_micro_local_POD_single_trajectory';
262 save(fullfile(rbmatlabresult,detailedfname),...
263 'detailed_data','model','params');
265 %
plot_sim_data(model, detailed_data, sim_data, plot_params);
266 % model.M = size(detailed_data.QM{1},2);
267 % model = model.set_mu(model, mu_default);
268 % ei_sim_data = detailed_ei_simulation(model, detailed_data);
269 % plot_params.title =
'detailed simulation with empirical interpolation';
270 %
plot_sim_data(model, detailed_data, ei_sim_data, plot_params);
272 case 5.2 % generate dummy RB by POD of one trajectory
274 tmp = load(fullfile(rbmatlabresult,CRBfname));
275 detailed_data = tmp.detailed_data;
279 model = model.set_mu(model, mu_default);
281 sim_data = detailed_simulation(model, detailed_data);
283 % the foloowing has not so good quality wrt initial data:
284 % RB = model.orthonormalize(model, model_data, sim_data.U);
285 % the following is better, but not well approximating riemann initdata:
286 % [RB,S,V] = svds(sim_data.U,300); % should be exact
287 % combination of inclusion of Riemann data components + svd of remainder
290 detailed_data.RB_info = [];
291 detailed_data.RB_info.M_train = [1 1;-3,-1];
292 RBinit = model.rb_init_data_basis(model, detailed_data);
293 % [RBinit,S,V] = svds(RBinit,203); % should be almost exact...
294 RBinit = orthonormalize_gram_schmidt(RBinit,1,eps);
295 Uorth = sim_data.U - RBinit * (RBinit
' * sim_data.U);
296 [RBorth,S,V] = svds(Uorth,200); % should be almost exact...
297 % [RB,S,V] = svds(RB); % should be almost exact...
298 i = find(diag(S)>=eps)
299 RB = [RBinit,RBorth(:,i)];
300 RB = orthonormalize_gram_schmidt(RB,1,eps);
302 % cut out lower triangle
303 Kerr = K - eye(size(K));
305 imat = (1:N)
'*ones(1,N);
307 triamask = imat<=jmat;
308 errvec = sum((Kerr.*triamask).^2);
309 i = find(errvec>1e-10);
311 % disp(
'discarding some vectors!');
313 % RB = RB(:,1:i(1)-1);
315 if max(max((K-eye(size(K))).^2))>1e-8
317 warning(
'problem in orthogonality of reduced basis, check!')
318 % add gram schmidt step at end!
320 disp([
'generated Nmax = ',num2str(size(RB,2)),
' RB vectors']);
322 plot_params.title =
'reduced basis by svd of single trajectory';
323 detailed_data.RB = RB;
324 plot_params.range_lim =
'auto';
325 plot_detailed_data(model,detailed_data,plot_params);
327 % detailedfname =
'hmm_micro_local_POD_single_trajectory';
328 save(fullfile(rbmatlabresult,detailedfname),...
329 'detailed_data',
'model',
'params');
331 %
plot_sim_data(model, detailed_data, sim_data, plot_params);
332 % model.M = size(detailed_data.QM{1},2);
333 % model = model.set_mu(model, mu_default);
334 % ei_sim_data = detailed_ei_simulation(model, detailed_data);
335 % plot_params.title =
'detailed simulation with empirical interpolation';
336 %
plot_sim_data(model, detailed_data, ei_sim_data, plot_params);
338 case 6 %compare rb_sim_data with original all should be identical
340 % compute detailed_ei_rb_proj_simulation and compare with
341 % detailed_ei_simulation
342 % and reconstructed rb_simulation
344 load(fullfile(rbmatlabresult,detailedfname));
345 %model = hmm_micro_local_model([]);
347 reduced_data = gen_reduced_data(model, detailed_data);
353 % automatically max is chosen?
355 model.M = length(detailed_data.TM{1});
357 model.N = size(detailed_data.RB,2);
360 sim_data= detailed_simulation(model,detailed_data);
363 % model.range_lim = [-6,6];
364 % model.enable_range_limiting = 0;
365 ei_global_sim_data= detailed_ei_simulation(model,detailed_data);
366 ei_local_sim_data= ...
367 nonlin_evol_detailed_local_ei_simulation(model,detailed_data);
368 ei_diff_sim_data = ei_local_sim_data;
369 ei_diff_sim_data.U = ei_diff_sim_data.U - ei_global_sim_data.U;
370 % rb_proj_sim_data= nonlin_evol_detailed_rb_proj_simulation(model,detailed_data,zeros(size(detailed_data.RB,1),1));
371 ei_rb_proj_sim_data = detailed_ei_rb_proj_simulation(model,detailed_data);
372 rb_sim_data = rb_simulation(model, reduced_data);
373 rb_sim_data = rb_reconstruction(model,detailed_data,rb_sim_data);
375 plot_params.title =
'detailed simulation without ei';
378 plot_params.title =
'detailed simulation with global empirical interpolation';
379 plot_sim_data(model, detailed_data, ei_global_sim_data, plot_params);
381 plot_params.title =
'detailed simulation with local empirical interpolation';
382 plot_sim_data(model, detailed_data, ei_local_sim_data, plot_params);
384 plot_params.title =
'error of ei with local and global empirical interpolation';
385 plot_sim_data(model, detailed_data, ei_diff_sim_data, plot_params);
389 % disp(
'todo: detailed_rb_proj sim_data');
390 % plot_params.title =
'detailed simulation with rb projection';
391 %
plot_sim_data(model, detailed_data, rb_proj_sim_data, plot_params);
393 plot_params.title =
'detailed simulation with rb-proj empirical interpolation';
394 plot_sim_data(model, detailed_data, ei_rb_proj_sim_data, plot_params);
396 plot_params.title =
'reduced simulation';
397 plot_sim_data(model, detailed_data, rb_sim_data, plot_params);
403 % detailedfname =
'hmm_micro_local_POD_single_trajectory';
406 if ~exist(
'detailedfname',
'var')
409 if isempty(detailedfname)
410 if isfield(model, 'model_type')
411 infix = model.model_type;
415 detailedfname = [model.name, infix, '_detailed.mat'];
420 load(fullfile(rbmatlabresult,detailedfname));
421 % model = hmm_micro_local_model([]);
423 %disp('reduced simulation:')
424 %reduced_data = gen_reduced_data(model, detailed_data);
425 % model.N = size(detailed_data.RB,2);
429 % model.M = size(detailed_data.QM{1},2);
433 %reduced_data = extract_reduced_data_subset(model, reduced_data);
435 %rb_sim_data = rb_simulation(model, reduced_data);
437 %disp([
'time for online phase: t = ',num2str(t)]);
439 %disp(
'full simulation:')
441 %sim_data = detailed_simulation(model, detailed_data);
443 %disp([
'time for detailed simulation: t = ',num2str(t)]);
445 % plot_params.plot = model.plot;
446 % plot_params.no_lines = 0;
447 % plot_params.axis_equal = 0;
449 demo_rb_gui(model, detailed_data, [], plot_params);
452 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
455 %time measurement of detailed versus reduced
457 % detailedfname =
'hmm_micro_local_POD_single_trajectory';
459 load(fullfile(rbmatlabresult,detailedfname));
461 % model = hmm_micro_local_model([]);
462 model = hmm_micro_local_model(params);
464 reduced_data_full = gen_reduced_data(model, detailed_data);
472 % Ns = [1,2,4,8,16,32,64,128,202];
473 % Ms = [1,3,6,12,24,48,96,192,296];
474 Ns = [1,2,4,8,16,32,64,128,256];
475 Ms = [1,2,4,8,16,32,64,128,256];
476 i = find(Ns<=size(detailed_data.RB,2));
479 i = find(Ms<=length(detailed_data.TM{1}));
487 disp(['reduced simulation for N = ',num2str(Ns(i)),...
488 ', M = ',num2str(Ms(i))]);
491 reduced_data = extract_reduced_data_subset(model, ...
494 rb_sim_data = rb_simulation(model, reduced_data);
496 disp([' time for online phase: t = ',num2str(t)]);
499 rb_sim_data = rb_reconstruction(model,detailed_data,rb_sim_data);
500 plot_params.title = 'reduced simulation';
501 plot_sim_data(model, detailed_data, rb_sim_data, plot_params);
504 disp('full simulation:')
506 sim_data = detailed_simulation(model, detailed_data);
508 disp(['time for detailed simulation: t = ',num2str(t)]);
510 % --------------------------------------------------------------
512 % --------------------------------------------------------------
514 % reduced simulation for N = 1, M = 1
515 % time for online phase: t = 1.0789
516 %reduced simulation for N = 2, M = 3
517 % time for online phase: t = 1.3841
518 %reduced simulation for N = 4, M = 6
519 % time for online phase: t = 1.6556
520 %reduced simulation for N = 8, M = 12
521 % time for online phase: t = 2.0574
522 %reduced simulation for N = 16, M = 24
523 % time for online phase: t = 3.1641
524 %reduced simulation for N = 32, M = 48
525 % time for online phase: t = 5.5709
526 %reduced simulation for N = 64, M = 96
527 % time for online phase: t = 8.7702
528 %reduced simulation for N = 128, M = 192
529 % time for online phase: t = 20.2892
530 %reduced simulation for N = 202, M = 296
531 % time for online phase: t = 51.7955
533 % full simulation: H = 300
534 %time for detailed simulation: t = 7.9688
537 %%%%%%%%%%%%%% for debugging only
539 case 9 % compare local and global evaluation of operator
541 load(fullfile(rbmatlabresult,detailedfname));
543 model.decomp_mode = 0;
545 U = model.init_values_algorithm(model, detailed_data);
548 nt = 150; % => nonzero!
550 QM = detailed_data.QM{1}(:,1:M);
551 BM = detailed_data.BM{1}(1:M,1:M);
552 TM = detailed_data.TM{1}(1:M);
554 dt = model.T/model.nt;
557 Uglob = zeros(length(U),nt+1);
561 inc = model.L_E_local_ptr(model, detailed_data, Unew, []);
563 ei_inc = QM * (BM \ inc_local);
564 Unew = Unew - dt * ei_inc;
569 grid = detailed_data.grid;
571 eind_ext = index_ext(grid,...
573 model.local_stencil_size);
576 tmp = zeros(1,grid.nelements);
577 tmp(eind_ext) = 1:length(eind_ext);
578 eind_local = tmp(eind);
580 if ~isfield(detailed_data,
'grid_local_ext');
581 detailed_data.grid_local_ext = {...
582 gridpart(grid,eind_ext)};
585 Uloc = zeros(length(U),nt+1);
590 U_local_ext = Unew2(eind_ext);
591 inc_local2 = model.L_E_local_ptr(model, detailed_data, ...
592 U_local_ext, eind_local);
593 ei_inc2 = QM * (BM \ inc_local2);
594 Unew2 = Unew2 - dt * ei_inc2;
598 plot([inc_local,inc_local2]);
599 disp([inc_local
'; inc_local2'])
600 % maxerr = max(abs(inc_local-inc_local2));
601 maxerr = max(abs(Unew-Unew2));
602 disp([
'max error: ',num2str(maxerr)])
604 % error('should be 0!!');
607 % eind_ext7 = index_ext(grid,...
609 % model.local_stencil_size);
610 % U_local_ext7 = U(eind_ext7);
611 % tmp = zeros(1,grid.nelements);
612 % tmp(eind_ext7) = 1:length(eind_ext7);
613 % eind_local7 = tmp(eind7);
615 % detailed_data7 = detailed_data;
616 % detailed_data7.grid_local_ext = {...
617 % gridpart(grid,eind_ext7)};
619 % inc_local = model.L_E_local_ptr(model, detailed_data7, ...
620 % U_local_ext7, eind_local7);
623 % u_new = Phasen_lokal(U,...
624 % model.epsilon,model.alpha,...
625 % [],model.dx,dt);%%%%%
628 % ind_ext = detailed_data.grid_local_ext{1}.global_eind;
629 % [u_new2,u_index2] = Phasen_lokal(U_local_ext,model.epsilon,...
630 % model.alpha,ind_ext,model.dx,dt);
635 inc = model.L_E_local_ptr(model, detailed_data, Unew, []);
638 U_local_ext = Unew(eind_ext);
639 inc_local2 = model.L_E_local_ptr(model, detailed_data, ...
640 U_local_ext, eind_local);
644 if isequal(inc_local,inc_local2);
645 disp(
'RBmatlab local model seems to be fine');
647 disp(
'RBmatlab local model seems to be buggy');
651 % direct call of Frederikes Model:
653 u_new = Phasen_lokal(Unew,...
654 model.epsilon,model.alpha,...
655 [],model.dx,dt);%%%%%
658 ind_ext = detailed_data.grid_local_ext{1}.global_eind;
659 U_local_ext = Unew(eind_ext);
660 [u_new2,u_index2] = Phasen_lokal(U_local_ext,model.epsilon,...
661 model.alpha,ind_ext,model.dx,dt);
666 if isequal(u_new2,u_new(u_index2));
667 disp(
'Frederikes local mode seems to be fine');
669 disp(
'Frederikes local mode seems to be buggy');
679 error(
'step number unknown')
function [ u , dx , fehler ] = test_lokal(Time, ul, ur, alpha, n_x)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % 1D, skalare Erhaltungsgl...
function demo_rb_gui(varargin)
reduced basis demo with sliders
function step6_demo_rb_gui()
script comparing time for a reduced and a detailed simulation and starting the demonstration GUI ...
function p = plot_sim_data(model, model_data, sim_data, plot_params)
function performing the plot of the simulation results as specified in model.
function step5_rb_generation()
script constructing a reduced basis space