1 function sim_data = lin_evol_detailed_simulation(model,model_data)
2 %
function sim_data = lin_evol_detailed_simulation(model,model_data)
4 % required fields of model:
6 % nt : number of time-intervals until T, i.e. nt+1
7 % solution slices are computed
8 % init_values_algorithm: name of
function for computing the
9 % initvalues-DOF with arguments (grid, params)
10 % example: init_values_cog
11 % operators_algorithm: name of
function for computing the
12 % L_E,L_I,b-operators with arguments (grid, params)
13 % example operators_conv_diff
14 % data_const_in_time :
if this optional field is 1, the time
15 % evolution is performed with constant operators,
16 % i.e. only the initial-time-operators are computed
17 % and used throughout the time simulation.
18 % compute_output_functional: flag indicating, whether output
19 % functional is to be computed
21 %
return fields of sim_data:
22 % U : sequence of DOF vectors
23 % y : sequence of output functional values (
if activated)
25 % optional fields of model:
26 % starting_time_step: starting time step
for simulation (
for example in
27 % t-partition). Default value is 0.
28 % stopping_tim_step: stopping time step
for simulation
30 % Bernard Haasdonk 27.8.2009
33 disp(
'entered detailed simulation ');
36 if isempty(model_data)
37 model_data = gen_model_data(model);
40 model.decomp_mode = 0; % == complete;
42 model.dt = model.T/model.nt;
43 %Fixing values for t-partition
44 if ~isfield(model,'starting_time_step')
47 t_ind_stop = model.nt;
49 t_ind_start = model.starting_time_step;
50 t_ind_stop = model.stopping_time_step;
51 model.t = (t_ind_start)*model.dt;
52 %model.nt = t_ind_stop-t_ind_start+1; %+1?
57 if isfield(model,'save_time_indices')
58 valid_save_time_indices1 = model.save_time_indices>=t_ind_start;
59 valid_save_time_indices2 = model.save_time_indices<=t_ind_stop;
60 valid_save_time_indices = valid_save_time_indices1.*valid_save_time_indices2;
61 valid_save_time_indices = find(valid_save_time_indices);
62 save_time_index = zeros(1,model.nt+1);
63 save_time_index(model.save_time_indices(valid_save_time_indices)+1) = 1;
64 save_time_indices = find(save_time_index==1)-1;
66 save_time_index = ones(1,model.nt+1);
67 save_time_indices = 1:(model.nt+1);
70 model.nt = t_ind_stop-t_ind_start+1;
72 % initial values by midpoint evaluation
73 Ut = model.init_values_algorithm(model,model_data);
76 U = zeros(length(Ut),length(save_time_indices));
77 time_indices = zeros(1,length(save_time_indices));
78 time = zeros(1,length(save_time_indices));
80 operators_required = 1;
81 %if ~isfield(model,'data_const_in_time')
82 % model.data_const_in_time = 0;
85 % get operator components
86 if model.affinely_decomposed
87 old_decomp_mode = model.decomp_mode;
88 model.decomp_mode = 1;
89 [L_I_comp, L_E_comp, b_comp] = model.operators_ptr(model, model_data);
90 model.decomp_mode = old_decomp_mode;
93 t_column = 1; % next column index to be filled in output
94 t_ind = t_ind_start; % t_ind between 0 and nt
95 t = model.t; % absolute time between 0 and T
98 if save_time_index(t_ind+1)
100 time_indices(t_column) = t_ind;
102 t_column = t_column + 1;
105 if ~isfield(model,'compute_output_functional')
106 model.compute_output_functional = 0;
109 %model.plot(U0,model_data.grid,params);
112 for t_ind = (t_ind_start):(t_ind_stop-1)
113 t = (t_ind)*model.dt;
115 disp(['entered time-loop step ',num2str(t_ind)]);
120 % get matrices and bias-vector
122 if operators_required
126 % new assembly in every iteration
127 if ~model.affinely_decomposed
128 [L_I, L_E, b] = model.operators_ptr(model, model_data);
130 % or simple assembly based on affine parameter decomposition
131 % => turns out to be slower for few components
133 % test affine decomposition
135 disp('test affine decomposition of operators')
136 test_affine_decomp(model.operators_ptr,3,1,model,model_data);
139 model.decomp_mode = 2;
140 [L_I_coeff, L_E_coeff, b_coeff] = model.operators_ptr(model, ...
142 model.decomp_mode = old_decomp_mode;
143 L_I = lincomb_sequence(L_I_comp, L_I_coeff);
144 L_E = lincomb_sequence(L_E_comp, L_E_coeff);
146 b = lincomb_sequence(b_comp, b_coeff);
148 b = zeros(size(L_I,1),1);
152 if model.data_const_in_time
153 operators_required = 0;
159 if isequal(L_I, speye(size(L_I)))
162 % solve linear system
163 % disp('check symmetry and choose solver accordingly!');
165 % nonsymmetric solvers:
166 % [U(:,t+1), flag] = bicgstab(L_I,rhs,[],1000);
167 % [U(:,t+1), flag] = cgs(L_I,rhs,[],1000);
169 % symmetric solver, non pd:
171 % see bug_symmlq.mat for a very strange bug: cannot solve identity system!
172 % reported to matlab central, but no solution up to now.
173 % [U(:,t+1), flag] = symmlq(L_I,rhs,[],1000);
175 % [U(:,t+1), flag] = minres(L_I,rhs,[],1000);
176 % symmetric solver, pd:
177 %[U(:,t+1), flag] = pcg(L_I,rhs,[],1000);
178 % bicgstab works also quite well:
179 %[U(:,t+1), flag] = bicgstab(L_I,rhs,[],1000);
183 % disp(['error in system solution, solver return flag = ', ...
190 if save_time_index(t_ind+2)
192 time_indices(t_column) = t_ind;
194 t_column = t_column + 1;
204 sim_data.time = time;
205 sim_data.time_indices = time_indices;
207 if model.compute_output_functional
209 v = model.operators_output(model,model_data);
210 sim_data.y = (v(:)') * U;
function r = verbose(level, message, messageId)
This function displays messages depending on a message-id and/or a level. Aditionally you can set/res...