rbmatlab  1.16.09
 All Classes Namespaces Files Functions Variables Modules Pages
basisgen_refined_tmp.m
1 function detailed_data = basisgen_refined_tmp(detailed_data, params)
2 %function detailed_data = basisgen_refined_tmp(detailed_data,params)
3 %
4 % script demonstrating the generation of a reduced basis with
5 % an refined parameter grid approach, either uniformly refined or
6 % adaptively.
7 %
8 % RB is started with initial data by varying all parameters as
9 % specified in the coarse parameter grid M_0.
10 % For current parameter mesh M_t, the basis is extended by
11 % greedy searching for the node with the maximum posterior error
12 % estimator, for this mu a single vector is chosen for basisextension.
13 % A validation error is determined for checking, whether overfitting is
14 % currently happening. If so, the grid is refined either uniformly
15 % or by an a-posteriori error estimation approach, and the
16 % basis-extension again started. The new training set is only the
17 % set of new nodes.
18 %
19 % produced fields of detailed_data:
20 % RB : collection of orthonormal reduced basis DOF
21 % vectors
22 % RB_info: datastructure with basis-generation
23 % information, see below
24 % grid : numerical grid, is constructed if not existent.
25 %
26 % generated fields in RB_info in addition to those by RB_greedy_extension:
27 % M_train_list : a cell array containing the training sets
28 % provided to RB_greedy_extension.
29 % MMesh_list : a cell array containing the sequence of
30 % increasingly refined meshes
31 % mesh_index_sequence : a list indicating, on which refinement-level
32 % each extension-mu-vector in
33 % mu_sequence was added.
34 % RB_stopped_on_max_refinement_level : indicates, whether
35 % the maximum level of refinements was reached
36 %
37 % required fields of params in addition to those of RB_greedy_extension:
38 % RB_refinement_mode: 'uniform' or 'adaptive'
39 % RB_max_refinement_level : level after which
40 % no further refinement is performed
41 % mu_ranges : cell array of intervals for the mu-components
42 % RB_numintervals : vector indicating the number of
43 % intervals for the coarse initial mesh
44 % RB_val_rand_seed: random seed for validation set
45 % creation, if not already given in RB_info.M_val
46 % RB_M_val_size: number of validation points for extension
47 % RB_refinement_theta: theta-value in [0,1] in case of
48 % adaptive refinement
49 
50 % Bernard Haasdonk 29.5.2007
51 
52 % generate grid
53 if ~isfield(detailed_data,'grid') | isempty(detailed_data.grid)
54  detailed_data.grid = construct_grid(params);
55 end;
56 
57 if ~isfield(params,'RB_max_refinement_level')
58  params.RB_max_refinement_level = inf;
59 end;
60 
61 %%%%%%%%%%%%%%% specific settings for adaptive algorithm
62 
63 par.range = params.mu_ranges;
64 par.numintervals = [params.RB_numintervals];
65 MMesh0 = cubegrid(par);
66 Mfine = get(MMesh0,'vertex')';
67 Mcoarse = zeros(size(Mfine,1),0); % i.e. empty, but row-number as Mfine
68 
69 % generate RB-start: all initial data constellations
70 detailed_data.RB = RB_init_data_basis(Mfine, ...
71  detailed_data.grid,...
72  params);
73 
74 %prepare params
75 params.Nmax = size(detailed_data.RB,2);
76 
77 disp(['Starting RB extension loop']);
78 
79 % start time measurement
80 tic;
81 
82 % generate RB_info if required
83 if ~isfield(detailed_data,'RB_info')
84  detailed_data.RB_info = [];
85 end;
86 
87 % generate validation data if required
88 if ~isfield(detailed_data.RB_info,'M_val')
89  rand('state',params.RB_val_rand_seed); % random but deterministic
90  detailed_data.RB_info.M_val = rand_uniform(params.RB_M_val_size,...
91  params.mu_ranges);
92 end;
93 
94 ref_loop = 0; % number of refinement loops
95 MMesh = MMesh0;
96 % monitor certain quantities during extension
97 MMesh_list = {};
98 M_train_list = {};
99 mesh_index_sequence = [];
100 
101 % loop over refinement levels
102 continue_loop = 1;
103 
104 while continue_loop
105 
106  % 1. Greedy search with current grid-vertices minus coarse grid-vertices
107 
108  % temporarily set M_train
109  % M_fine = get(MMesh, 'vertex')';
110 
111  % [M_train, i_fine_without_coarse] = vectorset_difference(Mfine,Mcoarse);
112  M_train = Mfine;
113  i_fine_without_coarse = 1:size(M_train,2);
114 
115  detailed_data.RB_info.M_train = M_train;
116 
117  disp(['Size of training set in next greedy: ',num2str(size(M_train,2))]);
118 % disp('check set difference!')
119 % keyboard;
120 
121  % extend basis
122  [detailed_data,offline_data] = RB_greedy_extension(detailed_data,params);
123 
124  % extend mesh-information if required
125  nRB_new = size(detailed_data.RB,2);
126  lmi = length(mesh_index_sequence);
127  if (nRB_new> lmi)
128  mesh_index_sequence = [mesh_index_sequence, ...
129  (ref_loop+1) * ones(1,nRB_new-lmi)];
130  MMesh_list = [MMesh_list, {MMesh}];
131  M_train_list = [M_train_list, {M_train}];
132  end;
133 
134  % determine eps max on Mfine cut Mcoarse
135  i_mask = zeros(1,size(Mfine,2));
136  i_mask(i_fine_without_coarse) = 1;
137  i_fine_cut_coarse = find(i_mask == 0);
138  epsmax_fine_cut_coarse = - inf;
139  if ~isempty(i_fine_cut_coarse)
140  Mfine_cut_coarse = Mfine(:,i_fine_cut_coarse);
141  tmpparams = params;
142  tmpparams.Nmax = size(detailed_data.RB,2);
143  tmpparams.N = params.Nmax;
144  errs_fine_cut_coarse = ...
145  rb_test_indicator(detailed_data, offline_data,...
146  Mfine_cut_coarse,...
147  [],tmpparams);
148  [max_errs, mu_inds] = max(errs_fine_cut_coarse);
149  epsmax_fine_cut_coarse = max_errs(1);
150  end;
151 
152  % construct eps_max on Mfine by both eps_max over (Mfine cut Mcoarse) and
153  % (Mfine \ Mcoarse))
154 
155  if isempty(detailed_data.RB_info.max_err_sequence)
156  % if no extension has been performed: get error from M_first_errs
157  epsmax_fine_without_coarse = ...
158  max(detailed_data.RB_info.M_first_errs);
159  else % an extension has been performed
160  epsmax_fine_without_coarse = ...
161  detailed_data.RB_info.max_err_sequence(end);
162  end;
163 
164  epsmax_fine = max(epsmax_fine_cut_coarse, ...
165  epsmax_fine_without_coarse);
166 
167  continue_loop = 1;
168 
169  % check if error on current grid is sufficient
170  stopped_on_epsilon = 0;
171  if epsmax_fine < params.RB_stop_epsilon
172  stopped_on_epsilon = 1;
173  continue_loop = 0;
174  end;
175 
176  stopped_on_max_refinement_level = 0;
177  levels = get(MMesh,'level');
178  maxlev = max(levels);
179  if maxlev >= params.RB_max_refinement_level
180  stopped_on_max_refinement_level = 1;
181  end;
182 
183  skip_refinement_and_break = stopped_on_epsilon | ...
184  detailed_data.RB_info.stopped_on_timeout | ...
185  detailed_data.RB_info.stopped_on_Nmax | ...
186  stopped_on_max_refinement_level;
187 
188  if skip_refinement_and_break
189  continue_loop = 0;
190  else
191  % 2. Grid refinement
192 
193  disp('detected grid refinement necessary!!');
194 
195  nleafelements = get(MMesh,'nleafelements');
196 
197  switch params.RB_refinement_mode
198  case 'uniform'
199 
200  MMesh = refine(MMesh,1:nleafelements);
201 
202  case 'adaptive'
203 
204  % compute error estimators
205  M_last_errs = -1 * ones(1,size(Mfine,2));
206  % last call of RB_greedy_extension gave fine_without_coarse
207  M_last_errs(i_fine_without_coarse) = detailed_data.RB_info.M_last_errs;
208  % call of rb_test_indicator gave fine_cut_coarse
209  if ~isempty(i_fine_cut_coarse)
210  M_last_errs(i_fine_cut_coarse) = errs_fine_cut_coarse;
211  end;
212 
213  if ~isempty(find(M_last_errs<0))
214  error('errors for refinement-estimators not all specified!!');
215  end;
216 
217  %for all leaf e in grid M_t compute
218  % eta(e) :=
219  % max_{mu in (V(e) and cog(e) Delta (Phi_t, mu)
220  %determine Theta most violating elements
221 
222  % eta: pol0-function on mu-grid-elements, only leaf elements
223  [eta, eta_info] = rb_mu_element_error_estimators ( ...
224  detailed_data, offline_data, MMesh, ...
225  M_last_errs, ...
226  params);
227 
228  disp(['saving before grid-refinement ',num2str(ref_loop)]);
229  fn = ['tmp_adapt',num2str(ref_loop),'.mat'];
230  save(fullfile(getenv('PDEMATLABTEMP'),fn));
231 
232  [dummy , ind] = sort(-eta);
233  nmax = ceil(nleafelements * params.RB_refinement_theta);
234  li = get(MMesh,'leafelements');
235  MMesh = refine(MMesh,ind(1:nmax));
236 
237  otherwise
238  error('refinement mode unknown!');
239  end;
240 
241  % compute new Mfine
242 
243  Mcoarse = Mfine;
244  Mfine = get(MMesh, 'vertex')';
245 
246  ref_loop = ref_loop + 1;
247 
248  disp('finished grid refinement!!');
249 
250  end % end skip_refinement_and_break select
251 
252 end % end loop over refinement levels
253 
254 % prepare output:
255 detailed_data.RB_info = rmfield(detailed_data.RB_info,'M_train');
256 
257 % the following field does also not make much sense in case of
258 % multiple greedy runs, so remove
259 detailed_data.RB_info = rmfield(detailed_data.RB_info,'M_first_errs');
260 
261 detailed_data.RB_info.mesh_index_sequence = mesh_index_sequence;
262 detailed_data.RB_info.M_train_list = M_train_list;
263 detailed_data.RB_info.MMesh_list = MMesh_list;
264 
265 detailed_data.RB_info.RB_stopped_on_max_refinement_level = ...
266  stopped_on_max_refinement_level;
267 
268 detailed_data.RB_info.RB_stopped_on_epsilon = ...
269  stopped_on_epsilon;
270 
271 disp(['Generated RB basis on refined grid with ',...
272  num2str(size(detailed_data.RB,2)), ...
273  ' basis vectors.']);
274 disp(['Number of MMesh-refinement loops = ',num2str(ref_loop)]);
275 t = toc;
276 disp(['Runtime = ',num2str(t)]);
277 
278 
279 
280 
281 
282 
283 
284 % TO BE ADJUSTED TO NEW SYNTAX
285 %| \docupdate
A hierarchical cubegrid of arbitrary dimension.
Definition: cubegrid.m:17
function [ detailed_data , reduced_data ] = RB_greedy_extension(model, detailed_data)
function performing a greedy search loop for reduced basis generation.
Customizable implementation of an abstract greedy algorithm.
Definition: DuneRBLeafNode.m:1