1 function sim_data = lin_evol_opt_detailed_der_simulation_t_part(model,model_data)
2 %
function sim_data = lin_evol_opt_detailed_der_simulation(model,model_data)
4 % Function performing a detailed simulation and calculating the derivative
5 % solution with respect to all parameters.
6 % required fields of model:
8 % nt : number of time-intervals until T, i.e. nt+1
9 % solution slices are computed
10 % init_values_algorithm: name of
function for computing the
11 % initvalues-DOF with arguments (grid, params)
12 % example: init_values_cog
13 % operators_algorithm: name of
function for computing the
14 % L_E,L_I,b-operators with arguments (grid, params)
15 % example operators_conv_diff
16 % data_const_in_time :
if this optional field is 1, the time
17 % evolution is performed with constant operators,
18 % i.e. only the initial-time-operators are computed
19 % and used throughout the time simulation.
20 % compute_output_functional: flag indicating, whether output
21 % functional is to be computed
23 % compute_derivative_info: flag
if derivative information shall be
26 % optional fields of model:
27 % starting_time_step: starting time step
for simulation (
for example in
28 % t-partition). Default value is 0.
29 % stopping_tim_step: stopping time step
for simulation
31 %
return fields of sim_data:
32 % U : sequence of DOF vectors
33 % y : sequence of output functional values (
if activated)
35 % Markus Dihlmann 15.10.10
38 disp(
'entered detailed simulation ');
41 if isempty(model_data)
42 model_data = gen_model_data(model);
45 if ~isfield(model,'compute_derivative_info')
46 model.compute_derivative_info = 0;
50 model.decomp_mode = 1; % == components;
51 model.dt = model.T/model.nt;
52 %Fixing values for t-partition
53 if ~isfield(model,'starting_time_step')
56 t_ind_stop = model.nt;
58 t_ind_start = model.starting_time_step;
59 t_ind_stop = model.stopping_time_step;
60 model.t = (t_ind_start)*model.dt;
61 %model.nt = t_ind_stop-t_ind_start+1; %+1?
65 if ~isfield(model,'optimization')
66 disp('Model has no field optimization.')
67 model.optimization =[];
71 if isfield(model.optimization,'params_to_optimize')
72 indx_mu_i = find(model.optimization.params_to_optimize); %indexes i of \mu_i for which derivative information is needed
74 indx_mu_i = 1:length(get_mu(model));
76 nr_indx_mu_i= length(indx_mu_i); %nr of parameters to be optimized
77 %nr_derivatives = size(model.optimization.params_to_optimize,2); %total nr of parameters in the model
81 if isfield(model,'save_time_indices')
82 valid_save_time_indices1 = model.save_time_indices>=t_ind_start;
83 valid_save_time_indices2 = model.save_time_indices<=t_ind_stop;
84 valid_save_time_indices = valid_save_time_indices1.*valid_save_time_indices2;
85 valid_save_time_indices = find(valid_save_time_indices);
86 save_time_index = zeros(1,model.nt+1);
87 save_time_index(model.save_time_indices(valid_save_time_indices)+1) = 1;
88 save_time_index(t_ind_start+1) = 1;
89 save_time_index(t_ind_stop+1) = 1;
90 save_time_indices = find(save_time_index==1)-1;
92 save_time_index = ones(1,model.nt+1);
93 save_time_indices = 1:(model.nt+1);
96 model.nt = t_ind_stop-t_ind_start+1;
98 %fixing initial values
100 % initial values by midpoint evaluation
101 U0_comp = model.init_values_algorithm(model,model_data);
103 model.decomp_mode=2; %==> coefficients
104 U0_coeff = model.init_values_algorithm(model, model_data);
105 if model.compute_derivative_info
106 U0_der_coeff = model.rb_derivative_init_values_coefficients(model);
109 %assembling initial values true solution
110 Ut = lincomb_sequence(U0_comp,U0_coeff);
112 %assembling initial values derivative solution
113 if model.compute_derivative_info
114 Ut_der = zeros(length(Ut),nr_indx_mu_i);
116 Ut_der(:,i) = lincomb_sequence(U0_comp, U0_der_coeff(indx_mu_i(i),:)'); %%% FRAGE: PASST DAS MIT Ut_der? Hat 3 dim statt 2?!
119 U_der=cell(nr_indx_mu_i,1);
121 U_der{i} = sparse(length(Ut), length(save_time_indices));
126 %take values from t-part foregoing solution
127 Ut= model.init_values_algorithm(model, model_data);
128 if model.compute_derivative_info
129 Ut_der = model.derivative_init_values_algorithm(model, model_data);
133 U = sparse(length(Ut),length(save_time_indices));
135 time_indices = zeros(1,length(save_time_indices));
136 time = zeros(1,length(save_time_indices));
138 %operators_required = 1;
139 if ~isfield(model,
'data_const_in_time')
140 model.data_const_in_time = 0;
143 % get operator components
144 if model.affinely_decomposed
145 old_decomp_mode = model.decomp_mode;
146 model.decomp_mode = 1;
147 [L_I_comp, L_E_comp, b_comp] = model.operators_ptr(model, model_data);
148 model.decomp_mode = old_decomp_mode;
151 t_column = 1; % next column index to be filled in output
152 t_ind = t_ind_start; % t_ind between 0 and nt
153 t = model.t; % absolute time between 0 and T
156 if save_time_index(t_ind+1)
158 if model.compute_derivative_info
160 U_der{i}(:,t_column) = Ut_der(:,i);
163 t_column = t_column + 1;
164 time_indices(t_column) = t_ind;
166 t_column = t_column + 1;
169 if ~isfield(model,
'compute_output_functional')
170 model.compute_output_functional = 0;
173 %model.plot(U0,model_data.grid,params);
175 %allocate memory for derivation matrices:
176 L_E_derived = cell(nr_indx_mu_i,1);
177 L_I_derived = cell(nr_indx_mu_i,1);
178 b_derived = cell(nr_indx_mu_i,1);
181 for t_ind = (t_ind_start):(t_ind_stop-1)
182 t = (t_ind)*model.dt;
184 disp(['entered time-loop step ',num2str(t_ind)]);
189 % get matrices and bias-vector
191 if (t_ind==t_ind_start)||(~model.data_const_in_time)%operators_required
195 % new assembly in every iteration
196 if ~model.affinely_decomposed
197 [L_I, L_E, b] = model.operators_ptr(model, model_data);
199 % or simple assembly based on affine parameter decomposition
200 % => turns out to be slower for few components
202 % test affine decomposition
204 disp('test affine decomposition of operators')
205 test_affine_decomp(model.operators_ptr,3,1,model,model_data);
208 model.decomp_mode = 2;
209 [L_I_coeff, L_E_coeff, b_coeff] = model.operators_ptr(model, ...
211 model.decomp_mode = old_decomp_mode;
212 L_I = lincomb_sequence(L_I_comp, L_I_coeff);
213 L_E = lincomb_sequence(L_E_comp, L_E_coeff);
214 b = lincomb_sequence(b_comp, b_coeff);
217 %if model.data_const_in_time
218 % operators_required = 0;
221 if model.compute_derivative_info
222 % Get Operators for derivative information:
223 model.decomp_mode = 2;
224 [sLi,sLe,sb] = model.rb_derivative_operators_coefficients_ptr(model);
227 if ~isempty(L_I_comp) %only do the linear combination if the coefficients sLi are not empty
230 L_I_derived{i}=sparse(size(L_I_comp{1},1),size(L_I_comp{1},2));
236 L_E_derived{i} = sparse(size(L_E_comp{1},1),size(L_E_comp{1},2));
237 b_derived{i} = sparse(length(b_comp{1}),1);
241 if ~isempty(sLi)%~isequal(sLi,[]) %only
do the linear combination
if the coefficients sLi are not empty
243 L_I_derived{i} = lincomb_sequence(L_I_comp,sLi(indx_mu_i(i)));
251 L_E_derived{i} = lincomb_sequence(L_E_comp,sLe(indx_mu_i(i),:));
252 b_derived{i} = lincomb_sequence(b_comp,sb(indx_mu_i(i),:));
260 rhs = (speye(size(L_E)) + model.dt*L_E) * Ut + model.dt*b;
263 if isempty(find(L_I, 1))%isequal(L_I, speye(size(L_I)))
266 % solve linear system
267 % disp(
'check symmetry and choose solver accordingly!');
269 % nonsymmetric solvers:
270 % [U(:,t+1), flag] = bicgstab(L_I,rhs,[],1000);
271 % [U(:,t+1), flag] = cgs(L_I,rhs,[],1000);
273 % symmetric solver, non pd:
275 % see bug_symmlq.mat
for a very strange bug: cannot solve identity system!
276 % reported to matlab central, but no solution up to now.
277 % [U(:,t+1), flag] = symmlq(L_I,rhs,[],1000);
279 % [U(:,t+1), flag] = minres(L_I,rhs,[],1000);
280 % symmetric solver, pd:
281 %[U(:,t+1), flag] = pcg(L_I,rhs,[],1000);
282 % bicgstab works also quite well:
283 %[U(:,t+1), flag] = bicgstab(L_I,rhs,[],1000);
284 Ut = (speye(size(L_I))-model.dt*L_I)\rhs;%Ut + model.dt*(L_I \ rhs);
287 % disp([
'error in system solution, solver return flag = ', ...
295 if save_time_index(t_ind+2)
296 U(:,t_column-1) = Ut;
297 time_indices(t_column-1) = t_ind+1;
298 time(t_column-1) = t+model.dt;
299 t_column = t_column + 1;
302 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
303 %calculation derivative solution
304 if model.compute_derivative_info
307 rhs_der = (speye(size(L_E)) + model.dt*L_E) * Ut_der(:,i) +...
308 model.dt*L_E_derived{i} * Ut_old + model.dt*L_I_derived{i}*Ut + model.dt*b_derived{i};
309 if isempty(find(L_I))%,zeros(size(L_I)))
310 Ut_der(:,i) = rhs_der;
312 Ut_der(:,i) = (speye(size(L_I))- model.dt*L_I)\rhs_der;
315 if save_time_index(t_ind+2)
316 U_der{i}(:,t_column-2) = Ut_der(:,i);
327 if model.verbose >= 9
333 if model.compute_derivative_info
334 sim_data.U_der = U_der;
337 sim_data.time = time;
338 sim_data.time_indices = time_indices;
343 if model.compute_output_functional
345 if isfield(model,
'get_output')
346 sim_data = model.get_output(model, sim_data,model_data);
349 v = model.operators_output(model,model_data);
350 sim_data.y = (v{1}
') * U;
352 if model.compute_derivative_info
353 for i = 1:nr_indx_mu_i
355 sim_data.y_der(i,:) = (v{1}')*U_der{i};
function r = verbose(level, message, messageId)
This function displays messages depending on a message-id and/or a level. Aditionally you can set/res...