1 classdef POD < 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 =
POD(generator,
id)
33 % function rbpcafe =
POD(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;
64 dmodel = rmodel.detailed_model;
65 nmus = size(M_train, 1);
66 RBinit = zeros(model_data.grid.nelements, nmus);
68 dmodel = set_mu(dmodel, M_train(i,:));
69 dmodel.decomp_mode = 0;
71 U0 = init_values_algorithm(dmodel, model_data);
75 RBinit = orthonormalize(dmodel, model_data, RBinit);
77 if size(RBinit, 2) == 0
78 RBinit = ones(size(RBinit,1),1);
80 if isstruct(model_data)
81 detailed_data = Greedy.DataTree.Detailed.RBLeafNode(model_data, this.
id);
82 set_fields(detailed_data, this, Greedy.Plugin.
POD.info_fields);
84 detailed_data.RB = RBinit;
85 disp(['found ',num2str(size(RBinit,2)),' basis functions for',...
86 ' init data variation.']);
89 function prepare_reduced_data(this, rmodel, detailed_data)
90 % function prepare_reduced_data(this, rmodel, detailed_data)
91 % @copybrief Greedy.Plugin::
Interfaceprepare_reduced_data()
93 % @copydetails Greedy.Plugin::
Interfaceprepare_reduced_data()
97 % detailed_data: of type .Greedy.DataTree.DetailedRBLeafNode
99 % if isempty(this.last_detailed_data_handle_for_prepare_reduced) ...
100 % || this.last_detailed_data_handle_for_prepare_reduced ~= detailed_data
101 this.reduced_data = gen_reduced_data(rmodel, detailed_data);
102 % this.last_detailed_data_handle_for_prepare_reduced = detailed_data;
106 function Uapprox = generate_reduced(this, rmodel, reduced_data, detailed_data, U)
107 % function Uapprox = generate_reduced(this, rmodel, reduced_data, detailed_data, U)
108 % @copybrief Greedy.Plugin::
Interfacegenerate_reduced()
110 % @copydetails Greedy.Plugin::
Interfacegenerate_reduced()
116 if isempty(reduced_data)
117 error('prepare_reduced_data() may need to be executed before generate_reduced()!');
119 rmodel.N = get_rb_size(detailed_data);
120 rmodel.M = get_ei_size(detailed_data) - rmodel.Mstrich;
122 rb_sim_data = rb_simulation(rmodel, reduced_data);
123 rb_sim_data = rb_reconstruction(rmodel, detailed_data, rb_sim_data);
124 Uapprox = rmodel.get_dofs_from_sim_data(rb_sim_data);
127 function [max_errs, max_err_sequence, max_mu_index] = ...
128 error_estimators(this, rmodel, detailed_data, M_train)
129 % function [max_errs, max_err_sequence, max_mu_index] = error_estimators(this, rmodel, detailed_data, M_train)
130 % @copybrief Greedy.Plugin::
Defaulterror_estimators()
132 % @copydetails Greedy.Plugin::
Defaulterror_estimators()
138 % The estimator `\eta_{N,M,M
'}` used here is described in @ref DHO11
139 % "[DHO11], Chapter 5".
141 nmus = size(M_train, 1);
143 tmp_max_errs = cell(nmus, 1);
144 tmp_max_err = zeros(nmus, 1);
146 reduced_data = this.reduced_data;
147 if rmodel.crb_enabled
148 if isa(reduced_data.M, 'DataTree.INode
')
149 rmodel.M = create_scalar_tree(reduced_data.M, @(x) get(x,1) - rmodel.Mstrich);
151 rmodel.M = reduced_data.M - rmodel.Mstrich;
155 if ~FakeMPI.MPI_IsActive
160 if trmodel.verbose >= 5
163 set_mu(trmodel, M_train(i,:));
165 sim_data = rb_simulation(trmodel, reduced_data);
166 tmp_max_errs{i} = trmodel.get_estimators_from_sim_data(sim_data);
167 tmp_max_err(i) = max(tmp_max_errs{i});
170 warning('RBmatlab:simulation_error
', ...
171 ['Method rb_simulation threw an error:
', tME.message]);
172 tmp_max_errs{i} = NaN;
173 tmp_max_err(i) = NaN;
178 buf.reduced_data = reduced_data;
180 buf.M_train = M_train;
181 FakeMPI.MPI_Deploy(buf);
185 nnodes = min(FakeMPI.MPI_Size, nmus-min_mu);
186 width = max((nmus-min_mu)/(2*nnodes),1);
189 bufs = cell(1,nnodes);
191 bufs{id}.st = round((id-1)*width)+1+min_mu;
192 bufs{id}.en = min(round(id*width)+min_mu, nmus);
193 % disp(['Send:
id =
', num2str(id), ' st =
', num2str(buf.st), ' en =
', num2str(buf.en)]);
195 FakeMPI.MPI_Send(ids, bufs);
196 min_mu = bufs{end}.en;
198 running=true(1,nnodes);
204 probed_ids = FakeMPI.MPI_Probe('*
');
205 id_mask = find(running);
206 for id=intersect(probed_ids, id_mask)
207 ret_struct = FakeMPI.MPI_Recv(id);
211 assert(bufs{id}.st == st)
212 assert(bufs{id}.en == en)
213 tmp_max_errs(st:en) = ret_struct.tmp_max_errs(:);
214 tmp_max_err(st:en) = ret_struct.tmp_max_err(:);
216 disp('received corrupt data
');
217 disp(exception.message);
218 disp('resend
this job
');
219 FakeMPI.MPI_ReDeploy(id);
220 FakeMPI.MPI_Send(id, bufs{id});
226 width = max((nmus-min_mu)/(2*nnodes), 1);
230 buf.en = min(min_mu + round(width), nmus);
232 % disp(['Send:
id =
', num2str(id), ' st =
', num2str(buf.st), ' en =
', num2str(buf.en)]);
233 FakeMPI.MPI_ReDeploy(id);
234 FakeMPI.MPI_Send(id, buf);
237 countdown = countdown - 1;
240 disp(['stopped parCompErrs #
', num2str(id)]);
245 if getenv('PBS_JOBID
')
251 [dummy, max_mu_index] = max(tmp_max_err);
252 max_errs = tmp_max_err;
253 max_err_sequence = tmp_max_errs{max_mu_index};
256 function [breakloop, reason] = pre_check_for_end(this, rmodel, detailed_data)
257 % function [breakloop, reason] = pre_check_for_end(this, rmodel, detailed_data)
258 % @copybrief Greedy.Plugin.Interface.pre_check_for_end()
260 % @copydetails Greedy.Plugin.Interface.pre_check_for_end()
263 % rmodel: of type .Greedy.User.IReducedModel
264 % detailed_data: of type .Greedy.DataTree.Detailed.RBLeafNode
266 dd_leaf = get_active_leaf(detailed_data, rmodel);
269 if stopped_on_active_child(detailed_data, 'stopped_on_empty_extension
', rmodel)
271 reason = 'RB extension stopped on empty extension
';
272 elseif size(detailed_data.RB, 2) >= get_field(dd_leaf, 'stop_Nmax
')
273 set_stop_flag(dd_leaf, 'stopped_on_max_basis_size
');
275 reason = 'RB extension reached maximum basis size
';
280 % breakloop = post_check_for_end(this, detailed_data);
282 function ret_detailed_data = basis_extension(this, rmodel, detailed_data, max_err_seq, mu)
283 % function ret_detailed_data = basis_extension(this, rmodel, detailed_data, max_err_seq, mu)
284 % @copybrief Greedy.Plugin.Interface.basis_extension()
286 % @copydetails Greedy.Plugin.Interface.basis_extension()
289 % rmodel: of type .Greedy.User.IReducedModel
290 % detailed_data: of type .Greedy.DataTree.Detailed.RBLeafNode
293 % detailed_data: of type .Greedy.DataTree.Detailed.RBLeafNode
295 % Implementation details:
297 % - This method computes first POD modes of this sequence of discrete
298 % functions and adds them to the reduced basis space.
299 % - The number of POD modes used can be specified by the #noof_pca_modes
301 % - The POD is computed by the PCA_fixspace() function.
303 % extend rb basis if you receive rbeidetailedata
304 if isa(detailed_data, 'Greedy.DataTree.Detailed.PODEILeafNode
')
305 ret_detailed_data = detailed_data;
306 detailed_data = ret_detailed_data.rb;
308 ret_detailed_data = detailed_data;
310 model_data = detailed_data.model_data;
311 max_err = max(max_err_seq);
312 disp(['Detected maximum error prediction
',num2str(max_err),...
316 dmodel = rmodel.detailed_model;
318 U = generate(this.generator, dmodel, model_data);
320 cur_stop_Nmax = get_field(detailed_data,
'stop_Nmax');
321 cur_noof_pca_modes = get_field(detailed_data,
'noof_pca_modes');
322 pca_modes = min(cur_noof_pca_modes, cur_stop_Nmax - size(detailed_data.RB,2));
323 W = dmodel.get_inner_product_matrix(model_data);
324 RBext = PCA_fixspace(U, detailed_data.RB,...
327 RB = [detailed_data.RB, RBext];
328 % RBo = model_orthonormalize_qr( dmodel.descr, model_data, RB, eps );
329 % K = inner_product(dmodel, model_data, RBo, RB);
330 % [dummy,i] = max(abs(K));
333 % K = rb_engine.inner_product(model,detailed_data, RB, RB);
334 % if max(max(abs(K-eye(size(K)))))>1e-10
335 % error('error in orthonormal basis extension!');
338 detailed_data.RB = RB;
339 N = size(detailed_data.RB, 2);
340 disp(['Extended RB to length ', num2str(N)]);
342 set_stop_flag(detailed_data, 'stopped_on_empty_extension');
343 %! \todo Remember to catch this one in the container above when
344 % combining eidetailed and pod-greedy!!!!