CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Wiki > PFV4 program script

PFV4 program script

From CFD-Wiki

Revision as of 19:03, 15 March 2013 by Jonas Holdeman (Talk | contribs)
(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)
Jump to: navigation, search
%LDC4W            LID-DRIVEN CAVITY 
% Finite element solution of the 2D Navier-Stokes equation using 
%   4-node, 6-DOF/node, C1,quartic-derived rectangular stream function 
%   elements for the Lid-Driven Cavity.
% 
% This could be characterized as a VELOCITY-STREAM FUNCTION or 
%   STREAM FUNCTION-VELOCITY method.
%
% Reference:  J. T. Holdeman, "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 NumNx nodes in the x-direction 
%   and NumNy nodes 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:
%    ELS4424r     - class with definitions of velocity shape functions
%    DMatW        - to evaluate element diffusion matrix 
%    CMatW        - to evaluate element convection matrix
%    GetPres3W    - to evaluate the pressure 
%    regrade      - to regrade the mesh 
% Uses
%   ilu          - incomplete LU preconditioner
%   gmres        - iterative solver
% Indirectly uses:
%    Gquad2      - Gauss integraion rules
%    ELG3412r    - class of pressure basis functions
%
% Jonas Holdeman    Updated August 2007, revised March 2013 

clear all;
eu = ELS4424r;       % class of functions for velocity

disp('Lid-driven cavity');
disp([' Four-node, 24 DOF, ' eu.name ' basis.']);

% -------------------------------------------------------------
  nnd = eu.nnodes;        % nnodes = number of element nodes
  ndof = eu.nndofs;       % nndofs = number of velocity dofs per node, (6);
  nd2=ndof*ndof; ND=1:ndof;  % Number of DOF per node - do not change!!
% -------------------------------------------------------------
ETstart=clock;

% Parameters for GMRES solver 
GM.Tol = 1.e-13;
GM.MIter = 20; 
GM.MRstrt = 6;
% parameters for ilu preconditioner
% Decrease su.droptol if ilu preconditioner fails
su.type='ilutp';
su.droptol=1.e-5;

% Suggested relaxation parameters for given Reynolds number 
% Re           100    1000    3200   5000   7500  10000  12500 
% ExpCR1                      .210 
% ExpCRO                      .402 
% RelxFacO:    1.0    1.11    .880   .84    .79    .76     .70
% CritFac:                   1.135

 % Define the problem geometry:
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;   % Set 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 = 18;  NumEy = NumEx;

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

%   Define problem parameters: 
 % Lid velocity (don't change, vary Re instead unless lid moves to left)
Vlid=1.;

 % Specify Reynolds number Re, set relaxation factor
Re=1000.; 

% factor for under/over-relaxation 
RelxFac = 1.;

% Start with simple non-linear iteration, then switch to Newton method 
%   when non-linear corrections are less than ItThreshold. 
% CAUTION! Large Re may require very small threshold. 
% If ItThreshold = 0 (or sufficiently small) method will never switch.
ItThreshold = 5.e-2; 

% NUMBER OF NONLINEAR ITERATIONS 
MaxNLit=12; %

%Selected and labeled contour levels, Ghia, Ghia & Shin for plotting
clSel=[-.1175,-.115,-.11,-.1,-.09,-.07,-.05,-.03,-.01,-.0001,1e-7,1e-5,...
        1e-4,5e-4,1e-3,1.5e-3,3e-3];
LclSel=[1,2,3,4,6,8,10,12,14,16];

%Contour and labeled levels, Ghia, Ghia & Shin, Re=100, 400, 1000, 3200, ...
clGGS=[-1.e-7,-1.e-5,-1.e-4,-.01,-.03,-.05,-.07,-.09,-.1,-.11,-.115,-.1175,...
    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];
LclGGS=[1,3,4,6,8,10,19,23];

% Select contour level option
%CL=clGGS;  CLL=LclGGS;   % GGS
CL=clSel;  CLL=LclSel;   % GGS (select)

if (Vlid<0), CL=-CL; end 

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

 % Viscosity for specified Reynolds number
 nu=Vlid*(Xmax-Xmin)/Re; 

if (xgrd~=1 || ygrd~=1), meshsp='graded'; else meshsp='uniform'; end
disp(['Reynolds number = ' num2str(Re) ',  ' num2str(NumEx) 'x' ...
      num2str(NumEy) ' element ' meshsp ' mesh' ]);
disp(['Maximum number of nonlinear iterations = ' num2str(MaxNLit)]);
pause(1);
  
% 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;
if ygrd ~= 1, YNc=regrade(YNc,ygrd,0); end;

[Xgrid,Ygrid]=meshgrid(XNc,YNc);% Generate the x- and y-coordinate meshes.

% Area-based nodal weights 
wx=zeros(1,NumNx);   wy=zeros(1,NumNy);
wx(1)=.5*(XNc(2)-XNc(1));
wx(2:NumNx-1)=.5*(XNc(3:NumNx)-XNc(1:NumNx-2));
wx(NumNx)=.5*(XNc(NumNx)-XNc(NumNx-1));
wy(1)=.5*(YNc(1)-YNc(2));
wy(2:NumNx-1)=.5*(YNc(1:NumNy-2)-YNc(3:NumNy));
wy(NumNx)=.5*(YNc(NumNy-1)-YNc(NumNy));
Wa=wy'*wx;
 
% Allocate storage for fields 
psi0=zeros(NumNy,NumNx);
u0=zeros(NumNy,NumNx);
v0=zeros(NumNy,NumNx);
pxx0=zeros(NumNy,NumNx);
pxy0=zeros(NumNy,NumNx);
pyy0=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']);

%-------------- End plotting Figure 1 ----------------------

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

NodXY=zeros(NumNod,2);
nn2nft=zeros(NumNod,2); % node number -> nf & nt
NodNdx=zeros(NumNod,2);
% Generate lists of active nodal indices, freedom number & type 
nn=0;  nf=-ndof+1;  nt=1;      %   _______
for nx=1:NumNx                 %  |       |
  for ny=1:NumNy               %  |       |
    nn=nn+1;                   %  |_______|
    NodNdx(nn,:)=[nx,ny];
    NodXY(nn,:)=[Xgrid(ny,nx),Ygrid(ny,nx)];
    nf=nf+ndof;               % all nodes have nd (6) dofs 
    nn2nft(nn,:)=[nf,nt];   % dof number & type (all nodes type 1)
  end;
end;

nf2nnt=zeros(MaxDof,2);  % (node, type) associated with dof
nd=0; dd=1:ndof;
for nn=1:NumNod
  for nf=1:ndof
    nf2nnt(nd+nf,:)=[nn,nf];
  end
  nd=nd+ndof;
end

NumEl=NumEx*NumEy;

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

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

nc=0;
for nf=1:MaxDof
  ni=nf2nnt(nf,1);
  nx=NodNdx(ni,1);  ny=NodNdx(ni,2);
  x=XNc(nx);        y=YNc(ny); 
  if(x==Xmin || x==Xmax)
    nt=nf2nnt(nf,2);
    switch nt;
    case {1, 2, 3, 5, 6}                                % psi,u,v,pxy,pyy 
      nc=nc+1; EBC.typ(nc)=nt; EBC.dof(nc)=nf; EBC.val(nc)=0; 
    end  % switch (type)
  elseif (y==Ymin)
    nt=nf2nnt(nf,2);
    switch nt;
    case {1, 2, 3, 4, 5}                                 % psi,u,v,pxx,pxy 
      nc=nc+1; EBC.typ(nc)=nt; EBC.dof(nc)=nf; EBC.val(nc)=0; 
    end  % switch (type)
  elseif (y==Ymax)
    nt=nf2nnt(nf,2);
    switch nt;
    case {1, 3, 4, 5}                                      % psi,v,pxx,pxy
      nc=nc+1; EBC.typ(nc)=nt; EBC.dof(nc)=nf; EBC.val(nc)=0; 
    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 correction for iteration

% Arrays for convergence norm info
MxNL=max(1,MaxNLit);
np0=zeros(1,MxNL);     % psi
nv0=zeros(1,MxNL);     % u & v 

Qs=[];

% Generate and assemble element matrices

Dmat = spalloc(MaxDof,MaxDof,56*MaxDof);   % to save the diffusion matrix
Vdof=zeros(ndof,nnd);
Xe=zeros(2,nnd);      % coordinates of element corners 

NLitr=0; 
ItType=0;   % Begin with simple iteration 

 % <<<<<<<< BEGIN NONLINEAR ITERATION >>>>>>>>>>>>
while (NLitr<MaxNLit), NLitr=NLitr+1;  
      
tclock=clock;   % Start assembly time <<<<<<<<<
% Generate and assemble element matrices
Mat=spalloc(MaxDof,MaxDof,56*MaxDof);  % (36)
RHS=spalloc(MaxDof,1,MaxDof);

if(ItType==0 && NLitr>1 && nv0(NLitr-1)<ItThreshold && np0(NLitr-1)<ItThreshold) 
  ItType=1; disp(' >>>>> Begin Newton method >>>>>');  % <<<<<<<<<<<<<<<
end   

% BEGIN GLOBAL MATRIX ASSEMBLY
for ne=1:NumEl   
 
  Xe(1:2,1:nnd)=NodXY(Elcon(ne,1:nnd),1:2)'; 
      
  if NLitr == 1     
%     Fluid element diffusion matrix, save on first iteration    
    [Emat,Rndx,Cndx] = DMatW(eu,Xe,Elcon(ne,:),nn2nft);
    Dmat=Dmat+sparse(Rndx,Cndx,Emat,MaxDof,MaxDof);  % Global diff matrix 
  end 

  if (NLitr>1)  
%     Get stream function and velocities, loop over local nodes of element
    for n=1:nnd  
      Vdof(ND,n)=Q((nn2nft(Elcon(ne,n),1)-1)+ND); 
    end
%    Fluid element convection matrix, first iteration uses Stokes equation. 
     if (ItType==0)
        % For simple iteration
       [Emat,Rndx,Cndx] = CMatW(eu,Xe,Elcon(ne,:),nn2nft,Vdof); 
     else
        % For Newton iteration
       [Emat,Rndx,Cndx,Rcm,RcNdx] = CMatW(eu,Xe,Elcon(ne,:),nn2nft,Vdof);
       RHS = RHS + sparse(RcNdx,1,Rcm,MaxDof,1); 
     end
    Mat=Mat+sparse(Rndx,Cndx,Emat,MaxDof,MaxDof);  % Assemble global 
  end   % end  if(NLitr>1)

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, Start solution time  <<<<<<<<<<<<<<
tclock=clock; 

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

[Lm,Um] = ilu(Matr,su);      % incomplete LU
Qr = gmres(Matr,RHSr,GM.MIter,GM.Tol,GM.MRstrt,Lm,Um,Qs);	% GMRES

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(k,1);
  nx=NodNdx(ni,1);  ny=NodNdx(ni,2);
  switch nf2nnt(k,2) % switch on dof type 
    case 1
      dsqp=dsqp+Wa(ny,nx)*(Q(k)-Q0(k))^2;
       errpsi(ny,nx)=Q0(k)-Q(k);  psi0(ny,nx)=Q(k);
    case 2
      dsqv=dsqv+Wa(ny,nx)*(Q(k)-Q0(k))^2;  u0(ny,nx)=Q(k);
    case 3
      dsqv=dsqv+Wa(ny,nx)*(Q(k)-Q0(k))^2;  v0(ny,nx)=Q(k);
    case 4
      pxx0(ny,nx)=Q(k);
    case 5
      pxy0(ny,nx)=Q(k);
    case 6
      pyy0(ny,nx)=Q(k);
  end  % switch on dof type 
end  % for 

np0(NLitr)=sqrt(dsqp); dP=np0(NLitr);
nv0(NLitr)=sqrt(dsqv); 

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

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

% Plot iteration 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');

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

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

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

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(k,1)==PBCnx && NodNdx(k,2)==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

[P,Px,Py] = GetPres3W(eu,NodXY,Elcon,nn2nft,Q,EBCp,0); % nu
P=reshape(P,NumNy,NumNx); 
Px = reshape(Px,NumNy,NumNx); 
Py = reshape(Py,NumNy,NumNx); 

% >>>>>>>>>>>>> End pressure recovery <<<<<<<<<<<<<

%--------------Plot of final result-------------
% ********************** CONTINUE FIGURE 1 *************************

figure(1); 
textid=['4424X' num2str(NumEx)];

% Stream function       (final)
subplot(2,2,3);
[CT,hn]=contour(Xgrid,Ygrid,psi0,CL,'k');  % Plot contours (trajectories)
clabel(CT,hn,CL(CLL));
hold on;
plot([Xmin,Xmin,Xmax,Xmax,Xmin],[Ymax,Ymin,Ymin,Ymax,Ymax],'k');
pcolor(Xgrid,Ygrid,psi0);
shading interp; %flat;  % or interp
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');
%pcolor(Xgrid,Ygrid,P);
%shading interp; %flat; 
hold off;
axis([Xmin,Xmax,Ymin,Ymax]);
axis equal;  axis image;
title(['Quartic pressure contours, Re=' num2str(Re)]);

% ********************* END FIGURE 1 *************************
% ------------ End final plot --------------

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