rbmatlab  1.16.09
 All Classes Namespaces Files Functions Variables Modules Pages
lin_evol_rb_simulation_primal_dual.m
Go to the documentation of this file.
1 function simulation_data = lin_evol_rb_simulation_primal_dual(model, reduced_data)
2 %simulation_data = lin_evol_rb_simulation_primal_dual(model, reduced_data)
3 %
4 % function, which performs a reduced basis online simulation for
5 % the parameter vector mu, which is assumed to be set in the model
6 % class, i.e. a previous model = model.set_mu(model,[...]) is assumed
7 % to have taken place.
8 %
9 % allowed dependency of data: Nmax, N, M, mu
10 % not allowed dependency of data: H
11 % allowed dependency of computation: Nmax, N, M, mu
12 % not allowed dependency of computation: H
13 % Unknown at this stage: ---
14 %
15 % Required fields of model as required by the numerical operators
16 % mu_names : the cell array of names of mu-components and
17 % for each of these stringt, a corresponding
18 % field in model is expected.
19 % T : end-time of simulation
20 % nt : number of timesteps to compute
21 % error_norm : 'l2'
22 % L_I_inv_norm_bound: upper bound on implicit operator inverse norm
23 % constant in time, a single scalar
24 % L_E_norm_bound: upper bounds on explicit operator
25 % constant in time, a single scalar
26 %
27 % optional fields of model:
28 % data_const_in_time : if this flag is set, then only operators
29 % for first time instance are computed
30 % name_output_functional: if this field is existent, then an
31 % output estimation is performed and error.estimtations
32 % starting_tim_step: in t-partition case this is the starting time step
33 % for the actual t-partitioni time step
34 % stopping_time_step: in t-partition this is the stopping time step for
35 % the actual t-partition
36 %
37 % optional fileds of reduced_data:
38 % Delta0: initial error
39 %
40 % generated fields of simulation_data:
41 % a: time sequence of solution coefficients, columns are
42 % 'a(:,k)' = `a^(k-1)`
43 % Delta:time sequence of `L^2`-posteriori error estimates
44 % 'Delta(k)' = `\Delta^{k-1}_N`
45 % or energy-norm-posterior error estimates
46 % 'Delta_energy(k)' = `\bar \Delta^{k-1}_N`
47 % depending on the field "error_norm" in model.
48 %
49 % if model.name_output_functional is set, then additionally, a
50 % sequence of output estimates 's(U(:,))' and error bound Delta_s is returned
51 %
52 % see the ***_gen_reduced_data() routine for specifications of the
53 % fields of reduced_data
54 
55 % Bernard Haasdonk 15.5.2007
56 % Markus Dihlmann Feb 2011
57 
58 % extended to primal/dual formulaton by Dominik Garmatter:
59 %
60 % There are now 3 controlfields (model.want_dual, model.want_impproved_output and model.use_scm).
61 % These enable the following cases:
62 %
63 
64 % - model.want_improved_output = 0 and model.want_dual = 1 or 0:
65 % => lin_evol_rb_simulation_data_primal_dual(model, reduced_data) with
66 % reduced_data including the primal or dual reduced_data. Simply performs a
67 % dual or primal rb_simulation. In the dual case a functional has to be
68 % specified via model.name_output_functional. In the dual case there is
69 % never an output. In the primal case one can choose to cpmoute the reduced
70 % output (and estimator) by specifying a functional via
71 % model.name_output_functional.
72 %
73 %- model.want_improved_output = 1 and model.want_dual = 0:
74 % => lin_evol_rb_simulation_primal_dual(model, reduced_data) with
75 % reduced_data including the primal and to some extend dual reduced_data!
79 %
80 % In all the cases it is possible to choose model.error_norm to be 'l2'
81 % (standard) or 'energy'. In the 'energy'-case a lower bound to the
82 % stability constant (coercivity or infsup) is required. Either it is set
83 % in model.constant_LB or one uses the SCM (model.use_scm = 1). In case you
84 % want to use SCM you have to set model.error_norm = 'energy' and
85 % model.use_scm = 1 before generating reduced_data. Also see scm_offline
86 % for other required quantities.
87 % In the above case (model.want_improved_output = 1 and model.want_dual = 0)
88 % the energy error estimator is always used.
89 %
90 % special Note: in the european_option_pricing model a functional named
91 % 'theta' includes a partial timederivative and therefore produces a
92 % slightly different dual problem. Therefore the 'theta'-case is sometimes
93 % treated differently. Don't care about this case if you don't use the
94 % european_option_pricing model.
95 
96 % Dominik Garmatter 03.08 2012
97 
98 if ~isfield(model, 'want_dual')
99  model.want_dual = 0;
100 end
101 if ~isfield(model, 'want_improved_output')
102  model.want_improved_output = 0;
103 end
104 if ~isfield(model, 'use_scm')
105  model.use_scm = 0;
106 end
107 
108 % if you want the improved output you need the dual rb_coefficients and the
109 % dual Delta for the estimator, so they are exracted here and the reduced_data
110 % have to include them if you want the improved output!
111 if model.want_improved_output && model.want_dual == 0
112  dual_data.a = reduced_data.dual_sim_data;
113  dual_data.Delta = reduced_data.dual_Delta;
114  if isequal(size(dual_data.Delta), [1,1]) == 0 % give a warning if the dual estimator isn't measured in the energy norm
115  warning('you might have chosen the wrong error estimator for your dual rb_simulation')
116  keyboard;
117  end
118 end
119 
120 if model.verbose >= 10
121  disp(['performing simulation for mu = ', mat2str(get_mu(model))]);
122 end;
123 
124 if ~isequal(model.error_norm,'l2') && ~isequal(model.error_norm,'energy')
125  error('error_norm unknown!!');
126 end;
127 
128 if isfield(model,'starting_time_step')
129  t_ind_start = model.starting_time_step;
130 else
131  t_ind_start = 0;
132 end
133 
134 if isfield(model, 'stopping_time_step')
135  t_ind_stop = model.stopping_time_step;
136 else
137  t_ind_stop = model.nt;
138 end
139 
140 % number of reduced Basis vectors
141 nRB = length(reduced_data.a0{1});
142 
143 a = zeros(nRB,t_ind_stop-t_ind_start+1); % matrix of RB-coefficients
144 Delta = zeros(1,t_ind_stop-t_ind_start+1); % vector of error estimators
145 model.dt = model.T/model.nt;
146 % initialise the correction term (only in primal case with improved output)
147 if model.want_improved_output && model.want_dual == 0
148  correction_term = 0;
149 end
150 
151 % initial data projection and linear combination
152 model.decomp_mode = 2;
153 sa0 = rb_init_values(model,[]); % rb_init_values have to include the differentiation between the primal and the dual case (see for example rb_init_values_primal_dual)
154 a0 = lincomb_sequence(reduced_data.a0,sa0);
155 a0 = a0(:);
156 a(:,1) = a0; % in the dual case this will be modified in the first timestep
157 
158 % initial error
159 if isfield(reduced_data,'Delta0')
160  Delta(1) = reduced_data.Delta0;
161 end
162 
163 LL_I_is_speye = -1; % unknown state here
164 
165 for t = (t_ind_start+1):t_ind_stop % loop over time steps: computation of a(t+1)==a^t
166  % set the time - relevant if your operators are dependant on time.
167  if model.want_dual == 0
168  model.t = (t-1)*model.dt; % forward time in primal case
169  else
170  model.t = model.T - (t-1)*model.dt; % the backwards time in dual case
171  end
172 
173  if model.verbose >= 20
174  disp(['entered time-loop step ',num2str(t)]);
175  end;
176 
177  % linear combination of all quantities - only once, if data is const in time
178  if (t==t_ind_start+1) || (~model.data_const_in_time)
179  if model.want_dual %dual
180  if strcmp(model.name_output_functional, 'theta') == 1 % in theta-case the M-1-timestep is a little different and thus the residualnorm in that timestep
181  [sLL_I, sLL_E, ~, sK_II, sK_IE, sK_EE, sm_I, sm_E, sm, sK_IdId, sK_IdE, ~, ~, ~] = ...
182  rb_operators(model,[]);
183  % additional quantities for the dual resiual in the theta-case
184  m_I = lincomb_sequence(reduced_data.m_I,sm_I);
185  m_E = lincomb_sequence(reduced_data.m_E,sm_E);
186  m = lincomb_sequence(reduced_data.m,sm);
187  else % non-theta dual case
188  [sLL_I, sLL_E, ~, sK_II, sK_IE, sK_EE, ~, ~, ~, sK_IdId, sK_IdE, ~, ~, ~] = ...
189  rb_operators(model,[]);
190  end
191  LL_I = lincomb_sequence(reduced_data.LL_I, sLL_I);
192  LL_E = lincomb_sequence(reduced_data.LL_E, sLL_E);
193  % for the dual residual
194  K_II = lincomb_sequence(reduced_data.K_II, sK_II);
195  K_IE = lincomb_sequence(reduced_data.K_IE, sK_IE);
196  K_EE = lincomb_sequence(reduced_data.K_EE, sK_EE);
197  K_IdId = lincomb_sequence(reduced_data.K_IdId, sK_IdId);
198  K_IdE = lincomb_sequence(reduced_data.K_IdE, sK_IdE);
199 
200  else % primal
201  [sLL_I, sLL_E, sbb, sK_II, sK_IE, sK_EE, sm_I, sm_E, sm, ~, ~, ~, ~, ~] = ...
202  rb_operators(model,[]);
203  LL_I = lincomb_sequence(reduced_data.LL_I,sLL_I);
204  LL_E = lincomb_sequence(reduced_data.LL_E,sLL_E);
205  bb = lincomb_sequence(reduced_data.bb,sbb);
206  K_II = lincomb_sequence(reduced_data.K_II,sK_II);
207  K_IE = lincomb_sequence(reduced_data.K_IE,sK_IE);
208  K_EE = lincomb_sequence(reduced_data.K_EE,sK_EE);
209  m_I = lincomb_sequence(reduced_data.m_I,sm_I);
210  m_E = lincomb_sequence(reduced_data.m_E,sm_E);
211  m = lincomb_sequence(reduced_data.m,sm);
212  if model.want_improved_output % imrproved output includes the quantities for the correction term
213  [sLL_I, sLL_E, sbb, sK_II, sK_IE, sK_EE, sm_I, sm_E, sm, ~, ~, sLL_I_correct, sLL_E_correct, sbb_correct] = ...
214  rb_operators(model,[]);
215  LL_I = lincomb_sequence(reduced_data.LL_I,sLL_I);
216  LL_E = lincomb_sequence(reduced_data.LL_E,sLL_E);
217  bb = lincomb_sequence(reduced_data.bb,sbb);
218  K_II = lincomb_sequence(reduced_data.K_II,sK_II);
219  K_IE = lincomb_sequence(reduced_data.K_IE,sK_IE);
220  K_EE = lincomb_sequence(reduced_data.K_EE,sK_EE);
221  m_I = lincomb_sequence(reduced_data.m_I,sm_I);
222  m_E = lincomb_sequence(reduced_data.m_E,sm_E);
223  m = lincomb_sequence(reduced_data.m,sm);
224  LL_I_correct = lincomb_sequence(reduced_data.LL_I_correct, sLL_I_correct);
225  LL_E_correct = lincomb_sequence(reduced_data.LL_E_correct, sLL_E_correct);
226  bb_correct = lincomb_sequence(reduced_data.bb_correct, sbb_correct);
227  end
228  end
229 
230 % fix for the dual endtimecondition: up to know it is only u^{du,M} = P_{X_N}(v_l) but we need u^{du,M} = -L_I^{du,-1}P_{X_N}(v_l)
231  if (t==t_ind_start+1) && model.want_dual
232  a(:,1) = - LL_I \ a0; % a0 was exactly P_{X_N}(v_l)
233  end
234 
235  if model.data_const_in_time
236  if model.want_dual % dual
237  % save the original L_E for the last timestep of the dual problem
238  % (it is a(:,t_ind_stop) = L_E*a(:,t_ind_stop+1), i.e. L_I =
239  % Id in the last timestep of the dual problem)
240  LL_E_save = LL_E;
241  LL_E = LL_I \ LL_E;
242  if strcmp(model.name_output_functional, 'theta') == 1 % prepare the inhomogeniety of the theta-case
243  a0 = LL_I \ a0;
244  end
245  LL_I = speye(size(LL_I));
246  LL_I_is_speye = 1;
247 
248  else % primal
249  LL_E = LL_I \ LL_E;
250  bb = LL_I \ bb;
251  LL_I = speye(size(LL_I));
252  LL_I_is_speye = 1;
253  end
254  elseif max(max(LL_I- speye(size(LL_I))))<1e-12
255  LL_I_is_speye = 1;
256  end
257  end
258 
259  if model.want_dual % solve LL_I * a(:,t_adj) = LL_E * a(:,t_adj+1), i.e. the dual problem
260  if (t==t_ind_start+1) && strcmp(model.name_output_functional, 'theta') == 1 % slightly different dual problem in case of the theta-functional
261  rhs = LL_E * a(:,t-t_ind_start) + a0;
262  elseif (t==t_ind_stop) % the last timestep of the dual problem
263  if model.data_const_in_time % use the saved L_E here (remember LL_I_is_speye = 1 anyway in this case)
264  rhs = LL_E_save * a(:,t-t_ind_start);
265  else % otherwise use the current L_E and set L_I = Id.
266  rhs = LL_E * a(:,t-t_ind_start);
267  LL_I = speye(size(LL_I));
268  LL_I_is_speye = 1;
269  end
270  else % the normal timesteps of the dual problem
271  rhs = LL_E * a(:,t-t_ind_start);
272  end
273 
274  if LL_I_is_speye == 1
275  a(:, t+1-t_ind_start) = rhs;
276  elseif LL_I_is_speye == -1 && isequal(LL_I, speye(size(LL_I)))
277  a(:, t+1-t_ind_start) = rhs;
278  else
279  [a(:, t+1-t_ind_start), flag] = bicgstab(LL_I,rhs,[],1000);
280  if flag>0
281  disp(['error in system solution, solver return flag = ', ...
282  num2str(flag)]);
283  keyboard;
284  end
285  end
286 % the norm of the dual residual
287  if strcmp(model.name_output_functional, 'theta') == 1 && (t==t_ind_start+1) % M-1-timestep in theta-case differs!
288  res_norm_sqr = ...
289  a(:,t+1-t_ind_start)' * K_II * a(:,t+1-t_ind_start) ...
290  -2* a(:,t+1-t_ind_start)' * K_IE * a(:,t-t_ind_start) ...
291  + a(:,t-t_ind_start)' * K_EE * a(:,t-t_ind_start) ...
292  + m ...
293  - 2* a(:,t+1-t_ind_start)' * m_I ...
294  + 2* a(:,t-t_ind_start)' * m_E;
295  elseif (t==t_ind_stop)
296  res_norm_sqr = ...
297  a(:,t+1-t_ind_start)' * K_IdId * a(:,t+1-t_ind_start) ...
298  -2* a(:,t+1-t_ind_start)' * K_IdE * a(:,t-t_ind_start) ...
299  + a(:,t-t_ind_start)' * K_EE * a(:,t-t_ind_start);
300  else
301  res_norm_sqr = ...
302  a(:,t+1-t_ind_start)' * K_II * a(:,t+1-t_ind_start) ...
303  -2* a(:,t+1-t_ind_start)' * K_IE * a(:,t-t_ind_start) ...
304  + a(:,t-t_ind_start)' * K_EE * a(:,t-t_ind_start);
305  end
306 
307  else % solve LL_I* a(:,t+1) = LL_E * a(:,t) + bb, i.e. the primal problem
308  rhs = LL_E * a(:,t-t_ind_start) + bb;
309  if (LL_I_is_speye == 1)
310  a(:,t+1-t_ind_start) = rhs;
311  elseif (LL_I_is_speye == -1) && isequal(LL_I, speye(size(LL_I)))
312  a(:,t+1-t_ind_start) = rhs;
313  else
314  [a(:,t+1-t_ind_start), flag] = bicgstab(LL_I,rhs,[],1000);
315  if flag>0
316  disp(['error in system solution, solver return flag = ', ...
317  num2str(flag)]);
318  keyboard;
319  end
320  end
321  res_norm_sqr = ...
322  a(:,t+1-t_ind_start)' * K_II * a(:,t+1-t_ind_start) ...
323  -2* a(:,t+1-t_ind_start)' * K_IE * a(:,t-t_ind_start) ...
324  + a(:,t-t_ind_start)' * K_EE * a(:,t-t_ind_start) ...
325  + m ...
326  - 2* a(:,t+1-t_ind_start)' * m_I ...
327  + 2* a(:,t-t_ind_start)' * m_E;
328  if model.want_improved_output && model.want_dual == 0 % sum up the correction term for the improved output
329  scalarproduct = (- a(:,t+1-t_ind_start)' * LL_I_correct * dual_data.a(:,t+1-t_ind_start)...
330  + a(:,t-t_ind_start)' * LL_E_correct * dual_data.a(:,t+1-t_ind_start)...
331  + bb_correct * dual_data.a(:,t+1-t_ind_start));
332  correction_term = correction_term + scalarproduct;
333  end
334  end
335 
336  % ensure real estimators (could be complex due to numerics)
337  if res_norm_sqr>=0
338  res_norm = sqrt(res_norm_sqr);
339  else
340  res_norm = 0;
341  end
342  % either use the estimator measured in l2 norm or energy norm. In the
343  % energy case just square the residuals - the rest will be done later.
344  if isequal(model.error_norm, 'l2')
345  Delta(t+1-t_ind_start) = model.L_I_inv_norm_bound * (res_norm + ...
346  model.L_E_norm_bound * Delta(t-t_ind_start));
347  elseif isequal(model.error_norm, 'energy') || (model.want_improved_output == 1 && model.want_dual == 0)
348  Delta(t+1-t_ind_start) = res_norm^2;
349  end
350 end %end of for-loop
351 
352 % Rest of energy error estimator
353 if isequal(model.error_norm, 'energy') || (model.want_improved_output == 1 && model.want_dual == 0)
354  % constant is a lower bound to alpha or beta (depending on your model)
355  % and can either be an empiric constant set in model.constant_LB or it
356  % can be computet via the SCM (see scm_offline and lin_evol_gen_reduced_data_primal_dual)
357  if model.use_scm
358  constant = scm_lower_bound(model, reduced_data);
359  else
360  constant = model.constant_LB;
361  end
362  Delta = 1/constant*sqrt(sum(Delta)); % sum up the Deltas and take the sqrt
363 end
364 
365 % return simulation result
366 simulation_data = [];
367 simulation_data.a = a;
368 simulation_data.Delta = Delta;
369 
370 
371 % in the primal case you can choose to have a reduced output
372 if isfield(model,'name_output_functional') && model.want_dual == 0
373  % ignore this 'theta'-case if you do not use the european option
374  % pricing model. See also the special note at the beginning of this function.
375  if strcmp(model.name_output_functional, 'theta') == 1
376  % the partial time-derivative in the theta-functional is approximated by \frac(P_N^k+1 - P_N^k}{model.deltaT}
377  % so we need a Matrix with the difference of the coefficients to form
378  % u_N^k+1 - u_N^k = RB * Matrix.
379  % This Matrix has the form:
380 
381  % a_1^1-a_1^0 a_1^2-a_1^1 ... a_1^M-a_1^M-1
382  % a_2^1-a_2^0 ... a_2^M-a_2^M-1
383  % . .
384  % . .
385  % . .
386  % a_N^1-a_N^0 a_N^2-a_N^1 ... a_N^M-a_N^M-1
387 
388  % for t^k und k = 0, ... ,M and N being the number of basis vectors.
389 
390  dP_dt = zeros(reduced_data.N , model.nt);
391  for i = 1 : size(dP_dt,2)
392  dP_dt(:,i) = a(:,i+1) - a(:,i);
393  end
394  if model.want_improved_output
395  % use both improved estimators (dual and primal) and multiply with
396  % the stability constant, since both estimators already include a
397  % 1/constant
398  s = reduced_data.s_RB * dP_dt(:,end) - correction_term; % in the improved output the functional is only evaluated at the last timestep
399  Delta_s = Delta*dual_data.Delta*constant;
400  else
401  s = reduced_data.s_RB * dP_dt;
402  if isequal(model.error_norm,'l2')
403  % normal sequence of l2-measured estimators
404  Delta_dt = zeros(1, model.nt);
405  for i = 1 : size(dP_dt,2)
406  Delta_dt(i) = simulation_data.Delta(1,i+1) + simulation_data.Delta(1,i); % this builds the standard theta error estimator, which was a sum of 2 l2-error estimators
407  end
408  Delta_s = reduced_data.s_l2norm .* Delta_dt; % continuousity constant * RB-Koeffs
409  elseif isequal(model.error_norm,'energy')
410  % estimator measured in energy norm
411  Delta_s = 2*reduced_data.s_l2norm * Delta; % 2*continuousity constant of the functional in theta case
412  end
413  end
414  simulation_data.s = s(:);
415  simulation_data.Delta_s = Delta_s(:);
416  else
417  if model.want_improved_output
418  % use both improved estimators (dual and primal) and multiply with
419  % the stability constant, since both estimators already include a
420  % 1/constant
421  s = reduced_data.s_RB * a(:,end) - correction_term; % in the improved output the functional is only evaluated at the last timestep
422  Delta_s = Delta.*dual_data.Delta*constant;
423  else
424  s = reduced_data.s_RB * a;
425  if isequal(model.error_norm,'l2')
426  % normal sequence of l2-measured estimators
427  Delta_s = reduced_data.s_l2norm .* Delta; % continuousity constant * RB-Koeffs
428  elseif isequal(model.error_norm,'energy')
429  % estimator measured in energy norm
430  Delta_s = reduced_data.s_l2norm * Delta; % remember: Delta is different now
431  end
432  end
433  simulation_data.s = s(:);
434  simulation_data.Delta_s = Delta_s(:);
435  end
436 end
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
function reduced_data = lin_evol_gen_reduced_data_primal_dual(model, detailed_data)
reduced_data = lin_evol_gen_reduced_data_primal_dual(model, detailed_data)
function reduced_data = lin_evol_primal_dual_gen_reduced_data(model, detailed_data)
reduced_data = lin_evol_primal_dual_gen_reduced_data(model, detailed_data)
function [ constant_LB , constant_UB ] = scm_lower_bound(model, reduced_data)
[constant_LB, constant_UB] = scm_lower_bound(model, reduced_data)
function a0 = rb_init_values_primal_dual(model, detailed_data)
a0 = rb_init_values_primal_dual(model, detailed_data)
function simulation_data = lin_evol_primal_dual_rb_simulation(model, reduced_data)
simulation_data = lin_evol_primal_dual_rb_simulation(model, reduced_data)
function scm_offline_data = scm_offline(model, detailed_data, M_train, D_train)
scm_offline_data = scm_offline(model, detailed_data, M_train, D_train)
Definition: scm_offline.m:17
function simulation_data = lin_evol_rb_simulation_primal_dual(model, reduced_data)
simulation_data = lin_evol_rb_simulation_primal_dual(model, reduced_data)
function detailed_data = lin_evol_primal_dual_gen_detailed_data(model, model_data)
detailed_data = lin_evol_primal_dual_gen_detailed_data(model, model_data)