3 %
function detailed_data =
6 %
function generating the reduced basis functions.
8 % Required fields of model:
9 % RB_generation_mode: there are two options: PCA_trajectories and
12 % Generated fields of detailed_data:
13 % RB: cell of the reduced basis on the subdomains
14 % I_0: cell of indizes into detailed_data.RB
15 % I_G: cell of indizes into detailed_data.RB
16 % RB_info: information about the generation process, depends on
21 detailed_data = model_data;
23 dir = model.dirichlet_side;
24 no_dir = mod(dir,2)+1;
26 [RB_0,RB_G] = model.get_rb_from_detailed_data(detailed_data);
29 switch model.RB_generation_mode
31 case 'PCA_trajectories'
33 % generate parameter sample
34 mu_params.numintervals = model.RB_numintervals;
35 mu_params.range = model.mu_ranges;
38 mu_sample =
get(mu_mesh,
'vertex');
42 for m = 1:size(mu_sample,1)
44 model = set_mu(model,mu_sample(m,:));
45 sim_data = detailed_simulation(model,detailed_data);
48 if sim_data.reached_maxiter
49 warning('detailed simulation may have diverged!');
51 % dof matrizes depend on model.sequence_mode
52 U = model.get_dofs_from_sim_data(model,sim_data);
62 [RB_0, RB_G] = model.orthonormalize(model, detailed_data, X, model.RB_generation_epsilon);
64 detailed_data = model.set_rb_in_detailed_data(detailed_data, RB_0, RB_G);
69 % check
if former greedy steps are available
70 if ~isfield(detailed_data,
'RB_info')
71 detailed_data.RB_info = [];
73 bases_sizes.N0 = cell(1,2);
74 bases_sizes.NG = cell(1,2);
79 step = detailed_data.RB_info.steps;
80 bases_sizes = detailed_data.RB_info.bases_sizes;
81 max_errs = detailed_data.RB_info.max_errs;
82 min_thetas = detailed_data.RB_info.min_thetas;
83 max_err_mus = detailed_data.RB_info.max_err_mus;
85 detailed_data.RB_info.stopped_on_tolerance = 0;
86 detailed_data.RB_info.stopped_on_empty_extension = 0;
88 % generate training parameter sample
89 par.numintervals = model.RB_numintervals;
90 par.range = model.base_model.mu_ranges;
93 detailed_data.RB_info.M_train = get(par_mesh,'vertex');
95 % step 0: initialize basis
96 [N0,NG] = model.get_rb_size(detailed_data);
101 switch model.RB_init_mode
103 case 'random' % random empirical mode
104 if ~all([N0{dir},NG{1},NG{2}])
105 model = set_mu(model,rand_uniform(1,model.mu_ranges)
');
107 old_ext_method = model.RB_extension_method;
108 model.RB_extension_method = 'A
';
109 [RBext_0,RBext_G] = model.RB_extension_PCA_fixspace(...
110 model,detailed_data);
111 model.RB_extension_method = old_ext_method;
114 case 'eigenbasis
' % eigenbasis on interface plus random empirical mode
116 model = set_mu(model,rand_uniform(1,model.mu_ranges)');
118 [RBext_0,RBext_G] = model.RB_extension_PCA_fixspace(...
119 model,detailed_data);
121 if (NG{1} < model.N_eigenbasis || NG{2} < model.N_eigenbasis)
122 RBext_G = model.RB_extension_eigenbasis(model,detailed_data);
124 RBext_G{i} = RBext_G{i}(:,1:(model.N_eigenbasis - NG{i}));
130 % prepare
for greedy algorithm
131 requireProjection1 = strcmp(model.RB_error_indicator, ...
132 'projection_error1');
133 requireProjection2 = strcmp(model.RB_error_indicator, ...
134 'projection_error2');
135 requireReducedData = ...
136 strcmp(model.RB_error_indicator,
'rb_error_estimate') | ...
137 strcmp(model.RB_error_indicator,
'rb_error');
138 requireBaseSimData = ...
139 strcmp(model.RB_error_indicator,
'projection_error1') | ...
140 strcmp(model.RB_error_indicator,
'projection_error2') | ...
141 strcmp(model.RB_error_indicator,
'rb_error');
143 % save detailed simulations of the base_model
for all
144 % parameters,
if needed
145 if (requireBaseSimData)
147 detailed_data.base_model_data, ...
148 detailed_data.RB_info.M_train
', ...
150 .detailed_train_savepath);
153 % inner product matrices for the computation of the projection
155 if (requireProjection1)
156 K = model.get_inner_product_matrices(detailed_data);
158 if (requireProjection2)
159 K = model.base_model.get_inner_product_matrix( ...
160 detailed_data.base_model_data);
167 while (continue_loop)
170 N_remain = model.N_max{1};
172 N_remain = min([N_remain, model.N_max{i} - N0{i} - NG{i}]);
175 N_extend = min([N_remain, size(RBext_G{i},2)]);
177 RB_G{i} = [RB_G{i}, RBext_G{i}(:, 1:N_extend)];
179 N_extend = min([model.N_max{i} - N0{i} - NG{i} - N_extend, ...
180 size(RBext_0{i},2)]);
182 RB_0{i} = [RB_0{i}, RBext_0{i}(:, 1:N_extend)];
186 detailed_data = model.set_rb_in_detailed_data(detailed_data, ...
190 [N0,NG] = model.get_rb_size(detailed_data);
192 % check if extension is empty
193 if (~any([N0{1}-N0_old{1}, N0{2}-N0_old{2}, ...
194 NG{1}-NG_old{1}, NG{2}-NG_old{2}]))
196 detailed_data.RB_info.stopped_on_empty_extension = 1;
202 % print new basis sizes
203 disp('new basis sizes:
');
204 disp(['(N1_0,N1_G,N2_0,N2_G) = (
',num2str(N0{dir}),',
', ...
205 num2str(NG{dir}),',
',num2str(N0{no_dir}),',
', ...
206 num2str(NG{no_dir}),')
']);
209 % save new basis sizes
211 bases_sizes.N0{i} = [bases_sizes.N0{i}; N0{i}];
212 bases_sizes.NG{i} = [bases_sizes.NG{i}; NG{i}];
218 % prepare error computation
219 errs = zeros(size(detailed_data.RB_info.M_train,1),1);
220 thetas = zeros(size(detailed_data.RB_info.M_train,1),1);
222 if (requireReducedData)
223 % generate required data for the reduced simulations
224 reduced_data = gen_reduced_data(model,detailed_data);
227 if (requireProjection1)
230 XN{i} = orthonormalize([RB_0{i},RB_G{i}],K{i},[],'qr
');
234 if (requireProjection2)
235 XN = zeros(detailed_data.base_model_data.df_info.ndofs,N0{1}+ ...
238 % insert XN_0 functions
239 XN(detailed_data.masks{1},1:N0{1}) = RB_0{1};
240 XN(detailed_data.masks{2},N0{1}+(1:N0{2})) = RB_0{2};
242 % insert XN_Gamma functions
243 XN(detailed_data.masks{1},(N0{1}+N0{2}+1):end) = RB_G{1};
244 XN(detailed_data.masks{2},(N0{1}+N0{2}+1):end) = RB_G{2};
246 XN = orthonormalize(XN,K,[],'qr
');
250 % testing on maximal error
251 fprintf('testing on maximal error
');
253 for m = 1:size(detailed_data.RB_info.M_train,1)
256 model = set_mu(model,detailed_data.RB_info.M_train(m,:));
258 if (requireBaseSimData)
259 % load detailed simulation of base model
260 rb_sim_data.base_sim_data = load_detailed_simulation( ...
261 m,model.base_model.detailed_train_savepath, ...
265 switch (model.RB_error_indicator)
267 case 'projection_error1
'
269 errs(m) = dom_dec_compute_projection_error(...
270 detailed_data,XN,K,rb_sim_data.base_sim_data);
272 case 'projection_error2
'
274 errs(m) = dom_dec_compute_projection_error(...
275 detailed_data,XN,K,rb_sim_data.base_sim_data);
277 case 'rb_error_estimate
'
279 rb_sim_data = rb_simulation(model,reduced_data);
280 if rb_sim_data.reached_maxiter
281 warning(['reached maximal iterations in reduced
', ...
284 errs(m) = rb_sim_data.Delta(end);
285 thetas(m) = rb_sim_data.RB_thetaN;
289 rb_sim_data = rb_simulation(model,reduced_data);
290 if rb_sim_data.reached_maxiter
291 warning(['reached maximal iterations in reduced
', ...
296 rb_sim_data = rb_reconstruction(model,detailed_data, ...
298 rb_sim_data = model.compute_error(model,detailed_data, ...
301 errs(m) = rb_sim_data.X_err(end);
302 thetas(m) = rb_sim_data.RB_thetaN;
306 if mod(m,ceil(size(detailed_data.RB_info.M_train,1)/45)) == 1
313 % print/save new information
315 [max_err, err_ind] = max(errs);
316 [min_theta, ~] = min(thetas);
317 mu_max_err = detailed_data.RB_info.M_train(err_ind,:);
319 disp([' maximal error:
',num2str(max_err)]);
320 fprintf([' maximal error mu-vector: (
', ...
321 num2str(mu_max_err(1))]);
322 for i = 2:length(model.mu_ranges)
323 fprintf([',
',num2str(mu_max_err(i))]);
326 disp([' minimal theta:
',num2str(min_theta)]);
328 max_errs = [max_errs; max_err];
329 min_thetas = [min_thetas; min_theta];
330 max_err_mus = [max_err_mus; mu_max_err];
334 if max_err < model.RB_greedy_tolerance
336 detailed_data.RB_info.stopped_on_tolerance = 1;
339 if continue_loop == 1
340 % compute basis extensions
341 model = set_mu(model,mu_max_err);
342 [RBext_0,RBext_G] = model.RB_extension_PCA_fixspace(model, ...
348 detailed_data.RB_info.steps = step;
349 detailed_data.RB_info.bases_sizes = bases_sizes;
350 detailed_data.RB_info.max_errs = max_errs;
351 detailed_data.RB_info.min_thetas = min_thetas;
352 detailed_data.RB_info.max_err_mus = max_err_mus;
353 if model.RB_sequence_mode
354 detailed_data.RB_info.extension_method = ...
355 [model.RB_extension_method, '-SEQ
'];
357 detailed_data.RB_info.extension_method = ...
358 [model.RB_extension_method, '-NSEQ
'];