1 function simulation_data = rb_burgers_fem_simulation(reduced_data,params)
2 %
function simulation_data = rb_burgers_fem_simulation(reduced_data,params)
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
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: ---
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
23 % generated fields of simulation_data:
24 % a: time sequence of solution coefficients, columns are
27 % see the rb_***_online_prep() routine for specifications of the
28 % fields of reduced_data
30 % Bernard Haasdonk 15.5.2007
33 disp(['performing simulation for mu = ',num2str(get_mu(params))]);
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;
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;
49 %(params.nt+1)*(maxNewtonSteps+1));
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
56 epssqr = params.newton_eps^2;
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);
66 %disp(
'Anfangswerte richtig setzen aus H muss a(:,1) berechnet werden!!')
72 disp(['startwert ', num2str(a(i,1))]);
76 %a0 = lincomb_sequence(reduced_data.H,theta_dir);
79 % loop over time steps: computation of a(t+1)==a^t
81 params.t = (t-1)*params.dt;
83 disp(['entered time-loop step ',num2str(t)]);
86 % linear combination of all quantities
87 % only once, if data is const in time
89 theta_eta = diffusivity([],[],params);
90 Q_eta = length(theta_eta);
92 theta_c = velocity([],[],params);
93 Q_c = length(theta_c);
97 B1 = lincomb_sequence(reduced_data.B1,theta_eta);
98 B2 = lincomb_sequence2(reduced_data.B2,theta_c, theta_dir);
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);
109 A = lincomb_sequence(reduced_data.A,theta_c);
117 newtonStep = newtonStep + 1;
118 disp(['newtonStep number ',num2str(newtonStep),' in time-loop ',num2str(t)]);
120 Au = zeros(size(A,2),size(A,3));
122 Au = Au + reshape(A(i,:,:), size(A,2),size(A,3)) * a_l_old(i);
125 rhs = -(M + params.dt * Au + params.dt * B) * a_l_old + params.dt * G + M * a(:,t);
127 S = M + 2*params.dt * Au + params.dt * B;
135 a_l_new = a_l_old + incr;
136 % b(:,t+newtonStep)=a_l_new;cd
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)
143 % elseif (LL_I_is_speye == -1) && isequal(LL_I, speye(size(LL_I)))
145 % else % solve linear system
146 % % disp('check symmetry and choose solver accordingly!');
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);
159 % disp(['error in system solution, solver return flag = ', ...
166 errsqr = (a_l_new-a_l_old)' * M * (a_l_new - a_l_old);
167 disp(['error ',num2str(errsqr)]);
169 if (errsqr <= epssqr) && (fixed_newton_steps==0)
173 if (newtonStep>=maxNewtonSteps) && (fixed_newton_steps==1)
177 %%% if newtonStep <=maxNewtonSteps
179 % b(:,(t*11)+newtonStep)=a_l_new;
180 disp(['write each
newton step', num2str(newtonStep),' value ']);
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));
199 % TO BE ADJUSTED TO NEW SYNTAX
function r = verbose(level, message, messageId)
This function displays messages depending on a message-id and/or a level. Aditionally you can set/res...
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 ...