1 function rb_sim_data = stokes_rb_simulation(model, reduced_data)
2 %
function rb_sim_data = stokes_rb_simulation(model, reduced_data)
4 % fixed point defect correction algorithm
10 model.decomp_mode = 2;
12 %
get parametric coefficients
13 [A_coeff, f_coeff] = model.operators(model, reduced_data);
15 % evaluate affine decomposition
16 AN = lincomb_sequence(reduced_data.AN_comp, A_coeff);
17 fN = lincomb_sequence(reduced_data.fN_comp, f_coeff);
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);
29 eval_nonlin = @(xN) lincomb_sequence(AN_nonlin, xN(1:reduced_data.N_nl));
35 [uN, defect, niterations] =
VecMat.fixed_point_algorithm(AN, fN, eval_nonlin, 1e-9, 10000, 1); %0.9
37 if isfield(model,
'RB_error_indicator') && ...
38 strcmp(model.RB_error_indicator,
'estimator')
40 % compute error estimate
41 neg_suN_coeff = -A_coeff * uN
';
43 %neg_suN_coeff = -uN * A_coeff';
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)';
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';
52 neg_suN_coeff = [neg_suN_coeff(:); neg_suN_nl_coeff(:)];
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;
58 res_norm_sqr(res_norm_sqr < 0) = 0;
59 res_norm = sqrt(res_norm_sqr);
61 % get stability constants:
63 if ~isfield(model, 'infsup_method
')
64 beta = model.coercivity_alpha(model);
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
')
71 mu(i) = (mu(i)-model.mu_ranges{i}(1))/ ...
72 (model.mu_ranges{i}(2)-model.mu_ranges{i}(1));
74 beta = reduced_data.infsup_constant(mu);
78 if model.has_nonlinearity
80 rho_sqr = reduced_data.rho_sqr;
81 tau = 4*sqrt(2)*rho_sqr * res_norm / beta^2;
83 rb_sim_data.Delta = beta / (2*sqrt(2)*rho_sqr) * (1-sqrt(1-tau));
85 rb_sim_data.Delta = tau;
87 warning('stokes:rb_simulation
', 'Delta is tau.
');
92 rb_sim_data.Delta = res_norm / beta;
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);