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)
9 %
function computing the time-dependent reduced basis operators and
12 % In contrast to lin_evol_rb_operators the time step (I+dt*L) is not
13 % included in the operators.
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
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
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
38 % In
'coefficients' mode the detailed data is empty.
40 % Markus Dihlmann 10.06.2010. Oliver Zeeb 14.10.10
42 % determine affine_decomposition_mode as integer
43 decomp_mode = model.decomp_mode;
47 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
49 %%% complete calculation matrices %%%
50 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
51 if decomp_mode == 0 % complete: simple projection on RB
53 %eval(
'u0 = ',params.init_value_algorithm,
'(grid,params)']);
54 [L_I, L_E, b] = model.operators_ptr(model,detailed_data);
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;
60 % fill output quantities
65 LL_I = detailed_data.RB
' * A * L_I * detailed_data.RB;
68 LL_E = detailed_data.RB' * A * L_E * detailed_data.RB;
71 bb = detailed_data.RB
' * A * b;
74 %Error estimation operators
76 % often used /Phi^T L_E^T A
77 RB_L_E = detailed_data.RB' * L_E
' * A;
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;
87 % often used /Phi^T L_I^T A
88 RB_L_I = detailed_data.RB' * L_I
' * A ;
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;
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);
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
104 indx_mu_i = 1:length(get_mu(model));
106 nr_indx_mu_i= length(indx_mu_i); %nr of parameters to be optimized
108 % create dummy files for derivative matrices --> memory
110 if ~isfield(model,'t
')
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
122 if isempty(dummy_M_EdEd)
123 M_EdEd = zeros(1,1,nr_indx_mu_i);
125 M_EdEd = zeros(nRB,nRB,nr_indx_mu_i);
128 if isempty(dummy_M_IdId)
129 M_IdId = zeros(1,1,nr_indx_mu_i);
131 M_IdId = zeros(nRB,nRB,nr_indx_mu_i);
134 if isempty(dummy_M_Ed)
135 M_Ed = zeros(1,1,nr_indx_mu_i);
137 M_Ed = zeros(nRB,nRB,nr_indx_mu_i);
140 if isempty(dummy_M_Id)
141 M_Id = zeros(1,1,nr_indx_mu_i);
143 M_Id = zeros(nRB,nRB,nr_indx_mu_i);
146 if isempty(dummy_M_IEd)
147 M_IEd = zeros(1,1,nr_indx_mu_i);
149 M_IEd = zeros(nRB,nRB,nr_indx_mu_i);
152 if isempty(dummy_M_IId)
153 M_IId = zeros(1,1,nr_indx_mu_i);
155 M_Id = zeros(nRB,nRB,nr_indx_mu_i);
158 if isempty(dummy_M_EEd)
159 M_EEd = zeros(1,1,nr_indx_mu_i);
161 M_EEd = zeros(nRB,nRB,nr_indx_mu_i);
164 if isempty(dummy_M_EId)
165 M_EId = zeros(1,1,nr_indx_mu_i);
167 M_EId = zeros(nRB,nRB,nr_indx_mu_i);
170 if isempty(dummy_M_EdId)
171 M_EdId = zeros(1,1,nr_indx_mu_i);
173 M_EdId = zeros(nRB,nRB,nr_indx_mu_i);
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
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);
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);
203 %do linear combination of components and coefficients
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));
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};
238 end %if model.compute_derivative_info
241 %set decomp_mode back to 0
242 model.decomp_mode = 0;
245 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
246 %%% decomp_mode 1 %%%
247 %%% calculation of components %%%
248 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
249 elseif decomp_mode == 1
251 %eval('u0 =
',params.init_value_algorithm,'(grid,params)
']);
252 [L_I, L_E, b] = model.operators_ptr(model, detailed_data);
259 Nmax = size(detailed_data.RB,2);
260 H = size(detailed_data.RB,1);
264 %reduced operator components:
266 LL_I = cell(Q_L_I,1);
268 LL_I(:) = {zeros(Nmax,Nmax)};
270 LL_I{q} = detailed_data.RB' * A * L_I{q} * detailed_data.RB;
274 LL_E = cell(Q_L_E,1);
275 LL_E(:) = {zeros(Nmax,Nmax)};
277 LL_E{q} = detailed_data.RB
' * A * L_E{q} * detailed_data.RB;
282 bb(:) = {zeros(Nmax,1)};
284 bb{q} = detailed_data.RB' * A * b{q};
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)};
300 M_Ib = cell(Q_L_I*Q_b,1); M_Ib(:) = {zeros(Nmax,1)};
303 % often used matrices:
304 RB_L_E = cell(Q_L_E,1);
305 RB_L_E(:) = {zeros(Nmax,H)};
307 RB_L_E{q}(:,:)= detailed_data.RB
' * L_E{q}' * A;
310 RB_L_I = cell(Q_L_I,1);
311 RB_L_I(:) = {zeros(Nmax,H)};
313 RB_L_I{q}(:,:)= detailed_data.RB
' * L_I{q}' * A;
316 L_E_RB = cell(Q_L_E,1);
317 L_E_RB(:) = {zeros(H,Nmax)};
319 L_E_RB{q} = L_E{q}*detailed_data.RB;
322 L_I_RB = cell(Q_L_I,1);
323 L_I_RB(:) = {zeros(H,Nmax)};
325 L_I_RB{q} = L_I{q}*detailed_data.RB;
329 % fill output quantities
331 M_E{q} = RB_L_E{q} * detailed_data.RB;
335 M_b{q} = b{q}
' * A * detailed_data.RB;
340 M_EE{q2+(q1-1)*Q_L_E} = RB_L_E{q1} * L_E_RB{q2};
346 M_Eb{q2+(q1-1)*Q_b} = RB_L_E{q1} * b{q2};
352 M_bb{q2+(q1-1)*Q_b} = b{q1}' * A * b{q2};
357 M_I{q} = RB_L_I{q} * detailed_data.RB;
362 M_II{q2+(q1-1)*Q_L_I} = RB_L_I{q1} * L_I_RB{q2};
368 M_IE{q2+(q1-1)*Q_L_E} = RB_L_I{q1} * L_E_RB{q2};
374 M_Ib{q2+(q1-1)*Q_b} = RB_L_I{q1} * b{q2};
379 if model.compute_derivative_info
382 M_EI{q2+(q1-1)*Q_L_I} = RB_L_E{q1} * L_I_RB{q2};
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
404 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
406 %%% calculation of coefficients
407 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
408 else % decomp_mode== 2 -> coefficients
409 [L_I, L_E, b] = model.operators_ptr(model,[]);
411 [sLi,sLe,sb] = model.rb_derivative_operators_coefficients_ptr(model); %gives coefficients to all the parameters in the model!
414 Q_L_I = length(L_I(:));
415 Q_L_E = length(L_E(:));
418 % rb-
operator coefficients can simply be forwarded
424 % rb error operators -coefficients
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);
435 % computation of derivative coefficients
436 if model.compute_derivative_info
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
441 indx_mu_i = 1:length(get_mu(model));
443 nr_indx_mu_i= length(indx_mu_i); %nr of parameters to be optimized
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!
452 sLi = sLi(indx_mu_i,:)';
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
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
473 end %
if model.compute_derivative_info