rbmatlab  1.16.09
 All Classes Namespaces Files Functions Variables Modules Pages
stokes_rb_simulation.m
1 function rb_sim_data = stokes_rb_simulation(model, reduced_data)
2 %function rb_sim_data = stokes_rb_simulation(model, reduced_data)
3 %
4 % fixed point defect correction algorithm
5 
6 
7 start_total = tic;
8 rb_sim_data = [];
9 
10 model.decomp_mode = 2;
11 
12 % get parametric coefficients
13 [A_coeff, f_coeff] = model.operators(model, reduced_data);
14 
15 % evaluate affine decomposition
16 AN = lincomb_sequence(reduced_data.AN_comp, A_coeff);
17 fN = lincomb_sequence(reduced_data.fN_comp, f_coeff);
18 
19 % prepare algorithm
20 if model.has_nonlinearity
21  AN_nonlin = cell(1, reduced_data.N_nl);
22  A_nl_coeff = model.operators(model, reduced_data);
23  A_nl_coeff = A_nl_coeff(model.matrix_nonlin_indices);
24  for i = 1:reduced_data.N_nl
25  AN_nonlin{i} = zeros(reduced_data.N);
26  AN_nonlin{i}(1:reduced_data.N_nl, 1:reduced_data.N_nl)...
27  = lincomb_sequence(reduced_data.AN_nl_comp{i}, A_nl_coeff);
28  end
29  eval_nonlin = @(xN) lincomb_sequence(AN_nonlin, xN(1:reduced_data.N_nl));
30 else
31  eval_nonlin = [];
32 end
33 
34 % execute algorithm
35 [uN, defect, niterations] = VecMat.fixed_point_algorithm(AN, fN, eval_nonlin, 1e-9, 10000, 1); %0.9
36 
37 if isfield(model, 'RB_error_indicator') && ...
38  strcmp(model.RB_error_indicator, 'estimator')
39 
40  % compute error estimate
41  neg_suN_coeff = -A_coeff * uN';
42  % enrich_rd
43  %neg_suN_coeff = -uN * A_coeff';
44 
45  if model.has_nonlinearity
46  neg_suN_nl_coeff = -A_nl_coeff * uN(1:reduced_data.N_nl)';
47  neg_suN_nl_coeff = neg_suN_nl_coeff(:) * uN(1:reduced_data.N_nl)';
48  % enrich_rd
49  %neg_suN_nl_coeff = uN(1:reduced_data.N_nl) * uN(1:reduced_data.N_nl)';
50  %neg_suN_nl_coeff = -neg_suN_nl_coeff(:) * A_nl_coeff';
51 
52  neg_suN_coeff = [neg_suN_coeff(:); neg_suN_nl_coeff(:)];
53  end
54  res_coeff = [f_coeff; neg_suN_coeff(:)];
55  res_norm_sqr = sum((reduced_data.G * res_coeff).^2);
56  %res_norm_sqr = res_coeff' * reduced_data.G * res_coeff;
57 
58  res_norm_sqr(res_norm_sqr < 0) = 0;
59  res_norm = sqrt(res_norm_sqr);
60 
61  % get stability constants:
62 
63  if ~isfield(model, 'infsup_method')
64  beta = model.coercivity_alpha(model);
65  else
66  if strcmp(model.infsup_method, 'scm')
67  beta = SCMonline(get_mu(model),model,reduced_data,reduced_data.infsup_constant,'lb');
68  elseif strcmp(model.infsup_method, 'rbf_interpolant')
69  mu = get_mu(model);
70  for i = 1:length(mu)
71  mu(i) = (mu(i)-model.mu_ranges{i}(1))/ ...
72  (model.mu_ranges{i}(2)-model.mu_ranges{i}(1));
73  end
74  beta = reduced_data.infsup_constant(mu);
75  end
76  end
77 
78  if model.has_nonlinearity
79 
80  rho_sqr = reduced_data.rho_sqr;
81  tau = 4*sqrt(2)*rho_sqr * res_norm / beta^2;
82  if tau < 1
83  rb_sim_data.Delta = beta / (2*sqrt(2)*rho_sqr) * (1-sqrt(1-tau));
84  else
85  rb_sim_data.Delta = tau;
86  if model.verbose > 2
87  warning('stokes:rb_simulation', 'Delta is tau.');
88  end
89  end
90  else
91 
92  rb_sim_data.Delta = res_norm / beta;
93  end
94 
95 end
96 
97 % save result
98 rb_sim_data.uN = uN;
99 rb_sim_data.defect = sqrt(defect' * reduced_data.KN * defect);
100 rb_sim_data.niterations = niterations;
101 rb_sim_data.total_time = toc(start_total);
102 
103 end