rbmatlab  1.16.09
 All Classes Namespaces Files Functions Variables Modules Pages
hmm_micro_local.m
1 function res = hmm_micro_local(step)
2 % function performing the multiscale micro local experiments
3 %
4 % step 1: simply call Frederikes test
5 % => 7.9897 sec.
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
9 % => 9,1 sec
10 %
11 % SMALL MODEL:
12 % step 3: generate approximation: empirical interpolation
13 % takes 10 minutes...
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.)
30 %
31 % LARGE MODEL:
32 % step 13: as step 3 but large model
33 % takes 10 minutes...
34 % very coarse grid: only corners of parameter square
35 % detailed one and time measurement
36 
37 % Fragen Frederike:
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?
61 %
62 % todo:
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
67 
68 % B. Haasdonk 26.5.2010
69 
70 if (nargin <1)
71  step = 1;
72 end;
73 
74 % for full model:
75 %params.Time = 2.24*10^(-5);
76 %CRBfname = 'hmm_micro_local_CRB';
77 %detailedfname = 'hmm_micro_local_POD_single_trajectory';
78 
79 if step <= 10
80  % for debugging:
81  params.Time = 2.24*10^(-5)/8400*200; % 100 timesteps
82  params.n_x = 100;
83  params.dt_factor = 0.05;
84  CRBfname = 'hmm_micro_local_CRB_small';
85  detailedfname = 'hmm_micro_local_POD_single_trajectory_small';
86 else
87  % for debugging:
88  params = [];
89 % params.Time = 2.24*10^(-5); % 100 timesteps
90 % params.n_x = 100;
91 % params.dt_factor = 0.05;
92  CRBfname = 'hmm_micro_local_CRB_large';
93  detailedfname = 'hmm_micro_local_POD_single_trajectory_large';
94 end;
95 
96 model = hmm_micro_local_model(params);
97 %keyboard;
98 model_data = gen_model_data(model);
99 plot_params = merge_model_plot_params(model,[]);
100 
101 switch step
102  case 1 % step 1: simply call Frederikes test
103  tic;
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);
106  t = toc;
107  disp(['computation time ',num2str(t),' sec.']);
108 % figure;
109 % [u,dx,fehler]=test_lokal(2.24*10^(-5),4,-2,4,700);
110 % t = toc;
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);
120  tic;
121  sim_data = detailed_simulation(model,model_data);
122  t = toc;
123  plot_sim_data(model, model_data, sim_data, plot_params);
124  disp(['computation time ',num2str(t),' sec.']);
125 
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';
130  tic;
131  detailed_data = gen_detailed_data(model, model_data);
132  t = toc;
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;
136  if isempty(CRBfname)
137  if isfield(model, 'model_type')
138  infix = model.model_type;
139  else
140  infix = '';
141  end
142  CRBfname = [model.name, infix, '_CRB_interpol.mat'];
143  end
144  save(fullfile(rbmatlabresult,CRBfname),...
145  'detailed_data','model');
146  hmm_micro_local(step+0.1); % plot_detailed_data;
147 
148  case {3.1,13.1}
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);
153  close(gcf-2);
154 
155  case 4 %check that ei_interpol_detailed_simulation is as good as full
156  disp(['comparison between detailed simulation with and without', ...
157  ' interpolation']);
158  tmp = load(fullfile(rbmatlabresult,CRBfname));
159  detailed_data = tmp.detailed_data;
160  mu_default = [4,-1];
161 
162  model = model.set_mu(model, mu_default);
163  sim_data = detailed_simulation(model, detailed_data);
164  plot_params.title = 'detailed simulation without interpolation';
165  plot_sim_data(model, detailed_data, sim_data, plot_params);
166  model.M = size(detailed_data.QM{1},2);
167 
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)))]);
180 
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)))]);
192 
193  case {5,15} % generate reduced basis by greedy
194 % params.M = 1;
195 % keyboard;
196  model.verbose = 5;
197  model.debug = 3;
198 % params.M = 296; % full model
200  % disp('todo...');
201  % later:
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);
207 
208  case 5.1 % generate dummy RB by POD of one trajectory
209 
210  tmp = load(fullfile(rbmatlabresult,CRBfname));
211  detailed_data = tmp.detailed_data;
212  mu_default = [1,-3];
213  model.debug = 3;
214 
215  model = model.set_mu(model, mu_default);
216 
217  sim_data = detailed_simulation(model, detailed_data);
218 
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
224  %...
225 
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)
236  keyboard;
237  RB = RB(:,i);
238  K = RB'*RB;
239  % cut out lower triangle
240  Kerr = K - eye(size(K));
241  N = size(K,2);
242  imat = (1:N)'*ones(1,N);
243  jmat = imat';
244  triamask = imat<=jmat;
245  errvec = sum((Kerr.*triamask).^2);
246  i = find(errvec>1e-10);
247  if ~isempty(i)
248  RB = RB(:,1:i(1)-1);
249  end;
250 % if max(max((K-eye(size(K))).^2)>1e-3
251 % pcolor(K);
252 % error('problem in orthogonality of reduced basis, check!')
253 % end;
254  disp(['generated Nmax = ',num2str(size(RB,2)),' RB vectors']);
255 
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);
260 
261 % detailedfname = 'hmm_micro_local_POD_single_trajectory';
262  save(fullfile(rbmatlabresult,detailedfname),...
263  'detailed_data','model','params');
264 
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);
271 
272  case 5.2 % generate dummy RB by POD of one trajectory
273 
274  tmp = load(fullfile(rbmatlabresult,CRBfname));
275  detailed_data = tmp.detailed_data;
276  mu_default = [1,-3];
277  model.debug = 3;
278 
279  model = model.set_mu(model, mu_default);
280 
281  sim_data = detailed_simulation(model, detailed_data);
282 
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
288  %...
289 
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);
301  K = RB'*RB;
302  % cut out lower triangle
303  Kerr = K - eye(size(K));
304  N = size(K,2);
305  imat = (1:N)'*ones(1,N);
306  jmat = imat';
307  triamask = imat<=jmat;
308  errvec = sum((Kerr.*triamask).^2);
309  i = find(errvec>1e-10);
310 % if ~isempty(i)
311 % disp('discarding some vectors!');
312 % keyboard;
313 % RB = RB(:,1:i(1)-1);
314 % end;
315  if max(max((K-eye(size(K))).^2))>1e-8
316  pcolor(K);
317  warning('problem in orthogonality of reduced basis, check!')
318  % add gram schmidt step at end!
319  end;
320  disp(['generated Nmax = ',num2str(size(RB,2)),' RB vectors']);
321 
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);
326 
327 % detailedfname = 'hmm_micro_local_POD_single_trajectory';
328  save(fullfile(rbmatlabresult,detailedfname),...
329  'detailed_data','model','params');
330 
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);
337 
338  case 6 %compare rb_sim_data with original all should be identical
339 
340  % compute detailed_ei_rb_proj_simulation and compare with
341  % detailed_ei_simulation
342  % and reconstructed rb_simulation
343 
344  load(fullfile(rbmatlabresult,detailedfname));
345  %model = hmm_micro_local_model([]);
346 
347  reduced_data = gen_reduced_data(model, detailed_data);
348 
349  model.Mstrich = 0;
350  model.verbose = 0;
351  model.debug = 0;
352 
353  % automatically max is chosen?
354 % keyboard;
355  model.M = length(detailed_data.TM{1});
356 % model.M = 7;
357  model.N = size(detailed_data.RB,2);
358 
359 
360  sim_data= detailed_simulation(model,detailed_data);
361 
362 % model.M = 1;
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);
374 
375  plot_params.title = 'detailed simulation without ei';
376  plot_sim_data(model, detailed_data, sim_data, plot_params);
377 
378  plot_params.title = 'detailed simulation with global empirical interpolation';
379  plot_sim_data(model, detailed_data, ei_global_sim_data, plot_params);
380 
381  plot_params.title = 'detailed simulation with local empirical interpolation';
382  plot_sim_data(model, detailed_data, ei_local_sim_data, plot_params);
383 
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);
386 
387  % keyboard;
388 
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);
392 
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);
395 
396  plot_params.title = 'reduced simulation';
397  plot_sim_data(model, detailed_data, rb_sim_data, plot_params);
398 
399 
400 
401  case 7 % demo_rb_gui
402 
403 % detailedfname = 'hmm_micro_local_POD_single_trajectory';
404 
405  % disp('todo...');
406  if ~exist('detailedfname','var')
407  detailedfname = [];
408  end;
409  if isempty(detailedfname)
410  if isfield(model, 'model_type')
411  infix = model.model_type;
412  else
413  infix = '';
414  end
415  detailedfname = [model.name, infix, '_detailed.mat'];
416  end
417 
419 
420  load(fullfile(rbmatlabresult,detailedfname));
421 % model = hmm_micro_local_model([]);
422 
423  %disp('reduced simulation:')
424  %reduced_data = gen_reduced_data(model, detailed_data);
425  % model.N = size(detailed_data.RB,2);
426 % model.N = 2;
427 % model.M = 40;
428  % model.N = 2;
429  % model.M = size(detailed_data.QM{1},2);
430  model.Mstrich = 0;
431  model.verbose = 5;
432  model.debug = 3;
433  %reduced_data = extract_reduced_data_subset(model, reduced_data);
434  %tic;
435  %rb_sim_data = rb_simulation(model, reduced_data);
436  %t = toc;
437  %disp(['time for online phase: t = ',num2str(t)]);
438 
439  %disp('full simulation:')
440  %tic;
441  %sim_data = detailed_simulation(model, detailed_data);
442  %t = toc;
443  %disp(['time for detailed simulation: t = ',num2str(t)]);
444 
445 % plot_params.plot = model.plot;
446 % plot_params.no_lines = 0;
447 % plot_params.axis_equal = 0;
448 
449  demo_rb_gui(model, detailed_data, [], plot_params);
450 
451 
452 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
453 
454 case 8
455  %time measurement of detailed versus reduced
456 
457 % detailedfname = 'hmm_micro_local_POD_single_trajectory';
458 
459  load(fullfile(rbmatlabresult,detailedfname));
460 
461 % model = hmm_micro_local_model([]);
462  model = hmm_micro_local_model(params);
463 
464  reduced_data_full = gen_reduced_data(model, detailed_data);
465 
466 % keyboard;
467 
468  model.Mstrich = 0;
469  model.verbose = 0;
470  model.debug = 0;
471 
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));
477  Ns = Ns(i);
478  Ms = Ms(i);
479  i = find(Ms<=length(detailed_data.TM{1}));
480  Ns = Ns(i)
481  Ms = Ms(i)
482 
483  % Ns = [];
484 % Ms = [];
485 
486  for i = 1:length(Ms)
487  disp(['reduced simulation for N = ',num2str(Ns(i)),...
488  ', M = ',num2str(Ms(i))]);
489  model.M = Ms(i);
490  model.N = Ns(i);
491  reduced_data = extract_reduced_data_subset(model, ...
492  reduced_data_full);
493  tic;
494  rb_sim_data = rb_simulation(model, reduced_data);
495  t = toc;
496  disp([' time for online phase: t = ',num2str(t)]);
497  end;
498 
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);
502 
503 
504  disp('full simulation:')
505  tic;
506  sim_data = detailed_simulation(model, detailed_data);
507  t = toc;
508  disp(['time for detailed simulation: t = ',num2str(t)]);
509 
510 % --------------------------------------------------------------
511 % Runtime results:
512 % --------------------------------------------------------------
513 %
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
532 %
533 % full simulation: H = 300
534 %time for detailed simulation: t = 7.9688
535 
536 
537 %%%%%%%%%%%%%% for debugging only
538 
539  case 9 % compare local and global evaluation of operator
540 
541  load(fullfile(rbmatlabresult,detailedfname));
542 
543  model.decomp_mode = 0;
544 
545  U = model.init_values_algorithm(model, detailed_data);
546  M = 1
547 % nt = 149; % => 0
548  nt = 150; % => nonzero!
549 
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);
553  BMinv = inv(BM);
554  dt = model.T/model.nt;
555 
556  % global evaluation:
557  Uglob = zeros(length(U),nt+1);
558  Uglob(:,1) = U;
559  Unew = U;
560  for t = 1:nt
561  inc = model.L_E_local_ptr(model, detailed_data, Unew, []);
562  inc_local = inc(TM);
563  ei_inc = QM * (BM \ inc_local);
564  Unew = Unew - dt * ei_inc;
565  Uglob(:,t+1) = Unew;
566  end;
567 
568  % local evaluation:
569  grid = detailed_data.grid;
570  eind = TM;
571  eind_ext = index_ext(grid,...
572  eind,...
573  model.local_stencil_size);
574 
575 
576  tmp = zeros(1,grid.nelements);
577  tmp(eind_ext) = 1:length(eind_ext);
578  eind_local = tmp(eind);
579 
580  if ~isfield(detailed_data,'grid_local_ext');
581  detailed_data.grid_local_ext = {...
582  gridpart(grid,eind_ext)};
583  end;
584 
585  Uloc = zeros(length(U),nt+1);
586  Uloc(:,1) = U;
587 
588  Unew2 = U;
589  for t =1:nt
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;
595  Uloc(:,t+1) = Unew2;
596  end;
597 
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)])
603  if (maxerr>1e-10)
604  % error('should be 0!!');
605 
606 % eind7 = TM(1:7);
607 % eind_ext7 = index_ext(grid,...
608 % eind7,...
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);
614 %
615 % detailed_data7 = detailed_data;
616 % detailed_data7.grid_local_ext = {...
617 % gridpart(grid,eind_ext7)};
618 %
619 % inc_local = model.L_E_local_ptr(model, detailed_data7, ...
620 % U_local_ext7, eind_local7);
621 %
622 % % global:
623 % u_new = Phasen_lokal(U,...
624 % model.epsilon,model.alpha,...
625 % [],model.dx,dt);%%%%%
626 %
627 % % local:
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);
631 %
632 % u_new2'
633 % u_new(u_index2)'
634 
635  inc = model.L_E_local_ptr(model, detailed_data, Unew, []);
636  inc_local = inc(TM);
637 
638  U_local_ext = Unew(eind_ext);
639  inc_local2 = model.L_E_local_ptr(model, detailed_data, ...
640  U_local_ext, eind_local);
641 
642  save debug_step9;
643 
644  if isequal(inc_local,inc_local2);
645  disp('RBmatlab local model seems to be fine');
646  else
647  disp('RBmatlab local model seems to be buggy');
648  keyboard;
649  end;
650 
651  % direct call of Frederikes Model:
652  % global:
653  u_new = Phasen_lokal(Unew,...
654  model.epsilon,model.alpha,...
655  [],model.dx,dt);%%%%%
656 
657  % local:
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);
662 
663  u_new2'
664  u_new(u_index2)'
665 
666  if isequal(u_new2,u_new(u_index2));
667  disp('Frederikes local mode seems to be fine');
668  else
669  disp('Frederikes local mode seems to be buggy');
670  end;
671 
672  keyboard;
673 
674  end;
675 
676  keyboard;
677 
678  otherwise
679  error('step number unknown')
680 end;
681 
function [ u , dx , fehler ] = test_lokal(Time, ul, ur, alpha, n_x)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % 1D, skalare Erhaltungsgl...
Definition: test_lokal.m:17
function demo_rb_gui(varargin)
reduced basis demo with sliders
Definition: demo_rb_gui.m:17
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.
Definition: plot_sim_data.m:17
function step5_rb_generation()
script constructing a reduced basis space