1 function simulation_data = rb_simulation_impl(rmodel, reduced_data)
2 %
function simulation_data = rb_simulation_impl(rmodel, reduced_data)
3 %
function, which performs a reduced basis online simulation
for
4 % the parameter vector `\mu \in {\cal P} \subset \mathbb{R}^p`, which is
7 % - allowed dependency of data: Nmax, N, M, mu
8 % - not allowed dependency of data: H
9 % - allowed dependency of computation: Nmax, N, M, mu
10 % - not allowed dependency of computation: H
11 % - Unknown at this stage: ---
13 % The behaviour of this simulation is controlled by the ModelDescr structure
17 % reduced_data: of type LinEvol.ReducedData
20 % simulation_data: struct holding the reduced simulation data.
22 % Required fields of descr:
23 % mu_names : the cell array of names of mu-components and
24 % for each of these stringt, a corresponding
25 % field in model is expected.
26 % T : end-time of simulation
27 % nt : number of timesteps to compute
28 % error_norm : string specifying the used error norm for residual
29 % computations. Possible values are 'l2' or 'energy'.
30 % L_I_inv_norm_bound: upper bound on implicit operator inverse norm
31 % constant in time, a single scalar
32 % L_E_norm_bound: upper bounds on explicit operator
33 % constant in time, a single scalar
34 % energy_norm_gamma: gamma >= 0 defining the weight in the energy norm
36 % optional fields of descr:
37 % data_const_in_time : if this flag is set, then only operators
38 % for first time instance are computed
39 % name_output_functional: if this field is existent, then an
40 % output estimation is performed and error.estimtations
41 % starting_time_step: in t-partition case this is the starting time step
42 % for the actual t-partitioni time step
43 % stopping_time_step: in t-partition this is the stopping time step for
44 % the actual t-partition
46 % optional fields of reduced_data:
47 % Delta0: initial error added to the total error
49 % generated fields of simulation_data:
50 % a: time sequence of solution coefficients, columns are 'a(:,k)' =
52 % Delta: time sequence of `L^2`-posteriori error estimates
53 % 'Delta(k)' = `\Delta^{k-1}_N` or energy-norm-posteriori error
54 % estimates 'Delta_energy(k)' = `\bar \Delta^{k-1}_N` depending on
55 % the field 'error_norm' in 'descr'.
56 % Delta_s: (optional) error bounds for the sequence of output
57 % estimates 's(U(:,))'. This is returned if
58 % 'model.name_output_functional' is set.
59 % s: (optional) sequence of output estimates 's(U(:,))'. This is
60 % returned if 'model.name_output_functional' is set.
61 % LL_I: @copybrief .LinEvol.ReducedData.LL_I
62 % LL_E: @copybrief .LinEvol.ReducedData.LL_E
63 % mu: parameter vector of simulation solution
65 % if 'model.name_output_functional' is set, then additionally, a sequence of
66 % output estimates 's(U(:,))' and error bound 'Delta_s' is returned
69 % Bernard Haasdonk 15.5.2007
70 % Markus Dihlmann Feb 2011
74 if model.verbose >= 10
75 disp([
'performing simulation for mu = ', mat2str(get_mu(model))]);
78 if (~isequal(model.error_norm,
'l2')) && ...
79 (~isequal(model.error_norm,
'energy'))
80 error(
'error_norm unknown!!');
83 %preparation
for t-partition
84 if ~isfield(model,'transition_model')
85 model.transition_model = 0;
88 if isfield(model,
'starting_time_step')
89 t_ind_start = model.starting_time_step;
94 if isfield(model, 'stopping_time_step')
95 t_ind_stop = model.stopping_time_step;
97 t_ind_stop = model.nt;
101 %cgrid = grid_geometry(params);
102 if ~model.transition_model
103 nRB = length(reduced_data.a0{1});
105 nRB = size(reduced_data.LL_E{1},2);
107 a = zeros(nRB,t_ind_stop-t_ind_start+1); % matrix of RB-coefficients
108 Delta = zeros(1,t_ind_stop-t_ind_start+1); % vector of error estimators
109 model.dt = model.T/model.nt;
112 % initial data projection and linear combination
113 model.decomp_mode = 2;
114 sa0 = rb_init_values(rmodel,[], 2);
115 a0 = lincomb_sequence(reduced_data.a0,sa0);
119 if isfield(reduced_data,
'Delta0')
120 Delta(1) = reduced_data.Delta0;
121 %disp(['initial error is ',num2str(Delta(1))]);
124 % loop over time steps: computation of a(t+1)==a^t
125 LL_I_is_speye = -1; % unknown state here
126 for t = (t_ind_start+1):t_ind_stop
127 model.t = (t-1)*model.dt;
129 disp(['entered time-loop step ',num2str(t)]);
132 % linear combination of all quantities
133 % only once, if data is const in time
135 if (t==t_ind_start+1) || (~model.data_const_in_time)
136 [sLL_I, sLL_E, sbb, sK_II, sK_IE, sK_EE, sm_I, sm_E, sm] = ...
137 reduced_data.rb_operators(rmodel,[], 2);
138 LL_I = lincomb_sequence(reduced_data.LL_I,sLL_I);
139 LL_E = lincomb_sequence(reduced_data.LL_E,sLL_E);
140 bb = lincomb_sequence(reduced_data.bb,sbb);
141 K_II = lincomb_sequence(reduced_data.K_II,sK_II);
142 K_IE = lincomb_sequence(reduced_data.K_IE,sK_IE);
143 K_EE = lincomb_sequence(reduced_data.K_EE,sK_EE);
144 m_I = lincomb_sequence(reduced_data.m_I,sm_I);
145 m_E = lincomb_sequence(reduced_data.m_E,sm_E);
146 m = lincomb_sequence(reduced_data.m,sm);
147 if model.data_const_in_time
148 % inverse of implicit matrix! L_I' = Id, L_E' = L_I^(-1) *
149 % L_E, bb = L_I^(-1) bb
153 LL_I = speye(size(LL_I));
155 elseif max(max(LL_I- speye(size(LL_I))))<1e-12
160 % solve LL_I* a(:,t+1) = LL_E * a(:,t) + bb
161 rhs = LL_E * a(:,t-t_ind_start) + bb;
162 % warning('this does not work for pure space operators (independent of delta t). In that case use the commented line below.');
163 %rhs = a(:,t) + model.dt * (LL_E * a(:,t) + bb);
164 % check whether pure explicit problem:
165 if (LL_I_is_speye == 1)
166 a(:,t+1-t_ind_start) = rhs;
167 elseif (LL_I_is_speye == -1) && isequal(LL_I, speye(size(LL_I)))
168 a(:,t+1-t_ind_start) = rhs;
169 else % solve linear system
170 % disp('check symmetry and choose solver accordingly!');
172 % nonsymmetric solvers:
173 % [U(:,t+1), flag] = bicgstab(Li,rhs,[],1000);
174 % [U(:,t+1), flag] = cgs(Li,rhs,[],1000);
175 % symmetric solver, non pd:
176 % see bug_symmlq for a very strange bug: cannot solve identity system!
177 % [U(:,t+1), flag] = symmlq(Li,rhs,[],1000);
178 % [U(:,t+1), flag] = minres(Li,rhs,[],1000);
179 % symmetric solver, pd:
180 % [a(:,t+1), flag] = pcg(LL_I,rhs,[],1000);
181 [a(:,t+1-t_ind_start), flag] = bicgstab(LL_I,rhs,[],1000);
183 disp(['error in system solution, solver return flag = ', ...
189 % compute error estimate recursively
190 % Delta^k := |L_I^(k-1)|^-1 * (res_norm(k-1)+
191 % |L_E^(k-1)|*Delta^(k-1))
195 a(:,t+1-t_ind_start)' * K_II * a(:,t+1-t_ind_start) ...
196 -2* a(:,t+1-t_ind_start)' * K_IE * a(:,t-t_ind_start) ...
197 + a(:,t-t_ind_start)' * K_EE * a(:,t-t_ind_start) ...
199 - 2* a(:,t+1-t_ind_start)' * m_I ...
200 + 2* a(:,t-t_ind_start)' * m_E;
201 % ensure real estimators (could be complex due to numerics)
203 res_norm = sqrt(res_norm_sqr);
208 if isequal(model.error_norm,'l2')
209 Delta(t+1-t_ind_start) = model.L_I_inv_norm_bound * ...
211 model.L_E_norm_bound * Delta(t-t_ind_start));
213 % only sum up the squares of the residuals, multiplication with coeff
214 % and squareroots are taken at the end.
215 Rkplus1sqr = max(res_norm_sqr/(params.dt.^2), 0);
216 Delta(t+1-t_ind_start) = Delta(t-t_ind_start)+Rkplus1sqr;
220 if isequal(model.error_norm,'energy')
221 % now perform scaling and squareroot of energy-error-estimate:
222 alpha = model.coercivity_bound_ptr(model);
223 if max(model.L_E_norm_bound>1) || (alpha == 0)
224 % no reasonable estimation possible
227 C = (sqrt(1-model.L_E_norm_bound.^2)+1) * 0.5;
228 Coeff = model.dt /(4*alpha*C*(1-model.energy_norm_gamma * C));
229 Delta = sqrt(Delta * Coeff);
233 % return simulation result
234 simulation_data = [];
235 simulation_data.a = a;
236 simulation_data.Delta = Delta;
237 simulation_data.LL_I = LL_I;
238 simulation_data.LL_E = LL_E;
239 simulation_data.mu = get_mu(rmodel);
241 if isfield(model,'output_function_ptr')
242 if isequal(model.error_norm,'l2')
243 s = reduced_data.s_RB * a;
244 Delta_s = reduced_data.s_l2norm * Delta;
245 simulation_data.s = s(:);
246 simulation_data.Delta_s = Delta_s(:);
248 error(['output estimation not implemented concerning energy-norm.' ...
249 ' choose l2 as error_norm!']);