rbmatlab  1.16.09
 All Classes Namespaces Files Functions Variables Modules Pages
lin_evol_opt_rb_derivative_simulation_separate_bases_tpart.m
1 function simulation_data = lin_evol_opt_rb_derivative_simulation_separate_bases_tpart(model, reduced_data)
2 %function simulation_data = lin_evol_opt_rb_derivative_simulation_separate_bases_tpart(model, reduced_data)
3 %
4 %performs in parallel a reduced simulation and calculation of the
5 %derivative @f$\partial_{\mu_i}u@f$
6 %
7 % This function is an extension of the \\c lin_evol_rb_simulation function.
8 % Except that in this method the time step is not included in the
9 % operators, but the time integration is conducted here by an explicit Euler.
10 % In parallel to the reduced basis online simulation a calculation of the
11 % derivative @f$\partial_{\mu_i}u@f$. The derivatives are approximated in
12 % different reduced spaces than the solution.
13 %
14 % It is assumed that the wanted mu is already set in the model.
15 % (setmu(model,mu))
16 %
17 % The calculation of the derivative is made in parallel by solving in each
18 % time step the ODE @f$\partial_t \partial_{\mu_i} a + L_E \partial_{\mu_i}
19 % a + \sum \partial_{\mu_i}(\Theta^q(\mu))L^q a -\partial_{mu_i} b(x,t,\mu) =
20 % 0@f$ plus derivative of the initial conditions.
21 % All information calculating the derivative is already contained in
22 % reduced_data or is already calculated for the normal rb_simulation except
23 % @f$\partial_{\mu_i}\Theta(\mu)@f$ and @f$\partial_{\mu_i}\Theta_b @f$
24 % which are obtained calling the field \\c model.rb_derivative_operators.
25 % The field \\c model.rb_derivative_init_values points to a method
26 % providing derivative initial values.
27 %
28 % allowed dependency of data: Nmax, N, M, mu
29 % not allowed dependency of data: H
30 % allowed dependency of computation: Nmax, N, M, mu
31 % not allowed dependency of computation: H
32 %
33 % Required fields of model
34 % mu_names : the cell array of names of mu-components and
35 % for each of these strings, a corresponding
36 % field in model is expected.
37 % T : end-time of simulation
38 % nt : number of timesteps to compute
39 % error_norm : 'l2' or 'energy'
40 % L_I_inv_norm_bound: upper bound on implicit operator inverse norm
41 % constant in time, a single scalar
42 % L_E_norm_bound: upper bounds on explicit operator
43 % constant in time, a single scalar
44 % L_E_der_norm_bound: upper bound on explicit derivative operator
45 % constant in time, a single scalar
46 % L_I_der_norm_bound: upper bound on implicit derivative operator
47 % constant in time, a single scalar
48 % energy_norm_gamma: gamma >= 0 defining the weight in the energy norm
49 % coercivity_bound_name : name of the function, which can
50 % computes a lower bound on the coercivity, e.g.
51 % fv_coercivity_bound
52 % compute_derivative_info: 0 if no derivative information is needed
53 % 1 if derivative information is needed
54 %
55 % plus fields required by coercivity_bound_name.m
56 %
57 % optional fields of model:
58 % data_const_in_time : if this flag is set, then only operators
59 % for first time instance are computed
60 % name_output_functional: if this field is existent, then an
61 % output estimation is performed and error.estimtations
62 %
63 % generated fields of simulation_data:
64 % a: time sequence of solution coefficients, columns are
65 % 'a(:,k)' = `a^(k-1)`
66 % c: matrix(3-dimensional) of time secquence representing the derivative of a for
67 % @f$\mu_i@f$
68 % Delta:time sequence of `L^2`-posteriori error estimates
69 % 'Delta(k)' = `\Delta^{k-1}_N`
70 % or energy-norm-posterior error estimates
71 % 'Delta_energy(k)' = `\bar \Delta^{k-1}_N`
72 % depending on the field "error_norm" in model.
73 %
74 % if model.name_output_functional is set, then additionally, a
75 % sequence of output estimates 's(U(:,))' and error bound Delta_s is returned
76 %
77 % see the ***_gen_reduced_data() routine for specifications of the
78 % fields of reduced_data
79 
80 % Markus Dihlmann 25.03.2011
81 
82 if model.verbose >= 10
83  if model.compute_derivative_info
84  printf('performing simulation and derivative calculation for mu = %s', mat2str(get_mu(model)));
85  else
86  printf('performing simulation');
87  end
88 end;
89 
90 if (~isequal(model.error_norm,'l2')) && ...
91  (~isequal(model.error_norm,'energy'))
92  error('error_norm unknown!!');
93 end;
94 
95 %t_part_preparation
96 if ~isfield(model,'transition_model')
97  model.transition_model = 0;
98 end
99 
100 if isfield(model,'starting_time_step')
101  t_ind_start = model.starting_time_step;
102 else
103  t_ind_start = 0;
104 end
105 
106 if isfield(model, 'stopping_time_step')
107  t_ind_stop = model.stopping_time_step;
108 else
109  t_ind_stop = model.nt;
110 end
111 if ~model.transition_model
112  nRB = length(reduced_data.a0{1});
113 else
114  nRB = size(reduced_data.LL_E{1},2);
115 end
116 
117 
118 
119 %needed for derivative information
120 if ~isfield(model,'optimization')
121  disp('Model has no field optimization.')
122  model.optimization =[];
123 end
124 if isfield(model.optimization,'params_to_optimize')
125  indx_mu_i = find(model.optimization.params_to_optimize); %indexes i of \mu_i for which derivative information is needed
126 else
127  indx_mu_i = 1:length(get_mu(model));
128 end
129 nr_indx_mu_i= length(indx_mu_i); %nr of parameters to be optimized
130 
131 if isfield(model,'compute_derivative_indices')
132  indx_derivatives = find(model.compute_derivative_indices);
133 else
134  indx_derivatives = 1:length(get_mu(model));
135 end
136 
137 %cgrid = grid_geometry(params);
138 %nRB = length(reduced_data.a0{1});
139 a = zeros(nRB,t_ind_stop-t_ind_start+1); % matrix of RB-coefficients
140 Delta = zeros(1,t_ind_stop-t_ind_start+1); % vector of error estimators
141 model.dt = model.T/model.nt;
142 
143 %debug
144 res_u = zeros(1,t_ind_stop-t_ind_start+1);
145 res_u_der = zeros(nr_indx_mu_i,t_ind_stop-t_ind_start+1);
146 %end debug
147 
148 
149 if model.compute_derivative_info
150  c = cell(nr_indx_mu_i,1);
151  c_0 = cell(nr_indx_mu_i,1);
152  nRB_der = zeros(nr_indx_mu_i,1);
153  Delta_der = zeros(nr_indx_mu_i, t_ind_stop-t_ind_start+1); % vector of derivative error estimators
154  for i=1:nr_indx_mu_i
155  nRB_der(i)=length(reduced_data.c0{indx_mu_i(i)}{1});
156  c{i} = zeros(nRB_der(i), t_ind_stop-t_ind_start+1); %cell-array of RB-derivative-coefficients
157  end
158 end
159 
160 %initial erro in t-part simulations
161 if isfield(reduced_data,'Delta0')
162  Delta(1) = reduced_data.Delta0;
163 end
164 if isfield(reduced_data,'Delta_der0')
165  Delta_der(:,1) = reduced_data.Delta_der0;
166 end
167 
168 % initial data projection and linear combination
169 model.decomp_mode = 2;
170 sa0 = rb_init_values(model,[]);
171 a0 = lincomb_sequence(reduced_data.a0,sa0);
172 a(:,1) = a0;
173 
174 if model.compute_derivative_info
175  sc_0 = model.rb_derivative_init_values_coefficients(model);
176  for i=1:nr_indx_mu_i
177  c_0{i} = lincomb_sequence(reduced_data.c0{indx_mu_i(i)},sc_0(indx_mu_i(i),:)');
178  c{i}(:,1) = c_0{i};
179  end
180 
181  %initialisation of error matrices
182  %initialize cells
183  L_E_dd = cell(nr_indx_mu_i,1);
184  L_I_dd = cell(nr_indx_mu_i,1);
185  dL_E_sd = cell(nr_indx_mu_i,1);
186  dL_I_sd = cell(nr_indx_mu_i,1);
187  db = cell(nr_indx_mu_i,1);
188  K_E = cell(nr_indx_mu_i,1);
189  K_I = cell(nr_indx_mu_i,1);
190  K_EE = cell(nr_indx_mu_i,1);
191  %K_Eb = cell(nr_indx_mu_i,1);
192  %K_bb = cell(nr_indx_mu_i,1);
193  K_II = cell(nr_indx_mu_i,1);
194  %K_Ib = cell(nr_indx_mu_i,1);
195  K_IE = cell(nr_indx_mu_i,1);
196  K_EdEd = cell(nr_indx_mu_i,1);
197  K_IdId = cell(nr_indx_mu_i,1);
198  K_bdbd = cell(nr_indx_mu_i,1);
199  K_Ed = cell(nr_indx_mu_i,1);
200  K_Id = cell(nr_indx_mu_i,1);
201  K_bd = cell(nr_indx_mu_i,1);
202  K_IEd = cell(nr_indx_mu_i,1);
203  K_IId = cell(nr_indx_mu_i,1);
204  K_Ibd = cell(nr_indx_mu_i,1);
205  K_EEd = cell(nr_indx_mu_i,1);
206  K_EId = cell(nr_indx_mu_i,1);
207  K_Ebd = cell(nr_indx_mu_i,1);
208  K_EdId = cell(nr_indx_mu_i,1);
209  K_Edbd = cell(nr_indx_mu_i,1);
210  K_Idbd = cell(nr_indx_mu_i,1);
211 
212  %fix size for all b-dependent matrices to prevent errors in
213  %approximation error calculation
214  for i=1:nr_indx_mu_i
215  K_bd{i} = zeros(1,nRB_der(i));
216  K_Ibd{i} = zeros(nRB_der(i),1);
217  K_Ebd{i} = zeros(nRB_der(i),1);
218  K_Edbd{i} = zeros(nRB_der(i),1);
219  K_Idbd{i} = zeros(nRB,1);
220  K_EId{i} = zeros(nRB_der(i),nRB);
221  K_EdId{i} = zeros(nRB,nRB);
222  K_I{i} = 0;
223  K_II{i} = 0;
224  K_IE{i} = 0;
225  K_IdId{i} = 0;
226  K_Id{i} = zeros(nRB_der(i),nRB);
227  K_IEd{i} = zeros(nRB_der(i),nRB);
228  K_IId{i} = zeros(nRB_der(i),nRB);
229 
230  end
231 
232 end
233 
234 
235 model.decomp_mode =2;
236 % loop over time steps: computation of a(t+1)==a^t
237 for t = (t_ind_start+1):t_ind_stop
238  model.t = (t-1)*model.dt;
239  if model.verbose >= 20
240  disp(['entered time-loop step ',num2str(t)]);
241  end;
242 
243  % linear combination of all quantities
244  % only once, if data is const in time
245  if (t==t_ind_start+1) || (~model.data_const_in_time)
246  if model.compute_derivative_info% NEU
247  [sLL_I, sLL_E, sbb, sM_E, sM_b, sM_EE, sM_Eb, sM_bb, sM_I, sM_II, sM_IE, sM_Ib, ...
248  sL_E_dd, sL_I_dd, sdL_E_sd, sdL_I_sd, sdb,...
249  sK_E, sK_I, sK_EE, sK_II, sK_IE,...
250  sK_EdEd, sK_IdId, sK_bdbd, sK_Ed, sK_Id, sK_bd, sK_IEd, sK_IId, sK_Ibd, sK_EEd, ...
251  sK_EId, sK_Ebd, sK_EdId, sK_Edbd, sK_Idbd] = rb_operators(model,[]);
252  else
253  [sLL_I, sLL_E, sbb, sM_E, sM_b, sM_EE, sM_Eb, sM_bb, sM_I, sM_II, sM_IE, sM_Ib] ...
254  = rb_operators(model,[]);
255  end
256 
257  LL_I = lincomb_sequence(reduced_data.LL_I,sLL_I);
258  LL_E = lincomb_sequence(reduced_data.LL_E,sLL_E);
259  bb = lincomb_sequence(reduced_data.bb,sbb);
260  M_E = lincomb_sequence(reduced_data.M_E,sM_E);
261  M_b = lincomb_sequence(reduced_data.M_b,sM_b);
262  M_EE = lincomb_sequence(reduced_data.M_EE,sM_EE);
263  M_Eb = lincomb_sequence(reduced_data.M_Eb,sM_Eb);
264  M_bb = lincomb_sequence(reduced_data.M_bb,sM_bb);
265  M_I = lincomb_sequence(reduced_data.M_I,sM_I);
266  M_II = lincomb_sequence(reduced_data.M_II,sM_II);
267  M_Ib = lincomb_sequence(reduced_data.M_Ib,sM_Ib);
268  M_IE = lincomb_sequence(reduced_data.M_IE,sM_IE);
269 
270 
271 %
272 % if model.data_const_in_time
273 % % inverse of implicit matrix! L_I' = Id, L_E' = L_I^(-1) *
274 % % L_E, bb = L_I^(-1) bb
275 % LL_E = (speye(size(LL_I)) - LL_I) \ LL_E;
276 % bb = (speye(size(LL_I)) - LL_I) \ bb;
277 % LL_I = zeros(size(LL_I));
278 % end;
279  end;
280 
281 
282  % solve (I-dt*LL_I) a(:,t+1) = (I+dt*LL_E) * a(:,t) + bb
283  rhs = (speye(size(LL_E)) + model.dt*LL_E) * a(:,t-t_ind_start) + model.dt*bb;
284 
285  % warning('this does not work for pure space operators (independent of delta t). In that case use the commented line below.');
286  %rhs = a(:,t) + model.dt * (LL_E * a(:,t) + bb);
287  % check whether pure explicit problem:
288  if isequal(LL_I, zeros(size(LL_I)))
289  a(:,t+1-t_ind_start) = rhs;
290  %debug 28.10.10: folgende zeile mit reingenommen, damit
291  %res_norm_times_dt_sqr ein skalar wird
292  M_Ib=zeros(nRB,1);
293  % end debug 28.10.10
294  else % solve linear system
295  a(:,t+1-t_ind_start) = (speye(size(LL_I)) - model.dt*LL_I)\rhs;
296  end;
297 
298 
299  % compute error estimate recursively
300  % Delta^k := |L_I^(k-1)|^-1 * (res_norm(k-1)+
301  % |L_E^(k-1)|*Delta^(k-1))
302  res_norm_times_dt_sqr = a(:,t+1-t_ind_start)'*a(:,t+1-t_ind_start) - 2*a(:,t+1-t_ind_start)' * a(:,t-t_ind_start) + a(:,t-t_ind_start)'*a(:,t-t_ind_start) ...
303  +2*model.dt*( a(:,t-t_ind_start)'*M_E*a(:,t-t_ind_start) - a(:,t-t_ind_start)'*M_E*a(:,t+1-t_ind_start) ...
304  +M_b*a(:,t-t_ind_start) - M_b*a(:,t+1-t_ind_start) ...
305  +a(:,t+1-t_ind_start)'*M_I*a(:,t-t_ind_start) - a(:,t+1-t_ind_start)'*M_I*a(:,t+1-t_ind_start)) ...
306  +((model.dt)^2)* ( a(:,t-t_ind_start)'*M_EE*a(:,t-t_ind_start) + a(:,t+1-t_ind_start)'*M_II*a(:,t+1-t_ind_start) + M_bb ...
307  +2*( a(:,t+1-t_ind_start)'*M_IE*a(:,t-t_ind_start) +a(:,t-t_ind_start)'*M_Eb + a(:,t+1-t_ind_start)'*M_Ib));
308 
309 
310 
311  % ensure real estimators (could be complex due to numerics)
312  if res_norm_times_dt_sqr>=0
313  res_norm_times_dt = sqrt(res_norm_times_dt_sqr);
314  else
315 % keyboard;
316  if model.verbose>=10
317  disp('residual < 0')
318  end
319  res_norm_times_dt = 0;
320  end;
321 
322  %debug
323  res_u(t-t_ind_start) = res_norm_times_dt;
324  %end debug
325 
326  if isequal(model.error_norm,'l2')
327  Delta(t+1-t_ind_start) = res_norm_times_dt + model.L_E_norm_bound * Delta(t-t_ind_start);
328  else
329  % only sum up the squares of the residuals, multiplication with coeff
330  % and squareroots are taken at the end.
331  keyboard;
332  Rkplus1sqr = max(res_norm_sqr/(params.dt.^2), 0);
333  Delta(t+1-t_ind_start) = Delta(t-t_ind_start)+Rkplus1sqr;
334  end;
335 
336 
337 
338 
339 
340 
341  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
342  %Calculation of derivative \partial_{\mu}a
343  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
344  if model.compute_derivative_info
345  model.decomp_mode=2;
346 
347  if (t==t_ind_start+1) || (~model.data_const_in_time)
348  %[sLi,sLe,sb] = model.rb_derivative_operators_coefficients_ptr(model);
349  %sLe = -sLe;
350 
351  %LL_E_derived=zeros(size(reduced_data.LL_E{1},1),size(reduced_data.LL_E{1},2),nr_indx_mu_i);
352  %bb_derived=zeros(length(reduced_data.bb{1}),nr_indx_mu_i);
353 
354  %LL_E stays the same
355  if model.compute_derivative_info
356  %[sdL_I_sd,sdL_E_sd,sdb] = model.rb_derivative_operators_coefficients_ptr(model);
357  for i=1:nr_indx_mu_i
358  if ~isempty(reduced_data.LL_I)
359  %dL_I_sd{i}=zeros(size(reduced_data.dL_I_sd{indx_mu_i(i)}{1},1),size(reduced_data.dL_I_sd{indx_mu_i(i)}{1},2),nr_indx_mu_i);
360  %L_I_dd{i}=zeros(size(reduced_data.L_I_dd{indx_mu_i(i)}{1},1),size(reduced_data.L_I_dd{indx_mu_i(i)}{1},2),nr_indx_mu_i);
361  else
362  dL_I_sd{i}=zeros(nRB_der(i),nRB);
363  end
364 
365  if ~isempty(sK_I) %only do the linear combination if the coefficients sLi are not empty
366  L_I_dd{i}(:,:) = lincomb_sequence(reduced_data.L_I_dd{indx_mu_i(i)}, sL_I_dd(:,:));
367  dL_I_sd{i}(:,:) = lincomb_sequence(reduced_data.dL_I_sd{indx_mu_i(i)}, sdL_I_sd(:,indx_mu_i(i)));
368  K_Ibd{i} = lincomb_sequence(reduced_data.K_Ibd{indx_mu_i(i)} ,sK_Ibd(:,indx_mu_i(i)));
369  K_EId{i} = lincomb_sequence(reduced_data.K_EId{indx_mu_i(i)} ,sK_EId(:,indx_mu_i(i)));
370  K_EdId{i} = lincomb_sequence(reduced_data.K_EdId{indx_mu_i(i)},sK_EdId(:,indx_mu_i(i)));
371  K_Idbd{i} = lincomb_sequence(reduced_data.K_Idbd{indx_mu_i(i)},sK_Idbd(:,indx_mu_i(i)));
372  K_I{i} = lincomb_sequence(reduced_data.K_I{indx_mu_i(i)},sK_I);
373  K_II{i} = lincomb_sequence(reduced_data.K_II{indx_mu_i(i)},sK_II);
374  K_IE{i} = lincomb_sequence(reduced_data.K_IE{indx_mu_i(i)},sK_IE);
375  K_IdId{i} = lincomb_sequence(reduced_data.K_IdId{indx_mu_i(i)},sK_IdId(:,indx_mu_i(i)));
376  K_Id{i} = lincomb_sequence(reduced_data.K_Id{indx_mu_i(i)} ,sK_Id(:,indx_mu_i(i)));
377  K_IEd{i} = lincomb_sequence(reduced_data.K_IEd{indx_mu_i(i)} ,sK_IEd(:,indx_mu_i(i)));
378  K_IId{i} = lincomb_sequence(reduced_data.K_IId{indx_mu_i(i)} ,sK_IId(:,indx_mu_i(i)));
379  end
380 
381  L_E_dd{i}(:,:) = lincomb_sequence(reduced_data.L_E_dd{indx_mu_i(i)}, sL_E_dd(:,:));
382  dL_E_sd{i}(:,:) = lincomb_sequence(reduced_data.dL_E_sd{indx_mu_i(i)}, sdL_E_sd(:,indx_mu_i(i)));
383  db{i} = lincomb_sequence(reduced_data.db{indx_mu_i(i)}, sdb(:,indx_mu_i(i)));
384  K_E{i} = lincomb_sequence(reduced_data.K_E{indx_mu_i(i)},sK_E);
385  %K_b = lincomb_sequence(reduced_data.K_b,sK_b);
386  K_EE{i} = lincomb_sequence(reduced_data.K_EE{indx_mu_i(i)},sK_EE);
387  %K_Eb = lincomb_sequence(reduced_data.K_Eb,sK_Eb);
388  %K_bb = lincomb_sequence(reduced_data.K_bb,sK_bb);
389  %K_Ib = lincomb_sequence(reduced_data.K_Ib,sK_Ib);
390  K_EdEd{i} = lincomb_sequence(reduced_data.K_EdEd{indx_mu_i(i)},sK_EdEd(:,indx_mu_i(i)));
391  K_bdbd{i} = lincomb_sequence(reduced_data.K_bdbd{indx_mu_i(i)},sK_bdbd(:,indx_mu_i(i)));
392  K_Ed{i} = lincomb_sequence(reduced_data.K_Ed{indx_mu_i(i)} ,sK_Ed(:,indx_mu_i(i)));
393  K_bd{i} = lincomb_sequence(reduced_data.K_bd{indx_mu_i(i)} ,sK_bd(:,indx_mu_i(i)));
394  K_EEd{i} = lincomb_sequence(reduced_data.K_EEd{indx_mu_i(i)} ,sK_EEd(:,indx_mu_i(i)));
395  K_Ebd{i} = lincomb_sequence(reduced_data.K_Ebd{indx_mu_i(i)} ,sK_Ebd(:,indx_mu_i(i)));
396  K_Edbd{i} = lincomb_sequence(reduced_data.K_Edbd{indx_mu_i(i)},sK_Edbd(:,indx_mu_i(i)));
397 
398  end
399  end
400 
401  end
402 
403  % solve d/dt c(:,:,i) = LL_E * c(:,:,i) + LL_E_derived a() +
404  % bb_derived
405  for i=1:nr_indx_mu_i
406  rhs_der = (speye(size(L_E_dd{i})) + model.dt*L_E_dd{i}) * c{i}(:,t-t_ind_start) + model.dt*dL_E_sd{i}* a(:,t-t_ind_start) ...
407  + model.dt*dL_I_sd{i}* a(:,t+1-t_ind_start) + model.dt*db{i};
408 
409  % check whether pure explicit problem:
410  if isequal(LL_I, zeros(size(LL_I)))
411  c{i}(:,t+1-t_ind_start) = rhs_der;
412  else %solve linear system
413  c{i}(:,t+1-t_ind_start) = (speye(size(L_I_dd{i})) - model.dt*L_I_dd{i})\rhs_der;
414  end
415 
416 
417  %compute error estrimator recuresivly
418  res_norm_der_times_dt_sqr = c{i}(:,t+1-t_ind_start)'*c{i}(:,t+1-t_ind_start) ...
419  + c{i}(:,t-t_ind_start)'*c{i}(:,t-t_ind_start) - 2*c{i}(:,t+1-t_ind_start)'*c{i}(:,t-t_ind_start) ...
420  + 2*model.dt*(-c{i}(:,t+1-t_ind_start)'*K_I{i}*c{i}(:,t+1-t_ind_start) - c{i}(:,t-t_ind_start)'*K_E{i}*c{i}(:,t+1-t_ind_start) ...
421  -a(:,t-t_ind_start)'*K_Ed{i}*c{i}(:,t+1-t_ind_start) - c{i}(:,t+1-t_ind_start)'*K_Id{i}*a(:,t+1-t_ind_start)...
422  -c{i}(:,t+1-t_ind_start)'*K_bd{i}' + c{i}(:,t+1-t_ind_start)'*K_I{i}*c{i}(:,t-t_ind_start) ...
423  +c{i}(:,t-t_ind_start)'*K_E{i}*c{i}(:,t-t_ind_start) + a(:,t-t_ind_start)'*K_Ed{i}*c{i}(:,t-t_ind_start) ...
424  +c{i}(:,t-t_ind_start)'*K_Id{i}*a(:,t+1-t_ind_start) + c{i}(:,t-t_ind_start)'*K_bd{i}' ) ... %13.10.10: bei M_bd aufpassen, nur ein Doppelpunkt?!
425  + ((model.dt)^2) * (c{i}(:,t+1-t_ind_start)'*K_II{i}*c{i}(:,t+1-t_ind_start) + c{i}(:,t-t_ind_start)'*K_EE{i}*c{i}(:,t-t_ind_start) ...
426  + a(:,t-t_ind_start)'*K_EdEd{i}*a(:,t-t_ind_start) + a(:,t+1-t_ind_start)'*K_IdId{i}*a(:,t+1-t_ind_start)...
427  + K_bdbd{i} + 2*c{i}(:,t+1-t_ind_start)'*K_IE{i}*c{i}(:,t-t_ind_start) ...
428  + 2*c{i}(:,t+1-t_ind_start)'*K_IEd{i}*a(:,t-t_ind_start) + 2*c{i}(:,t+1-t_ind_start)'*K_IId{i}*a(:,t+1-t_ind_start)...
429  + 2*c{i}(:,t+1-t_ind_start)'*K_Ibd{i} + 2*c{i}(:,t-t_ind_start)'*K_EEd{i}*a(:,t-t_ind_start) ...
430  + 2*c{i}(:,t-t_ind_start)'*K_EId{i}*a(:,t+1-t_ind_start) + 2*c{i}(:,t-t_ind_start)'*K_Ebd{i} ...
431  + 2*a(:,t-t_ind_start)'*K_EdId{i}*a(:,t+1-t_ind_start) + 2*a(:,t-t_ind_start)'*K_Edbd{i}...
432  + 2*a(:,t+1-t_ind_start)'*K_Idbd{i} );
433 
434 
435 
436 
437  % ensure real estimators (could be complex due to numerics)
438  if res_norm_der_times_dt_sqr>=0
439  res_norm_der_times_dt = sqrt(res_norm_der_times_dt_sqr);
440  else
441  res_norm_der_times_dt = 0;
442  end;
443 
444  %debug
445  res_u_der(t-t_ind_start,i) = res_norm_der_times_dt;
446  %end debug
447 
448  if isequal(model.error_norm,'l2')
449  if model.time_const
450  t_old = model.t;
451  model.t=0;
452  end
453  Delta_der(i,t+1-t_ind_start) = model.L_E_norm_bound*Delta_der(i,t-t_ind_start) + model.L_E_der_norm_bound(model,indx_derivatives(i))*Delta(t-t_ind_start) ...
454  + model.L_I_der_norm_bound(model,indx_derivatives(i))*Delta(t+1-t_ind_start) + res_norm_der_times_dt ;
455  if model.time_const
456  model.t=t_old;
457  end
458  end
459 
460 
461 
462  end %for i=1:nr_indx_mu
463  end %if model.compute_derivative_info
464 end;
465 
466 if isequal(model.error_norm,'energy')
467  % now perform scaling and squareroot of energy-error-estimate:
468  alpha = model.coercivity_bound_ptr(model);
469  if max(model.L_E_norm_bound>1) || (alpha == 0)
470  % no reasonable estimation possible
471  Delta(:) = nan;
472  else
473  C = (sqrt(1-model.L_E_norm_bound.^2)+1) * 0.5;
474  Coeff = model.dt /(4*alpha*C*(1-model.energy_norm_gamma * C));
475  Delta = sqrt(Delta * Coeff);
476  end;
477 end;
478 
479 % return simulation result
480 simulation_data = [];
481 simulation_data.a = a;
482 simulation_data.Delta = Delta;
483 %simulation_data.Delta_rel = Delta./max(u_N_norm);
484 simulation_data.LL_I = LL_I;
485 simulation_data.LL_E = LL_E;
486 simulation_data.nRB = nRB;
487 
488 %debug
489 simulation_data.res_u = res_u;
490 simulation_data.res_u_der = res_u_der;
491 %end debug
492 
493 if model.compute_derivative_info
494  simulation_data.Delta_der = Delta_der;
495  simulation_data.c = c;
496  simulation_data.nRB_der = nRB_der;
497  %for i=1:size(Delta_der,1)
498  % simulation_data.Delta_der_rel(i,:) = Delta_der(i,:)./max(u_N_der_norm(i,:));
499  %end
500 end
501 
502 if isfield(model,'name_output_functional')
503 
504  if isfield(model,'get_output')
505  simulation_data = model.get_output(model, simulation_data, reduced_data);
506  else
507  s = reduced_data.s_RB'*a;
508  if isfield(model,'s_l2norm')
509  Delta_s = model.s_l2norm * Delta; %war ursprünglich in reduced_data.s_l2norm
510  simulation_data.Delta_s = Delta_s;
511  else
512  if (model.verbose >=10)
513  disp('warning: no L2_norm for output defined in lin_evol_rb_simulation...')
514  end
515  end
516  simulation_data.s = s;
517 
518  if model.compute_derivative_info
519  %derivative output and error estimation of gradient
520  %sum=0;%for l2-norm
521  for i=1:length(c)
522  s_der{i} = reduced_data.s_RB_der{indx_mu_i(i)}'*c{i}(:,:);
523  %sum = sum + max(Delta_der(i,:))^2;
524  end
525  simulation_data.s_der = s_der;
526  simulation_data.Delta_s_grad = model.s_l2norm * max(Delta_der,[],2);%sqrt(sum); %Gradient error bound
527 
528  end
529 
530 
531  end;
532 
533 end;
534 
535 
function r = verbose(level, message, messageId)
This function displays messages depending on a message-id and/or a level. Aditionally you can set/res...
Definition: verbose.m:17