rbmatlab  1.16.09
 All Classes Namespaces Files Functions Variables Modules Pages
t_part_gen_detailed_data.m
1 function detailed_data = t_part_gen_detailed_data(model, model_data)
2 %function detailed_data = t_part_gen_detailed_data(model, model_data)
3 %
4 % function generating the detailed data (the RB) for a t-partition model
5 %
6 % required fields of model:
7 % basis_vector_overlap: number of basis vectors from the foregoing t-part RB
8 % used as initial basis in the next one.
9 % RB_stop_epsilon_global: This is equal to the base model's 'RB_stop_epsilon'
10 % field, which defines the maximum error estimator
11 % tolerance for the generated basis.
12 %
13 % there are different generation modes: for fixed t-partition and adaptive.
14 %
15 
16 % Markus Dihlmann 15.02.2011
17 
18 model.t_part_in_basis_generation =1;
19 model.reduced_data_subset = @(model, reduced_data)reduced_data; %temp replace to circumvent problems in rb_test_indicator
20 
21 detailed_data = [];
22 
23 detailed_data = structcpy(detailed_data,model_data);
24 
25 detailed_data.t_part_range = model.t_part_range;
26 
27 base_model = model.base_model;
28 
29 if ~isfield(model,'t_part_rb_generation_mode')
30  error('For t-partition please set the field t_part_rb_generation_mode in t_part_model');
31 end
32 
33 switch model.t_part_rb_generation_mode
34  case 'fixed'
35  nr_t_parts = length(model.t_part_range);
36  %t_part_detailed_data = cell(1,nr_t_parts);
37 
38  for ind = 1:nr_t_parts
39  filecache_clear;
40  model.t_part_for_simulation = ind;
41  if ind>1
42  %overlap from foregoing t-part RB as initial basis in next
43  %t-part RB
44  model.init_values_algorithm = @(model, detailed_data)t_part_detailed_data{ind-1}.RB(:,1:model.basis_vector_overlap);
45 
46  %try with POD of some final time step solutions of the
47  %foregoing partition
48  %model.init_values_algorithm = @t_part_initial_cond_POD;
49  %model.rb_init_data_basis = @t_part_initial_cond_POD;
50 
51  end
52  model.RB_stop_epsilon = model.RB_stop_epsilon_global/model.nt*(model.t_part_range{ind}(2)-model.t_part_range{ind}(1));
53  t_part_detailed_data{ind} = rb_basis_generation(model, detailed_data);
54  % prevent recursive detailed data structures
55  t_part_detailed_data{ind} = rmfield(t_part_detailed_data{ind}, 't_part_detailed_data');
56  detailed_data.t_part_detailed_data = t_part_detailed_data;
57  end
58 
59  case 'adaptive'
60  %make a refinement of the ranges which have the biggest error gain
61  %(model.t_part_theta controls the partition of total refinement)
62 
63  case 'adaptive_error_restriction'
64  %make refinement so that the error gain is equally repartitioned
65  nr_t_parts = length(model.t_part_range);
66 
67  if ~isfield(model,'max_t_part_refinement_level')
68  model.max_t_part_refinement_level = 4;
69  end
70  loop_nr =1;
71  continue_loop = 1;
72 
73  t_parts_for_basis_generation = 1:nr_t_parts;
74 
75  while continue_loop
76 
77 
78  refine_p=[]; % array of partition indices marked for refinement
79  nr_t_parts = length(t_parts_for_basis_generation);
80  if isfield(detailed_data,'t_part_detailed_data')
81  t_part_detailed_data = detailed_data.t_part_detailed_data;
82  end
83 
84  for ind = 1:nr_t_parts
85  filecache_clear;
86  t_part_ind = t_parts_for_basis_generation(ind);
87  model.t_part_for_simulation = t_part_ind;
88  model.t_part_id_in_basisgen = t_part_ind;
89  if t_part_ind>1
90  %overlap from foregoing t-part RB as initial basis in next
91  %t-part RB
92  %
93  %model.init_values_algorithm = @(model, detailed_data)...
94  % t_part_detailed_data{t_part_ind-1}.RB(:,1:model.basis_vector_overlap);
95  %
96  %model.init_values_algorithm = @t_part_initial_cond_POD;
97  model.rb_init_data_basis = @t_part_initial_cond_POD;
98  end
99  model.RB_stop_epsilon = model.RB_stop_epsilon_global ...
100  /model.nt...
101  *(model.t_part_range{t_part_ind}(2)-model.t_part_range{t_part_ind}(1));
102 
103  % call the usual reduced basis generation POD-greedy algorithm for the current t partition
104  t_part_detailed_data{t_part_ind} = rb_basis_generation(model, detailed_data);
105 
106  % prevent recursive detailed data structures
107  if isfield(t_part_detailed_data{t_part_ind},'t_part_detailed_data')
108  t_part_detailed_data{t_part_ind} = rmfield(t_part_detailed_data{t_part_ind}, 't_part_detailed_data');
109  end
110  detailed_data.t_part_detailed_data = t_part_detailed_data;
111 
112  end
113 
114  t_parts_for_basis_generation=[];
115 
116  % create maximumm error curve over randon test sample
117  model.t_part_for_simulation = length(model.t_part_range);
118  reduced_data = gen_reduced_data(model, detailed_data);
119  sim_data = get_max_sim_data(model, reduced_data);
120 
121 
122  for p=1:length(sim_data.t_part_sim_data)
123  % compute maximum epsilon tolerance for this partition
124  epsilon_tol_p = model.RB_stop_epsilon_global/ model.nt* (model.t_part_range{p}(2)-model.t_part_range{p}(1));
125 
126  if (sim_data.t_part_sim_data{p}.Delta(end) > epsilon_tol_p)||(t_part_detailed_data{p}.RB_info.stopped_on_Nmax)
127  disp(['refinement of partition ',num2str(p),' is neccessary']);
128  % mark this partition for refinement
129  refine_p=[refine_p,p];
130  end
131  end
132 
133  if isempty(refine_p)||(loop_nr>=model.max_t_part_refinement_level)
134  continue_loop = 0;
135  else
136  % refine the partition
137  [model, detailed_data] = refine_t_part(model, detailed_data, refine_p, sim_data);
138 
139  %delete temp files for mesh refinement
140  delete(fullfile(rbmatlabtemp,'tmp*'));
141  if isfield(model,'t_parts_for_basis_generation')
142  t_parts_for_basis_generation = model.t_parts_for_basis_generation;
143  end
144 
145  % check for termination conditions: no refinement or maximum number of refinements
146 
147 
148  %check consistency
149  %number of t_part_ranges shoud be the same as number of RBs
150  if length(model.t_part_range) ~= length(detailed_data.t_part_detailed_data)
151  error('different number of t-partitions in model and detailed_data. Probably problem with refinement')
152  end
153 
154 
155  end
156  loop_nr=loop_nr+1;
157 
158  % do not stop on Nmax for last loop
159  if (loop_nr==model.max_t_part_refinement_level)
160  model.RB_stop_Nmax = 100;
161  end
162 
163  end
164 end
165 
166 
167 detailed_data.t_part_detailed_data = t_part_detailed_data;
168 
169 end
170 
171 function max_sim_data = get_max_sim_data(model, reduced_data)
172  rand('state',model.random_seed);
173  M_val=rand_uniform(20, model.mu_ranges);
174 
175  model = set_mu(model,M_val(:,1));
176 
177  max_sim_data = rb_simulation(model, reduced_data);
178 
179  for i=1:size(M_val,2)
180  model = set_mu(model, M_val(:,i));
181  sim_data = rb_simulation(model, reduced_data);
182 
183  for p=1:length(sim_data.t_part_sim_data)
184  for k=1:length(sim_data.t_part_sim_data{p}.Delta)
185  if (max_sim_data.t_part_sim_data{p}.Delta(k)<sim_data.t_part_sim_data{p}.Delta(k))
186  max_sim_data.t_part_sim_data{p}.Delta(k) = sim_data.t_part_sim_data{p}.Delta(k);
187  end
188  end
189  end
190 
191  end
192 
193  if isfield(max_sim_data,'Delta')
194  max_sim_data = rmfield(max_sim_data,'Delta');
195  end
196 
197  %test_error = rb_test_estimator(model, reduced_data,M_val);
198  %[error_max,mu_ind] = max(test_error);
199 
200  %model = set_mu(model,M_val(:,mu_ind));
201  %sim_data = rb_simulation(model, reduced_data);
202 
203 end