Code: Lid driven cavity using pressure free velocity form

From CFD-Wiki

Revision as of 04:19, 29 June 2011

Lid-driven cavity using pressure-free velocity formulation

%LDCW            LID-DRIVEN CAVITY 
% Finite element solution of the 2D Navier-Stokes equation using 4-node, 12 DOF,
%  (3-DOF/node), simple-cubic-derived rectangular Hermite basis for 
%   the Lid-Driven Cavity problem.
%
% This could also be characterized as a VELOCITY-STREAM FUNCTION or 
%   STREAM FUNCTION-VELOCITY method.
%
% Reference:  "A Hermite finite element method for incompressible fluid flow", 
%    Int. J. Numer. Meth. Fluids, 64, P376-408 (2010). 
%
% Simplified Wiki version 
% The rectangular problem domain is defined between Cartesian 
%   coordinates Xmin & Xmax and Ymin & Ymax.
% The computational grid has NumEx elements in the x-direction 
%   and NumEy elements in the y-direction. 
% The nodes and elements are numbered column-wise from the  
%   upper left corner to the lower right corner. 
%
%This script calls the user-defined functions:
% regrade      - to regrade the mesh 
% DMatW        - to evaluate element diffusion matrix 
% CMatW        - to evaluate element convection matrix
% GetPresW     - to evaluate the pressure 
% ilu_gmres_with_EBC - to solve the system with essential/Dirichlet BCs 
%
% Jonas Holdeman   August 2007, revised June 2011

clear all;
disp('Lid-driven cavity');
disp(' Four-node, 12 DOF, simple-cubic stream function basis.');

% -------------------------------------------------------------
  nd = 3; nd2=nd*nd;  % Number of DOF per node - do not change!!
% -------------------------------------------------------------
ETstart=clock;

% Parameters for GMRES solver 
GMRES.Tolerance=1.e-14;
GMRES.MaxIterates=20; 
GMRES.MaxRestarts=6;

% Optimal relaxation parameters for given Reynolds number
% (see IJNMF reference)
% Re          100   1000   3200   5000   7500  10000  12500 
% RelxFac:  1.04    1.11   .860   .830   .780   .778   .730 
% ExpCR1    1.488   .524   .192   .0378   --     --     -- 
% ExpCRO    1.624   .596   .390   .331   .243   .163   .133
% CritFac:  1.82    1.49   1.14  1.027   .942   .877   .804 

% Define the problem geometry, set mesh bounds:
Xmin = 0.0; Xmax = 1.0; Ymin = 0.0; Ymax = 1.0; 

% Set mesh grading parameters (set to 1 if no grading).
% See below for explanation of use of parameters. 
xgrd = .75; ygrd=.75;   % (xgrd = 1, ygrd=1 for uniform mesh) 

% Set " RefineBoundary=1 " for additional refinement at boundary, 
%  i.e., split first element along boundary into two. 
RefineBoundary=1; 

%     DEFINE THE MESH  
% Set number of elements in each direction
NumEx = 16;   NumEy = NumEx;

% PLEASE CHANGE OR SET NUMBER OF ELEMENTS TO CHANGE/SET NUMBER OF NODES!
NumNx=NumEx+1;  NumNy=NumEy+1;

%   Define problem parameters: 
 % Lid velocity
Vlid=1.;

 % Reynolds number
Re=1000.; 

% factor for under/over-relaxation starting at iteration RelxStrt 
RelxFac = 1.;  % 

% Number of nonlinear iterations
MaxNLit=10; %

%--------------------------------------------------------

 % Viscosity for specified Reynolds number
 nu=Vlid*(Xmax-Xmin)/Re; 
 
% Grade the mesh spacing if desired, call regrade(x,agrd,e). 
% if e=0: refine both sides, 1: refine upper, 2: refine lower
% if agrd=xgrd|ygrd is the parameter which controls grading, then
%   if agrd=1 then leave array unaltered.
%   if agrd<1 then refine (make finer) towards the ends
%   if agrd>1 then refine (make finer) towards the center.
% 
%  Generate equally-spaced nodal coordinates and refine if desired.
if (RefineBoundary==1)
  XNc=linspace(Xmin,Xmax,NumNx-2); 
  XNc=[XNc(1),(.62*XNc(1)+.38*XNc(2)),XNc(2:end-1),(.38*XNc(end-1)+.62*XNc(end)),XNc(end)];
  YNc=linspace(Ymax,Ymin,NumNy-2); 
  YNc=[YNc(1),(.62*YNc(1)+.38*YNc(2)),YNc(2:end-1),(.38*YNc(end-1)+.62*YNc(end)),YNc(end)];
else
  XNc=linspace(Xmin,Xmax,NumNx); 
  YNc=linspace(Ymax,Ymin,NumNy); 
end
if xgrd ~= 1 XNc=regrade(XNc,xgrd,0); end;  % Refine mesh if desired
if ygrd ~= 1 YNc=regrade(YNc,ygrd,0); end;
[Xgrid,Ygrid]=meshgrid(XNc,YNc);% Generate the x- and y-coordinate meshes.

% Allocate storage for fields 
psi0=zeros(NumNy,NumNx);
u0=zeros(NumNy,NumNx);
v0=zeros(NumNy,NumNx);

%--------------------Begin grid plot-----------------------
% ********************** FIGURE 1 *************************
% Plot the grid 
figure(1);
clf;
orient portrait;  orient tall; 
subplot(2,2,1);
hold on;
plot([Xmax;Xmin],[YNc;YNc],'k');
plot([XNc;XNc],[Ymax;Ymin],'k');
hold off;
axis([Xmin,Xmax,Ymin,Ymax]);
axis equal;
axis image;
title([num2str(NumNx) 'x' num2str(NumNy) ...
      ' node mesh for Lid-driven cavity']);
pause(.1);
%-------------- End plotting Figure 1 ----------------------


%Contour levels, Ghia, Ghia & Shin, Re=100, 400, 1000, 3200, ...
clGGS=[-.1175,-.1150,-.11,-.1,-.09,-.07,-.05,-.03,-.01,-1.e-4,-1.e-5,-1.e-7,-1.e-10,...
      1.e-8,1.e-7,1.e-6,1.e-5,5.e-5,1.e-4,2.5e-4,5.e-4,1.e-3,1.5e-3,3.e-3];
CL=clGGS;   % Select contour level option
if (Vlid<0) CL=-CL; end

NumNod=NumNx*NumNy;     % total number of nodes
MaxDof=nd*NumNod;        % maximum number of degrees of freedom
EBC.Mxdof=nd*NumNod;        % maximum number of degrees of freedom

nn2nft=zeros(2,NumNod); % node number -> nf & nt
NodNdx=zeros(2,NumNod);
% Generate lists of active nodal indices, freedom number & type 
ni=0;  nf=-nd+1;  nt=1;          %   ________
for nx=1:NumNx                   %  |        |
   for ny=1:NumNy                %  |        |
      ni=ni+1;                   %  |________|
      NodNdx(:,ni)=[nx;ny];
      nf=nf+nd;               % all nodes have 4 dofs 
      nn2nft(:,ni)=[nf;nt];   % dof number & type (all nodes type 1)
   end;
end;
%NumNod=ni;     % total number of nodes
nf2nnt=zeros(2,MaxDof);  % (node, type) associated with dof
ndof=0; dd=[1:nd];
for n=1:NumNod
  for k=1:nd
    nf2nnt(:,ndof+k)=[n;k];
  end
  ndof=ndof+nd;
end

NumEl=NumEx*NumEy;

% Generate element connectivity, from upper left to lower right. 
Elcon=zeros(4,NumEl);
ne=0;  LY=NumNy;
for nx=1:NumEx
  for ny=1:NumEy
    ne=ne+1;
    Elcon(1,ne)=1+ny+(nx-1)*LY; 
    Elcon(2,ne)=1+ny+nx*LY;
    Elcon(3,ne)=1+(ny-1)+nx*LY;
    Elcon(4,ne)=1+(ny-1)+(nx-1)*LY;
  end  % loop on ny
end  % loop on nx

% Begin essential boundary conditions, allocate space 
MaxEBC = nd*2*(NumNx+NumNy-2);
EBC.dof=zeros(MaxEBC,1);  % Degree-of-freedom index  
EBC.typ=zeros(MaxEBC,1);  % Dof type (1,2,3)
EBC.val=zeros(MaxEBC,1);  % Dof value 

 X1=XNc(2);  X2=XNc(NumNx-1);
nc=0;
for nf=1:MaxDof
  ni=nf2nnt(1,nf);
  nx=NodNdx(1,ni);
  ny=NodNdx(2,ni);
  x=XNc(nx);
  y=YNc(ny); 
  if(x==Xmin | x==Xmax | y==Ymin)
    nt=nf2nnt(2,nf);
    switch nt;
    case {1, 2, 3}
      nc=nc+1; EBC.typ(nc)=nt; EBC.dof(nc)=nf; EBC.val(nc)=0;  % psi, u, v 
    end  % switch (type)
  elseif (y==Ymax)
    nt=nf2nnt(2,nf);
    switch nt;
    case {1, 3}
      nc=nc+1; EBC.typ(nc)=nt; EBC.dof(nc)=nf; EBC.val(nc)=0;   % psi, v 
    case 2
      nc=nc+1; EBC.typ(nc)=nt; EBC.dof(nc)=nf; EBC.val(nc)=Vlid;   % u
    end  % switch (type) 
  end  % if (boundary)
end  % for nf 
EBC.num=nc; 
  
if (size(EBC.typ,1)>nc)   % Truncate arrays if necessary 
   EBC.typ=EBC.typ(1:nc);
   EBC.dof=EBC.dof(1:nc);
   EBC.val=EBC.val(1:nc);
end     % End ESSENTIAL (Dirichlet) boundary conditions

% partion out essential (Dirichlet) dofs
p_vec = [1:EBC.Mxdof]';         % List of all dofs
EBC.p_vec_undo = zeros(1,EBC.Mxdof);
% form a list of non-diri dofs
EBC.ndro = p_vec(~ismember(p_vec, EBC.dof));	% list of non-diri dofs
% calculate p_vec_undo to restore Q to the original dof ordering
EBC.p_vec_undo([EBC.ndro;EBC.dof]) = [1:EBC.Mxdof]; %p_vec';

Q=zeros(MaxDof,1); % Allocate space for solution (dof) vector

% Initialize fields to boundary conditions
for k=1:EBC.num
   Q(EBC.dof(k))=EBC.val(k); 
end;

errpsi=zeros(NumNy,NumNx);  % error correct for iteration

MxNL=max(1,MaxNLit);
np0=zeros(1,MxNL);     % Arrays for convergence info
nv0=zeros(1,MxNL);

Qs=[];
   
Dmat = spalloc(MaxDof,MaxDof,36*MaxDof);   % to save the diffusion matrix
Vdof=zeros(nd,4);
Xe=zeros(2,4);      % coordinates of element corners 

NLitr=0; ND=1:nd;
while (NLitr<MaxNLit), NLitr=NLitr+1;   % <<< BEGIN NONLINEAR ITERATION 
      
tclock=clock;   % Start assembly time <<<<<<<<<
% Generate and assemble element matrices
Mat=spalloc(MaxDof,MaxDof,36*MaxDof);
RHS=spalloc(MaxDof,1,MaxDof);
%RHS = zeros(MaxDof,1);
Emat=zeros(1,16*nd2);         % Values 144=4*4*(nd*nd) 

% BEGIN GLOBAL MATRIX ASSEMBLY
for ne=1:NumEl   
  
  for k=1:4
     ki=NodNdx(:,Elcon(k,ne));
     Xe(:,k)=[XNc(ki(1));YNc(ki(2))];   
  end
   
  if NLitr == 1    
%     Fluid element diffusion matrix, save on first iteration    
     [DEmat,Rndx,Cndx] = DMatW(Xe,Elcon(:,ne),nn2nft);
     Dmat=Dmat+sparse(Rndx,Cndx,DEmat,MaxDof,MaxDof);  % Global diffusion mat 
   end 
   
   if (NLitr>1) 
%    Get stream function and velocities
     for n=1:4  
       Vdof(ND,n)=Q((nn2nft(1,Elcon(n,ne))-1)+ND); % Loop over local element nodes
     end
%     Fluid element convection matrix, first iteration uses Stokes equation. 
       [Emat,Rndx,Cndx] = CMatW(Xe,Elcon(:,ne),nn2nft,Vdof);  
      Mat=Mat+sparse(Rndx,Cndx,-Emat,MaxDof,MaxDof);  % Global convection assembly 
   end

end;  % loop ne over elements 
% END GLOBAL MATRIX ASSEMBLY

Mat = Mat -nu*Dmat;    % Add in cached/saved global diffusion matrix 

disp(['(' num2str(NLitr) ') Matrix assembly complete, elapsed time = '...
      num2str(etime(clock,tclock)) ' sec']);  % Assembly time <<<<<<<<<<<
pause(1);

Q0 = Q;  % Save dof values 

% Solve system
tclock=clock; %disp('start solution'); % Start solution time  <<<<<<<<<<<<<<

RHSr=RHS(EBC.ndro)-Mat(EBC.ndro,EBC.dof)*EBC.val;
Matr=Mat(EBC.ndro,EBC.ndro);
Qs=Q(EBC.ndro);

Qr=ilu_gmres_with_EBC(Matr,RHSr,[],GMRES,Qs);

Q=[Qr;EBC.val];        % Augment active dofs with esential (Dirichlet) dofs
Q=Q(EBC.p_vec_undo);   % Restore natural order
   
stime=etime(clock,tclock); % Solution time <<<<<<<<<<<<<<

% ****** APPLY RELAXATION FACTOR *********************
if(NLitr>1) Q=RelxFac*Q+(1-RelxFac)*Q0; end
% ****************************************************

% Compute change and copy dofs to field arrays
dsqp=0; dsqv=0;
for k=1:MaxDof
  ni=nf2nnt(1,k); nx=NodNdx(1,ni); ny=NodNdx(2,ni);
  switch nf2nnt(2,k) % switch on dof type 
    case 1
      dsqp=dsqp+(Q(k)-Q0(k))^2; psi0(ny,nx)=Q(k);
      errpsi(ny,nx)=Q0(k)-Q(k);  
    case 2
      dsqv=dsqv+(Q(k)-Q0(k))^2; u0(ny,nx)=Q(k);
    case 3
      dsqv=dsqv+(Q(k)-Q0(k))^2; v0(ny,nx)=Q(k);
  end  % switch on dof type 
end  % for 
np0(NLitr)=sqrt(dsqp); 
nv0(NLitr)=sqrt(dsqv); 

if (np0(NLitr)<=1e-15|nv0(NLitr)<=1e-15) 
  MaxNLit=NLitr; np0=np0(1:MaxNLit); nv0=nv0(1:MaxNLit);   end;
disp(['(' num2str(NLitr) ') Solution time for linear system = '...
     num2str(etime(clock,tclock)) ' sec']); % Solution time <<<<<<<<<<<<
 
%---------- Begin plot of intermediate results ----------
% ********************** FIGURE 2 *************************
figure(1);

% Stream function (intermediate) 
subplot(2,2,3);
contour(Xgrid,Ygrid,psi0,8,'k');  % Plot contours (trajectories)
axis([Xmin,Xmax,Ymin,Ymax]);
title(['Lid-driven cavity,  Re=' num2str(Re)]);
axis equal; axis image;

% Plot convergence info 
subplot(2,2,2);
semilogy(1:NLitr,nv0(1:NLitr),'k-+',1:NLitr,np0(1:NLitr),'k-o');
xlabel('Nonlinear iteration number');
ylabel('Nonlinear correction');
axis square; 
title(['Iteration conv.,  Re=' num2str(Re)]);
legend('U','Psi');

% Plot nonlinear iteration correction contours 
subplot(2,2,4);
contour(Xgrid,Ygrid,errpsi,8,'k');  % Plot contours (trajectories)
axis([Xmin,Xmax,Ymin,Ymax]);
axis equal; axis image;
title(['Iteration correction']);
pause(1);
% ********************** END FIGURE 2 *************************
%----------  End plot of intermediate results  ---------

if (nv0(NLitr)<1e-15) break; end  % Terminate iteration if non-significant 

end;   % <<< (while) END NONLINEAR ITERATION

format short g;
disp('Convergence results by iteration: velocity, stream function');
disp(['nv0:  ' num2str(nv0)]); disp(['np0:  ' num2str(np0)]); 

% >>>>>>>>>>>>>> BEGIN PRESSURE RECOVERY <<<<<<<<<<<<<<<<<<
% Essential pressure boundary condition 
% Index of node to apply pressure BC, value at node
PBCnx=fix((NumNx+1)/2);   % Apply at center of mesh
PBCny=fix((NumNy+1)/2);
PBCnod=0;
for k=1:NumNod
  if (NodNdx(1,k)==PBCnx & NodNdx(2,k)==PBCny) PBCnod=k; break; end
end
if (PBCnod==0) error('Pressure BC node not found'); 
else
  EBCp.nodn = [PBCnod];  % Pressure BC node number
  EBCp.val = [0];  % set P = 0.
end
% Cubic pressure 
[P,Px,Py] = GetPresW(NumNod,NodNdx,Elcon,nn2nft,Xgrid,Ygrid,Q,EBCp,nu);
% ******************** END PRESSURE RECOVERY *********************

% ********************** CONTINUE FIGURE 1 *************************
figure(1);

% Stream function    (final)
subplot(2,2,3);
[CT,hn]=contour(Xgrid,Ygrid,psi0,CL,'k');  % Plot contours (trajectories)
clabel(CT,hn,CL([1,3,5,7,9,10,11,19,23]));
hold on;
plot([Xmin,Xmin,Xmax,Xmax,Xmin],[Ymax,Ymin,Ymin,Ymax,Ymax],'k');
hold off;
axis([Xmin,Xmax,Ymin,Ymax]);
axis equal;  axis image;
title(['Stream lines, ' num2str(NumNx) 'x' num2str(NumNy) ...
    ' mesh, Re=' num2str(Re)]);

% Plot pressure contours   (final)
subplot(2,2,4);
CPL=[-.002,0,.02,.05,.07,.09,.11,.12,.17,.3];
[CT,hn]=contour(Xgrid,Ygrid,P,CPL,'k');  % Plot pressure contours
clabel(CT,hn,CPL([3,5,7,10]));
hold on;
plot([Xmin,Xmin,Xmax,Xmax,Xmin],[Ymax,Ymin,Ymin,Ymax,Ymax],'k');
hold off;
axis([Xmin,Xmax,Ymin,Ymax]);
axis equal;  axis image;
title(['Simple cubic pressure contours, Re=' num2str(Re)]);
% ********************* END FIGURE 1 *************************

disp(['Total elapsed time = '...
   num2str(etime(clock,ETstart)/60) ' min']); % Elapsed time from start <<<

Diffusion matrix for pressure-free velocity method

Convection matrix for pressure-free velocity method

Consistent pressure for pressure-free velocity method

GMRES solver with ILU preconditioning and Essential BC for pressure-free velocity method

grade node spacing for pressure-free velocity method