dom_dec_gen_detailed_data.m

Go to the documentation of this file.

function detailed_data = dom_dec_gen_detailed_data(model, ...
                                                  model_data)
% function detailed_data = dom_dec_gen_detailed_data(model,model_data)
% function generating the reduced basis functions.
%
% Required fields of model:
%  RB_generation_mode: there are two options: PCA_trajectories and
%  POD_greedy
%
% Generated fields of detailed_data:
%  RB: cell of the reduced basis on the subdomains
%  I_0: cell of indizes into detailed_data.RB
%  I_G: cell of indizes into detailed_data.RB
%  RB_info: information about the generation process, depends on
%  generation mode

% I.Maier, 19.07.2011

detailed_data = model_data;

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

[RB_0,RB_G] = model.get_rb_from_detailed_data(detailed_data);


switch model.RB_generation_mode
  
 case 'PCA_trajectories'
  
  % generate parameter sample
  mu_params.numintervals = model.RB_numintervals;
  mu_params.range = model.mu_ranges;
  
  mu_mesh = cubegrid(mu_params);
  mu_sample = get(mu_mesh,'vertex');
  
  X = cell(1, 2);
  
  for m = 1:size(mu_sample,1)
    
    model = set_mu(model,mu_sample(m,:));
    sim_data = detailed_simulation(model,detailed_data);
    fprintf('.');
    
    if sim_data.reached_maxiter
        warning('detailed simulation may have diverged!');
    end;
    % dof matrizes depend on model.sequence_mode
    U = model.get_dofs_from_sim_data(model,sim_data);
    
    for i = 1:2
        X{i} = [X{i}, U{i}];
    end
    
  end;
  
  fprintf('\n');
  
  [RB_0, RB_G] = model.orthonormalize(model, detailed_data, X, model.RB_generation_epsilon);
  
  detailed_data = model.set_rb_in_detailed_data(detailed_data, RB_0, RB_G);
  
  
 case 'POD_greedy'
    
  % check if former greedy steps are available
  if ~isfield(detailed_data,'RB_info')
    detailed_data.RB_info = [];
    step = 0;
    bases_sizes.N0 = cell(1,2);
    bases_sizes.NG = cell(1,2);
    max_errs = [];  
    min_thetas = [];
    max_err_mus = [];
  else
    step = detailed_data.RB_info.steps;
    bases_sizes = detailed_data.RB_info.bases_sizes;
    max_errs = detailed_data.RB_info.max_errs;
    min_thetas = detailed_data.RB_info.min_thetas;
    max_err_mus = detailed_data.RB_info.max_err_mus;
  end;
  detailed_data.RB_info.stopped_on_tolerance = 0;
  detailed_data.RB_info.stopped_on_empty_extension = 0;
  
  % generate training parameter sample
  par.numintervals = model.RB_numintervals;
  par.range = model.base_model.mu_ranges;
  
  par_mesh = cubegrid(par);
  detailed_data.RB_info.M_train = get(par_mesh,'vertex');
    
  % step 0: initialize basis
  [N0,NG] = model.get_rb_size(detailed_data);

  RBext_0 = cell(1,2);
  RBext_G = cell(1,2);
  
  switch model.RB_init_mode
    
   case 'random' % random empirical mode
    if ~all([N0{dir},NG{1},NG{2}])
      model = set_mu(model,rand_uniform(1,model.mu_ranges)');
      
      old_ext_method = model.RB_extension_method;
      model.RB_extension_method = 'A';
      [RBext_0,RBext_G] = model.RB_extension_PCA_fixspace(...
          model,detailed_data);
      model.RB_extension_method = old_ext_method;
    end;
    
   case 'eigenbasis' % eigenbasis on interface plus random empirical mode
    if (~N0{dir})
      model = set_mu(model,rand_uniform(1,model.mu_ranges)');
      
      [RBext_0,RBext_G] = model.RB_extension_PCA_fixspace(...
          model,detailed_data);
    end;
    if (NG{1} < model.N_eigenbasis || NG{2} < model.N_eigenbasis)
      RBext_G = model.RB_extension_eigenbasis(model,detailed_data);
      for i = 1:2
        RBext_G{i} = RBext_G{i}(:,1:(model.N_eigenbasis - NG{i}));
      end;
    end;

  end;
  
  % prepare for greedy algorithm
  requireProjection1 = strcmp(model.RB_error_indicator, ...
                              'projection_error1');
  requireProjection2 = strcmp(model.RB_error_indicator, ...
                              'projection_error2');
  requireReducedData = ...
      strcmp(model.RB_error_indicator,'rb_error_estimate') | ...
      strcmp(model.RB_error_indicator,'rb_error');
  requireBaseSimData = ...
      strcmp(model.RB_error_indicator,'projection_error1') | ...
      strcmp(model.RB_error_indicator,'projection_error2') | ...
      strcmp(model.RB_error_indicator,'rb_error');
  
  % save detailed simulations of the base_model for all
  % parameters, if needed
  if (requireBaseSimData)
    save_detailed_simulations(model.base_model, ...
                              detailed_data.base_model_data, ...
                              detailed_data.RB_info.M_train', ...
                              model.base_model ...
                              .detailed_train_savepath);
  end;
  
  % inner product matrices for the computation of the projection
  % errors
  if (requireProjection1)
    K = model.get_inner_product_matrices(detailed_data);
  end;
  if (requireProjection2)
    K = model.base_model.get_inner_product_matrix( ...
        detailed_data.base_model_data);
  end;
  
  
  % start loop
  continue_loop = 1;
  
  while (continue_loop)
    
    % extend bases
    N_remain = model.N_max{1};
    for i = 1:2
      N_remain = min([N_remain, model.N_max{i} - N0{i} - NG{i}]);
    end;
    for i = 1:2
      N_extend = min([N_remain, size(RBext_G{i},2)]);
      if (N_extend > 0)
        RB_G{i} = [RB_G{i}, RBext_G{i}(:, 1:N_extend)];
      end;
      N_extend = min([model.N_max{i} - N0{i} - NG{i} - N_extend, ...
                      size(RBext_0{i},2)]);
      if (N_extend > 0)
        RB_0{i} = [RB_0{i}, RBext_0{i}(:, 1:N_extend)];
      end;
    end;
    
    detailed_data = model.set_rb_in_detailed_data(detailed_data, ...
                                                  RB_0,RB_G);
    N0_old = N0;
    NG_old = NG;
    [N0,NG] = model.get_rb_size(detailed_data);
    
    % check if extension is empty
    if (~any([N0{1}-N0_old{1}, N0{2}-N0_old{2}, ...
              NG{1}-NG_old{1}, NG{2}-NG_old{2}]))       
      continue_loop = 0;
      detailed_data.RB_info.stopped_on_empty_extension = 1;
      
      break;
      
    end;
    
    % print new basis sizes
    disp('new basis sizes:');
    disp(['(N1_0,N1_G,N2_0,N2_G) = (',num2str(N0{dir}),',', ...
          num2str(NG{dir}),',',num2str(N0{no_dir}),',', ...
          num2str(NG{no_dir}),')']);
    fprintf('\n');
    
    % save new basis sizes
    for i = 1:2
      bases_sizes.N0{i} = [bases_sizes.N0{i}; N0{i}];
      bases_sizes.NG{i} = [bases_sizes.NG{i}; NG{i}];
    end;
    
    % next step
    step = step + 1;
    
    % prepare error computation
    errs = zeros(size(detailed_data.RB_info.M_train,1),1);
    thetas = zeros(size(detailed_data.RB_info.M_train,1),1);
    
    if (requireReducedData)
      % generate required data for the reduced simulations
      reduced_data = gen_reduced_data(model,detailed_data);    
    end;
    
    if (requireProjection1)
      XN = cell(1,2);
      for i = 1:2
          XN{i} = orthonormalize([RB_0{i},RB_G{i}],K{i},[],'qr');
      end;
    end;
    
    if (requireProjection2)
      XN = zeros(detailed_data.base_model_data.df_info.ndofs,N0{1}+ ...
                 N0{2}+NG{1});
      
      % insert XN_0 functions
      XN(detailed_data.masks{1},1:N0{1}) = RB_0{1};
      XN(detailed_data.masks{2},N0{1}+(1:N0{2})) = RB_0{2};
      
      % insert XN_Gamma functions
      XN(detailed_data.masks{1},(N0{1}+N0{2}+1):end) = RB_G{1};
      XN(detailed_data.masks{2},(N0{1}+N0{2}+1):end) = RB_G{2};
      
      XN = orthonormalize(XN,K,[],'qr');
    end;
    
    
    % testing on maximal error
    fprintf('testing on maximal error ');
    
    for m = 1:size(detailed_data.RB_info.M_train,1)
      
      rb_sim_data = [];
      model = set_mu(model,detailed_data.RB_info.M_train(m,:));
      
      if (requireBaseSimData)
        % load detailed simulation of base model
        rb_sim_data.base_sim_data = load_detailed_simulation( ...
            m,model.base_model.detailed_train_savepath, ...
            model.base_model);  
      end;
      
      switch (model.RB_error_indicator)
        
       case 'projection_error1'

        errs(m) = dom_dec_compute_projection_error(...
            detailed_data,XN,K,rb_sim_data.base_sim_data);
        
       case 'projection_error2'

        errs(m) = dom_dec_compute_projection_error(...
            detailed_data,XN,K,rb_sim_data.base_sim_data);
        
       case 'rb_error_estimate'

        rb_sim_data = rb_simulation(model,reduced_data);
        if rb_sim_data.reached_maxiter
          warning(['reached maximal iterations in reduced', ...
                   ' simulation!']);      
        end;
        errs(m) = rb_sim_data.Delta(end);
        thetas(m) = rb_sim_data.RB_thetaN;
        
       case 'rb_error'
        
        rb_sim_data = rb_simulation(model,reduced_data);
        if rb_sim_data.reached_maxiter
          warning(['reached maximal iterations in reduced', ...
                   ' simulation!']);
        end;
          
        % compute error
        rb_sim_data = rb_reconstruction(model,detailed_data, ...
                                        rb_sim_data);
        rb_sim_data = model.compute_error(model,detailed_data, ...
                                          rb_sim_data);
        
        errs(m) = rb_sim_data.X_err(end);
        thetas(m) = rb_sim_data.RB_thetaN;
        
      end;

      if mod(m,ceil(size(detailed_data.RB_info.M_train,1)/45)) == 1
        fprintf('.');
      end;
      
    end;

    
    % print/save new information
    fprintf('\n');
    [max_err, err_ind] = max(errs);
    [min_theta, ~] = min(thetas);    
    mu_max_err = detailed_data.RB_info.M_train(err_ind,:);
    
    disp(['    maximal error: ',num2str(max_err)]);
    fprintf(['    maximal error mu-vector: (', ...
             num2str(mu_max_err(1))]);    
    for i = 2:length(model.mu_ranges)
      fprintf([',',num2str(mu_max_err(i))]);
    end;
    fprintf(')\n');
    disp(['    minimal theta: ',num2str(min_theta)]);
    
    max_errs = [max_errs; max_err];
    min_thetas = [min_thetas; min_theta];
    max_err_mus = [max_err_mus; mu_max_err];
    
    
    % test tolerance
    if max_err < model.RB_greedy_tolerance
      continue_loop = 0;
      detailed_data.RB_info.stopped_on_tolerance = 1;
    end;
        
    if continue_loop == 1      
      % compute basis extensions
      model = set_mu(model,mu_max_err);
      [RBext_0,RBext_G] = model.RB_extension_PCA_fixspace(model, ...
                                                  detailed_data);
    end;
    
  end;
  
  detailed_data.RB_info.steps = step;
  detailed_data.RB_info.bases_sizes = bases_sizes;
  detailed_data.RB_info.max_errs = max_errs;
  detailed_data.RB_info.min_thetas = min_thetas;
  detailed_data.RB_info.max_err_mus = max_err_mus;
  if model.RB_sequence_mode
    detailed_data.RB_info.extension_method = ...
        [model.RB_extension_method, '-SEQ'];
  else
    detailed_data.RB_info.extension_method = ...
        [model.RB_extension_method, '-NSEQ'];
  end;
  
end;