rbmatlab  1.16.09
 All Classes Namespaces Files Functions Variables Modules Pages
lin_evol_opt_rb_operators.m
1 function [LL_I, LL_E, bb, M_E, M_b, M_EE, M_Eb, M_bb, M_I, M_II, M_IE, M_Ib, ...
2  varargout] = lin_evol_opt_rb_operators(model, detailed_data)
3 %function [LL_I, LL_E, bb, M_E, M_b, M_EE, M_Eb, M_bb, M_I, M_II, M_IE, M_Ib, ...
4 % {M_EdEd, M_IdId, M_bdbd, M_Ed, M_Id, M_bd, M_IEd, M_IId, M_Ibd, M_EEd, ...
5 % M_EId, M_Ebd, M_EdId, M_Edbd, M_Idbd]} = ...
6 % lin_evol_opt_rb_operators(model, detailed_data)
7 %
8 %
9 % function computing the time-dependent reduced basis operators and
10 % vectors.
11 %
12 % In contrast to lin_evol_rb_operators the time step (I+dt*L) is not
13 % included in the operators.
14 %
15 % required fields of model
16 % operators_algorithm: name of function for computing the
17 % L_E,L_I,b-operators with arguments (grid, params)
18 % example fv_operators_implicit
19 % compute_derivative_info: 1 if derivative information is wanted
20 % 0 if derivative information is not wanted
21 %
22 % Function supports affine decomposition, i.e. different operation modes
23 % guided by optional field decomp_mode in params. See also the
24 % contents.txt for general explanation
25 %
26 % optional fields of model:
27 % mu_names : names of fields to be regarded as parameters in vector mu
28 % decomp_mode: operation mode of the function
29 % -0= 'complete' (default): no parameter dependence or decomposition is
30 % performed. output is as described above.
31 % -1= 'components': For each output argument a cell array of output
32 % arguments is returned representing the q-th component
33 % independent of the parameters given in mu_names
34 % -2= 'coefficients': For each output argument a cell array of output
35 % arguments is returned representing the q-th coefficient
36 % dependent of the parameters given in mu_names
37 %
38 % In 'coefficients' mode the detailed data is empty.
39 
40 % Markus Dihlmann 10.06.2010. Oliver Zeeb 14.10.10
41 
42 % determine affine_decomposition_mode as integer
43 decomp_mode = model.decomp_mode;
44 
45 
46 
47 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
48 %%% decomp_mode 0 %%%
49 %%% complete calculation matrices %%%
50 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
51 if decomp_mode == 0 % complete: simple projection on RB
52 
53  %eval('u0 = ',params.init_value_algorithm,'(grid,params)']);
54  [L_I, L_E, b] = model.operators_ptr(model,detailed_data);
55 
56  % apply all operator contributions to RB:
57  %$L_E_RB = L_E * detailed_data.RB;
58  %$L_I_RB = L_I * detailed_data.RB;
59 
60  % fill output quantities
61  % LL_I : matrix
62  A = detailed_data.W;
63 
64  %reduced Operators
65  LL_I = detailed_data.RB' * A * L_I * detailed_data.RB;
66 
67  % LL_E : matrix
68  LL_E = detailed_data.RB' * A * L_E * detailed_data.RB;
69 
70  % bb : vector
71  bb = detailed_data.RB' * A * b;
72 
73 
74  %Error estimation operators
75 
76  % often used /Phi^T L_E^T A
77  RB_L_E = detailed_data.RB' * L_E' * A;
78 
79  M_E = RB_L_E * detailed_data.RB;
80  M_b = b' * A * detailed_data.RB;
81  M_EE = RB_L_E * L_E * detailed_data.RB;
82  M_Eb = RB_L_E * b;
83  M_bb = b' * A * b;
84 
85 
86  % Implicit operators
87  % often used /Phi^T L_I^T A
88  RB_L_I = detailed_data.RB' * L_I' * A ;
89 
90  M_I = RB_L_I * detailed_data.RB;
91  M_II = RB_L_I * L_I * detailed_data.RB;
92  M_IE = RB_L_I * L_E * detailed_data.RB;
93  M_Ib = RB_L_I * b;
94 
95 
96  % creation of derivative matrices
97  % get components and coefficients and do linear combination
98  if model.compute_derivative_info
99  nRB = size(detailed_data.RB,2);
100 
101  if isfield(model,'compute_derivative_indices')
102  indx_mu_i = find(model.compute_derivative_indices); %indexes i of \mu_i for which derivative information is needed
103  else
104  indx_mu_i = 1:length(get_mu(model));
105  end
106  nr_indx_mu_i= length(indx_mu_i); %nr of parameters to be optimized
107 
108  % create dummy files for derivative matrices --> memory
109  % preallocation
110  if ~isfield(model,'t')
111  model.t=1;
112  end
113 
114  model.decomp_mode = 1;
115  [dummy_LL_I, dummy_LL_E, dummy_bb, dummy_M_E, dummy_M_b, dummy_M_EE, ...
116  dummy_M_Eb, dummy_M_bb, dummy_M_I, dummy_M_II, dummy_M_IE, ...
117  dummy_M_Ib, dummy_M_EdEd, dummy_M_IdId, dummy_M_bdbd, ...
118  dummy_M_Ed, dummy_M_Id, dummy_M_bd, dummy_M_IEd, dummy_M_IId, ...
119  dummy_M_Ibd, dummy_M_EEd, dummy_M_EId, dummy_M_Ebd, dummy_M_EdId,...
120  ]=rb_operators(model,detailed_data); %dummy_M_Edbd, dummy_M_Idbd] = rb_operators(model,detailed_data); --> MEdbd, M_Idbd not needed
121 
122  if isempty(dummy_M_EdEd)
123  M_EdEd = zeros(1,1,nr_indx_mu_i);
124  else
125  M_EdEd = zeros(nRB,nRB,nr_indx_mu_i);
126  end
127 
128  if isempty(dummy_M_IdId)
129  M_IdId = zeros(1,1,nr_indx_mu_i);
130  else
131  M_IdId = zeros(nRB,nRB,nr_indx_mu_i);
132  end
133 
134  if isempty(dummy_M_Ed)
135  M_Ed = zeros(1,1,nr_indx_mu_i);
136  else
137  M_Ed = zeros(nRB,nRB,nr_indx_mu_i);
138  end
139 
140  if isempty(dummy_M_Id)
141  M_Id = zeros(1,1,nr_indx_mu_i);
142  else
143  M_Id = zeros(nRB,nRB,nr_indx_mu_i);
144  end
145 
146  if isempty(dummy_M_IEd)
147  M_IEd = zeros(1,1,nr_indx_mu_i);
148  else
149  M_IEd = zeros(nRB,nRB,nr_indx_mu_i);
150  end
151 
152  if isempty(dummy_M_IId)
153  M_IId = zeros(1,1,nr_indx_mu_i);
154  else
155  M_Id = zeros(nRB,nRB,nr_indx_mu_i);
156  end
157 
158  if isempty(dummy_M_EEd)
159  M_EEd = zeros(1,1,nr_indx_mu_i);
160  else
161  M_EEd = zeros(nRB,nRB,nr_indx_mu_i);
162  end
163 
164  if isempty(dummy_M_EId)
165  M_EId = zeros(1,1,nr_indx_mu_i);
166  else
167  M_EId = zeros(nRB,nRB,nr_indx_mu_i);
168  end
169 
170  if isempty(dummy_M_EdId)
171  M_EdId = zeros(1,1,nr_indx_mu_i);
172  else
173  M_EdId = zeros(nRB,nRB,nr_indx_mu_i);
174  end
175 
176  M_bdbd = zeros(1 ,1 ,nr_indx_mu_i);
177  M_bd = zeros(1 ,nRB,nr_indx_mu_i);
178  M_Ibd = zeros(nRB,1 ,nr_indx_mu_i);
179  M_Ebd = zeros(nRB,1 ,nr_indx_mu_i);
180  M_Edbd = zeros(nRB,1 ,nr_indx_mu_i);
181  M_Idbd = zeros(nRB,1 ,nr_indx_mu_i);
182  % end of: create dummy files for derivative matrices --> memory
183  % preallocation
184 
185 
186  %get components
187  model.decomp_mode = 1;
188  [compLL_I, compLL_E, compbb, compM_E, compM_b, compM_EE, compM_Eb,...
189  compM_bb, compM_I, compM_II, compM_IE, compM_Ib, compM_EdEd,...
190  compM_IdId, compM_bdbd, compM_Ed, compM_Id, compM_bd, compM_IEd,...
191  compM_IId, compM_Ibd, compM_EEd, compM_EId, compM_Ebd,...
192  compM_EdId, compM_Edbd, compM_Idbd] = rb_operators(model,detailed_data);
193 
194 
195  %get coefficients
196  model.decomp_mode = 2;
197  [sLL_I, sLL_E, sbb, sM_E, sM_b, sM_EE, sM_Eb, sM_bb, sM_I, sM_II,...
198  sM_IE, sM_Ib, sM_EdEd, sM_IdId, sM_bdbd, sM_Ed, sM_Id, sM_bd,...
199  sM_IEd, sM_IId, sM_Ibd, sM_EEd, sM_EId, sM_Ebd, sM_EdId,...
200  sM_Edbd, sM_Idbd] = rb_operators(model,detailed_data);
201 
202 
203  %do linear combination of components and coefficients
204  for i=1:nr_indx_mu_i
205  M_EdEd(:,:,i) = lincomb_sequence(compM_EdEd,sM_EdEd(:,i));
206  M_IdId(:,:,i) = lincomb_sequence(compM_IdId,sM_IdId(:,i));
207  M_bdbd(:,:,i) = lincomb_sequence(compM_bdbd,sM_bdbd(:,i));
208  M_Ed(:,:,i) = lincomb_sequence(compM_Ed ,sM_Ed(:,i));
209  M_Id(:,:,i) = lincomb_sequence(compM_Id ,sM_Id(:,i));
210  M_bd(:,:,i) = lincomb_sequence(compM_bd ,sM_bd(:,i));
211  M_IEd(:,:,i) = lincomb_sequence(compM_IEd ,sM_IEd(:,i));
212  M_IId(:,:,i) = lincomb_sequence(compM_IId ,sM_IId(:,i));
213  M_Ibd(:,:,i) = lincomb_sequence(compM_Ibd ,sM_Ibd(:,i));
214  M_EEd(:,:,i) = lincomb_sequence(compM_EEd ,sM_EEd(:,i));
215  M_EId(:,:,i) = lincomb_sequence(compM_EId ,sM_EId(:,i));
216  M_Ebd(:,:,i) = lincomb_sequence(compM_Ebd ,sM_Ebd(:,i));
217  M_EdId(:,:,i) = lincomb_sequence(compM_EdId,sM_EdId(:,i));
218  M_Edbd(:,:,i) = lincomb_sequence(compM_Edbd,sM_Edbd(:,i));
219  M_Idbd(:,:,i) = lincomb_sequence(compM_Idbd,sM_Idbd(:,i));
220  end
221 
222  varargout( 1) = {M_EdEd};
223  varargout( 2) = {M_IdId};
224  varargout( 3) = {M_bdbd};
225  varargout( 4) = {M_Ed};
226  varargout( 5) = {M_Id};
227  varargout( 6) = {M_bd};
228  varargout( 7) = {M_IEd};
229  varargout( 8) = {M_IId};
230  varargout( 9) = {M_Ibd};
231  varargout(10) = {M_EEd};
232  varargout(11) = {M_EId};
233  varargout(12) = {M_Ebd};
234  varargout(13) = {M_EdId};
235  varargout(14) = {M_Edbd};
236  varargout(15) = {M_Idbd};
237 
238  end %if model.compute_derivative_info
239 
240 
241 %set decomp_mode back to 0
242 model.decomp_mode = 0;
243 
244 
245 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
246 %%% decomp_mode 1 %%%
247 %%% calculation of components %%%
248 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
249 elseif decomp_mode == 1
250 
251  %eval('u0 = ',params.init_value_algorithm,'(grid,params)']);
252  [L_I, L_E, b] = model.operators_ptr(model, detailed_data);
253 
254  Q_L_I = length(L_I);
255  Q_L_E = length(L_E);
256  Q_b = length(b);
257 
258 
259  Nmax = size(detailed_data.RB,2);
260  H = size(detailed_data.RB,1);
261  A = detailed_data.W;
262 
263 
264  %reduced operator components:
265  % LL_I : matrices
266  LL_I = cell(Q_L_I,1);
267  if ~isempty(L_I)
268  LL_I(:) = {zeros(Nmax,Nmax)};
269  for q = 1:Q_L_I;
270  LL_I{q} = detailed_data.RB' * A * L_I{q} * detailed_data.RB;
271  end;
272  end
273  % LL_E : matrices
274  LL_E = cell(Q_L_E,1);
275  LL_E(:) = {zeros(Nmax,Nmax)};
276  for q = 1:Q_L_E;
277  LL_E{q} = detailed_data.RB' * A * L_E{q} * detailed_data.RB;
278  end;
279 
280  % bb : vectors
281  bb = cell(Q_b,1);
282  bb(:) = {zeros(Nmax,1)};
283  for q = 1:Q_b;
284  bb{q} = detailed_data.RB' * A * b{q};
285  end;
286 
287 
288  % init output quantities
289  M_E = cell(Q_L_E,1); M_E(:) = {zeros(Nmax,Nmax)};
290  M_b = cell(Q_b,1); M_b(:) = {zeros(Nmax,1)};
291  M_EE = cell(Q_L_E*Q_L_E,1); M_EE(:) = {zeros(Nmax,Nmax)};
292  M_Eb = cell(Q_L_E*Q_b,1); M_Eb(:) = {zeros(Nmax,1)};
293  M_bb = cell(Q_b*Q_b,1); M_bb(:) = {zeros(1,1)};
294  M_I = cell(Q_L_I,1); M_I(:) = {zeros(Nmax,Nmax)};
295  M_II = cell(Q_L_I*Q_L_I,1); M_II(:) = {zeros(Nmax,Nmax)};
296  M_IE = cell(Q_L_I*Q_L_E,1); M_IE(:) = {zeros(Nmax,Nmax)};
297  if model.compute_derivative_info
298  M_EI = cell(Q_L_I*Q_L_E,1); M_EI(:) = {zeros(Nmax,Nmax)};
299  end
300  M_Ib = cell(Q_L_I*Q_b,1); M_Ib(:) = {zeros(Nmax,1)};
301 
302 
303  % often used matrices:
304  RB_L_E = cell(Q_L_E,1);
305  RB_L_E(:) = {zeros(Nmax,H)};
306  for q=1:Q_L_E
307  RB_L_E{q}(:,:)= detailed_data.RB' * L_E{q}' * A;
308  end;
309 
310  RB_L_I = cell(Q_L_I,1);
311  RB_L_I(:) = {zeros(Nmax,H)};
312  for q=1:Q_L_I
313  RB_L_I{q}(:,:)= detailed_data.RB' * L_I{q}' * A;
314  end;
315 
316  L_E_RB = cell(Q_L_E,1);
317  L_E_RB(:) = {zeros(H,Nmax)};
318  for q=1:Q_L_E
319  L_E_RB{q} = L_E{q}*detailed_data.RB;
320  end
321 
322  L_I_RB = cell(Q_L_I,1);
323  L_I_RB(:) = {zeros(H,Nmax)};
324  for q=1:Q_L_I
325  L_I_RB{q} = L_I{q}*detailed_data.RB;
326  end
327 
328 
329  % fill output quantities
330  for q=1:Q_L_E
331  M_E{q} = RB_L_E{q} * detailed_data.RB;
332  end
333 
334  for q=1:Q_b
335  M_b{q} = b{q}' * A * detailed_data.RB;
336  end
337 
338  for q1=1:Q_L_E
339  for q2= 1:Q_L_E
340  M_EE{q2+(q1-1)*Q_L_E} = RB_L_E{q1} * L_E_RB{q2};
341  end
342  end
343 
344  for q1=1:Q_L_E
345  for q2=1:Q_b
346  M_Eb{q2+(q1-1)*Q_b} = RB_L_E{q1} * b{q2};
347  end;
348  end
349 
350  for q1=1:Q_b
351  for q2=1:Q_b
352  M_bb{q2+(q1-1)*Q_b} = b{q1}' * A * b{q2};
353  end
354  end
355 
356  for q=1:Q_L_I
357  M_I{q} = RB_L_I{q} * detailed_data.RB;
358  end
359 
360  for q1=1:Q_L_I
361  for q2= 1:Q_L_I
362  M_II{q2+(q1-1)*Q_L_I} = RB_L_I{q1} * L_I_RB{q2};
363  end
364  end
365 
366  for q1=1:Q_L_I
367  for q2= 1:Q_L_E
368  M_IE{q2+(q1-1)*Q_L_E} = RB_L_I{q1} * L_E_RB{q2};
369  end
370  end
371 
372  for q1=1:Q_L_I
373  for q2=1:Q_b
374  M_Ib{q2+(q1-1)*Q_b} = RB_L_I{q1} * b{q2};
375  end;
376  end
377 
378 
379  if model.compute_derivative_info
380  for q1=1:Q_L_E
381  for q2= 1:Q_L_I
382  M_EI{q2+(q1-1)*Q_L_I} = RB_L_E{q1} * L_I_RB{q2};
383  end
384  end
385  varargout( 1) = {M_EE}; % M_EdEd
386  varargout( 2) = {M_II}; % M_IdId
387  varargout( 3) = {M_bb}; % M_bdbd
388  varargout( 4) = {M_E} ; % M_Ed
389  varargout( 5) = {M_I} ; % M_Id
390  varargout( 6) = {M_b} ; % M_bd
391  varargout( 7) = {M_IE}; % M_IEd
392  varargout( 8) = {M_II}; % M_IId
393  varargout( 9) = {M_Ib}; % M_Ibd
394  varargout(10) = {M_EE}; % M_EEd
395  varargout(11) = {M_EI}; % M_EId
396  varargout(12) = {M_Eb}; % M_Ebd
397  varargout(13) = {M_EI}; % M_EdId
398  varargout(14) = {M_Eb}; % M_Edbd
399  varargout(15) = {M_Ib}; % M_Idbd
400  end
401 
402 
403 
404 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
405 %%% decomp_mode 2
406 %%% calculation of coefficients
407 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
408 else % decomp_mode== 2 -> coefficients
409  [L_I, L_E, b] = model.operators_ptr(model,[]);
410 
411  [sLi,sLe,sb] = model.rb_derivative_operators_coefficients_ptr(model); %gives coefficients to all the parameters in the model!
412 
413 
414  Q_L_I = length(L_I(:));
415  Q_L_E = length(L_E(:));
416  Q_b = length(b(:));
417 
418  % rb-operator coefficients can simply be forwarded
419  LL_E = L_E;
420  LL_I = L_I;
421  bb = b;
422 
423 
424  % rb error operators -coefficients
425  M_E = L_E;
426  M_b = b;
427  M_I = L_I;
428  M_EE = repmatrows(L_E(:),Q_L_E) .* repmat(L_E(:),Q_L_E,1);
429  M_Eb = repmatrows(L_E, Q_b) .* repmat(b,Q_L_E,1);
430  M_bb = repmatrows(b, Q_b) .* repmat(b,Q_b,1);
431  M_II = repmatrows(L_I(:),Q_L_I) .* repmat(L_I(:),Q_L_I,1);
432  M_Ib = repmatrows(L_I, Q_b) .* repmat(b,Q_L_I,1);
433  M_IE = repmatrows(L_I(:),Q_L_E) .* repmat(L_E(:),Q_L_I,1);
434 
435  % computation of derivative coefficients
436  if model.compute_derivative_info
437 
438  if isfield(model.optimization,'params_to_optimize')
439  indx_mu_i = find(model.optimization.params_to_optimize); %indexes i of \mu_i for which derivative information is needed
440  else
441  indx_mu_i = 1:length(get_mu(model));
442  end
443  nr_indx_mu_i= length(indx_mu_i); %nr of parameters to be optimized
444 
445 
446 
447  % use only those coefficients of parameters which are to be optimimzed
448  sLe = sLe(indx_mu_i,:)';
449  sb = sb(indx_mu_i,:)';
450  % in case of no implicit operator (e.g. sLi=[]) repmat does not work!
451  if ~isempty(sLi)
452  sLi = sLi(indx_mu_i,:)';
453  else
454  sLi=repmat(repmatrows(L_I(:),Q_L_I),1,nr_indx_mu_i); %sLi ist a 0xn empty matrix, not only a 0x1 empty matrix
455  end
456 
457  varargout( 1) = {repmatrows(sLe,Q_L_E) .* repmat(sLe,Q_L_E,1)}; % M_EdEd
458  varargout( 2) = {repmatrows(sLi,Q_L_I) .* repmat(sLi,Q_L_I,1)}; % M_IdId
459  varargout( 3) = {repmatrows(sb,Q_b) .* repmat(sb,Q_b,1)}; % M_bdbd
460  varargout( 4) = {sLe}; % M_Ed
461  varargout( 5) = {sLi}; % M_Id
462  varargout( 6) = {sb}; % M_bd
463  varargout( 7) = {repmat(repmatrows(L_I(:),Q_L_E),1,nr_indx_mu_i) .* repmat(sLe,Q_L_I,1)}; % M_IEd
464  varargout( 8) = {repmat(repmatrows(L_I(:),Q_L_I),1,nr_indx_mu_i) .* repmat(sLi,Q_L_I,1)}; % M_IId
465  varargout( 9) = {repmat(repmatrows(L_I(:),Q_b ),1,nr_indx_mu_i) .* repmat(sb,Q_L_I,1)}; % M_Ibd
466  varargout(10) = {repmat(repmatrows(L_E(:),Q_L_E),1,nr_indx_mu_i) .* repmat(sLe,Q_L_E,1)}; % M_EEd
467  varargout(11) = {repmat(repmatrows(L_E(:),Q_L_I),1,nr_indx_mu_i) .* repmat(sLi,Q_L_E,1)}; % M_EId
468  varargout(12) = {repmat(repmatrows(L_E(:),Q_b ),1,nr_indx_mu_i) .* repmat(sb,Q_L_E,1)}; % M_Ebd
469  varargout(13) = {repmatrows(sLe,Q_L_I) .* repmat(sLi,Q_L_E,1)}; % M_EdId
470  varargout(14) = {repmatrows(sLe,Q_b) .* repmat(sb,Q_L_E,1)}; % M_Edbd
471  varargout(15) = {repmatrows(sLi,Q_b) .* repmat(sb,Q_L_I,1)}; % M_Idbd
472 
473  end %if model.compute_derivative_info
474 
475 end