PFV GMRES solver

From CFD-Wiki

GMRES solver ilu_gmres_with_EBC.m with ILU preconditioner and essential BC for pressure-free velocity method.

function Q = ilu_gmres_with_EBC(mat,rhs,EBC,GMRES,Q0,DropTol)
% ILU_GMRES_WITH_EBCx  Solves matrix equation mat*Q = rhs.
% 
% Solves the matrix equation mat*Q = rhs, optionally constrained 
% by Dirichlet boundary conditions described in diri_list, using 
% Matlab's preconditioned gmres sparse solver.  When Dirichlet 
% boundary conditions are provided, the routine enforces them by 
% reordering to partition out Dirichlet degrees of freedom. 
% usage:  Q = ilu_gmres_with_EBC(mat,rhs,EBC,GMRES,Q0,DropTol)
%    or:  Q = ilu_gmres_with_EBC(mat,rhs,EBC,GMRES,Q0)
%    or:  Q = ilu_gmres_with_EBC(mat,rhs,EBC,GMRES)
%    or:  Q = ilu_gmres_with_EBC(mat,rhs,EBC)
%    or:  Q = ilu_gmres_with_EBC(mat,rhs)
%
%  mat  - matrix of linear system to be solved. 
%  rhs  - right hand side of linear system.
%  EBC.dof, EBC.val - (optional) list of Dirichlet boundary 
%          conditions (constraints). May be empty ([]).
%  GMRES - Structure specifying tolerance, max iterations and
%          restarts. Use [] for default values.
%  Q0   - (optional) initial approximation to solution for restart, 
%          may be empty ([]).
%  DropTol - drop tolerance for luinc preconditioner (default=1e-6).
%
% The solution Q is reported with the original ordering restored. 
%
% If specified and not empty, diri_list should have two columns 
%   total and one row for each diri degree of freedom.  The first 
%   column of each row must contain global index of the degree of 
%   freedom.  The second column contains the actual Dirichlet value. 				 
%   If there are no Dirichlet nodes, omit this parameter if not using
%   the restart capabilities, or supply empty matrix ([]) as diri_list.
%
% Nominal values for GMRES for a large difficult problem might be: 
%  GMRES.Tolerance  = 1.e-12, 
%  GMRES.MaxIterates = 75, 
%  GMRES.MaxRestarts = 14.
%
% Jonas Holdeman,   revised February, 2009. 

% Drop tolerance for luinc preconditioner, nominal value - 1.e-6
if (nargin<6 | isempty(DropTol) | DropTol<=0) 
  droptol = 1.e-6;        % default
else droptol = DropTol;   % assigned
end

if nargin<=3 | isempty(GMRES)
   Tol = 1.e-12;  MaxIter = 75;  MaxRstrt = 14;
else
% Tol=tolerance for residual, increase it if solution takes too long.
   Tol=GMRES.Tolerance;
% MaxIter = maximum number of iterations before restart (MT used 5)
   MaxIter=GMRES.MaxIterates;
% MaxRstrt = maximum number of restarts before giving up (MT used 10)
   MaxRstrt=GMRES.MaxRestarts;
end

% check arguments for reasonableness
tdof = size(mat,1);	 % total number of degrees of freedom

% good mat is a square tdof by tdof matrix
if (size(mat,2)~=size(mat,1) | size(size(mat),2)~=2)
  error('mat must be a square matrix')
end

% valid rhs has the dimensions [tdof, 1]
if (size(rhs,1)~=tdof | size(rhs,2)~=1)
  error('rhs must be a column matrix with the same number of rows as mat')
end

% valid dimensions for optional diri_list
if nargin<=2
   EBC=[];
elseif ~isempty(EBC) & (size(EBC.val,1)>=tdof ...
      | size(EBC.dof,1)>=tdof)
   error('check dimensions of EBC')
end
% (optional) valid Q0 is empty or has the dimensions [tdof, 1]
if nargin<=4
   Q0=[];
elseif ~isempty(Q0) & (size(Q0,1)~=tdof | size(rhs,2)~=1)
  error('Q0 must be a column matrix with the same number of rows as mat')
end

% handle the case of no Dirichlet dofs separately
if isempty(EBC)
% skip diri partitioning, solve the system
  [L,U] = luinc(mat,droptol);      % incomplete LU
  Q = gmres(mat,rhs,MaxIter,Tol,MaxRstrt,L,U,Q0);	% GMRES
  
else
% Form list of all DOFs   
  p_vec = [1:tdof]';
% partion out diri dofs
  EBCdofs = EBC.dof(:,1);	 % list of dofs which are Dirichlet
  EBCvals = EBC.val(:,1);  % Dirichlet dof values
  
% form a list of non-diri dofs
  ndro = p_vec(~ismember(p_vec, EBCdofs));	% list of non-diri dofs
  
% Move Dirichlet DOFs to right side
  rhs_reduced = rhs(ndro) - mat(ndro, EBCdofs) * EBCvals;
  
% solve the reduced system (preconditioned gmres)
  A = mat(ndro,ndro);
   
% Compute incomplete LU preconditioner
  [L,U] = luinc(A,droptol);      % incomplete LU
   
% Remove Dirichlet DOFs from initial estimate
  if ~isempty(Q0)  Q0=Q0(ndro);  end
 
% solve the reduced system (preconditioned gmres)
  Q_reduced = gmres(A,rhs_reduced,MaxIter,Tol,MaxRstrt,L,U,Q0);
  
% insert the Dirichlet values into the solution
  Q = [Q_reduced; EBCvals];

% calculate p_vec_undo to restore Q to the original dof ordering
  p_vec_undo = zeros(1,tdof);
  p_vec_undo([ndro;EBCdofs]) = [1:tdof];
  
% restore the original ordering
  Q = Q(p_vec_undo);
end