dom_dec_detailed_simulation.m

Go to the documentation of this file.

function sim_data = dom_dec_detailed_simulation(model,model_data)
% function sim_data = dom_dec_detailed_simulation(model,model_data)
% the detailed solution is computed via the iterative
% Dirichlet-Neumann schema. 
%
% Return values:
%  sim_data: struct containing the solution and further information
%
% Generated fields of sim_data:
%  uh: cell of discrete functions, which contains the last
%  iteratives
%  all_dofs: cell of matrizes containing the dofs of all iteratives
%  theta: the last relaxation-parameter
%  numiter: number of iterations
%  reached_maxiter: indicates the cause of termination
%  reached_tolerance: indicates the cause of termination

% I.Maier, 31.10.2011

model.decomp_mode = 0;
dir = model.dirichlet_side;
no_dir = mod(dir,2)+1;
  
sim_data = [];
sim_data.reached_maxiter = 0;
sim_data.reached_tolerance = 0;

if model.verbose >= 1
  fprintf('start detailed simulation ');
end;

% compute required matrizes and vectors
[A, r, A_gamma, r_gamma] = model.operators(model,model_data);

% initialize iterative schema
uh = cell(1,2);
uh_old = cell(1,2);
uh_tilde = cell(1,2);
zh = cell(1,2);
for i = 1:2
  uh{i} = femdiscfunc([],model_data.df_infos{i});
  uh_old{i} = zeros(uh{i}.ndofs,1);
  uh_tilde{i} = zeros(uh{i}.ndofs,1);
  zh{i} = zeros(uh{i}.ndofs,1);
end;
a = [];
sigma = 0;
tau = 0;
theta = [];

% initialize output
sim_data.uh = cell(1,2);
for i = 1:2
  sim_data.uh{i} = uh{i};
end;
sim_data.all_dofs = cell(1,2);
for i = 1:2
  sim_data.all_dofs{i} = zeros(uh{i}.ndofs,model.maxiter);
end;

% step 0
uh{dir}.dofs = A_gamma{dir} \ r_gamma{dir};
uh{no_dir}.dofs = A{no_dir} \ r{no_dir};
uh_0_no_dir = uh{no_dir}.dofs;

for i = 1:2
  r_gamma{i} = zeros(size(r_gamma{i}));
  r_gamma{i}(model_data.gamma_dofs{i}) = ...
      uh{i}.dofs(model_data.gamma_dofs{i});
  uh_tilde{i} = A_gamma{i} \ r_gamma{i};
end;

% delete matrizes not needed any longer
A_gamma_dir = A_gamma{dir};
clear('A_gamma');
r_dir = r{dir};
clear('r');

% step >=1
continue_loop = 1;
numiter = 0;
diff_norm = zeros(1,2);
while continue_loop
  
  % ...
  numiter = numiter+1;
  for i = 1:2
    uh_old{i} = uh{i}.dofs;
  end;
  
  % ...
  r_gamma{dir}(model_data.gamma_dofs{dir}) = ...
      uh_tilde{no_dir}(model_data.gamma_dofs{no_dir});
  zh{dir} = A_gamma_dir \ r_gamma{dir};
  
  % ...
  a = (zh{dir}'*A{dir}*zh{dir}) / (uh_tilde{no_dir}'*A{no_dir}* ...
                                   uh_tilde{no_dir});
  sigma = max(sigma,a);
  tau = max(tau,1/a);
  theta = (tau+1)/(sigma^2*tau+tau+2)/2*model.relaxation_factor;
  
  % ...
  uh{dir}.dofs = uh{dir}.dofs + theta*(zh{dir}-uh_tilde{dir});
  uh_tilde{dir} = uh_tilde{dir} + theta*(zh{dir}-uh_tilde{dir});
  
  % ...
  zh{dir} = r_dir - A{dir}*uh{dir}.dofs;
  r_gamma{no_dir}(model_data.gamma_dofs{no_dir}) = ...
      zh{dir}(model_data.gamma_dofs{dir});
  zh{no_dir} = A{no_dir} \ r_gamma{no_dir};
  
  % ...
  uh{no_dir}.dofs = uh_0_no_dir + zh{no_dir};
  uh_tilde{no_dir} = uh_tilde{no_dir} + uh{no_dir}.dofs - ...
      uh_old{no_dir};
      
  % save dofs
  for i = 1:2
    sim_data.all_dofs{i}(:,numiter) = uh{i}.dofs;
  end;
  
  % stopping criteria
  for i = 1:2
    diff_norm(i) = sqrt((uh{i}.dofs-uh_old{i})' * A{i} * ...
                        (uh{i}.dofs-uh_old{i}));
  end;
  if sum(diff_norm) < model.det_sim_tol
    continue_loop = 0;
    sim_data.reached_tolerance = 1;
  end;
  if numiter == model.maxiter
    continue_loop = 0;
    sim_data.reached_maxiter = 1;
  end;

  if (mod(numiter,100) == 1) && (model.verbose >= 1)
    fprintf('.');
  end;
  
end;

% further information
sim_data.theta = theta;
sim_data.numiter = numiter;
for i = 1:2
  sim_data.all_dofs{i} = sim_data.all_dofs{i}(:,1:numiter);
end;

% compute output
if model.base_model.compute_output_functional
    s = 0;
    for i = 1:2
        model.base_model.decomp_mode = model.decomp_mode;
        
        model_data.df_info = model_data.df_infos{i};
        v = model.operators_output(model.base_model,model_data);
        s = s + (v(:)') * sim_data.uh{i}.dofs;
    end
    sim_data.s = s;
end

if model.verbose >= 1
  fprintf('\n');
end;