rbmatlab  1.16.09
 All Classes Namespaces Files Functions Variables Modules Pages
detailed_simulation.m
1 function sim_data = detailed_simulation(dmodel,model_data)
2 % function sim_data = detailed_simulation(this, model_data);
3 % executes a detailed simulation for a given parameter
4 %
5 % This function computes a numerical scheme defined by the descr for the
6 % parameter set via the set_mu() method.
7 %
8 % Return values:
9 % sim_data: structure holding the `H`-dimensional simulation data.
10 %
11 % Generated fields of sim_data:
12 % U: Dof vectors for the solution snapshots `u_h(\cdot, t^{k_i})`
13 % where `(k_i)_{i=1}^{\hat{K}}` are the 'time_indices'.
14 % time_indices: This is a vector of stored time indices
15 % `(k_i)_{i=1}^{\hat{K}}`
16 % time: This is a vector of stored time instances
17 % `(t^{k_i})_{i=1}^{\hat{K}}`
18 % y: Vector of output functional values (if activated)
19 
20 % Required fields of dmodel:
21 % T : final time
22 % nt : number of time-intervals until T, i.e. nt+1
23 % solution slices are computed
24 % init_values_algorithm: name of function for computing the
25 % initvalues-DOF with arguments (grid, params)
26 % example: init_values_cog
27 % operators_algorithm: name of function for computing the
28 % L_E,L_I,b-operators with arguments (grid, params)
29 % example operators_conv_diff
30 % data_const_in_time : if this optional field is 1, the time
31 % evolution is performed with constant operators,
32 % i.e. only the initial-time-operators are computed
33 % and used throughout the time simulation.
34 % compute_output_functional: flag indicating, whether output
35 % functional is to be computed
36 %
37 % Optional fields of dmodel:
38 % starting_time_step: starting time step for simulation (for example in
39 % t-partition). Default value is 0.
40 % stopping_tim_step: stopping time step for simulation
41 %
42 % Bernard Haasdonk 27.8.2009
43 
44 if dmodel.verbose >= 9
45  disp('entered detailed simulation ');
46 end;
47 
48 if isempty(model_data)
49  model_data = gen_model_data(dmodel);
50 end;
51 
52 dmodel.decomp_mode = 0; % == complete;
53 
54 model = dmodel.descr;
55 
56 model.dt = model.T/model.nt;
57 %Fixing values for t-partition
58 if ~isfield(model,'starting_time_step')
59  model.t = 0;
60  t_ind_start =0;
61  t_ind_stop = model.nt;
62 else
63  t_ind_start = model.starting_time_step;
64  t_ind_stop = model.stopping_time_step;
65  model.t = (t_ind_start)*model.dt;
66  %model.nt = t_ind_stop-t_ind_start+1; %+1?
67 end
68 
69 
70 
71 if isfield(model,'save_time_indices')
72  valid_save_time_indices1 = model.save_time_indices>=t_ind_start;
73  valid_save_time_indices2 = model.save_time_indices<=t_ind_stop;
74  valid_save_time_indices = logical(valid_save_time_indices1.*valid_save_time_indices2);
75  save_time_index = zeros(1,model.nt+1);
76  save_time_index(model.save_time_indices(valid_save_time_indices)+1) = 1;
77  save_time_indices = find(save_time_index==1)-1;
78 else
79  save_time_index = ones(1,model.nt+1);
80  save_time_indices = 1:(model.nt+1);
81 end;
82 
83 model.nt = t_ind_stop-t_ind_start+1;
84 
85 % initial values by midpoint evaluation
86 Ut = model.init_values_algorithm(model,model_data);
87 
88 % preallocate U
89 U = zeros(length(Ut),length(save_time_indices));
90 time_indices = zeros(1,length(save_time_indices));
91 time = zeros(1,length(save_time_indices));
92 
93 operators_required = 1;
94 %if ~isfield(model,'data_const_in_time')
95 % model.data_const_in_time = 0;
96 %end,
97 
98 % get operator components
99 if model.affinely_decomposed
100  old_decomp_mode = model.decomp_mode;
101  model.decomp_mode = 1;
102  [L_I_comp, L_E_comp, b_comp] = model.operators_ptr(model, model_data);
103  model.decomp_mode = old_decomp_mode;
104 end;
105 
106 t_column = 1; % next column index to be filled in output
107 t_ind = t_ind_start; % t_ind between 0 and nt
108 t = model.t; % absolute time between 0 and T
109 
110 % store init data
111 if save_time_index(t_ind+1)
112  U(:,t_column) = Ut;
113  time_indices(t_column) = t_ind;
114  time(t_column) = t;
115  t_column = t_column + 1;
116 end;
117 
118 if ~isfield(model,'compute_output_functional')
119  model.compute_output_functional = 0;
120 end,
121 
122 %model.plot(U0,model_data.grid,params);
123 
124 % loop over fv-steps
125 for t_ind = (t_ind_start):(t_ind_stop-1)
126  t = (t_ind)*model.dt;
127  if model.verbose >= 10
128  disp(['entered time-loop step ',num2str(t_ind)]);
129  end;
130  if model.verbose == 9
131  fprintf('.');
132  end;
133  % get matrices and bias-vector
134 
135  if operators_required
136 
137  model.t = t;
138 
139  % new assembly in every iteration
140  if ~model.affinely_decomposed
141  [L_I, L_E, b] = model.operators_ptr(model, model_data);
142  else
143  % or simple assembly based on affine parameter decomposition
144  % => turns out to be slower for few components
145 
146  % test affine decomposition
147  if model.debug
148  disp('test affine decomposition of operators')
149  test_affine_decomp(model.operators_ptr,3,1,model,model_data);
150  disp('OK')
151  end;
152  model.decomp_mode = 2;
153  [L_I_coeff, L_E_coeff, b_coeff] = model.operators_ptr(model, ...
154  model_data);
155  model.decomp_mode = old_decomp_mode;
156  L_I = lincomb_sequence(L_I_comp, L_I_coeff);
157  L_E = lincomb_sequence(L_E_comp, L_E_coeff);
158  b = lincomb_sequence(b_comp, b_coeff);
159  end;
160 
161  if model.data_const_in_time
162  operators_required = 0;
163  end;
164 
165  end;
166 
167  rhs = L_E * Ut + b;
168  if isequal(L_I, speye(size(L_I)))
169  Ut = rhs;
170  else
171  % solve linear system
172  % disp('check symmetry and choose solver accordingly!');
173  % keyboard;
174  % nonsymmetric solvers:
175  % [U(:,t+1), flag] = bicgstab(L_I,rhs,[],1000);
176  % [U(:,t+1), flag] = cgs(L_I,rhs,[],1000);
177  %
178  % symmetric solver, non pd:
179  % omit symmlq:
180  % see bug_symmlq.mat for a very strange bug: cannot solve identity system!
181  % reported to matlab central, but no solution up to now.
182  % [U(:,t+1), flag] = symmlq(L_I,rhs,[],1000);
183  %
184  % [U(:,t+1), flag] = minres(L_I,rhs,[],1000);
185  % symmetric solver, pd:
186  %[U(:,t+1), flag] = pcg(L_I,rhs,[],1000);
187  % bicgstab works also quite well:
188  %[U(:,t+1), flag] = bicgstab(L_I,rhs,[],1000);
189  Ut = L_I \ rhs;
190 
191  % if flag>0
192  % disp(['error in system solution, solver return flag = ', ...
193  % num2str(flag)]);
194  % keyboard;
195  % end;
196 
197  end;
198 
199  if save_time_index(t_ind+2)
200  U(:,t_column) = Ut;
201  time_indices(t_column) = t_ind;
202  time(t_column) = t;
203  t_column = t_column + 1;
204  end;
205 
206 end;
207 
208 if model.verbose == 9
209  fprintf('\n');
210 end;
211 
212 sim_data.U = U;
213 sim_data.time = time;
214 sim_data.time_indices = time_indices;
215 sim_data.mu = get_mu(dmodel);
216 
217 if model.compute_output_functional
218  % get operators
219  v = model.operators_output(model,model_data);
220  sim_data.y = (v(:)') * U;
221 end;
222 
function r = verbose(level, message, messageId)
This function displays messages depending on a message-id and/or a level. Aditionally you can set/res...
Definition: verbose.m:17