dictionary_rb_simulation.m

Go to the documentation of this file.

function rb_sim_data = dictionary_rb_simulation(model, ...
                                                reduced_data);
%function rb_sim_data = dictionary_rb_simulation(model, ...
%                                               reduced_data);
%  Online rb simulation by online greedy minimizing redidual norm
%  over dictionary, or normal RB simulation:
%
%  model.rb_online_mode: 0 standard rb, no adaptivity
%  model.rb_online_mode: 1 Delta-extended indicator for online extension
%
% B. Haasdonk 5.3.2015

% prepare mu-dependent quantities

rb_sim_data =[];

old_mode = model.decomp_mode;
model.decomp_mode = 2; % coefficients

tic; % measure linear combination time

[a_coeff,f_coeff] = ...
    model.operators(model,[]);

AN = lincomb_sequence(reduced_data.AN_comp, a_coeff);
fN = lincomb_sequence(reduced_data.fN_comp, f_coeff);
Qf = length(f_coeff);
Qa = length(a_coeff);
N = reduced_data.N;

% transform G matrix, such that later simple mult:
% G_mu matrix has <vf,vf>, <vf,va_phi_1(mu)> , <vf,va_phi_N> 
% as first row and column, and  <va_phi_i,va_phi_j(mu)> as
% remaining entries 
% the following is the correct code, if lin_stat_reduced_data is
% passed as argument.
%if ~isfield(reduced_data,'used_dictionary_gen_reduced_data');
%  Saa = kron(speye(N),sparse(1:Qa,1,a_coeff(:),Qa,1)); 
%  S = [f_coeff(:), sparse([],[],[],Qf, N); ...
%       sparse([],[],[], Qa * N, 1), Saa];
%  G = S' * reduced_data.G * S;
%else
Hff = f_coeff' * reduced_data.bar_Hff * f_coeff;
af_coeff = a_coeff * f_coeff';
Haf = reduced_data.bar_Haf * af_coeff(:);
aa_coeff = a_coeff * a_coeff';
Haa_vec = reduced_data.bar_Haa * aa_coeff(:);
Haa = reshape(Haa_vec,[N,N]);
G = [Hff,Haf'; Haf,Haa];  
%end;
G = 0.5*(G+G'); % G is positive definite.
diag_G = diag(G);

t_lincomb = toc;

%if 0
%  disp('check, whether G has correct entries');
%  % get full operators
%  model.decomp_mode = 0; % complete  
%  [A ,f] = ...
%      model.operators(model,reduced_data);
%  K = model.get_inner_product_matrix(reduced_data);
%  vf = K\f;
%  G(1,1)
%  vf'*K * vf
%  va_phi1 = K\(A*reduced_data.RB(:,1)); 
%  G(2,1)
%  vf'*K*va_phi1
%  G(2,2)
%  va_phi1'*K*va_phi1 
%  keyboard;
%end;

tic;

if model.rb_online_mode == 0 
  % standard rb simulation, no online enrichment  
  uN = AN\fN;  
  % the following works with standard lin_stat_gen_reduced_data:
  %  Q_r = size(reduced_data.G,1);
  %  neg_auN_coeff = -A_coeff * uN'; 
  %  res_coeff = [f_coeff; neg_auN_coeff(:)];
  %  res_norm_sqr = res_coeff' * reduced_data.G * res_coeff;
  % this is required for dictionary_gen_reduced_data:
  res_coeff = [1;-uN];
  res_norm_sqr = res_coeff' * G * res_coeff;
  res_norm_sqr = max(res_norm_sqr,0);    
  res_norm = sqrt(res_norm_sqr);
  rb_sim_data.Delta = ...
      res_norm/model.coercivity_alpha(model); 
  
elseif model.rb_online_mode == 1 
  % online enrichment by eta(v) = - Delta_N(v; Phi \cup v)
  diag_AN = diag(AN);
  
  %uN = ones(size(AN,1),1);
  
  % run iterative loop of online_greedy by minimum residual
  
  eps_sequence = [];
  current_eps = 1e10; % do at least one extension
  basis_indices = [];
  non_basis_indices = 1:N; 
  iter = 0;
  ANinv = [];
  UN = zeros(0,N);
  uN = zeros(0,1);
  
  %  if model.RB_online_stop_Nmax == 2 
  %    keyboard;
  %  end;
  
  if N<model.RB_online_stop_Nmax
    model.RB_online_stop_Nmax = N;
    disp(['warning: dictionary too small, Now decreasing RB_online_stop_Nmax']);
  end;

  while (current_eps > model.RB_stop_epsilon) && (iter < model.RB_online_stop_Nmax)
    % compute all new uN vectors at once:
    a = diag_AN(non_basis_indices)';
    bar_a = AN(basis_indices,non_basis_indices);
    bar_bar_a = AN(non_basis_indices,basis_indices)';
    %  tilde_f = fN(non_basis_indices)' - sum(bar_bar_a.*UN,1);
    tilde_f = fN(non_basis_indices)' - uN' * bar_bar_a;
    ANinv_bar_a = (ANinv*bar_a);
    tilde_a = a - sum(bar_bar_a .* ANinv_bar_a,1); 
    %  UN_old = UN;
    UN = uN * ones(1,length(non_basis_indices));
    uN_old = uN;
    n = length(tilde_f);
    UN = [UN; zeros(1,length(non_basis_indices))] ...
         + [-ANinv_bar_a;ones(1,length(non_basis_indices))] * ...
         sparse(1:n,1:n,tilde_f./tilde_a,n,n);
%        spdiag((tilde_f./tilde_a),length(tilde_f));
    UN_new = UN;
    % compute all residualnormssquare at once:
    G11 = G([1,1+basis_indices],[1,1+basis_indices]);
    UN1 = [ones(1,length(non_basis_indices)); - UN(1:end-1,:)];
    UN2 = -UN(end,:);
    G12 = G([1,1+basis_indices],1+non_basis_indices);
    n = length(UN2);
    %    g = [G11 * UN1 + G12*spdiag(UN2(:),length(UN2)); ...
    g = [G11 * UN1 + G12*sparse(1:n,1:n,UN2,n,n); ...
         sum(G12 .* UN1,1) + diag_G(1+non_basis_indices)'.*UN2];
    resnormsqr = sum([ones(1,length(non_basis_indices)); - UN] .* g,1);
    if min(resnormsqr<-1e-2)
      disp('error: negative residual, inspect workspace!');
      % test for first vector:
      keyboard;
    end;
    %  if resnormsqr(3)>1e-10
    %    disp('error: residual for func 3 should be 0!!!');  
    %    plot(resnormsqr);
    %    title('residual norms');
    %    % the following test was successful
    %    %    G1 = [G11, G12(:,1); G12(:,1)', diag_G(1+1)];
    %    %    g1 = G1 * [1; -UN(:,1)]
    %    %    g(:,1)    
    %    %    resnormsqr1 = [1; -UN(:,1)]' * G1 * [1; -UN(:,1)]
    %    %    resnormsqr(1) 
    %    %    G3 = [G11, G12(:,3); G12(:,3)', diag_G(1+3)];
    %    %    resnormsqr3 = [1; -UN(:,3)]' * G3 * [1; -UN(:,3)]
    %    resnormsqr(3) 
    %    keyboard;
    %  end;
    %  end;
    % set NaN values to very high, in order not set to 0 by clipping:
    p = find(isnan(resnormsqr));
    resnormsqr(p) = 1e10;
    % take positive part to avoid numerical instabilities  
    resnormsqr = max(resnormsqr,0);    
    % search maximum  
    [minresnormsqr, new_index] = min(resnormsqr);
    % select one if multiple minima
    current_eps = sqrt(minresnormsqr(1));
    eps_sequence = [eps_sequence, current_eps];
    new_index = new_index(1);
    % remove new basis vector from search set
    basis_indices = [basis_indices, non_basis_indices(new_index)];
    non_basis_indices = [non_basis_indices(1:new_index-1),...
                    non_basis_indices(new_index+1:end)];
    ANinv = inv(AN(basis_indices,basis_indices));
    uN = UN(:,new_index);
    UN = [UN(:,1:new_index-1),UN(:,new_index+1:end)];
    iter = iter + 1;
    if iter > 100
      disp('investigate iteration')
      % true reduced solution with selected basis vectors of previous
      % step:
      detailed_data = reduced_data;
      detailed_data.RB = detailed_data.RB(:,basis_indices(1:end-1));
      reduced_data_tmp = gen_reduced_data(model,detailed_data);
      rb_sim_data_tmp = rb_simulation(model,reduced_data_tmp);
      disp('correct UN old:')
      rb_sim_data_tmp.uN
      disp('current UN_old:')
      uN_old
      %    UN_new(:,new_index)
      
      % true reduced solution with selected basis vectors:
      detailed_data = reduced_data;
      detailed_data.RB = detailed_data.RB(:,basis_indices);
      reduced_data_tmp = gen_reduced_data(model,detailed_data);
      rb_sim_data_tmp = rb_simulation(model,reduced_data_tmp);
      disp('correct UN:')
      rb_sim_data_tmp.uN
      disp('current UN:')
      uN
      %    UN_new(:,new_index)
      
      % recompute by hand:
      %disp('recomputed current uN')
      %uN_new = [uN_old; 0] + tilde_f(new_index)/tilde_a(new_index)*[-ANinv_bar_a(:,new_index);1]
      
      keyboard; 
    end;
    %  keyboard;
  end;
  
  % return outputs

  rb_sim_data.eps_sequence = eps_sequence;
  rb_sim_data.basis_indices = basis_indices;
  rb_sim_data.iter = iter;
  rb_sim_data.current_eps = current_eps;
  rb_sim_data.Delta = ...
      current_eps/model.coercivity_alpha(model); 
end; % rb_online_mode

t_solution = toc;

rb_sim_data.uN = uN;
rb_sim_data.t_lincomb = t_lincomb;
rb_sim_data.t_solution = t_solution;

model.decomp_mode = old_mode;