1 function simulation_data = lin_evol_rb_simulation(model, reduced_data)
2 %
function simulation_data = lin_evol_rb_simulation(model, reduced_data)
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
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: ---
15 % Required fields of model as required by the numerical operators and
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' or 'energy'
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 % energy_norm_gamma: gamma >= 0 defining the weight in the energy norm
27 % coercivity_bound_name : name of the function, which can
28 % computes a lower bound on the coercivity, e.g.
31 % plus fields required by coercivity_bound_name.m
33 % optional fields of model:
34 % data_const_in_time : if this flag is set, then only operators
35 % for first time instance are computed
36 % name_output_functional: if this field is existent, then an
37 % output estimation is performed and error.estimtations
38 % starting_tim_step: in t-partition case this is the starting time step
39 % for the actual t-partitioni time step
40 % stopping_time_step: in t-partition this is the stopping time step for
41 % the actual t-partition
43 % optional fileds of reduced_data:
44 % Delta0: initial error
46 % generated fields of simulation_data:
47 % a: time sequence of solution coefficients, columns are
48 % 'a(:,k)' = `a^(k-1)`
49 % Delta:time sequence of `L^2`-posteriori error estimates
50 % 'Delta(k)' = `\Delta^{k-1}_N`
51 % or energy-norm-posterior error estimates
52 % 'Delta_energy(k)' = `\bar \Delta^{k-1}_N`
53 % depending on the field "error_norm" in model.
55 % if model.name_output_functional is set, then additionally, a
56 % sequence of output estimates 's(U(:,))' and error bound Delta_s is returned
58 % see the ***_gen_reduced_data() routine for specifications of the
59 % fields of reduced_data
61 % Bernard Haasdonk 15.5.2007
62 % Markus Dihlmann Feb 2011
64 if model.verbose >= 10
65 disp(['performing simulation for mu = ', mat2str(get_mu(model))]);
68 if (~isequal(model.error_norm,
'l2')) && ...
69 (~isequal(model.error_norm,
'energy')) && ...
70 (~isequal(model.error_norm,
'l2l2'))
71 error(
'error_norm unknown!!');
74 %preparation
for t-partition
75 if ~isfield(model,'transition_model')
76 model.transition_model = 0;
79 if isfield(model,
'starting_time_step')
80 t_ind_start = model.starting_time_step;
85 if isfield(model, 'stopping_time_step')
86 t_ind_stop = model.stopping_time_step;
88 t_ind_stop = model.nt;
92 %cgrid = grid_geometry(params);
93 if ~model.transition_model
94 nRB = length(reduced_data.a0{1});
96 nRB = size(reduced_data.LL_E{1},2);
98 a = zeros(nRB,t_ind_stop-t_ind_start+1); % matrix of RB-coefficients
99 Delta = zeros(1,t_ind_stop-t_ind_start+1); % vector of error estimators
100 model.dt = model.T/model.nt;
103 % initial data projection and linear combination
104 model.decomp_mode = 2;
105 sa0 = rb_init_values(model,[]);
106 a0 = lincomb_sequence(reduced_data.a0,sa0);
110 if isfield(reduced_data,
'Delta0')
111 Delta(1) = reduced_data.Delta0;
112 %disp(['initial error is ',num2str(Delta(1))]);
115 % loop over time steps: computation of a(t+1)==a^t
116 LL_I_is_speye = -1; % unknown state here
117 for t = (t_ind_start+1):t_ind_stop
118 model.t = (t-1)*model.dt;
120 disp(['entered time-loop step ',num2str(t)]);
123 % linear combination of all quantities
124 % only once, if data is const in time
126 if (t==t_ind_start+1) || (~model.data_const_in_time)
127 [sLL_I, sLL_E, sbb, sK_II, sK_IE, sK_EE, sm_I, sm_E, sm] = ...
128 rb_operators(model,[]);
129 LL_I = lincomb_sequence(reduced_data.LL_I,sLL_I);
130 LL_E = lincomb_sequence(reduced_data.LL_E,sLL_E);
131 bb = lincomb_sequence(reduced_data.bb,sbb);
132 K_II = lincomb_sequence(reduced_data.K_II,sK_II);
133 K_IE = lincomb_sequence(reduced_data.K_IE,sK_IE);
134 K_EE = lincomb_sequence(reduced_data.K_EE,sK_EE);
135 m_I = lincomb_sequence(reduced_data.m_I,sm_I);
136 m_E = lincomb_sequence(reduced_data.m_E,sm_E);
137 m = lincomb_sequence(reduced_data.m,sm);
138 if model.data_const_in_time
139 % inverse of implicit matrix! L_I' = Id, L_E' = L_I^(-1) *
140 % L_E, bb = L_I^(-1) bb
144 LL_I = speye(size(LL_I));
146 elseif max(max(LL_I- speye(size(LL_I))))<1e-12
151 % solve LL_I* a(:,t+1) = LL_E * a(:,t) + bb
152 rhs = LL_E * a(:,t-t_ind_start) + bb;
153 % warning('this does not work for pure space operators (independent of delta t). In that case use the commented line below.');
154 %rhs = a(:,t) + model.dt * (LL_E * a(:,t) + bb);
155 % check whether pure explicit problem:
156 if (LL_I_is_speye == 1)
157 a(:,t+1-t_ind_start) = rhs;
158 elseif (LL_I_is_speye == -1) && isequal(LL_I, speye(size(LL_I)))
159 a(:,t+1-t_ind_start) = rhs;
160 else % solve linear system
161 % disp('check symmetry and choose solver accordingly!');
163 % nonsymmetric solvers:
164 % [U(:,t+1), flag] = bicgstab(Li,rhs,[],1000);
165 % [U(:,t+1), flag] = cgs(Li,rhs,[],1000);
166 % symmetric solver, non pd:
167 % see bug_symmlq for a very strange bug: cannot solve identity system!
168 % [U(:,t+1), flag] = symmlq(Li,rhs,[],1000);
169 % [U(:,t+1), flag] = minres(Li,rhs,[],1000);
170 % symmetric solver, pd:
171 % [a(:,t+1), flag] = pcg(LL_I,rhs,[],1000);
172 [a(:,t+1-t_ind_start), flag] = bicgstab(LL_I,rhs,[],1000);
174 disp(['error in system solution, solver return flag = ', ...
180 % compute error estimate recursively
181 % Delta^k := |L_I^(k-1)|^-1 * (res_norm(k-1)+
182 % |L_E^(k-1)|*Delta^(k-1))
186 a(:,t+1-t_ind_start)' * K_II * a(:,t+1-t_ind_start) ...
187 -2* a(:,t+1-t_ind_start)' * K_IE * a(:,t-t_ind_start) ...
188 + a(:,t-t_ind_start)' * K_EE * a(:,t-t_ind_start) ...
190 - 2* a(:,t+1-t_ind_start)' * m_I ...
191 + 2* a(:,t-t_ind_start)' * m_E;
192 % ensure real estimators (could be complex due to numerics)
194 res_norm = sqrt(res_norm_sqr);
199 if isequal(model.error_norm,'l2')
200 Delta(t+1-t_ind_start) = model.L_I_inv_norm_bound * ...
202 model.L_E_norm_bound * Delta(t-t_ind_start));
204 % only sum up the squares of the residuals, multiplication with coeff
205 % and squareroots are taken at the end.
206 Rkplus1sqr = max(res_norm_sqr/(model.dt.^2), 0);
207 Delta(t+1-t_ind_start) = Delta(t-t_ind_start)+Rkplus1sqr;
211 if isequal(model.error_norm,'energy')
212 % now perform scaling and squareroot of energy-error-estimate:
213 alpha = model.coercivity_bound_ptr(model);
214 if max(model.L_E_norm_bound>1) || (alpha == 0)
215 % no reasonable estimation possible
218 C = (sqrt(1-model.L_E_norm_bound.^2)+1) * 0.5;
219 Coeff = model.dt /(4*alpha*C*(1-model.energy_norm_gamma * C));
220 Delta = sqrt(Delta * Coeff);
224 if isequal(model.error_norm,'l2l2')
225 % Delta_space_time(mu) = 1/beta * sqrt(sum delta t * |Rk|^2);
226 % now perform scaling and squareroot of error-estimate:
227 beta = model.inf_sup_beta_lower_bound(model);
228 Delta = sqrt(model.dt^3 * Delta)/beta;
231 % return simulation result
232 simulation_data = [];
233 simulation_data.a = a;
234 simulation_data.Delta = Delta;
235 simulation_data.LL_I = LL_I;
236 simulation_data.LL_E = LL_E;
238 if isfield(model,'name_output_functional')
239 if isequal(model.error_norm,'l2')
240 s = a' * reduced_data.s_RB;
241 Delta_s = reduced_data.s_l2norm * Delta;
242 simulation_data.s = s(:);
243 simulation_data.Delta_s = Delta_s(:);
245 error(['output estimation not implemented concerning energy-norm.' ...
246 ' choose l2 as error_norm!']);