1 function ei_sim_data = nonlin_evol_detailed_ei_simulation(model, detailed_data)
2 %
function ei_sim_data = nonlin_evol_detailed_ei_simulation(model, detailed_data)
4 %
function performing a general time evolution and computing the
5 % solution DOF-vector ei_sim_data.U(:,t)
for all times t=1,...,nt+1
6 % instead of the real
explicit operator, the empirical
7 % interpolation
operator given in detailed_data is used.
8 % So
this function can be used
for testing the stability of the
9 % empirical interpolation
operator. The coefficient vectors of the
10 % empirical interpolation of the space
operator are returned as
11 % columns of ei_sim_data.c
13 % required fields of model:
15 % nt : number of time-intervals until T, i.e. nt+1
16 % solution slices are computed
17 % init_values_algorithm: name of
function for computing the
18 % initvalues-DOF with arguments (grid)
19 % example: init_values_cog
20 % implicit_operators_algorithm: name of
function for computing the
21 % L_I, b_I-operators with arguments (grid)
22 % example fv_operators_implicit
23 % data_const_in_time :
if this optional field is 1, the time
24 % evolution is performed with constant operators,
25 % i.e. only the initial-time-operators are computed
26 % and used throughout the time simulation.
27 % L_E_local_ptr :
function pointer to the local space-discretization
28 %
operator evaluation with syntax
30 % INC_local = L_local(U_local_ext, ind_local,
32 % where *_local indicates results on a set of
33 % elements, *_local_ext includes information including
34 % their neighbours, i.e. the grid must be known
35 % also
for the neighbours, the values of the
36 % previous timstep must be known on the
37 % neighbours. The subset of elements, which are
38 % to be computed are given in ind_local and the
39 % values produced are INC_local.
40 % A time-step can then be realized by
41 % NU_local = U_local_ext(ind_local) - dt * INC_local
42 % M : number of colateral basis vectors, that will be used
for
43 % empirical inteprolation of the
explicit operator
45 % Bernard Haasdonk 24.6.2006
48 error(
'wrong number of parameters!');
52 disp([
'entered detailed simulation with interpolation ',...
53 'of the explicit operator']);
56 if ~iscell(detailed_data.BM)
57 detailed_data.BM = { detailed_data.BM, detailed_data.BM, detailed_data.BM };
58 detailed_data.QM = { detailed_data.QM, detailed_data.QM, detailed_data.QM };
59 detailed_data.TM = { detailed_data.TM, detailed_data.TM, detailed_data.TM };
62 M = min(model.M, max(cellfun(@(X)(size(X,2)), detailed_data.BM)));
64 if model.separate_CRBs
65 Mexpl = min(model.M, size(detailed_data.BM{detailed_data.explicit_crb_index},2));
69 QM = cell(1,length(detailed_data.BM));
70 BM = cell(1,length(detailed_data.BM));
71 TM = cell(1,length(detailed_data.BM));
72 BMinv = cell(1,length(detailed_data.BM));
74 for oi=1:length(detailed_data.BM)
75 M = min(model.M, max(size(detailed_data.BM{oi},2)));
76 % prepare online-data
for online-ei-components (e.g. varying M)
77 QM{oi} = detailed_data.QM{oi}(:,1:M);
78 BM{oi} = detailed_data.BM{oi}(1:M,1:M);
79 TM{oi} = detailed_data.TM{oi}(1:M);
80 BMinv{oi} = inv(BM{oi});
83 % test of the effect of magic-point exchange
84 %TM = TM(randperm(model.M));
86 model.decomp_mode = 0;
88 % initial values by midpoint evaluation
89 U0 = model.init_values_algorithm(model, detailed_data);
91 U = zeros(length(U0(:)),model.nt+1);
94 c = zeros(Mexpl,model.nt);
96 model.dt = model.T/model.nt;
97 operators_required = 1;
98 if isempty(model.data_const_in_time
')
99 model.data_const_in_time = 0;
102 impl_CRB_ind = detailed_data.implicit_crb_index;
103 expl_CRB_ind = detailed_data.explicit_crb_index;
107 if model.verbose > 19
108 disp(['entered time-loop step
',num2str(t)]);
109 elseif model.verbose > 9
112 % get matrices and bias-vector
113 if operators_required
114 model.t = (t-1)*model.dt;
116 if all([model.fv_impl_conv_weight, model.fv_impl_diff_weight, model.fv_impl_react_weight] == [0, 0, 0])
117 % purely explicit scheme
119 tmpmodel.decomp_mode = 0;
121 model.implicit_operators_algorithm(tmpmodel, detailed_data, []);
122 elseif ~model.newton_solver % implicit scheme without newton steps
124 tmpmodel.decomp_mode = 0;
126 [L_I_red, b_I_red] = ...
127 model.implicit_operators_algorithm(tmpmodel, detailed_data, []);
129 LL_I = QM{impl_CRB_ind} * BMinv{impl_CRB_ind} ...
130 * L_I_red(TM{impl_CRB_ind},:);
131 bb_I = QM{impl_CRB_ind} * BMinv{impl_CRB_ind} ...
132 * b_I_red(TM{impl_CRB_ind});
135 if model.data_const_in_time
136 operators_required = 0;
141 % compute real explicit operator result, but only use the magic-point entries
142 % and perform an empirical interpolation instead
143 inc = model.L_E_local_ptr(model, detailed_data, U(:,t), []);
145 inc_local = inc(TM{expl_CRB_ind});
146 %ei_coefficients = BM \ inc_local;
147 %ei_inc = detailed_data.QM(:,1:model.M)*ei_coefficients;
149 c(:,t) = BM{expl_CRB_ind} \ inc_local;
151 ei_inc = QM{expl_CRB_ind} * (BM{expl_CRB_ind} \ inc_local);
153 if model.newton_solver
154 % dTM = detailed_data.TM{impl_CRB_ind};
155 % mask = zeros(1, detailed_data.grid.nelements);
156 % nbi = detailed_data.grid.NBI(dTM,:);
158 % % get indices of TM's neighbour-neighbours
159 % nnbi = detailed_data.grid.NBI(unique(nbi(i)),:);
160 % ni = find(nnbi > 0);
161 % [nnbuniq,tally] = unique(sort(nnbi(ni)),
'first');
162 % tally = [tally(2:end);length(nnbi(ni))+1] - tally;
163 % mask(nnbuniq) = tally;
166 % mask(mask < 2) = 0;
169 % eind_speye = sparse(1:length(eind), eind, ones(1,length(eind)), ...
171 % detailed_data.grid.nelements);
173 % reduced_data.grid = gridpart(detailed_data.grid,eind);
175 % glob2loc = zeros(detailed_data.grid.nelements,1);
176 % glob2loc(eind) = 1:length(eind);
177 % reduced_data.TM_local = glob2loc(dTM);
179 updelta = zeros(length(U0),model.newton_steps+1);
180 exprhs = model.dt * inc;
181 updelta(:,1) = U(:,t) - exprhs;
188 for n=1:model.newton_steps;
189 if model.verbose > 10
190 disp([
'Computing time step no. ', num2str(t), ...
191 ', Newton step no. ', num2str(n), ...
192 '... Residual norm: ']);
194 model.t = model.t + model.dt;
195 Urhs = model.L_I_local_ptr(model, detailed_data, updelta(:,n), []);
197 % project result of L_I onto CRB
198 Urhs_local = Urhs(TM{impl_CRB_ind});
199 Urhs_ei_coeffs = BM{impl_CRB_ind} \ Urhs_local;
200 Urhs_ei = QM{impl_CRB_ind}(:,1:model.M) * Urhs_ei_coeffs;
202 %
get derivative DL_I (interpolation is done only implicitly)
203 [gradL_I_local, gL_I_offset_local] = ...
204 model.implicit_gradient_operators_algorithm(model, detailed_data, ...
207 model.t = model.t - model.dt;
208 % gL_I_offset_ei = QM{impl_CRB_ind} * (BM{impl_CRB_ind} \ ...
209 % gL_I_offset_local(reduced_data.TM_local));
211 % gradL_I_ei = QM{impl_CRB_ind} * (BM{impl_CRB_ind} \ gradL_I_local);
212 gradL_I_ei = gradL_I_local;
214 gradL_I = gradL_I_ei * model.dt + speye(size(gradL_I_ei));
216 rhs = - updelta(:,n) + U(:,t) - exprhs - model.dt * Urhs_ei;
218 delta = gradL_I \ rhs;
220 % project delta on CRB
221 delta_local = delta(TM{impl_CRB_ind});
222 delta_ei_coeffs = BM{impl_CRB_ind} \ delta_local;
223 delta_ei = QM{impl_CRB_ind}(:,1:model.M) * delta_ei_coeffs;
225 residual = sqrt(model.inner_product(model, detailed_data, delta_ei, delta_ei));
226 if model.verbose > 10
227 disp(num2str(residual));
229 updelta(:,n+1) = updelta(:,n) + delta_ei;
230 iltz = find(updelta(:,n+1) < 0);
232 updelta(iltz,n+1) = 0.001;
234 if residual < model.newton_epsilon
239 U(:,t+1) = updelta(:,n+1);
240 Unewton = [ Unewton, updelta(:,1:n+1) ];
242 L_I = LL_I * model.dt + speye(size(LL_I));
244 rhs = U(:,t) - model.dt * ei_inc - model.dt * bb_I;
246 % disp('check rhs and inc!!');
249 if isequal(L_I, speye(size(L_I)))
252 % solve linear system
253 % disp('check symmetry and choose solver accordingly!');
255 % nonsymmetric solvers:
256 % [U(:,t+1), flag] = bicgstab(L_I,rhs,[],1000);
257 % [U(:,t+1), flag] = cgs(L_I,rhs,[],1000);
259 % symmetric solver, non pd:
261 % see bug_symmlq.mat for a very strange bug: cannot solve identity system!
262 % reported to matlab central, but no solution up to now.
263 % [U(:,t+1), flag] = symmlq(L_I,rhs,[],1000);
265 % [U(:,t+1), flag] = minres(L_I,rhs,[],1000);
266 % symmetric solver, pd:
267 U(:,t+1) = L_I \ rhs;
268 % [U(:,t+1), flag] = bicgstab(L_I,rhs,[],1000);
270 disp(['error in system solution, solver return flag = ', ...
276 % plot_element_data(grid,U(:,t));
280 if model.newton_solver && model.
verbose > 0
281 disp(['computation ready: We needed at most ', num2str(nmax), ' Newton steps']);
290 if model.newton_solver
291 ei_sim_data.Unewton = Unewton;
function r = verbose(level, message, messageId)
This function displays messages depending on a message-id and/or a level. Aditionally you can set/res...