dom_dec_rb_simulation.m

function rb_sim_data = dom_dec_rb_simulation(model,reduced_data)
% function rb_sim_data = dom_dec_rb_simulation(model,reduced_data)
%
% realization of the iterative schema for RB methods according to
% my diploma thesis.
%
% Generated fields of rb_sim_data:
%  uN: cell, coefficient vectors of all iteratives
%  RB_thetaN: relaxation-parameter
%  numiter: number of iterations
%  reached_maxiter: indicates the cause of termination
%  reached_tolerance: indicates the cause of termination
%  Delta: vector of error estimates for all iteratives in uN

% I. Maier, 19.07.2011

rb_sim_data = [];
rb_sim_data.reached_maxiter = 0;
rb_sim_data.reached_tolerance = 0;

model.decomp_mode = 2;

dir = model.dirichlet_side;
no_dir = mod(dir,2)+1;

% get coefficients
[A_coeff, f_coeff, dummy1, dummy2] = model.operators(model,[]);

AN = cell(1,2);
fN = cell(1,2);
for i = 1:2
  AN{i} = lincomb_sequence(reduced_data.AN_comp{i},A_coeff{i});
  fN{i} = lincomb_sequence(reduced_data.fN_comp{i},f_coeff{i});
end;

% generate matrizes for iterative schema
AN_00 = cell(1,2);
AN_0G = cell(1,2);
for i = 1:2
  AN_00{i} = AN{i}(reduced_data.I_0{i},reduced_data.I_0{i});
  AN_0G{i} = AN{i}(reduced_data.I_0{i},reduced_data.I_G{i});
end;


if model.RB_approximate_thetaN
  %%%%%%%%%%%% approximation of optimal theta
  
  % generate random test vectors
  ntest_vecs = 10000;
  intervals = cell(1,length(reduced_data.I_G{no_dir}));
  for j = 1:length(intervals)
    intervals{j} = [0,1];
  end;
  Rtest_G = rand_uniform(ntest_vecs,intervals);
  
  % compute harmonic extensions:
  Rtest_norm = cell(1,2);
  for i = 1:2
    Rtest_0 = AN_00{i} \ (-AN_0G{i}*Rtest_G);
    Rtest = [Rtest_0; Rtest_G];

    Rtest_norm{i} = sum((AN{i}' * Rtest) .* Rtest,1);
  end;
  
  sigma = max( Rtest_norm{dir} ./ Rtest_norm{no_dir} );
  tau = max( Rtest_norm{no_dir} ./ Rtest_norm{dir} );
  
  RB_thetaN = (tau+1)/(sigma^2*tau+tau+2)/2;

end;

restart = 0;
rb_sim_data.reached_maxiter = 1;
while rb_sim_data.reached_maxiter
  
  rb_sim_data.reached_maxiter = 0;
  
  % initialize iteration
  uN = cell(1,2);
  uN_old = cell(1,2);
  uN_tilde = cell(1,2);
  zN = cell(1,2);
  BN = zeros(reduced_data.N{no_dir},1);
  
  rb_sim_data.uN = cell(1,2);
  for i = 1:2
    rb_sim_data.uN{i} = zeros(reduced_data.N{i},model.RB_maxiterN);
  end;
  
  for i = 1:2
    uN{i} = zeros(reduced_data.N{i},1);
    uN_tilde{i} = zeros(reduced_data.N{i},1);
    zN{i} = zeros(reduced_data.N{i},1);
  end;
  
  if model.RB_approximate_thetaN == 0
    a = [];
    sigma = 0;
    tau = 0;
    RB_thetaN = [];
  end;
  
  uN{dir}(reduced_data.I_0{dir}) = AN_00{dir} \ ...
      fN{dir}(reduced_data.I_0{dir});
  uN{no_dir} = AN{no_dir} \ fN{no_dir};
  uN_0_no_dir = uN{no_dir};
  
  for i = 1:2
    uN_tilde{i}(reduced_data.I_0{i}) = ...
        AN_00{i} \ (-AN_0G{i}*uN{i}(reduced_data.I_G{i}));
    uN_tilde{i}(reduced_data.I_G{i}) = uN{i}(reduced_data.I_G{i});  
  end;
  
  AN_dir = AN_00{dir}^(-1) * (-AN_0G{dir});
  fN_dir_G = fN{dir}(reduced_data.I_G{dir});
  AN_dir_G = AN{dir}(reduced_data.I_G{dir},:);
  AN_no_dir_inv = AN{no_dir}^(-1);
    
  if (restart == 2 || ~model.RB_approximate_thetaN)
    clear('AN_00');
    clear('AN_0G');
    clear('fN');
  end;
  
  %%%%%%%%%%%% iterative schema %%%%%%%%%%%%
  continue_loop = 1;
  numiter = 0;
  diff_norm = zeros(1,2);
  while continue_loop
    
    % ...
    numiter = numiter+1;
    for i = 1:2
      uN_old{i} = uN{i};
    end;
    
    % ...
    zN{dir}(reduced_data.I_0{dir}) = AN_dir * ...
        uN_tilde{no_dir}(reduced_data.I_G{no_dir});
    zN{dir}(reduced_data.I_G{dir}) = ...
        uN_tilde{no_dir}(reduced_data.I_G{no_dir});
    
    % ...
    if model.RB_approximate_thetaN == 0
      a = (zN{dir}'*AN{dir}*zN{dir}) / (uN_tilde{no_dir}'*AN{no_dir}* ...
                                        uN_tilde{no_dir});
      sigma = max(sigma,a);
      tau = max(tau,1/a);
      RB_thetaN = (tau+1)/(sigma^2*tau+tau+2)/2;
    end;
    
    % ...
    uN{dir} = uN{dir} + RB_thetaN*(zN{dir}-uN_tilde{dir});
    uN_tilde{dir} = uN_tilde{dir} + RB_thetaN*(zN{dir}- ...
                                               uN_tilde{dir});
    
    % ...
    BN(reduced_data.I_G{no_dir}) = fN_dir_G - AN_dir_G*uN{dir};
    zN{no_dir} = AN_no_dir_inv * BN;
    
    % ...
    uN{no_dir} = uN_0_no_dir + zN{no_dir};
    uN_tilde{no_dir} = uN_tilde{no_dir} + uN{no_dir} - ...
        uN_old{no_dir};
    
    % save uN
    for i = 1:2
      rb_sim_data.uN{i}(:,numiter) = uN{i};
    end;
    
    % stopping criteria
    for i = 1:2
      diff_norm(i) = sqrt((uN{i}-uN_old{i})' * AN{i} * ...
                          (uN{i}-uN_old{i}));
    end;
    
    if sum(diff_norm) < model.RB_sim_tol
      continue_loop = 0;
      rb_sim_data.reached_tolerance = 1;
    end;
    if numiter == model.RB_maxiterN
      continue_loop = 0;
      rb_sim_data.reached_maxiter = 1;
    end;
    
  end;
  
  rb_sim_data.RB_thetaN = RB_thetaN;
  rb_sim_data.numiter = numiter;
  for i = 1:2
    rb_sim_data.uN{i} = rb_sim_data.uN{i}(:,1:numiter);
  end;
  
  % if reached_maxiter procedure is started again 2 times
  % with new relaxation parameter
  if (restart == 2 || ~model.RB_approximate_thetaN)
    break;
  end;
  RB_thetaN = RB_thetaN/2;
  restart = restart+1;

end;

%%%%%%%%%%%% error estimation %%%%%%%%%%%%

% jump over interface
jumpex_norm_sqr = zeros(numiter,1);
for j = 1:numiter
  jumpex_norm_sqr(j) = rb_sim_data.uN{dir}(:,j)' ...
      * reduced_data.R_11 * rb_sim_data.uN{dir}(:,j) ...
      - 2*rb_sim_data.uN{dir}(:,j)' ...
      * reduced_data.R_12 * rb_sim_data.uN{no_dir}(:,j) ...
      + rb_sim_data.uN{no_dir}(:,j)' ...
      * reduced_data.R_22 * rb_sim_data.uN{no_dir}(:,j);
end;
jumpex_norm_sqr(find(jumpex_norm_sqr < 0)) = 0;

% residual
res_norm_sqr = cell(1,2);
for i = 1:2
  res_norm_sqr{i} = zeros(numiter,1);
  for j = 1:numiter
    neg_a1uN_coeff = -A_coeff{1} * rb_sim_data.uN{1}(:,j)';
    neg_a2uN_coeff = -A_coeff{2} * rb_sim_data.uN{2}(:,j)';
    res_coeff = [f_coeff{1}; neg_a1uN_coeff(:); neg_a2uN_coeff(:)];
    res_norm_sqr{i}(j) = res_coeff' * reduced_data.G{i} * res_coeff;
  end;
  res_norm_sqr{i}(find(res_norm_sqr{i} < 0)) = 0;
  
end;

% error estimator parts
Delta_part1 = sqrt(res_norm_sqr{1} + res_norm_sqr{2})/ ...
    model.coercivity_alpha(model);
Delta_part2 = (1+sqrt(model.continuity_gamma(model)/ ...
                      model.coercivity_alpha(model))) * ...
    sqrt(jumpex_norm_sqr);

% complete error estimator
rb_sim_data.Delta = Delta_part1 + Delta_part2;