1 function rb_sim_data = lin_ds_rb_simulation(model,reduced_data)
2 %
function rb_sim_data = lin_ds_rb_simulation(model,reduced_data)
4 %
function performing a reduced basis simulation with a
5 % theta-scheme time discretization.
7 % Bernard Haasdonk 2.4.2009
10 disp(
'entered reduced simulation ');
15 % possibly shrink reduced_data detected by
this routine
16 reduced_data_t = lin_ds_reduced_data_subset(model,reduced_data);
18 %model.affine_decomp_mode =
'complete';
19 model.decomp_mode = 2; % coefficients
21 model = model.set_time(model,0);
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;
29 save_time_index = ones(1,model.nt+1);
30 save_time_indices = 1:(model.nt+1);
33 x0coeff = model.x0_function_ptr(model,[]);
34 xr = lincomb_sequence(reduced_data_t.x0r,x0coeff);
37 % instead of large dense matrix allocation equivalent:
39 % a' G b - (a' W) (V' G b) - (a' G V) (W' b) + (a' W) (V' G V) (W' b)
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:
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
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));
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));
65 if model.enable_error_estimator
66 DeltaX = zeros(1,length(save_time_indices));
69 time_indices = zeros(1,length(save_time_indices));
70 time = zeros(1,length(save_time_indices));
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
76 if model.enable_error_estimator
77 DeltaXt = model.state_bound_constant_C1*sqrt(max(e0sqr,0));
80 if save_time_index(t_ind+1)
81 Xr(:,t_column) = xr(:);
82 time_indices(t_column) = t_ind;
84 %u = model.u_function_ptr(model);
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;
91 if model.enable_error_estimator
92 DeltaX(t_column) = DeltaXt;
94 t_column = t_column + 1;
99 model.dt = model.T/model.nt;
103 operators_required = 1;
104 if ~isfield(model,'data_const_in_time')
105 model.data_const_in_time = 0;
107 previous_implicit_operators_available = 0;
110 for t_ind = 1:model.nt
113 disp(['entered time-loop step ',num2str(t_ind)]);
120 model = model.set_time(model,t);
123 if operators_required
125 if previous_implicit_operators_available
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);
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!!
144 model = model.set_time(model,t+dt);
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;
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);
158 if model.data_const_in_time
159 operators_required = 0;
163 u0 = model.u_function_ptr(model);
165 ddt_xr = Ar0 * xr + Br0 * u0;
167 xr = xr + dt * ddt_xr;
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);
180 if model.enable_error_estimator
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));
195 if save_time_index(t_ind+1)
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);
204 Y(:,t_column) = Cr1 * xr + D1 * u;
207 if model.enable_error_estimator
208 DeltaX(t_column) = DeltaXt;
209 Rnorms(t_column) = sqrt(max(Rnormsqr,0));
212 t_column = t_column + 1;
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;
229 rb_sim_data.time = time;
230 rb_sim_data.time_indices = time_indices;
233 % if max(rb_sim_data.Y)>1e10
234 % error('stability problem, please inspect!');
241 % disp('
short before leavin lin_ds_rb_simulation');
function r = verbose(level, message, messageId)
This function displays messages depending on a message-id and/or a level. Aditionally you can set/res...