vi_rb_simulation.m

function rb_sim_data = vi_rb_simulation(model,reduced_data)
%function rb_sim_data = vi_rb_simulation(model,reduced_data)
%
% function performing a reduced simulation of a vi-model.
%
% Generated fields of rb_sim_data:
%  UN: reduced solution vector
%  LN: reduced Lagrange multiplier vector
%  Delta_lambda: error estimator for lambda variable
%  Delta_u: error estimator for u

% I.Maier 31.05.2011

rb_sim_data = [];

model.decomp_mode = 2;

[A_coeff,dummy,f_coeff,g_coeff] = model.operators(model,[]);

% compute required matrices
AN = lincomb_sequence(reduced_data.AN_comp,A_coeff);
BN = reduced_data.BN;
fN = lincomb_sequence(reduced_data.fN_comp,f_coeff);
gN = lincomb_sequence(reduced_data.gN_comp,g_coeff);

warning off
options = optimset;
%options.TolFun = 0;
options.MaxIter = 1000; % required for large bases...
options.Display = 'off';
rb_sim_data.UN = quadprog(AN,-fN,BN',gN,[],[],[],[],[],options);
warning on

% BN may have low rank
rb_sim_data.LN = pinv(BN)*(fN - AN*rb_sim_data.UN);

% a-posteriori error estimators
if model.enable_error_estimator
  % currently: no offline-online, high.dim evaluation:
  model.decomp_mode = 0;  
  [A,B,f,g] = model.operators(model,reduced_data);
  K = model.get_inner_product_matrix(model,reduced_data);
  K_v_r = (f-A*(reduced_data.RB_U*rb_sim_data.UN)...
           -B*reduced_data.RB_L*rb_sim_data.LN);
  v_r = K\K_v_r;
  delta_r = sqrt(max(sum(v_r .* K_v_r),0));

  %  Kinv_eta_s = B*reduced_data.RB_U*rb_sim_data.UN-g;
  %  eta_s = K * Kinv_eta_s;  
  %  pi_eta_s = max(eta_s,0); % projection pi!!
  %  Kinv_pi_eta_s = K\pi_eta_s;
  %  delta_s1 = pi_eta_s'* Kinv_pi_eta_s;
  %  delta_s2 = (reduced_data.RB_L*rb_sim_data.LN)' * Kinv_pi_eta_s;
  
  Kinv_eta_s = B*reduced_data.RB_U*rb_sim_data.UN-g;
  Kinv_pi_eta_s = max(Kinv_eta_s,0); % projection pi!!  
  delta_s1 = sqrt(Kinv_pi_eta_s'* K * Kinv_pi_eta_s);
  delta_s2 = (reduced_data.RB_L*rb_sim_data.LN)' * Kinv_pi_eta_s;
  
  %disp(['delta_s1: ',num2str(delta_s1),...
  %     ', delta_s2: ',num2str(delta_s2),...
  %     ', delta_r: ',num2str(delta_r)]);

  % compare with true residual
  %sim_data = detailed_simulation(model,reduced_data);
  %tilde_Kinv_eta_s = B*sim_data.U-g;
  %plot([Kinv_eta_s, tilde_Kinv_eta_s]);
  
  gamma_a = model.gamma_a(model);
%  gamma_a = svds(K^(-0.5)*A*K^(-0.5),1);
  alpha_a = model.alpha_a(model);
%  alpha_a = eigs(K^(-0.5)*A*K^(-0.5),1,0);
  beta =  model.beta_b;
%  keyboard;
  %  beta =  svds(B,1,'s')
%  disp(['mu: ',num2str(transpose(model.get_mu(model)))]);
%  disp(['beta: ',num2str(beta),', alpha:',num2str(alpha_a),...
%       ', gamma: ',num2str(gamma_a)]);
  c1 = (delta_r + delta_s1 * gamma_a / beta)/(2*alpha_a);
  c2 = (delta_s1*delta_r/beta + delta_s2)/alpha_a;
  Delta_u = c1 + sqrt(c1^2 + c2);
  rb_sim_data.Delta_u = Delta_u;
  rb_sim_data.Delta_lambda = (delta_r + gamma_a * Delta_u)/beta;
%  keyboard;
end;