dom_dec_rb_simulation.m

Go to the documentation of this file.

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, ~, ~] = 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;

if model.base_model.compute_output_functional && ~model.base_model.dual_mode
    
    model.base_model.dual_mode = 1;
    rb_dual_sim_data = rb_simulation(model, reduced_data.dual_reduced_data);
    model.base_model.dual_mode = 0;
    
    model.base_model.decomp_mode = model.decomp_mode;
    
    sN = 0;
    sN1 = 0;
    v_coeff = model.operators_output(model.base_model, []);
    
    % l(u_N)
    for i = 1:2
        vN = lincomb_sequence(reduced_data.vN_comp{i},v_coeff);
        sN = sN + (vN(:)') * rb_sim_data.uN{i}(:,numiter);
    end
    
    % l(R_1gu_{N,i})
    for i = 1:2
        vN = lincomb_sequence(reduced_data.vNext_comp{i},v_coeff);
        sN = sN + (-1)^i * (vN(:)') * rb_sim_data.uN{i}(:,numiter);
    end
    
    % f(\psi_N)
    for i = 1:2
        vN = lincomb_sequence(reduced_data.fN_d_comp{i},f_coeff{i});
        sN1 = sN1 + (vN(:)') * rb_dual_sim_data.uN{i}(:,rb_dual_sim_data.numiter);
    end
    
    % f(R_2g\psi_{N,i})
    for i = 1:2
        vN = lincomb_sequence(reduced_data.vNfR_comp{i},f_coeff{dir});
        sN1 = sN1 - (-1)^i * (vN(:)') * rb_dual_sim_data.uN{i}(:,rb_dual_sim_data.numiter);
    end
    
    % A(u_N,\psi_N)
    for i = 1:2
        vN = lincomb_sequence(reduced_data.AN_pd_comp{i},A_coeff{i}) * ...
            rb_sim_data.uN{i}(:,numiter);
        sN1 = sN1 - (vN(:)') * rb_dual_sim_data.uN{i}(:,rb_dual_sim_data.numiter);
    end
    
    % A(u_N,R_2g\psi_{N,i})
    for i = 1:2
        vN = lincomb_sequence(reduced_data.vNAR1_comp{i},A_coeff{no_dir}) * ...
            rb_sim_data.uN{no_dir}(:,numiter);
        sN1 = sN1 + (-1)^i * (vN(:)') * rb_dual_sim_data.uN{i}(:,rb_dual_sim_data.numiter);
    end
    
    % A(R_1gu_{N,i},\psi_N)
    for i = 1:2
        vN = lincomb_sequence(reduced_data.vNAR2_comp{i},A_coeff{dir}) * ...
            rb_sim_data.uN{i}(:,numiter);
        sN1 = sN1 - (-1)^i * (vN(:)') * rb_dual_sim_data.uN{dir}(:,rb_dual_sim_data.numiter);
    end
    
    rb_sim_data.sN = sN - sN1;
    
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(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}(res_norm_sqr{i} < 0) = 0;
  
end;

% error estimator parts
if model.base_model.compute_output_functional || model.base_model.dual_mode
    % here: use \bar{u}_N^n approximation
    
    Delta_part1 = sqrt(res_norm_sqr{1} + res_norm_sqr{2})/ ...
        sqrt(model.coercivity_alpha(model));
    Delta_part2 = sqrt(model.continuity_gamma(model)) * sqrt(jumpex_norm_sqr);
else
    
    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);
end

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

% output error bound
if model.base_model.compute_output_functional && ~model.base_model.dual_mode
    
    rb_sim_data.Delta_s = rb_sim_data.Delta(end) * rb_dual_sim_data.Delta(end);
    rb_sim_data.rb_dual_sim_data = rb_dual_sim_data;
end