rbmatlab  1.16.09
 All Classes Namespaces Files Functions Variables Modules Pages
t_part_script.m
1 function t_part_script(step,eps,Nmax)
2 %function t_part_script(step)
3 %
4 %script testing t-partition function with the linear advection problem
5 %
6 % step 1: testing fixed t-partition detailed simulations
7 % step 2: reduced simulation in t-partition mit selbstgebastelter RB
8 % step 3: generate t-partition RB with a fixed t-partition
9 % sollte immer gesetzt sein: t_model.t_part_rb_generation_cut_error = 1
10 % evtl. Fragen am Telefon klären :-)
11 % step 4: Test: genauigkeitsvergleich zwischen t-partition und normal
12 % ist aber obsolet und funktioniert auch nicht mehr
13 % step 5: Vergleich ohne t-partition und mit
14 % Basen mit und ohne t-partition werden geladen und reduzierte
15 % Simulationen durchgeführt. Dann Vergleich der Fehlerverläufe.
16 % step 6: Basisgenerierung mit adaptiver t-partitionszerlegung
17 %
18 % Markus Dihlmann 11.02.11
19 %
20 
21 
22 switch(step)
23 
24  case 1
25 
26  mu=[0;1;0];
27 
28  params.coarse_factor = 4;
29 
30  model = advection_fv_output_model(params);
31  model = set_mu(model, mu);
32  %model.save_time_indices = 0:model.nt;
33  model_data1 = gen_model_data(model);
34 
35  sim_data1 = detailed_simulation(model, model_data1);
36 
37  plot_params.plot_function = @plot_element_data;
38  plot_data_sequence(model, model_data1.grid, sim_data1.U, plot_params);
39  title('original solution')
40 
41  %wrap a t_part_model
42 
43  %t_part_range{1} = [0,floor(model.nt/3)];
44  %t_part_range{2} = [floor(model.nt/3),floor(2*model.nt/3)];
45  %t_part_range{3} = [floor(2*model.nt/3),model.nt];
46  t_part_range{1} = [0,floor(model.nt/2)];
47  t_part_range{2} = [floor(model.nt/2),model.nt];
48 
49  t_params.t_part_range= t_part_range;
50 
51  model = t_part_model(model,t_params);
52  model = set_mu(model,mu);
53  model_data = gen_model_data(model);
54 
55  %detailed_simulation for whole t-domain
56  sim_data_all = detailed_simulation(model, model_data);
57  plot_data_sequence(model, model_data.grid, sim_data_all.U, plot_params);
58  title('t-part whole traject.')
59 
60  U1 = sim_data1.U;
61  U2 = sim_data_all.U;
62  diff = fv_l2_error(U1,U2,model_data1.W);
63  disp(diff(end));
64 
65  %detailed_simulation for domain 2
66  model.t_part_for_simulation=2;
67  sim_data = detailed_simulation(model, model_data);
68  plot_data_sequence(model, model_data.grid, sim_data.U, plot_params);
69  title('t-part 2')
70  %U_end1=U1(:,floor(2*model.nt/3));
71  %U_end2=sim_data.U(:,end);
72 
73  %err = fv_l2_error(U_end1,U_end2,model_data1.grid,[]);
74  %disp(err);
75 
76  keyboard;
77 
78 
79  case 2
80 
81  load('detailed_data_coarse4_30_adapt.mat');
82 
83  %rb simulation
84  params.coarse_factor = 4;
85 
86  model = advection_fv_output_model(params);
87  model_data = gen_model_data(model);
88  reduced_data = gen_reduced_data(model, detailed_data);
89 
90  rb_sim_data1 = rb_simulation(model, reduced_data);
91  rb_sim_data1 = rb_reconstruction(model, detailed_data, rb_sim_data1);
92 
93 
94  %only one partition
95  model.rb_simulation = @lin_evol_rb_simulation_t_part;
96  t_model = t_part_model(model);
97  model_data = gen_model_data(t_model);
98 
99  %construct detailed_data
100  t_part_detailed_data{1}=detailed_data;
101  detailed_data =[];
102  detailed_data.t_part_detailed_data = t_part_detailed_data;
103 
104  %gen reduced_data
105  t_reduced_data = gen_reduced_data(t_model, detailed_data);
106  rb_sim_data2 = rb_simulation(t_model, t_reduced_data);
107  rb_sim_data2 = rb_reconstruction(t_model, detailed_data, rb_sim_data2);
108 
109 
110  error = fv_l2_error(rb_sim_data1.U, rb_sim_data2.U, model_data.W);
111  plot(error)
112 
113  %error is 0...
114  %-->two partitions with same RB
115  detailed_data.t_part_detailed_data{2} = detailed_data.t_part_detailed_data{1};
116  detailed_data.W = model_data.W;
117  detailed_data.grid = model_data.grid;
118  t_part_range{1} = [0,floor(model.nt/2)];
119  t_part_range{2} = [floor(model.nt/2),model.nt];
120  %t_part_range{3} = [floor(2*model.nt/3),model.nt];
121  t_params.t_part_range= t_part_range;
122  t_params.transition_model =1;
123 
124  t_model = t_part_model(model,t_params);
125  model_data = gen_model_data(t_model);
126  reduced_data = gen_reduced_data(t_model, detailed_data);
127 
128  rb_sim_data3 = rb_simulation(t_model, reduced_data);
129  rb_sim_data3 = rb_reconstruction(t_model, detailed_data, rb_sim_data3);
130 
131  error = fv_l2_error(rb_sim_data1.U, rb_sim_data3.U, model_data.W);
132  plot(error)
133 
134  keyboard;
135 
136 
137  case 3
138  %%%%%%%%%%%%%%%%%
139  %generate t-part reduced basis for fixed t-partition
140  %%%%%%%%%%%%%%%%%%%%%%%%%%5
141 
142 
143  % fix base model
144  params.coarse_factor = 4;
145  model = advection_fv_output_model(params);
146 
147  model.random_seed=1234;
148  model.RB_numintervals = [2];
149  model.RB_generation_mode = 'greedy_refined';%'greedy_uniform_fixed'; %
150  model.RB_refinement_mode = 'adaptive';
151  model.RB_error_indicator = 'estimator';%'error';%'estimator';
152  model.RB_detailed_train_savepath = fullfile(RBMATLABTEMP,'dump');
153  model.RB_extension_algorithm = @RB_extension_PCA_fixspace;
154  model.RB_stop_epsilon = 1e-2;%1e-1;
155  model.RB_stop_timeout = 24*60*60; % 1 minute
156  model.RB_stop_Nmax = 30;
157  model.RB_stop_max_val_train_ratio = 1;
158  model.RB_refinement_theta = 0.2;
159  model.RB_val_rand_seed = 543;
160  model.RB_M_val_size = 10;
161  model.RB_stop_max_refinement_level = 5;
162  model.RB_element_indicator_mode='nodes_cogs_skippedrefs';
163  model.RB_element_indicator_s_max = 1;
164 
165 
166  %rb_simulation without t-partition basis
167  model_data = gen_model_data(model);
168  detailed_data = gen_detailed_data(model, model_data);
169  %reduced_data = gen_reduced_data(model, detailed_data);
170  %sim_data = rb_simulation(model, reduced_data);
171  %sim_data_original = rb_reconstruction(model, detailed_data, sim_data);
172 
173  save('Greedybase_1p_30_coarse4_adapt.mat','detailed_data');
174 
175  model.rb_simulation = @lin_evol_rb_simulation_t_part;
176  % fix t-partitions
177  %t_part_range{1} = [0,floor(model.nt/2)];
178  %t_part_range{2} = [floor(model.nt/2),model.nt];
179  t_part_range{1} = [0,floor(model.nt/3)];
180  t_part_range{2} = [floor(model.nt/3),floor(2*model.nt/3)];
181  t_part_range{3} = [floor(2*model.nt/3),model.nt];
182  t_params.t_part_range= t_part_range;
183  %generate t-partition model
184  t_model = t_part_model(model,t_params);
185  t_model.basis_vector_overlap = 5; %Wieviele Basisvektoren aus der alten basis werden als initial Basis für die nächste Partition verwendet
186 
187  model_data = gen_model_data(t_model);
188 
189  %generate detailed_data without error cut
190  %t_model.t_part_rb_generation_cut_error = 0;
191  %detailed_data = gen_detailed_data(t_model, model_data);
192 
193  %simulate
194  %reduced_data = gen_reduced_data(t_model, detailed_data);
195  %sim_data = rb_simulation(t_model, reduced_data);
196  %sim_data = rb_reconstruction(t_model, detailed_data, sim_data);
197 
198  %save('Greedybase_50_t_part_coarse4_222.mat','detailed_data');
199 
200 
201  %generate detailed_data with error cut
202  t_model.t_part_rb_generation_cut_error = 1;
203  detailed_data = gen_detailed_data(t_model, model_data);
204 
205  %simulate
206  %reduced_data = gen_reduced_data(t_model, detailed_data);
207  %sim_data2 = rb_simulation(t_model, reduced_data);
208  %sim_data2 = rb_reconstruction(t_model, detailed_data, sim_data2);
209  save('compare_proj_Greedybase_30_t_part_3parts_coarse4_adapt_cutoff_PODinit.mat','detailed_data');
210 
211 
212  keyboard;
213 
214 
215  case 4
216  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
217  % Test: wieviele Basisvektoren werden benötigt um eine bestimme
218  % Genauigkeit zu erreichen?
219  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
220 
221  load('Greedybase_30_coarse4_222.mat')
222  Nmax=50;
223 
224  params.coarse_factor = 4;
225  model = advection_fv_output_model(params);
226  model.random_seed = 1234;
227 
228  model_data = gen_model_data(model);
229 
230  reduced_data= gen_reduced_data(model, detailed_data);
231 
232 
233  %loop over number of basis vectors
234  N=10:10:30;
235  av_error=zeros(1,size(N,2));
236  error_max=zeros(1,size(N,2));
237  rb_sim_av_time = zeros(1,size(N,2));
238  for i=1:size(N,2)
239  model.N=N(i);
240  %reduced_data.N=N(i);
241  disp(['Caculating error with ',num2str(N(i)),' basis vectors']);
242 
243  reduced_data_sub = model.reduced_data_subset(model, reduced_data);
244 
245  [av_error(i), error_max(i), rb_sim_av_time(i), std_dev]=calculate_test_estimator(model, reduced_data_sub);
246  end
247 
248 
249  %----------------T-part model simulations-----------------
250  model.rb_simulation = @lin_evol_rb_simulation_t_part;
251  % fix t-partitions
252  %t_part_range{1} = [0,floor(model.nt/2)];
253  %t_part_range{2} = [floor(model.nt/2),model.nt];
254  t_part_range{1} = [0,floor(model.nt/3)];
255  t_part_range{2} = [floor(model.nt/3),floor(2*model.nt/3)];
256  t_part_range{3} = [floor(2*model.nt/3),model.nt];
257  t_params.t_part_range= t_part_range;
258  %generate t-partition model
259  t_model = t_part_model(model,t_params);
260  model_data = gen_model_data(t_model);
261 
262  %generate detailed_data without error cut
263  t_model.t_part_rb_generation_cut_error = 0;
264  load('Greedybase_30_t_part_3parts_coarse4_222_cutoff.mat')
265  t_reduced_data = gen_reduced_data(t_model, detailed_data);
266 
267  av_error_t_part=zeros(1,size(N,2));
268  error_max_t_part= zeros(1,size(N,2));
269  rb_sim_av_time_t_part = zeros(1,size(N,2));
270 
271 
272  for i=1:size(N,2)
273  t_model.N=N(i);
274  %reduced_data.N=N(i);
275  disp(['Caculating error with ',num2str(N(i)),' basis vectors']);
276 
277  t_reduced_data_sub = t_model.reduced_data_subset(t_model, t_reduced_data);
278 
279  [av_error_t_part(i), error_max_t_part(i), rb_sim_av_time_t_part(i), std_dev]=calculate_test_estimator(t_model, t_reduced_data_sub);
280  end
281 
282 
283  figure
284  plot(N,av_error)
285  hold on
286  plot(N,error_max,'b-.')
287  plot(N,av_error_t_part,'r');
288  plot(N,error_max_t_part,'r-.')
289  title('Errors over number of basisvectors')
290  legend('average error','maximal error','t-part av error','t_part max error')
291 
292  keyboard;
293 
294  case 5
295  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
296  % einfach mal eine simulation anschauen und vergleichen
297  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
298 
299 
300  %ursprüngliches modell laden und generieren
301  load('Greedybase_1p_coarse4_adapt_0.001.mat')
302 
303  params.coarse_factor = 4;
304  model = advection_fv_output_model(params);
305  model_data = gen_model_data(model);
306  model = set_mu(model, [0.712]);
307  reduced_data = gen_reduced_data(model, detailed_data);
308 
309  %real detailed_sim
310  sim_data_det = detailed_simulation(model, model_data);
311 
312  disp(['Simulationen für ',num2str(get_mu(model)')]);
313  sim_data = rb_simulation(model, reduced_data);
314  sim_data = rb_reconstruction(model, detailed_data, sim_data);
315 
316 
317  %t-model laden und simulieren
318  load('Greedybase_adaptive_1p_t_part_coarse4_adapt_cutoff_eps0.001_N50.mat');
319 
320  model.rb_simulation = @lin_evol_rb_simulation_t_part;
321  % fix t-partitions
322  if isfield(detailed_data,'t_part_range')
323  t_part_range = detailed_data.t_part_range;
324  else
325  %t_part_range{1} = [0,floor(model.nt/2)];
326  %t_part_range{2} = [floor(model.nt/2),model.nt];
327  t_part_range{1} = [0,floor(model.nt/3)];
328  t_part_range{2} = [floor(model.nt/3),floor(2*model.nt/3)];
329  t_part_range{3} = [floor(2*model.nt/3),model.nt];
330  end
331  t_params.t_part_range= t_part_range;
332  t_params.transition_model = 0;
333  %generate t-partition model
334  t_model = t_part_model(model,t_params);
335 
336 
337  t_reduced_data = gen_reduced_data(t_model, detailed_data);
338  disp(['Simulationen für ',num2str(get_mu(t_model)')]);
339  sim_data_t = rb_simulation(t_model, t_reduced_data);
340  sim_data_t = rb_reconstruction(t_model, detailed_data,sim_data_t);
341 
342  %calculate real error
343  err = fv_l2_error(sim_data.U(:,model.save_time_indices+1), sim_data_det.U, model_data.grid,[]);
344  err_t = fv_l2_error(sim_data_t.U(:,model.save_time_indices+1), sim_data_det.U, model_data.grid,[]);
345 
346  %plot
347  figure
348  plot(sim_data.Delta)
349  hold on
350  plot(sim_data_t.Delta,'r')
351  plot(err,'b-.');
352  plot(err_t, 'r-.');
353  title('Fehlerplot')
354  legend('Error estimator normal model','Error estimator t-part model', 'real error normal model', 'real error t-part model')
355 
356  plot_params.plot_function = @plot_element_data;
357  plot_data_sequence(model, model_data.grid,sim_data.U,plot_params)
358  title('original solution')
359 
360  plot_data_sequence(t_model, model_data.grid,sim_data_t.U,plot_params)
361  title('t-part_model')
362 
363 
364 
365 
366  keyboard;
367 
368 
369  case 6
370  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
371  % adaptive basis generation
372  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
373 
374  if nargin<2||isempty(eps)
375  eps = 1e-2;
376  end
377 
378  if nargin<3||isempty(Nmax)
379  Nmax = 35;
380  end
381 
382  params.coarse_factor = 4;
383  model = advection_fv_output_model(params);
384 
385  model.random_seed = 1234;
386  model.RB_numintervals = [2];
387  model.RB_stop_epsilon = eps;%1e-1;
388  model.RB_stop_Nmax = Nmax;
389  model.RB_refinement_mode = 'adaptive';
390  model.RB_error_indicator = 'estimator';%'error';%'estimator';
391  model.RB_detailed_train_savepath = fullfile(rbmatlabtemp,'dump');
392  model.RB_extension_algorithm = @RB_extension_PCA_fixspace;
393  model.RB_stop_epsilon = eps;%1e-1;
394  model.RB_stop_timeout = 24*60*60; %
395  model.RB_stop_max_val_train_ratio = 1;
396  model.RB_refinement_theta = 0.2;
397  model.RB_val_rand_seed = 543;
398  model.RB_M_val_size = 10;
399  model.RB_stop_max_refinement_level = 100;
400  model.RB_element_indicator_mode='nodes_cogs_skippedrefs';
401  model.RB_element_indicator_s_max = 1;
402 
403  model_data = gen_model_data(model);
404 
405 
406  model.rb_simulation = @lin_evol_rb_simulation_t_part;
407  % fix t-partitions
408  t_params=[];
409  %t_part_range{1} = [0,floor(model.nt/2)];
410  %t_part_range{2} = [floor(model.nt/2),model.nt];
411  %t_params.t_part_range= t_part_range;
412  %generate t-partition model
413  t_model = t_part_model(model,t_params);
414 
415 
416  t_model.t_part_rb_generation_mode ='adaptive_error_restriction';
417  t_model.t_part_rb_generation_cut_error = 1;
418  t_model.max_t_part_refinement_level = 5; %fix maximum refinement level
419  %t_model.RB_stop_max_val_train_ratio = 1;
420  %middle point bisection!!!!!!!!!!!!!!!!!
421  %t_model.t_part_middle_point_bisection=1;
422  %!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
423 
424  tstart = tic;
425  detailed_data = gen_detailed_data(t_model, model_data);
426  detailed_data.total_elapsed_time = toc(tstart);
427 
428  %simulate
429  %reduced_data = gen_reduced_data(t_model, detailed_data);
430  %sim_data2 = rb_simulation(t_model, reduced_data);
431  %sim_data2 = rb_reconstruction(t_model, detailed_data, sim_data2);
432  %save(['Greedybase_adaptive_1p_t_part_coarse4_adapt_cutoff_eps',num2str(eps),'_N',num2str(Nmax),'.mat'],'detailed_data');
433  save(fullfile(rbmatlabresult,['Greedybase_adaptive_1p_t_part_coarse',num2str(params.coarse_factor), '_adapt_cutoff_eps',num2str(eps),'_N',num2str(Nmax),'.mat']),'detailed_data','t_model');
434  %load(fullfile(rbmatlabresult,['Greedybase_adaptive_1p_t_part_coarse4_adapt_cutoff_eps',num2str(eps),'_N',num2str(Nmax),'.mat']));%,'detailed_data')
435 
436  %keyboard;
437 
438  t_reduced_data = gen_reduced_data(t_model, detailed_data);
439  t_model.t_part_range = detailed_data.t_part_range;
440 
441  [max_sim_data_t, rtimes] = get_max_sim_data(t_model, t_reduced_data);
442 
443  max_sim_data_t.Delta = max_sim_data_t.t_part_sim_data{1}.Delta;
444  for p=2:length(max_sim_data_t.t_part_sim_data)
445  max_sim_data_t.Delta = [max_sim_data_t.Delta, max_sim_data_t.t_part_sim_data{p}.Delta(2:end)];
446  end
447 
448  offline_time = detailed_data.total_elapsed_time;
449 
450  columntitles = {'minslicesize', 'offline_time', 'noofslices', 'averageN', 'minN' ,'maxN', 'averagedtime'};
451  values = zeros(7,1);
452  values(1) = min(cellfun(@(X) X(2) - X(1), detailed_data.t_part_range, 'UniformOutput', true));
453  values(2) = offline_time;
454  values(3) = length(detailed_data.t_part_range);
455  Nsizes = cellfun(@(X) size(X.RB,2), detailed_data.t_part_detailed_data);
456  values(4) = sum(Nsizes) / values(3);
457  values(5) = min(Nsizes);
458  values(6) = max(Nsizes);
459  values(7) = mean(rtimes);
460  save(fullfile(rbmatlabresult, 'data'), 'max_sim_data_t', 'values', 'columntitles', 'rtimes');
461  print_datatable(fullfile(rbmatlabresult,'tableline'),columntitles,values);
462 
463  case 7
464  %%%%%%%%%%%%%%%%%
465  %generate normal reduced basis
466  %%%%%%%%%%%%%%%%%%%%%%%%%%5
467 
468 
469  % fix base model
470  params.coarse_factor = 4;
471  model = advection_fv_output_model(params);
472  model.rb_simulation = @lin_evol_rb_simulation_t_part;
473 
474  model.random_seed=1234;
475  model.RB_numintervals = [2];
476  model.RB_generation_mode = 'greedy_refined';%'greedy_uniform_fixed'; %
477  model.RB_refinement_mode = 'adaptive';
478  model.RB_error_indicator = 'estimator';%'error';%'estimator';
479  model.RB_detailed_train_savepath = fullfile(RBMATLABTEMP,'dump');
480  model.RB_extension_algorithm = @RB_extension_PCA_fixspace;
481  model.RB_stop_epsilon = eps;%1e-1;
482  model.RB_stop_timeout = 24*60*60; %
483  model.RB_stop_Nmax = 500;
484  model.RB_stop_max_val_train_ratio = 1;
485  model.RB_refinement_theta = 0.2;
486  model.RB_val_rand_seed = 543;
487  model.RB_M_val_size = 10;
488  model.RB_stop_max_refinement_level = 5;
489  model.RB_element_indicator_mode='nodes_cogs_skippedrefs';
490  model.RB_element_indicator_s_max = 1;
491 
492 
493  %rb_simulation without t-partition basis
494  model_data = gen_model_data(model);
495  detailed_data = gen_detailed_data(model, model_data);
496  %reduced_data = gen_reduced_data(model, detailed_data);
497  %sim_data = rb_simulation(model, reduced_data);
498  %sim_data_original = rb_reconstruction(model, detailed_data, sim_data);
499 
500  save(['Greedybase_1p_coarse4_adapt_',num2str(eps),'.mat'],'detailed_data');
501 
502  %keyboard
503 
504 
505  case 8
506  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
507  % generate the plot average sim time vs. epsilon_tol
508  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
509 
510  %Part 1:loading and caclculating erro ans sim time data
511 
512  %The bases should be saved as (example) :Greedybase_1p_coarse4_adapt_0.05.mat
513  %and:
514  %Greedybase_adaptive_1p_t_part_coarse4_adapt_cutoff_eps0.06_N30.mat
515 
516  %be sure to be in the right directory
517 
518  params.coarse_factor = 4;
519  model = advection_fv_output_model(params);
520  model.rb_simulation = @lin_evol_rb_simulation_t_part;
521 
522  model_data = gen_model_data(model);
523 
524  s=what;
525  s=s.mat;
526 
527  eps_tol = [];
528  av_error = [];
529  max_error = [];
530  av_sim_time = [];
531  N_normal=[];
532 
533  s_norm_base = 'Greedybase_1p_coarse4_adapt_';
534  %go through all normal bases and calculate average error/max
535  %error/av simtime
536  for i=1:length(s)
537  if strncmp(s{i},s_norm_base,length(s_norm_base))
538  v=strfind(s{i},'.');
539  indx=v(end);
540  eps_tol = [eps_tol,str2double(s{i}(length(s_norm_base)+1:indx-1))];
541  load(s{i});
542  N_normal = [N_normal,detailed_data.N];
543  reduced_data = gen_reduced_data(model, detailed_data);
544  %start testing
545  [av_error_new, max_error_new, av_sim_time_new, std_dev]=calculate_test_estimator(model, reduced_data);
546  %append
547  av_error = [av_error, av_error_new];
548  max_error = [max_error, max_error_new];
549  av_sim_time = [av_sim_time, av_sim_time_new];
550  end
551  end
552 
553  av_error_normal = av_error;
554  eps_tol_normal = eps_tol;
555  max_error_normal = max_error;
556  av_sim_time_normal = av_sim_time;
557 
558 
559  %go through all t-part bases
560  eps_tol = [];
561  av_error = [];
562  max_error = [];
563  av_sim_time = [];
564  Nmax = [];
565  av_N = [];
566 
567  s_norm_base = 'Greedybase_adaptive_middle_1p_t_part_coarse4_adapt_cutoff_eps';
568  %s_norm_base = 'Greedybase_adaptive_1p_t_part_coarse4_adapt_cutoff_eps';
569  %go through all normal bases and calculate average error/max
570  %error/av simtime
571  for i=1:length(s);
572  if strncmp(s{i},s_norm_base,length(s_norm_base))
573  v=strfind(s{i},'_');
574  indx=v(10);%indx=v(9);%not nice but works here!! Change when needed!!
575  v=strfind(s{i},'.');
576  indx_dot = v(end);
577  eps_tol = [eps_tol,str2double(s{i}(length(s_norm_base)+1:indx-1))];
578  Nmax= [Nmax,str2double(s{i}(indx+2:indx_dot))];
579  load(s{i});
580  %calculate average number of basis vectors
581  xsum=0;
582  for j=1:length(detailed_data.t_part_detailed_data)
583  xsum=xsum+detailed_data.t_part_detailed_data{j}.N;
584  end
585  av_N = [av_N,xsum/length(detailed_data.t_part_detailed_data)];
586 
587  t_params.t_part_range = detailed_data.t_part_range;
588  t_model = t_part_model(model, t_params);
589  reduced_data = gen_reduced_data(t_model, detailed_data);
590  %start testing
591  [av_error_new, max_error_new, av_sim_time_new, std_dev]=calculate_test_estimator(t_model, reduced_data);
592  %append
593  av_error = [av_error, av_error_new];
594  max_error = [max_error, max_error_new];
595  av_sim_time = [av_sim_time, av_sim_time_new];
596  end
597  end
598 
599  av_error_t = av_error;
600  eps_tol_t = eps_tol;
601  max_error_t = max_error;
602  av_sim_time_t = av_sim_time;
603  N_max_t = Nmax;
604  av_N_t = av_N;
605 
606  save('error_sim_time_data.mat','av_error_t','av_error_normal','eps_tol_t','eps_tol_normal','max_error_normal','max_error_t','N_max_t','av_sim_time_t','av_sim_time_normal');
607 
608 
609  %figure
610  %semilogy(max_error_normal,av_sim_time_normal,'x-')
611  %hold on
612  %semilogy(max_error_t,av_sim_time_t,'rx-');
613  %xlabel('max error')
614  %ylabel('av_sim_time')
615  %legend('normal','t-part')
616 
617  %sort for epsion_tol
618  [eps_tol_t,av_error_t,max_error_t,av_sim_time_t,N_max_t] = sort_low_to_high(eps_tol_t, ...
619  av_error_t,max_error_t,av_sim_time_t, N_max_t);
620  [eps_tol_normal,av_error_normal,max_error_normal,av_sim_time_normal] = sort_low_to_high(eps_tol_normal, ...
621  av_error_normal,max_error_normal,av_sim_time_normal, []);
622 
623 
624  %sort out for N_max settings
625  N_look_up=[];
626  T_data={};
627  for i = 1:length(N_max_t)
628  %is this N_max already in N_look_up?
629  is_element=0;
630  for j=1:length(N_look_up)
631  if N_look_up(j) == N_max_t(i)
632  is_element=1;
633  pos=j;
634  end
635  end
636  if ~is_element
637  pos=length(N_look_up);
638  N_look_up(pos+1) = N_max_t(i);
639  T_data{pos+1}.av_sim_time_t=av_sim_time_t(i);
640  T_data{pos+1}.max_error_t = max_error_t(i);
641  else
642  T_data{pos}.av_sim_time_t(end+1) = av_sim_time_t(i);
643  T_data{pos}.max_error_t(end+1) = max_error_t(i);
644  end
645  end
646 
647  if length(N_look_up) ==1
648  figure
649  semilogy(max_error_normal,av_sim_time_normal,'bx-','LineWidth',2);%,'Color',[1,0.6,0])
650  hold on
651  semilogy(max_error_t,av_sim_time_t,'rx-','LineWidth',2);
652  title('Approximation error vs. average simulation time')
653  xlabel('maximum estimated error')
654  ylabel('average simulation time in sec.')
655  legend('no T-partition','adaptive T-partition Nmax40')
656  set(gca,'ytick',[0.5,0.55,0.575,0.6,0.65,0.7,0.75])
657  else
658  colorstyles={'rx-','gx-','mx-','yx-'};
659  legend_text={'normal reduced model'};
660  figure
661  semilogy(max_error_normal,av_sim_time_normal,'bx-','LineWidth',2);%,'Color',[1,0.6,0])
662  hold on
663  title('Approximation error vs. average simulation time')
664  xlabel('max error')
665  ylabel('average simulation time')
666  for i=1:length(N_look_up)
667  semilogy(T_data{i}.max_error_t,T_data{i}.av_sim_time_t,colorstyles{i},'LineWidth',2);
668  legend_text{i+1}=num2str(N_look_up(i));
669  end
670  legend(legend_text);
671  set(gca,'ytick',[0.5,0.55,0.6,0.65,0.7,0.75])
672  end
673  keyboard;
674 
675 
676 
677  case 9
678  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
679  % fixed partition maximal error plot
680  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
681 
682  %ursprüngliches modell laden und generieren
683  load('Greedybase_1p_50_coarse4_adapt.mat')
684 
685  params.coarse_factor = 4;
686  model = advection_fv_output_model(params);
687  model.rb_simulation = @lin_evol_rb_simulation_t_part;
688  model_data = gen_model_data(model);
689 
690  reduced_data = gen_reduced_data(model, detailed_data);
691 
692  max_sim_data_normal = get_max_sim_data(model, reduced_data);
693 
694  %t-model laden und simulieren
695  load('Greedybase_1p_50_t_part_3parts_coarse4_adapt_cutoff.mat');
696 
697  model.rb_simulation = @lin_evol_rb_simulation_t_part;
698  % fix t-partitions
699  if isfield(detailed_data,'t_part_range')
700  t_part_range = detailed_data.t_part_range;
701  else
702  %t_part_range{1} = [0,floor(model.nt/2)];
703  %t_part_range{2} = [floor(model.nt/2),model.nt];
704  t_part_range{1} = [0,floor(model.nt/3)];
705  t_part_range{2} = [floor(model.nt/3),floor(2*model.nt/3)];
706  t_part_range{3} = [floor(2*model.nt/3),model.nt];
707  end
708  t_params.t_part_range= t_part_range;
709  t_params.transition_model = 0;
710  %generate t-partition model
711  t_model = t_part_model(model,t_params);
712 
713 
714  t_reduced_data = gen_reduced_data(t_model, detailed_data);
715 
716  max_sim_data_t = get_max_sim_data(t_model, t_reduced_data);
717 
718  max_sim_data_t.Delta = max_sim_data_t.t_part_sim_data{1}.Delta;
719  for p=2:length(max_sim_data_t.t_part_sim_data)
720  max_sim_data_t.Delta = [max_sim_data_t.Delta, max_sim_data_t.t_part_sim_data{p}.Delta(2:end)];
721  end
722 
723  figure
724  plot(max_sim_data_normal.Delta,'LineWidth',2)
725  hold on
726  plot(max_sim_data_t.Delta,'r','LineWidth',2)
727  xlabel('time step')
728  ylabel('estimated maximal error')
729  title('Evolution of maximal error of a validation set')
730  legend('normal','t-part')
731 
732 
733  %plot partition borders
734  limy = ylim;
735  for p=1:length(t_model.t_part_range)
736  x_pos = t_model.t_part_range{p}(2);
737 
738  plot([x_pos,x_pos],limy,'-.','Color',[1,0.6,0]);
739  end
740 
741  legend('normal','t-part','interval-border')
742 
743  keyboard;
744 
745  case 10
746 
747  %%%%%%%%%%%%%%%%%
748  % generate reduced bases for time intervals of length 'granularity'
749  %%%%%%%%%%%%%%%%%%%%%%%%%%5
750 
751  if ~exist('eps', 'var') || isempty(eps)
752  eps = 1e-2;
753  end
754 
755  if ~exist('Nmax', 'var') || isempty(Nmax)
756  Nmax = 30;
757  end
758 
759 
760  params.coarse_factor = 8;
761  model = advection_fv_output_model(params);
762  model.rb_simulation = @lin_evol_rb_simulation_t_part;
763 
764  model.random_seed = 1234;
765  model.RB_numintervals = [2];
766  model.RB_generation_mode = 'greedy_refined';%'greedy_uniform_fixed'; %
767  model.RB_refinement_mode = 'adaptive';
768  model.RB_error_indicator = 'estimator';%'error';%'estimator';
769  model.RB_detailed_train_savepath = fullfile(RBMATLABTEMP,'dump');
770  model.RB_extension_algorithm = @RB_extension_PCA_fixspace;
771  model.RB_stop_epsilon = eps;%1e-1;
772  model.RB_stop_timeout = 24*60*60; % 24 stunden
773  model.RB_stop_Nmax = Nmax;
774  model.RB_stop_max_val_train_ratio = 1;
775  model.RB_refinement_theta = 0.2;
776  model.RB_val_rand_seed = 543;
777  model.RB_M_val_size = 10;
778  model.RB_stop_max_refinement_level = 5;
779  model.RB_element_indicator_mode = 'nodes_cogs_skippedrefs';
780  model.RB_element_indicator_s_max = 1;
781  model.rb_simulation = @lin_evol_rb_simulation_t_part;
782 
783  model_data = gen_model_data(model);
784 
785 
786  % fix t-partitions
787  t_params=[];
788  t_model = t_part_model(model,t_params);
789 
790  t_model.t_part_rb_generation_mode ='fixed';
791 
792  granularity = 2;
793 
794  t_model.t_part_range = cell(1,ceil(model.nt/granularity));
795  for i=1:length(t_model.t_part_range)-1
796  t_model.t_part_range{i} = [granularity*(i-1),granularity*i];
797  end
798  t_model.t_part_range{end} = [granularity*i,model.nt];
799 
800 
801 % t_part_range{1} = [0,4];%floor(model.nt/2)];
802 % t_part_range{2} = [4,8];%floor(model.nt/2),model.nt];
803 % t_model.t_part_range = t_part_range;
804 
805  t_model.t_part_rb_generation_cut_error = 1;
806 
807  tstart = tic;
808  detailed_data = gen_detailed_data(t_model, model_data);
809  t=toc(tstart);
810  detailed_data.total_elapsed_time = t;
811  infix = ['eps', num2str(eps) , '_Nmax', num2str(Nmax), '_coarse', num2str(params.coarse_factor), '_gran', num2str(granularity)];
812  save(fullfile(rbmatlabresult, ['Greedybase_',infix]),'detailed_data');
813 
814  %rb_simulation without t-partition basis
815 % model_data = gen_model_data(model);
816 % detailed_data = gen_detailed_data(model, model_data);
817  %reduced_data = gen_reduced_data(model, detailed_data);
818  %sim_data = rb_simulation(model, reduced_data);
819  %sim_data_original = rb_reconstruction(model, detailed_data, sim_data);
820 
821  t_reduced_data = gen_reduced_data(t_model, detailed_data);
822 
823  [max_sim_data_t, rtimes] = get_max_sim_data(t_model, t_reduced_data);
824 
825  max_sim_data_t.Delta = max_sim_data_t.t_part_sim_data{1}.Delta;
826  for p=2:length(max_sim_data_t.t_part_sim_data)
827  max_sim_data_t.Delta = [max_sim_data_t.Delta, max_sim_data_t.t_part_sim_data{p}.Delta(2:end)];
828  end
829 
830  offline_time = detailed_data.total_elapsed_time;
831 
832  columntitles = {'minslicesize', 'offline_time', 'noofslices', 'averageN', 'minN' ,'maxN', 'averagedtime'};
833  values = zeros(7,1);
834  values(1) = granularity;
835  values(2) = offline_time;
836  values(3) = length(t_model.t_part_range);
837  Nsizes = cellfun(@(X) size(X.RB,2), detailed_data.t_part_detailed_data);
838  values(4) = sum(Nsizes) / values(3);
839  values(5) = min(Nsizes);
840  values(6) = max(Nsizes);
841  values(7) = mean(rtimes);
842  save(fullfile(rbmatlabresult, ['data',infix]), 'max_sim_data_t', 'values', 'columntitles', 'rtimes');
843  print_datatable(fullfile(rbmatlabresult,['tableline',infix]),columntitles,values);
844 
845 
846  case 11
847  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
848  % adaptive basis generation with middelpoint bisection
849  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
850 
851  if nargin<2||isempty(eps)
852  eps = 1e-2;
853  end
854 
855  if nargin<3||isempty(Nmax)
856  Nmax = 35;
857  end
858 
859  params.coarse_factor = 4;
860  model = advection_fv_output_model(params);
861 
862  model.random_seed = 1234;
863  model.RB_numintervals = [2];
864  model.RB_stop_epsilon = eps;%1e-1;
865  model.RB_stop_Nmax = Nmax;
866  model.RB_refinement_mode = 'adaptive';
867  model.RB_error_indicator = 'estimator';%'error';%'estimator';
868  model.RB_detailed_train_savepath = fullfile(rbmatlabtemp,'dump');
869  model.RB_extension_algorithm = @RB_extension_PCA_fixspace;
870  model.RB_stop_epsilon = eps;%1e-1;
871  model.RB_stop_timeout = 24*60*60; %
872  model.RB_stop_max_val_train_ratio = 1;
873  model.RB_refinement_theta = 0.2;
874  model.RB_val_rand_seed = 543;
875  model.RB_M_val_size = 10;
876  model.RB_max_refinement_level = 50;
877  model.RB_element_indicator_mode='nodes_cogs_skippedrefs';
878  model.RB_element_indicator_s_max = 1;
879 
880  model_data = gen_model_data(model);
881 
882 
883  model.rb_simulation = @lin_evol_rb_simulation_t_part;
884  % fix t-partitions
885  t_params=[];
886  %t_part_range{1} = [0,floor(model.nt/2)];
887  %t_part_range{2} = [floor(model.nt/2),model.nt];
888  %t_params.t_part_range= t_part_range;
889  %generate t-partition model
890  t_model = t_part_model(model,t_params);
891 
892 
893  t_model.t_part_rb_generation_mode ='adaptive_error_restriction';
894  t_model.t_part_rb_generation_cut_error = 1;
895  t_model.max_t_part_refinement_level = 5; %fix maximum refinement level
896  %t_model.RB_stop_max_val_train_ratio = 1;
897  %middle point bisection!!!!!!!!!!!!!!!!!
898  t_model.t_part_middle_point_bisection=1;
899  %!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
900 
901  tstart = tic;
902  detailed_data = gen_detailed_data(t_model, model_data);
903  detailed_data.total_elapsed_time = toc(tstart);
904 
905  %simulate
906  %reduced_data = gen_reduced_data(t_model, detailed_data);
907  %sim_data2 = rb_simulation(t_model, reduced_data);
908  %sim_data2 = rb_reconstruction(t_model, detailed_data, sim_data2);
909  %save(['Greedybase_adaptive_1p_t_part_coarse4_adapt_cutoff_eps',num2str(eps),'_N',num2str(Nmax),'.mat'],'detailed_data');
910  save(fullfile(rbmatlabresult,['Greedybase_adaptive_middle_1p_t_part_coarse',num2str(params.coarse_factor), '_adapt_cutoff_eps',num2str(eps),'_N',num2str(Nmax),'.mat']),'detailed_data','t_model');
911  %load(fullfile(rbmatlabresult,['Greedybase_adaptive_1p_t_part_coarse4_adapt_cutoff_eps',num2str(eps),'_N',num2str(Nmax),'.mat']));%,'detailed_data')
912 
913  %keyboard;
914 
915  t_reduced_data = gen_reduced_data(t_model, detailed_data);
916  t_model.t_part_range = detailed_data.t_part_range;
917 
918  [max_sim_data_t, rtimes] = get_max_sim_data(t_model, t_reduced_data);
919 
920  max_sim_data_t.Delta = max_sim_data_t.t_part_sim_data{1}.Delta;
921  for p=2:length(max_sim_data_t.t_part_sim_data)
922  max_sim_data_t.Delta = [max_sim_data_t.Delta, max_sim_data_t.t_part_sim_data{p}.Delta(2:end)];
923  end
924 
925  offline_time = detailed_data.total_elapsed_time;
926 
927  columntitles = {'minslicesize', 'offline_time', 'noofslices', 'averageN', 'minN' ,'maxN', 'averagedtime'};
928  values = zeros(7,1);
929  values(1) = min(cellfun(@(X) X(2) - X(1), detailed_data.t_part_range, 'UniformOutput', true));
930  values(2) = offline_time;
931  values(3) = length(detailed_data.t_part_range);
932  Nsizes = cellfun(@(X) size(X.RB,2), detailed_data.t_part_detailed_data);
933  values(4) = sum(Nsizes) / values(3);
934  values(5) = min(Nsizes);
935  values(6) = max(Nsizes);
936  values(7) = mean(rtimes);
937  save(fullfile(rbmatlabresult, 'data'), 'max_sim_data_t', 'values', 'columntitles', 'rtimes');
938  print_datatable(fullfile(rbmatlabresult,'tableline'),columntitles,values);
939 
940 
941  case 12
942  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
943  % generate entire data and figure for average_sim_time vs.
944  % max_error
945  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
946 
947  %generate normal bases
948  t_part_script(7,0.1,40);
949  t_part_script(7,0.05,40);
950  t_part_script(7,0.03,40);
951  t_part_script(7,0.01,40);
952  t_part_script(7,0.007,40);
953  t_part_script(7,0.005,40);
954  t_part_script(7,0.003,40);
955  t_part_script(7,0.001,40);
956  t_part_script(11,0.1,40);
957  t_part_script(11,0.08,40);
958  t_part_script(11,0.05,40);
959  t_part_script(11,0.01,40);
960  t_part_script(11,0.008,40);
961  t_part_script(11,0.005,40);
962  t_part_script(11,0.003,40);
963 
964  t_part_script(8);
965 
966 
967 end
968 
969 
970 
971 end
972 
973 function [max_sim_data, rtimes] = get_max_sim_data(model, reduced_data)
974  rand('state',model.random_seed);
975  M_val=rand_uniform(20, model.mu_ranges);
976 
977  model = set_mu(model,M_val(:,1));
978 
979  rtimes = zeros(1,100);
980 
981 
982  if isfield(model,'t_part_range')
983  %tpart model...
984 
985  max_sim_data = rb_simulation(model, reduced_data);
986 
987  if ~isfield(max_sim_data, 't_part_sim_data')
988  max_sim_data.t_part_sim_data{1} = max_sim_data;
989  end
990 
991 
992  for i=1:size(M_val,2)
993  model = set_mu(model, M_val(:,i));
994  tic
995  sim_data = rb_simulation(model, reduced_data);
996  rtimes(i) = toc;
997 
998 
999  if ~isfield(sim_data, 't_part_sim_data')
1000  sim_data.t_part_sim_data{1} = sim_data;
1001  end
1002 
1003 
1004  for p=1:length(sim_data.t_part_sim_data)
1005  for k=1:length(sim_data.t_part_sim_data{p}.Delta)
1006  if (max_sim_data.t_part_sim_data{p}.Delta(k)<sim_data.t_part_sim_data{p}.Delta(k))
1007  max_sim_data.t_part_sim_data{p}.Delta(k) = sim_data.t_part_sim_data{p}.Delta(k);
1008  end
1009  end
1010  end
1011 
1012  end
1013 
1014  if isfield(max_sim_data,'Delta')
1015  max_sim_data = rmfield(max_sim_data,'Delta');
1016  end
1017 
1018  %test_error = rb_test_estimator(model, reduced_data,M_val);
1019  %[error_max,mu_ind] = max(test_error);
1020 
1021  %model = set_mu(model,M_val(:,mu_ind));
1022  %sim_data = rb_simulation(model, reduced_data);
1023  else
1024  max_sim_data = rb_simulation(model, reduced_data);
1025 
1026  for i=1:size(M_val,2)
1027  model = set_mu(model, M_val(:,i));
1028  tic;
1029  sim_data = rb_simulation(model, reduced_data);
1030  rtimes(i) = toc;
1031 
1032  for k=1:length(sim_data.Delta)
1033  if(max_sim_data.Delta(k)<sim_data.Delta(k))
1034  max_sim_data.Delta(k) = sim_data.Delta(k);
1035  end
1036  end
1037  end
1038 
1039  end
1040 
1041 end
1042 
1043 
1044 function draw_t_part_lines(model)
1045 
1046  limy = ylim;
1047 
1048  x_pos = [];
1049  for p=1:length(model.t_part_range)
1050  x_pos = model.t_part_range{p}(2);
1051 
1052  plot([x_pos,x_pos],limy,'r-.');
1053  end
1054 
1055 
1056 
1057 
1058 
1059 end
1060 
1061 
1062 
1063 function [s1, s2, s3 ,s4, s5] = sort_low_to_high(v1, v2, v3, v4, v5)
1064 %function sorts the entries of v1 from low to high and changes the order of
1065 %the other vectors in the same manner as v1.
1066 %!!!!!!!!!!!works only for vectors with non negative entries
1067 
1068 s1=zeros(size(v1,1),size(v1,2));
1069 s2=zeros(size(v2,1),size(v2,2));
1070 s3=zeros(size(v3,1),size(v3,2));
1071 s4=zeros(size(v4,1),size(v4,2));
1072 s5=zeros(size(v5,1),size(v5,2));
1073 
1074 
1075 for i = 1:size(s1,2)
1076  %find starting value
1077  continue_loop=1;
1078  indx=0;
1079  while (continue_loop)&&(indx<=size(v1,2))
1080  indx=indx+1;
1081  if v1(indx)~=-1
1082  continue_loop=0;
1083  end
1084  end
1085  minim=v1(indx);
1086  %find miminum
1087  for j=1:size(v1,2)
1088  if (v1(j)~=-1)&&(v1(j)<minim)
1089  minim = v1(j);
1090  indx=j;
1091  end
1092  end
1093  %fill new vectors
1094  s1(i)=v1(indx);
1095  if ~isempty(v2)
1096  s2(i)=v2(indx);
1097  end
1098  v1(indx)=-1;
1099  if ~isempty(v3)
1100  s3(i)=v3(indx);
1101  end
1102  if ~isempty(v4)
1103  s4(i)=v4(indx);
1104  end
1105  if ~isempty(v5)
1106  s5(i)=v5(indx);
1107  end
1108 
1109 end
1110 
1111 end
1112 
1113 
Test
Test reduced basis implementation.
Definition: DetailedData.m:1
fv_l2_error
function l2_error = fv_l2_error(U1, U2, W)
function computing the l2-error between the two fv-functions or function sequences in U1...
Definition: fv_l2_error.m:17
Generated on Fri Sep 30 2016 10:51:25 for rbmatlab by   doxygen 1.8.8