1 function detailed_data = basisgen_refined_tmp(detailed_data, params)
2 %
function detailed_data = basisgen_refined_tmp(detailed_data,params)
4 % script demonstrating the generation of a reduced basis with
5 % an refined parameter grid approach, either uniformly refined or
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
19 % produced fields of detailed_data:
20 % RB : collection of orthonormal reduced basis DOF
22 % RB_info: datastructure with basis-generation
23 % information, see below
24 % grid : numerical grid, is constructed if not existent.
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
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
50 % Bernard Haasdonk 29.5.2007
53 if ~isfield(detailed_data,'grid') | isempty(detailed_data.grid)
54 detailed_data.grid = construct_grid(params);
57 if ~isfield(params,
'RB_max_refinement_level')
58 params.RB_max_refinement_level = inf;
61 %%%%%%%%%%%%%%% specific settings for adaptive algorithm
63 par.range = params.mu_ranges;
64 par.numintervals = [params.RB_numintervals];
66 Mfine = get(MMesh0,'vertex')';
67 Mcoarse = zeros(size(Mfine,1),0); % i.e. empty, but row-number as Mfine
69 % generate RB-start: all initial data constellations
70 detailed_data.RB = RB_init_data_basis(Mfine, ...
71 detailed_data.grid,...
75 params.Nmax = size(detailed_data.RB,2);
77 disp(['Starting RB extension loop']);
79 % start time measurement
82 % generate RB_info if required
83 if ~isfield(detailed_data,'RB_info')
84 detailed_data.RB_info = [];
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,...
94 ref_loop = 0; % number of refinement loops
96 % monitor certain quantities during extension
99 mesh_index_sequence = [];
101 % loop over refinement levels
106 % 1. Greedy search with current grid-vertices minus coarse grid-vertices
108 % temporarily set M_train
109 % M_fine = get(MMesh, 'vertex')';
111 % [M_train, i_fine_without_coarse] = vectorset_difference(Mfine,Mcoarse);
113 i_fine_without_coarse = 1:size(M_train,2);
115 detailed_data.RB_info.M_train = M_train;
117 disp([
'Size of training set in next greedy: ',num2str(size(M_train,2))]);
118 % disp(
'check set difference!')
124 % extend mesh-information if required
125 nRB_new = size(detailed_data.RB,2);
126 lmi = length(mesh_index_sequence);
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}];
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);
142 tmpparams.Nmax = size(detailed_data.RB,2);
143 tmpparams.N = params.Nmax;
144 errs_fine_cut_coarse = ...
148 [max_errs, mu_inds] = max(errs_fine_cut_coarse);
149 epsmax_fine_cut_coarse = max_errs(1);
152 % construct eps_max on Mfine by both eps_max over (Mfine cut Mcoarse) and
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);
164 epsmax_fine = max(epsmax_fine_cut_coarse, ...
165 epsmax_fine_without_coarse);
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;
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;
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;
188 if skip_refinement_and_break
193 disp('detected grid refinement necessary!!');
195 nleafelements = get(MMesh,'nleafelements');
197 switch params.RB_refinement_mode
200 MMesh = refine(MMesh,1:nleafelements);
204 % compute error estimators
205 M_last_errs = -1 * ones(1,size(Mfine,2));
207 M_last_errs(i_fine_without_coarse) = detailed_data.RB_info.M_last_errs;
209 if ~isempty(i_fine_cut_coarse)
210 M_last_errs(i_fine_cut_coarse) = errs_fine_cut_coarse;
213 if ~isempty(find(M_last_errs<0))
214 error('errors for refinement-estimators not all specified!!');
217 %for all
leaf e in grid M_t compute
219 % max_{mu in (V(e) and cog(e) Delta (Phi_t, mu)
220 %determine Theta most violating elements
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, ...
228 disp([
'saving before grid-refinement ',num2str(ref_loop)]);
229 fn = ['tmp_adapt',num2str(ref_loop),'.mat'];
230 save(fullfile(getenv(
'PDEMATLABTEMP'),fn));
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));
238 error(
'refinement mode unknown!');
244 Mfine =
get(MMesh,
'vertex')
';
246 ref_loop = ref_loop + 1;
248 disp('finished grid refinement!!
');
250 end % end skip_refinement_and_break select
252 end % end loop over refinement levels
255 detailed_data.RB_info = rmfield(detailed_data.RB_info,'M_train
');
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
');
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;
265 detailed_data.RB_info.RB_stopped_on_max_refinement_level = ...
266 stopped_on_max_refinement_level;
268 detailed_data.RB_info.RB_stopped_on_epsilon = ...
271 disp(['Generated RB basis on refined grid with
',...
272 num2str(size(detailed_data.RB,2)), ...
274 disp(['Number of MMesh-refinement loops =
',num2str(ref_loop)]);
276 disp(['Runtime =
',num2str(t)]);
284 % TO BE ADJUSTED TO NEW SYNTAX