rbmatlab  1.16.09
 All Classes Namespaces Files Functions Variables Modules Pages
lin_ds_rb_simulation.m
1 function rb_sim_data = lin_ds_rb_simulation(model,reduced_data)
2 %function rb_sim_data = lin_ds_rb_simulation(model,reduced_data)
3 %
4 % function performing a reduced basis simulation with a
5 % theta-scheme time discretization.
6 
7 % Bernard Haasdonk 2.4.2009
8 
9 if model.verbose >= 9
10  disp('entered reduced simulation ');
11 end;
12 
13 %keyboard;
14 
15 % possibly shrink reduced_data detected by this routine
16 reduced_data_t = lin_ds_reduced_data_subset(model,reduced_data);
17 
18 %model.affine_decomp_mode = 'complete';
19 model.decomp_mode = 2; % coefficients
20 
21 model = model.set_time(model,0);
22 %model.t = 0;
23 
24 if isfield(model,'save_time_indices')
25  save_time_index = zeros(1,model.nt+1);
26  save_time_index(model.save_time_indices(:)+1) = 1;
27  save_time_indices = find(save_time_index==1)-1;
28 else
29  save_time_index = ones(1,model.nt+1);
30  save_time_indices = 1:(model.nt+1);
31 end;
32 
33 x0coeff = model.x0_function_ptr(model,[]);
34 xr = lincomb_sequence(reduced_data_t.x0r,x0coeff);
35 
36 
37 % instead of large dense matrix allocation equivalent:
38 % a' Mtemp b =
39 % a' G b - (a' W) (V' G b) - (a' G V) (W' b) + (a' W) (V' G V) (W' b)
40 %
41 % these components can not be summed up here,
42 % as we must be able to extract online data for other reduced dimensions!!!
43 % hence no longer computation of m here
44 % but decomposition as follows:
45 %
46 % m00: 2d cell array of values x0{q}'G x0{q'}
47 % VtGx0: 1d cell array of of vectors (V' G x0{q})
48 % Wtx0: 1d cell array of of vectors (W' x0{q})
49 % VtGV : matrix V' G V
50 
51 if model.enable_error_estimator
52  %e0sqr = lincomb_sequence2(reduced_data_t.m,x0coeff,x0coeff);
53  x0tGx0 = lincomb_sequence2(reduced_data_t.m00,x0coeff,x0coeff);
54  Wtx0 = lincomb_sequence(reduced_data_t.Wtx0,x0coeff);
55  VtGx0 = lincomb_sequence(reduced_data_t.VtGx0,x0coeff);
56  VtGV = reduced_data_t.VtGV;
57  e0sqr = x0tGx0 - (Wtx0)'*VtGx0 - (VtGx0)' * Wtx0 + (Wtx0)' * VtGV * Wtx0;
58  Rnormsqrs = zeros(1,length(save_time_indices));
59 end;
60 
61 D = model.D_function_ptr(model,[]);
62 Xr = zeros(length(xr(:)),length(save_time_indices));
63 Y = zeros(size(reduced_data_t.Cr{1},1),size(Xr,2));
64 
65 if model.enable_error_estimator
66  DeltaX = zeros(1,length(save_time_indices));
67 end;
68 
69 time_indices = zeros(1,length(save_time_indices));
70 time = zeros(1,length(save_time_indices));
71 
72 t_column = 1; % next column index to be filled in output
73 t_ind = 0; % current t_ind between 0 and nt
74 t = 0; % absolute time between 0 and T
75 
76 if model.enable_error_estimator
77  DeltaXt = model.state_bound_constant_C1*sqrt(max(e0sqr,0));
78 end;
79 
80 if save_time_index(t_ind+1)
81  Xr(:,t_column) = xr(:);
82  time_indices(t_column) = t_ind;
83  time(t_column) = t;
84  %u = model.u_function_ptr(model);
85 
86  Ccoeff = model.C_function_ptr(model,[]);
87  Cr = lincomb_sequence(reduced_data_t.Cr,Ccoeff);
88  u = model.u_function_ptr(model);
89  Y(:,t_column) = Cr * xr + D * u;
90 
91  if model.enable_error_estimator
92  DeltaX(t_column) = DeltaXt;
93  end;
94  t_column = t_column + 1;
95 end;
96 
97 %keyboard;
98 
99 model.dt = model.T/model.nt;
100 dt = model.dt;
101 theta = model.theta;
102 
103 operators_required = 1;
104 if ~isfield(model,'data_const_in_time')
105  model.data_const_in_time = 0;
106 end;
107 previous_implicit_operators_available = 0;
108 
109 % loop over time
110 for t_ind = 1:model.nt
111 % keyboard;
112  if model.verbose >= 10
113  disp(['entered time-loop step ',num2str(t_ind)]);
114  end;
115  if model.verbose == 9
116  fprintf('.');
117  end;
118 
119  t = (t_ind-1)*dt;
120  model = model.set_time(model,t);
121 % model.t = t;
122 
123  if operators_required
124 
125  if previous_implicit_operators_available
126  Ar0 = Ar1;
127  Br0 = Br1;
128  else
129  A0coeff = model.A_function_ptr(model,[]);
130  Ar0 = lincomb_sequence(reduced_data_t.Ar,A0coeff);
131  B0coeff = model.B_function_ptr(model,[]);
132  Br0 = lincomb_sequence(reduced_data_t.Br,B0coeff);
133  end;
134 
135  if model.enable_error_estimator
136  M1 = lincomb_sequence2(reduced_data_t.M1,A0coeff,A0coeff);
137  M2 = lincomb_sequence2(reduced_data_t.M2,B0coeff,B0coeff);
138  M3 = reduced_data_t.M3;
139  M4 = lincomb_sequence2(reduced_data_t.M4,B0coeff,A0coeff);
140  M5 = lincomb_sequence(reduced_data_t.M5,A0coeff);
141  M6 = lincomb_sequence(reduced_data_t.M6,B0coeff);
142  % output matrices required at next time!!
143  end;
144  model = model.set_time(model,t+dt);
145  if theta > 0
146  A1coeff = model.A_function_ptr(model,[]);
147  Ar1 = lincomb_sequence(reduced_data_t.Ar,A1coeff);
148  B1coeff = model.B_function_ptr(model,[]);
149  Br1 = lincomb_sequence(reduced_data_t.Br,B1coeff);
150  previous_implicit_operators_available = 1;
151  end;
152  % model.t = t+model.dt;
153  C1coeff = model.C_function_ptr(model,[]);
154  Cr1 = lincomb_sequence(reduced_data_t.Cr,C1coeff);
155  D1 = model.D_function_ptr(model,[]);
156  model = model.set_time(model,t);
157 % model.t = t;
158  if model.data_const_in_time
159  operators_required = 0;
160  end;
161  end;
162 
163  u0 = model.u_function_ptr(model);
164  xr_old = xr;
165  ddt_xr = Ar0 * xr + Br0 * u0;
166  if theta == 0
167  xr = xr + dt * ddt_xr;
168  else
169  model = model.set_time(model,t+dt);
170  u1 = model.u_function_ptr(model);
171  model = model.set_time(model,t+dt);
172  % solve M1 x1= rhs := M0 x0 + b
173  Ma1 = eye(size(Ar0))-dt * theta * Ar1;
174  Ma0 = eye(size(Ar0))+dt * (1-theta) * Ar0;
175  b = dt * ( theta * Br1 * u1 + (1-theta)*Br0 *u0);
176  rhs = Ma0 * xr + b;
177  xr = Ma1 \ rhs;
178  end;
179 
180  if model.enable_error_estimator
181 
182  % residual at previous timestep => "rectangular" integration
183  Rnormsqr = xr_old' * M1 * xr_old + u0' * M2 * u0 + ...
184  ddt_xr' * M3 * ddt_xr + ...
185  2*u0' * M4 * xr_old - 2*ddt_xr' * M5 * xr_old - ...
186  2 * ddt_xr' * M6 * u0;
187  % error estimator at next time
188  DeltaXt = DeltaXt + ...
189  model.state_bound_constant_C1 * model.dt ...
190  * sqrt(max(Rnormsqr,0));
191 % keyboard;
192  end;
193 
194  % store result
195  if save_time_index(t_ind+1)
196  Xr(:,t_column) = xr;
197  time_indices(t_column) = t_ind;
198  % for output u required at next timestep!
199  model = model.set_time(model,t+dt);
200 % model.t = t+model.dt;
201  u = model.u_function_ptr(model);
202  model = model.set_time(model,t);
203 % model.t = t;
204  Y(:,t_column) = Cr1 * xr + D1 * u;
205  time(t_column) = t;
206 
207  if model.enable_error_estimator
208  DeltaX(t_column) = DeltaXt;
209  Rnorms(t_column) = sqrt(max(Rnormsqr,0));
210  end;
211 
212  t_column = t_column + 1;
213  end;
214 
215 end;
216 
217 if model.verbose == 9
218  fprintf('\n');
219 end;
220 
221 rb_sim_data = [];
222 rb_sim_data.Xr = Xr;
223 rb_sim_data.Y = Y;
224 if model.enable_error_estimator
225  rb_sim_data.DeltaX = DeltaX;
226  rb_sim_data.DeltaY = DeltaX * model.output_bound_constant_C2;
227  rb_sim_data.Rnorms = Rnorms;
228 end;
229 rb_sim_data.time = time;
230 rb_sim_data.time_indices = time_indices;
231 
232 %if model.debug
233 % if max(rb_sim_data.Y)>1e10
234 % error('stability problem, please inspect!');
235 % end;
236 %end;
237 
238 %keyboard;
239 
240 %if model.debug
241 % disp('short before leavin lin_ds_rb_simulation');
242 %end;
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