dom_dec_final.m

Go to the documentation of this file.

% script generating figures of the paper "An Iterative Domain
% Decomposition Procedure for The Reduced Basis Method"

% IM, 11.06.2012


clear all;
close all;

fprintf('\n');
disp('Script generating the figures of the paper');
disp('==========================================');
disp(['"An Iterative Domain Decomposition Procedure for The Reduced' ...
      ' Basis Method"']);
disp(['============================================================' ...
      '==============']);
fprintf('\n\n');


% initialize
disp('MODEL');
disp('--------------------------------------------------------');
fprintf('Initialization ...');

params = [];
params.numintervals = 84;
params.RB_numintervals = 4;
params.detailed_train_savepath = 'detailed_train_84_rb4';

base_model = thermalblock_dd_model(params);
base_model.compute_output_functional = 0;
model = dom_dec_model(base_model, params);

clear('params');

model_data = gen_model_data(model);
mu_test = rand_uniform(100, model.mu_ranges)';

input_string = [];
input_string_filter = {'y', 'n'};

fprintf('done!\n\n');


% display discretization settings:
disp('Finite Element discretization: ');
disp(['mesh size: ', ...
      num2str(model_data.base_model_data.df_info.grid.hmin)]);
disp(['number of dofs: ', ...
      num2str(model_data.base_model_data.df_info.ndofs)]);
disp(['polynomial degree: ', num2str(base_model.pdeg)]);
fprintf('\n');


% visualize source function
% Figure 6.1
while ~any(strcmp(input_string, input_string_filter))
  input_string = input('-> Visualize source function? (y/n)', 's');
end;

if strcmp(input_string, 'y')
  source = femdiscfunc([], model_data.base_model_data.df_info);
  source_mu_list = {[0.1 0.1 0.1 0.35], [0.1 0.1 0.1 1]};
  
  for i = 1:2
    base_model = set_mu(base_model, source_mu_list{i});
    source = fem_interpol_local(base_model.source, source, base_model);
    
    disp(['plot ', num2str(i), ': mu4 = ', ...
          num2str(source_mu_list{i}(4))]);
    figure,
    plot(source, []);
    
    shading flat;
    axis equal;
    axis tight;
    set(gca, 'FontSize', 15);
    
  end;
  
  clear('source');
  clear('source_mu_list');
  
end;

fprintf('\n\n');


disp('BASIS GENERATION');
disp('--------------------------------------------------------');

% display basis generation settings:
disp('Basis generation settings:');
disp(['size of training set: ', ...
      num2str(model.RB_numintervals(1)+1), '^', ...
      num2str(length(model.RB_numintervals))]);
disp(['error tolerance: ', num2str(model.RB_greedy_tolerance)]);
disp(['maximal total basis size: {', num2str(model.N_max{1}), ...
      ', ', num2str(model.N_max{2}), '}']);
fprintf('\n');

method_strings = {'A-SEQ', 'B-SEQ', 'A-NSEQ', 'B-NSEQ'};

% basis generation
input_string = [];
while ~any(strcmp(input_string, input_string_filter))
  input_string = input('-> Generate bases? (unrecommended)(y/n)', 's');
end;

if strcmp(input_string, 'y')
  
  input_string_filter = {'A', 'B'};
  while ~any(strcmp(input_string, input_string_filter))
    input_string = input('Extension method? (A/B)', 's');
  end;
  
  model.RB_extension_method = input_string;
  
  input_string_filter = {'y', 'n'};
  while ~any(strcmp(input_string, input_string_filter))
    input_string = input('Sequence mode? (y/n)', 's');
  end;
  
  model.RB_sequence_mode = strcmp(input_string, 'y');
  
  method_index = 1 + strcmp(model.RB_extension_method, 'B') + ...
      2 * (1 - model.RB_sequence_mode);
  
  bases_filename = ['dom_dec_bases_', ...
                    method_strings{method_index}, '.mat'];
  
  fprintf('\nstart basis generation ...\n\n');
  
  tic;
  detailed_data = gen_detailed_data(model, model_data);
  reduced_data = gen_reduced_data(model, detailed_data);
  basis_gen_time = toc;
  
  detailed_data = rmfield(detailed_data, 'reduced_data');
  
  disp(['save data in file "',bases_filename,'"']);
  save(bases_filename, 'detailed_data', 'reduced_data', ...
       'basis_gen_time');
  
  fprintf('\nbasis generation done!\n');
  
  clear('input_string');
  clear('input_string_filter');
  clear('method_index');
  clear('method_strings');
  clear('bases_filename');
  return;

end;

fprintf('\n');


% compare generation methods
% Figure 6.3, Figure 6.4
input_string = [];
while ~any(strcmp(input_string, input_string_filter))
  input_string = input('-> Compare basis generation methods? (y/n)', ...
                       's');
end;

if strcmp(input_string, 'y')

  bases_data = cell(4,1);
  for i = 1:4
    bases_data{i} = load(['dom_dec_bases_',method_strings{i},'.mat']);
  end;
  
  % display computation times
  hours = zeros(4,1);
  minutes = zeros(4,1);
  seconds = zeros(4,1);
  for i = 1:4
    hours(i) = floor(bases_data{i}.basis_gen_time / 3600);
    minutes(i) = floor((bases_data{i}.basis_gen_time - ...
                        hours(i) * 3600) / 60);
    seconds(i) = round(bases_data{i}.basis_gen_time - ...
                       hours(i) * 3600 - minutes(i) * 60);
  end;
  disp('Computation times:');
  for i = 1:4
    disp(['method ', method_strings{i}, ': ', ...
          num2str(hours(i)), 'h ', num2str(minutes(i)), ...
          'm ', num2str(seconds(i)), 's']);
  end;
  fprintf('\n');
  
  % compare error decay
  % Figure 6.3 left
  input_string = [];
  while ~any(strcmp(input_string, input_string_filter))
    input_string = input('--> Compare error decays? (y/n)', 's');
  end;
  
  if strcmp(input_string, 'y')
    figure,
    obj1 = semilogy(...
        bases_data{1}.detailed_data.RB_info.max_errs, ...
        'k', 'LineWidth', 2);
    hold on;
    obj2 = semilogy(...
        bases_data{2}.detailed_data.RB_info.max_errs, ...
        'r', 'LineWidth', 2);
    obj3 = semilogy(...
        bases_data{3}.detailed_data.RB_info.max_errs, ...
        '--k', 'LineWidth', 2);
    obj4 = semilogy(...
        bases_data{4}.detailed_data.RB_info.max_errs, ...
        '--r', 'LineWidth', 2);
    hold off;
    
    set(gca, 'FontSize', 14);
    legend([obj1 obj2 obj3 obj4], ...
           'method A-SEQ', 'method B-SEQ', ...
           'method A-NSEQ', 'method B-NSEQ');
    xlabel('greedy-step');
    ylabel('maximal error estimate');    
    
  end;
  
  % compare final bases sizes
  % Figure 6.3 right
  input_string = [];
  while ~any(strcmp(input_string, input_string_filter))
    input_string = input('--> Compare final bases sizes? (y/n)', 's');
  end;
  
  if strcmp(input_string, 'y')
    final_bases_sizes = zeros(4,1);
    for i = 1:4
      final_bases_sizes(i) = bases_data{i}.reduced_data.N{1} + ...
          bases_data{i}.reduced_data.N{2};
    end;
    
    figure,
    bar(final_bases_sizes);
    
    set(gca, 'FontSize', 14);
    set(gca, 'xTickLabel', method_strings);
    xlabel('method');
    ylabel('final total basis size');
    
  end;
  
  % compare condition numbers
  % Figure 6.4
  input_string = [];
  while ~any(strcmp(input_string, input_string_filter))
    input_string = input(['--> Compare condition numbers of ', ...
                    'Neumann-side system matrices? (y/n)'], 's');
  end;
  
  if strcmp(input_string, 'y')
    old_decomp_mode = model.decomp_mode;
    model.decomp_mode = 2;
    
    dir = model.dirichlet_side;
    no_dir = mod(dir,2)+1;
    
    final_min_mean_max_cond = zeros(4,3);
    max_cond = cell(4,1);
    
    for i = 1:4
      all_cond = ...
          zeros(size(mu_test, 1), ...
                bases_data{i}.detailed_data.RB_info.steps);
      
      for mu_ind = 1:size(mu_test, 1)
        model = set_mu(model, mu_test(mu_ind, :));
        
        [A_coeff, dummy1, dummy2, dummy3] = ...
            model.operators(model, []);
        A = lincomb_sequence(...
            bases_data{i}.reduced_data.AN_comp{no_dir}, ...
            A_coeff{no_dir});
        
        for j = 1:size(all_cond, 2)
          indices = ...
              [bases_data{i}.reduced_data.I_0{no_dir}(...
                  1:bases_data{i}.detailed_data.RB_info ...
                  .bases_sizes.N0{no_dir}(j)) ...
               bases_data{i}.reduced_data.I_G{no_dir}(...
                   1:bases_data{i}.detailed_data.RB_info ...
                   .bases_sizes.NG{no_dir}(j))];
          
          all_cond(mu_ind, j) = cond(A(indices, indices));
          
        end;
        max_cond{i} = max(all_cond);
        
      end;
      
      final_min_mean_max_cond(i,1) = min(all_cond(:, end));
      final_min_mean_max_cond(i,2) = mean(all_cond(:, end));
      final_min_mean_max_cond(i,3) = max(all_cond(:, end));
      
    end;
    
    model.decomp_mode = old_decomp_mode;
    
    % Figure 6.4 left
    figure,
    bar(final_min_mean_max_cond);
    
    set(gca, 'yScale', 'log');
    set(gca, 'FontSize', 14);
    set(gca, 'xTickLabel', method_strings);
    xlabel('method');
    ylabel('cond(A_{N,2}(mu))');
    
    % Figure 6.4 right
    figure,
    obj1 = semilogy(max_cond{1}, 'k', 'LineWidth', 2);
    hold on;
    obj2 = semilogy(max_cond{2}, 'r', 'LineWidth', 2);
    obj3 = semilogy(max_cond{3}, '--k', 'LineWidth', 2);
    obj4 = semilogy(max_cond{4}, '--r', 'LineWidth', 2);
    hold off;
    
    set(gca, 'FontSize', 14);
    legend([obj1 obj2 obj3 obj4], ...
           'method A-SEQ', 'method B-SEQ', ...
           'method A-NSEQ', 'method B-NSEQ', ...
           'Location', 'SouthEast');
    xlabel('greedy-step');
    ylabel('max_mu cond(A_{N,2}(mu))');
    
  end;
    
end;

fprintf('\n\n');


disp('REDUCED SIMULATION');
disp('--------------------------------------------------------');

% settings for simulations
model.det_sim_tol = 1e-10;
model.RB_sim_tol = 1e-10;
model.RB_approximate_thetaN = 0;

% display simulation settings:
disp('Settings for the detailed simulations:');
disp(['tolerance: ', num2str(model.det_sim_tol)]);
disp(['maximal iteration-number: ', num2str(model.maxiter)]);
disp('Settings for the reduced simulations:');
disp(['tolerance: ', num2str(model.RB_sim_tol)]);
disp(['maximal iteration-number: ', ...
      num2str(model.RB_maxiterN)]);
fprintf('\n');


% visualize reduced iteration
% Figure 6.5, Figure 6.6 left
input_string = [];
while ~any(strcmp(input_string, input_string_filter))
  input_string = input('-> Visualize reduced iteration? (y/n)', 's');
end;

if strcmp(input_string, 'y')
  load('dom_dec_bases_B-NSEQ.mat');
  
  model = set_mu(model, [1.3 2.6 2.4 0.56]);
  
  sim_data = detailed_simulation(model, model_data);
  sim_data = model.compute_error(model, model_data, sim_data);
  
  rb_sim_data = rb_simulation(model, reduced_data);
  rb_sim_data = rb_reconstruction(model, detailed_data, ...
                                  rb_sim_data);
  rb_sim_data = model.compute_error(model, detailed_data, ...
                                    rb_sim_data);
  
  plot_params.clim = [0 0.02];
  iteration_indices = [1 5 rb_sim_data.numiter];
  
  % Figure 6.5
  for i = 1:length(iteration_indices)
    for j = 1:2
      rb_sim_data.uh{j}.dofs = ...
          rb_sim_data.all_dofs{j}(:, iteration_indices(i));
    end;
    
    if (iteration_indices(i) == rb_sim_data.numiter)
      disp(['plot ', num2str(i), ': n = n_N(mu) = ', ...
            num2str(iteration_indices(i))]);
    else
      disp(['plot ', num2str(i), ': n = ', ...
            num2str(iteration_indices(i))]);
    end;
    plot_sim_data(model, detailed_data, rb_sim_data, plot_params);
    
    shading flat;
    axis equal;
    axis tight;
    set(gca, 'FontSize', 15);
    
  end;
  
  % Figure 6.6 left
  disp(['plot ', num2str(i+1), ...
        ': comparison with detailed iteration']);
  figure,
  obj1 = semilogy(sim_data.X_err, 'b', 'LineWidth', 2);
  hold on;
  obj2 = semilogy(1:2:rb_sim_data.numiter, ...
                  rb_sim_data.X_err(1:2:end), 'bo', 'LineWidth', 2);
  obj3 = semilogy(1:2:rb_sim_data.numiter, ...
                  rb_sim_data.Delta(1:2:end), 'b.', 'LineWidth', 2);
  hold off;
  
  set(gca, 'FontSize', 13);
  legend([obj1 obj2 obj3], 'error in detailed simulation', ...
         'error in reduced simulation', ...
         'error estimate in reduced simulation');
  xlabel('iteration number');
  
end;

fprintf('\n');


% investigate effectivity of error estimator
% Figure 6.6 right
input_string = [];
while ~any(strcmp(input_string, input_string_filter))
  input_string = input(['-> Investigate effectivity ', ...
                    'of error estimator? (y/n)'], 's');
end;

if strcmp(input_string, 'y')
  load('dom_dec_bases_B-NSEQ.mat');
  model.RB_sim_tol = 1e-06;
  eff = zeros(size(mu_test, 1), 1);
  
  fprintf(['Start reduced simulations of ', ...
        num2str(size(mu_test, 1)), ' test parameters ']);
  
  for mu_ind = 1:size(mu_test, 1)
    model = set_mu(model, mu_test(mu_ind, :));
    
    rb_sim_data = rb_simulation(model, reduced_data);
    rb_sim_data = rb_reconstruction(model, detailed_data, ...
                                    rb_sim_data);
    rb_sim_data = model.compute_error(model, detailed_data, ...
                                      rb_sim_data);
    
    eff(mu_ind) = rb_sim_data.Delta(end) / rb_sim_data.X_err(end);
    fprintf('.');
    
  end;
  
  fprintf('done!\n');
  
  figure,
  hist(eff);
  
  set(gca, 'FontSize', 13);
  xlabel('effectivity of error estimate');
  ylabel('frequency');
  
end;

fprintf('\n');


% compare simulation times
input_string = [];
while ~any(strcmp(input_string, input_string_filter))
  input_string = input(['-> Compare simulation time ', ...
                    'of detailed and reduced simulations? (y/n)'], 's');
end;

if strcmp(input_string, 'y')
  load('dom_dec_bases_B-NSEQ.mat');
  
  det_sim_times = zeros(size(mu_test, 1), 1);
  rb_sim_times = zeros(size(mu_test, 1), 1);
  
  for mu_ind = 1:size(mu_test, 1)
    model = set_mu(model, mu_test(mu_ind, :));
    
    tic;
    sim_data = detailed_simulation(model, model_data);
    det_sim_times(mu_ind) = toc;
    
    tic;
    rb_sim_data = rb_simulation(model, reduced_data);
    rb_sim_data = rb_reconstruction(model, detailed_data, ...
                                    rb_sim_data);
    rb_sim_times(mu_ind) = toc;
    
  end;
  
  disp(['Detailed simulation times range from ', ...
        num2str(min(det_sim_times)), 's to ', ...
        num2str(max(det_sim_times)), 's.']);
  disp(['Reduced simulation times range from ', ...
        num2str(min(rb_sim_times)), 's to ', ...
        num2str(max(rb_sim_times)), 's.']);
  disp(['Mean speed up factor: ', ...
        num2str(mean(det_sim_times./rb_sim_times))]);  
  
end;

fprintf('\n');