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
5 % This
function computes a numerical scheme defined by the descr
for the
6 % parameter set via the set_mu() method.
9 % sim_data: structure holding the `H`-dimensional simulation data.
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)
20 % Required fields of dmodel:
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
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
42 % Bernard Haasdonk 27.8.2009
44 if dmodel.verbose >= 9
45 disp(
'entered detailed simulation ');
48 if isempty(model_data)
49 model_data = gen_model_data(dmodel);
52 dmodel.decomp_mode = 0; % == complete;
56 model.dt = model.T/model.nt;
57 %Fixing values for t-partition
58 if ~isfield(model,'starting_time_step')
61 t_ind_stop = model.nt;
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?
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;
79 save_time_index = ones(1,model.nt+1);
80 save_time_indices = 1:(model.nt+1);
83 model.nt = t_ind_stop-t_ind_start+1;
85 % initial values by midpoint evaluation
86 Ut = model.init_values_algorithm(model,model_data);
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));
93 operators_required = 1;
94 %if ~isfield(model,'data_const_in_time')
95 % model.data_const_in_time = 0;
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;
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
111 if save_time_index(t_ind+1)
113 time_indices(t_column) = t_ind;
115 t_column = t_column + 1;
118 if ~isfield(model,'compute_output_functional')
119 model.compute_output_functional = 0;
122 %model.plot(U0,model_data.grid,params);
125 for t_ind = (t_ind_start):(t_ind_stop-1)
126 t = (t_ind)*model.dt;
128 disp(['entered time-loop step ',num2str(t_ind)]);
133 % get matrices and bias-vector
135 if operators_required
139 % new assembly in every iteration
140 if ~model.affinely_decomposed
141 [L_I, L_E, b] = model.operators_ptr(model, model_data);
143 % or simple assembly based on affine parameter decomposition
144 % => turns out to be slower for few components
146 % test affine decomposition
148 disp('test affine decomposition of operators')
149 test_affine_decomp(model.operators_ptr,3,1,model,model_data);
152 model.decomp_mode = 2;
153 [L_I_coeff, L_E_coeff, b_coeff] = model.operators_ptr(model, ...
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);
161 if model.data_const_in_time
162 operators_required = 0;
168 if isequal(L_I, speye(size(L_I)))
171 % solve linear system
172 % disp('check symmetry and choose solver accordingly!');
174 % nonsymmetric solvers:
175 % [U(:,t+1), flag] = bicgstab(L_I,rhs,[],1000);
176 % [U(:,t+1), flag] = cgs(L_I,rhs,[],1000);
178 % symmetric solver, non pd:
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);
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);
192 % disp(['error in system solution, solver return flag = ', ...
199 if save_time_index(t_ind+2)
201 time_indices(t_column) = t_ind;
203 t_column = t_column + 1;
213 sim_data.time = time;
214 sim_data.time_indices = time_indices;
215 sim_data.mu = get_mu(dmodel);
217 if model.compute_output_functional
219 v = model.operators_output(model,model_data);
220 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...