1 classdef PODDune < Greedy.Plugin.PODCommon
2 % extension
class implementing the POD-greedy algorithm extension routines
6 % cell array of field names which shall be copied to the detailed data node
7 % structure on creation in init_basis() method.
8 info_fields = {
'stop_Nmax',
'noof_pca_modes'};
12 % maximum number of generated reduced basis vectors.
15 % number of most significant POD modes to be added to the reduced basis in
16 % each extension step.
18 % See also: basis_extension() for more information.
23 % @copybrief .Greedy::
Interfacegenerated_basis_type
25 % @copydetails .Greedy::
Interfacegenerated_basis_type
27 generated_basis_type = 'rb';
32 function rbpcafe =
PODDune(generator,
id)
33 % function rbpcafe =
PODDune(generator)
34 % constructor for the
POD-greedy extension algorithms.
37 % generator:
object of type SnapshotsGenerator.Trajectories generating
38 % the (high dimensional) basis functions, i.e. solutions of
39 % the numerical scheme for the parametrized partial
40 % differential equation.
42 rbpcafe = rbpcafe@Greedy.Plugin.
PODCommon(generator);
45 function detailed_data = init_basis(this, rmodel, model_data, M_train)
46 % function detailed_data = init_basis(this, rmodel, model_data, M_train)
47 % @copybrief .Greedy.Plugin::
Interfaceinit_basis()
49 % @copydetails Greedy.Plugin::
Interfaceinit_basis()
55 % detailed_data: of type .Greedy.DataTree.Detailed.RBLeafNode
57 % Computes the initial data functions for all parameters from the
58 % training set 'M_train' and initializes the reduced basis with an
59 % orthonormalization of these DOF vectors. The result is stored in a
60 % detailed data tree
leaf node
object.
61 if isa(model_data, 'Greedy.DataTree.Detailed.INode')
62 model_data = model_data.model_data;
65 RBinit = rmodel.descr.mexptr('init_data_basis');
67 if isstruct(model_data)
68 detailed_data = Greedy.DataTree.Detailed.DuneRBLeafNode(model_data, this.
id);
69 set_fields(detailed_data, this, Greedy.Plugin.
PODDune.info_fields);
72 detailed_data.RB = RBinit;
73 disp(['found ',num2str(size(RBinit,2)),' basis functions for',...
74 ' init data variation.']);
77 function prepare_reduced_data(this, rmodel, detailed_data)
78 % function prepare_reduced_data(this, rmodel, detailed_data)
79 % @copybrief Greedy.Plugin::
Interfaceprepare_reduced_data()
81 % @copydetails Greedy.Plugin::
Interfaceprepare_reduced_data()
85 % detailed_data: of type .Greedy.DataTree.DetailedRBLeafNode
87 % if isempty(this.last_detailed_data_handle_for_prepare_reduced) ...
88 % || this.last_detailed_data_handle_for_prepare_reduced ~= detailed_data
89 this.reduced_data = gen_reduced_data(rmodel, detailed_data);
90 % this.last_detailed_data_handle_for_prepare_reduced = detailed_data;
94 function Uapprox = generate_reduced(this, rmodel, reduced_data, detailed_data, U)
95 % function Uapprox = generate_reduced(this, rmodel, reduced_data, detailed_data, U)
96 % @copybrief Greedy.Plugin::
Interfacegenerate_reduced()
98 % @copydetails Greedy.Plugin::
Interfacegenerate_reduced()
104 if isempty(reduced_data)
105 error('prepare_reduced_data() may need to be executed before generate_reduced()!');
107 rmodel.N = get_rb_size(detailed_data);
108 rmodel.M = get_ei_size(detailed_data) - rmodel.Mstrich;
110 rb_sim_data = rb_simulation(rmodel, reduced_data);
111 rb_sim_data = rb_reconstruction(rmodel, detailed_data, rb_sim_data);
112 Uapprox = rmodel.get_dofs_from_sim_data(rb_sim_data);
115 function [max_errs, max_err_sequence, max_mu_index] = ...
116 error_estimators(this, rmodel, detailed_data, M_train)
117 % function [max_errs, max_err_sequence, max_mu_index] = error_estimators(this, rmodel, detailed_data, M_train)
118 % @copybrief Greedy.Plugin::
Defaulterror_estimators()
120 % @copydetails Greedy.Plugin::
Defaulterror_estimators()
126 % The estimator `\eta_{N,M,M
'}` used here is described in @ref DHO11
127 % "[DHO11], Chapter 5".
129 nmus = size(M_train, 1);
131 tmp_max_errs = cell(nmus, 1);
132 tmp_max_err = zeros(nmus, 1);
134 reduced_data = this.reduced_data;
135 if rmodel.crb_enabled
136 if isa(reduced_data.M, 'DataTree.INode
')
137 rmodel.M = create_scalar_tree(reduced_data.M, @(x) get(x,1) - rmodel.Mstrich);
139 rmodel.M = reduced_data.M - rmodel.Mstrich;
143 if ~FakeMPI.MPI_IsActive
148 if trmodel.verbose >= 5
151 set_mu(trmodel, M_train(i,:));
153 sim_data = rb_simulation(trmodel, reduced_data);
154 tmp_max_errs{i} = trmodel.get_estimators_from_sim_data(sim_data);
155 tmp_max_err(i) = max(tmp_max_errs{i});
158 warning('RBmatlab:simulation_error
', ...
159 ['Method rb_simulation threw an error:
', tME.message]);
160 tmp_max_errs{i} = NaN;
161 tmp_max_err(i) = NaN;
166 buf.reduced_data = reduced_data;
168 buf.M_train = M_train;
169 FakeMPI.MPI_Deploy(buf);
173 nnodes = min(FakeMPI.MPI_Size, nmus-min_mu);
174 width = max((nmus-min_mu)/(2*nnodes),1);
177 bufs = cell(1,nnodes);
179 bufs{id}.st = round((id-1)*width)+1+min_mu;
180 bufs{id}.en = min(round(id*width)+min_mu, nmus);
181 % disp(['Send:
id =
', num2str(id), ' st =
', num2str(buf.st), ' en =
', num2str(buf.en)]);
183 FakeMPI.MPI_Send(ids, bufs);
184 min_mu = bufs{end}.en;
186 running=true(1,nnodes);
192 probed_ids = FakeMPI.MPI_Probe('*
');
193 id_mask = find(running);
194 for id=intersect(probed_ids, id_mask)
195 ret_struct = FakeMPI.MPI_Recv(id);
199 assert(bufs{id}.st == st)
200 assert(bufs{id}.en == en)
201 tmp_max_errs(st:en) = ret_struct.tmp_max_errs(:);
202 tmp_max_err(st:en) = ret_struct.tmp_max_err(:);
204 disp('received corrupt data
');
205 disp(exception.message);
206 disp('resend
this job
');
207 FakeMPI.MPI_ReDeploy(id);
208 FakeMPI.MPI_Send(id, bufs{id});
214 width = max((nmus-min_mu)/(2*nnodes), 1);
218 buf.en = min(min_mu + round(width), nmus);
220 % disp(['Send:
id =
', num2str(id), ' st =
', num2str(buf.st), ' en =
', num2str(buf.en)]);
221 FakeMPI.MPI_ReDeploy(id);
222 FakeMPI.MPI_Send(id, buf);
225 countdown = countdown - 1;
228 disp(['stopped parCompErrs #
', num2str(id)]);
233 if getenv('PBS_JOBID
')
239 [dummy, max_mu_index] = max(tmp_max_err);
240 max_errs = tmp_max_err;
241 max_err_sequence = tmp_max_errs{max_mu_index};
244 function [breakloop, reason] = pre_check_for_end(this, rmodel, detailed_data)
245 % function [breakloop, reason] = pre_check_for_end(this, rmodel, detailed_data)
246 % @copybrief Greedy.Plugin.Interface.pre_check_for_end()
248 % @copydetails Greedy.Plugin.Interface.pre_check_for_end()
251 % rmodel: of type .Greedy.User.IReducedModel
252 % detailed_data: of type .Greedy.DataTree.Detailed.RBLeafNode
254 dd_leaf = get_active_leaf(detailed_data, rmodel);
257 if stopped_on_active_child(detailed_data, 'stopped_on_empty_extension
', rmodel)
259 reason = 'RB extension stopped on empty extension
';
260 elseif size(detailed_data.RB, 2) >= get_field(dd_leaf, 'stop_Nmax
')
261 set_stop_flag(dd_leaf, 'stopped_on_max_basis_size
');
263 reason = 'RB extension reached maximum basis size
';
268 % breakloop = post_check_for_end(this, detailed_data);
270 function detailed_data = basis_extension(this, rmodel, detailed_data, max_err_seq, mu)
271 % function ret_detailed_data = basis_extension(this, rmodel, detailed_data, max_err_seq, mu)
272 % @copybrief Greedy.Plugin.Interface.basis_extension()
274 % @copydetails Greedy.Plugin.Interface.basis_extension()
277 % rmodel: of type .Greedy.User.IReducedModel
278 % detailed_data: of type .Greedy.DataTree.Detailed.RBLeafNode
281 % detailed_data: of type .Greedy.DataTree.Detailed.RBLeafNode
283 % Implementation details:
285 % - This method computes first POD modes of this sequence of discrete
286 % functions and adds them to the reduced basis space.
287 % - The number of POD modes used can be specified by the #noof_pca_modes
289 % - The POD is computed by the PCA_fixspace() function.
291 % extend rb basis if you receive rbeidetailedata
293 max_err = max(max_err_seq);
294 disp(['Detected maximum error prediction
',num2str(max_err),...
298 dmodel = rmodel.detailed_model;
301 cur_stop_Nmax = get_field(detailed_data,
'stop_Nmax');
302 cur_noof_pca_modes = get_field(detailed_data,
'noof_pca_modes');
303 pca_modes = min(cur_noof_pca_modes, cur_stop_Nmax - size(detailed_data.RB,2));
305 RBext = rmodel.descr.mexptr(
'rb_extension_PCA', mu, pca_modes);
308 RB = [detailed_data.RB, RBext];
309 % RBo = model_orthonormalize_qr( dmodel.descr, model_data, RB, eps );
310 % K = inner_product(dmodel, model_data, RBo, RB);
311 % [dummy,i] = max(abs(K));
314 % K = rb_engine.inner_product(model,detailed_data, RB, RB);
315 % if max(max(abs(K-eye(size(K)))))>1e-10
316 % error('error in orthonormal basis extension!');
319 detailed_data.RB = RB;
320 N = size(detailed_data.RB, 2);
321 disp(['Extended RB to length ', num2str(N)]);
323 set_stop_flag(detailed_data, 'stopped_on_empty_extension');
324 %! \todo Remember to catch this one in the container above when
325 % combining eidetailed and pod-greedy!!!!