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.
27 % M_train_list : a cell array containing the training sets
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
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!
')
122 [detailed_data,offline_data] = RB_greedy_extension(detailed_data,params);
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 = ...
145 rb_test_indicator(detailed_data, offline_data,...
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));
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;
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
A hierarchical cubegrid of arbitrary dimension.
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.