rbmatlab  1.16.09
 All Classes Namespaces Files Functions Variables Modules Pages
rb_burgers_fem_simulation.m
1 function simulation_data = rb_burgers_fem_simulation(reduced_data,params)
2 %function simulation_data = rb_burgers_fem_simulation(reduced_data,params)
3 %
4 % function, which performs a reduced basis online simulation for
5 % the parameter vector mu, which is assumed to be set in the params
6 % struct, i.e. a previous params = params.set_mu([...],params) is assumed
7 % to have taken place.
8 %
9 % allowed dependency of data: Nmax, N, M, mu
10 % not allowed dependency of data: H
11 % allowed dependency of computation: Nmax, N, M, mu
12 % not allowed dependency of computation: H
13 % Unknown at this stage: ---
14 %
15 % Required fields of params as required by the numerical operators and
16 % mu_names : the cell array of names of mu-components and
17 % for each of these stringt, a corresponding
18 % field in params is expected.
19 % T : end-time of simulation
20 % nt : number of timesteps to compute
21 % newton_eps : epsilon for newton loop
22 %
23 % generated fields of simulation_data:
24 % a: time sequence of solution coefficients, columns are
25 % a(:,k) = a^(k-1)
26 %
27 % see the rb_***_online_prep() routine for specifications of the
28 % fields of reduced_data
29 
30 % Bernard Haasdonk 15.5.2007
31 
32 if params.verbose >= 10
33  disp(['performing simulation for mu = ',num2str(get_mu(params))]);
34 end;
35 
36 % set this value, if fixed number of newton steps is desired
37 fixed_newton_steps = params.fixed_newton_steps;
38 %fixed_newton_steps = 0;
39 % only relevant in case of fixed_newton_steps:
40 maxNewtonSteps = params.max_newton_steps;
41 
42 %grid = grid_geometry(params);
43 nRB = length(reduced_data.H{1});
44 a = zeros(nRB,params.nt+1); % matrix of RB-coefficients
45 params.dt = params.T/params.nt;
46 
47 
48 b = zeros(nRB, 1);
49 %(params.nt+1)*(maxNewtonSteps+1));
50 
51 % initial data projection and linear combination
52 params.affine_decomp_mode = 'coefficients';
53 % init data is given as the dirichlet data at time 0
54 
55 % Newton-eps squared:
56 epssqr = params.newton_eps^2;
57 
58  theta_dir = dirichlet_values([],[],params);
59 
60 M = reduced_data.M;
61 H = lincomb_sequence(reduced_data.H,theta_dir);
62 H1 = lincomb_sequence(reduced_data.H1,theta_dir);
63 H2 = lincomb_sequence(reduced_data.H2,theta_dir);
64 Hr = H1 + H2;
65 
66 %disp('Anfangswerte richtig setzen aus H muss a(:,1) berechnet werden!!')
67 
68 a(:,1) = M \ Hr;
69 b(:,1)= a(:,1);
70 
71 for i=1:size(a,1)
72  disp(['startwert ', num2str(a(i,1))]);
73  end;
74 
75 
76 %a0 = lincomb_sequence(reduced_data.H,theta_dir);
77 %a(:,1) = a0;
78 
79 % loop over time steps: computation of a(t+1)==a^t
80 for t = 1:params.nt
81  params.t = (t-1)*params.dt;
82  if params.verbose >= 10
83  disp(['entered time-loop step ',num2str(t)]);
84  end;
85 
86  % linear combination of all quantities
87  % only once, if data is const in time
88 
89  theta_eta = diffusivity([],[],params);
90  Q_eta = length(theta_eta);
91 
92  theta_c = velocity([],[],params);
93  Q_c = length(theta_c);
94 
95  % M and H are known
96 
97  B1 = lincomb_sequence(reduced_data.B1,theta_eta);
98  B2 = lincomb_sequence2(reduced_data.B2,theta_c, theta_dir);
99  % B2 = transpose(B2);
100  B = B1+B2;
101 
102  G1 = lincomb_sequence2(reduced_data.G1,theta_eta, theta_dir);
103  G2 = zeros(size(G1));
104  G3 = lincomb_sequence3(reduced_data.G3,theta_c,theta_dir, theta_dir);
105  G4 = lincomb_sequence2(reduced_data.G4,theta_eta, theta_dir);
106 
107  G = G1+G2+G3+G4;
108 
109  A = lincomb_sequence(reduced_data.A,theta_c);
110 
111  continue_loop = 1;
112 
113  a_l_old = a(:,t);
114  newtonStep = 0;
115 % newton-iteration...
116  while continue_loop
117  newtonStep = newtonStep + 1;
118  disp(['newtonStep number ',num2str(newtonStep),' in time-loop ',num2str(t)]);
119 
120  Au = zeros(size(A,2),size(A,3));
121  for i = 1:size(A,1)
122  Au = Au + reshape(A(i,:,:), size(A,2),size(A,3)) * a_l_old(i);
123  end;
124 
125  rhs = -(M + params.dt * Au + params.dt * B) * a_l_old + params.dt * G + M * a(:,t);
126 
127  S = M + 2*params.dt * Au + params.dt * B;
128 
129  %if (t==0.1)
130  % keyboard;
131  % end;
132 %%% s * incr = rhs
133  incr = S \ rhs;
134 
135  a_l_new = a_l_old + incr;
136 % b(:,t+newtonStep)=a_l_new;cd
137 % keyboard;
138 % % solve LL_I* a(:,t+1) = LL_E * a(:,t) + bb
139 % rhs= LL_E * a(:,t) + bb;
140 % % check whether pure explicit problem:
141 % if (LL_I_is_speye == 1)
142 % a(:,t+1) = rhs;
143 % elseif (LL_I_is_speye == -1) && isequal(LL_I, speye(size(LL_I)))
144 % a(:,t+1) = rhs;
145 % else % solve linear system
146 % % disp('check symmetry and choose solver accordingly!');
147 % % keyboard;
148  % nonsymmetric solvers:
149  % [U(:,t+1), flag] = bicgstab(Li,rhs,[],1000);
150  % [U(:,t+1), flag] = cgs(Li,rhs,[],1000);
151  % symmetric solver, non pd:
152  % see bug_symmlq for a very strange bug: cannot solve identity system!
153  % [U(:,t+1), flag] = symmlq(Li,rhs,[],1000);
154  % [U(:,t+1), flag] = minres(Li,rhs,[],1000);
155  % symmetric solver, pd:
156  % [a(:,t+1), flag] = pcg(LL_I,rhs,[],1000);
157 % [a(:,t+1), flag] = bicgstab(LL_I,rhs,[],1000);
158 % if flag>0
159 % disp(['error in system solution, solver return flag = ', ...
160 % num2str(flag)]);
161 % keyboard;
162 % end;
163 % end;
164 
165  % fehler checken:
166  errsqr = (a_l_new-a_l_old)' * M * (a_l_new - a_l_old);
167  disp(['error ',num2str(errsqr)]);
168 
169  if (errsqr <= epssqr) && (fixed_newton_steps==0)
170  continue_loop = 0;
171  end;
172 
173  if (newtonStep>=maxNewtonSteps) && (fixed_newton_steps==1)
174  continue_loop = 0;
175  end;
176 
177 %%% if newtonStep <=maxNewtonSteps
178  b = [b, a_l_new];
179 % b(:,(t*11)+newtonStep)=a_l_new;
180  disp(['write each newton step', num2str(newtonStep),' value ']);
181 % end;
182 
183  a_l_old = a_l_new;
184  end;
185  a(:,t+1) = a_l_new;
186  %a(:,9)=a_l_new;
187 
188 end;
189 
190 %a(:,10)=M\H1;
191 %a(:,9)=M\H2;
192 
193 % return simulation result
194 simulation_data = [];
195 simulation_data.a = a;
196 simulation_data.b = b;
197 %simulation_data.b = b(:,1:(1+newtonStep));
198 
199 % TO BE ADJUSTED TO NEW SYNTAX
200 %| \docupdate
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
function Udirichlet = dirichlet_values(model, X, Y)
UDIRICHLET = DIRICHLET_VALUES([X],[Y], MODEL) Examples dirichlet_values([0,1,2],[1,1,1],struct(name_dirichlet_values, homogeneous, ... c_dir, 1)) dirichlet_values([0:0.1:1],[0],struct(name_dirichlet_values, xstripes, ... c_dir, [0 1 2], ... dir_borders, [0.3 0.6])) dirichlet_values([0:0.1:1],[0],struct(name_dirichlet_values, box, ... c_dir, 1, ... dir_box_xrange, [0.3 0.6], ... dir_box_yrange, [-0.1 0.1]))
function [ descr , rdescr , dmodel , rmodel ] = newton(steps, combined, M_by_N_ratio, noof_ei_extensions, use_laplacian, model_size, random, num_cpus, Mstrich)
small script demonstrating a buckley leverett problem discretization and RB model reduction ...
Definition: newton.m:17