From CFD-Wiki
%LDC3W 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:
% ELS3412r - class of velocity basis functions
% DMatW - to evaluate element diffusion matrix
% CMatW - to evaluate element convection matrix
% GetPresW - to evaluate the pressure
% regrade - to regrade the mesh
% Uses
% ilu - incomplete LU preconditioner
% gmres - iterative solver
% Indirectly uses:
% Gquad2 - Gauss integraion rules for rectangle
% ELG3412r - class of pressure basis functions
%
% Jonas Holdeman August 2007, revised March 2013
clear all;
eu = ELS3412r; % class of functions for velocity
disp('Lid-driven cavity');
disp([' Four-node, 12 DOF, ' eu.name ' basis.']);
% ------------------------------------------------------------------------
ndf = eu.nndofs; % nndofs = number of velocity dofs per node, (3);
nnd = eu.nnodes; % nnodes = number of element nodes
nd2=ndf*ndf; % Number of DOF per node - do not change!!
% ------------------------------------------------------------------------
ETstart=clock;
% Parameters for GMRES solver
GM.Tol=1.e-12;
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;
% 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 = 18; 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.; %
% 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 = .1; %1.e-2;
% Number of nonlinear iterations
MaxNLit=12; %20; %
%--------------------------------------------------------
% 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']);
%-------------- 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=ndf*NumNod; % maximum number of degrees of freedom
EBC.Mxdof=ndf*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
ni=0; nf=-ndf+1; nt=1; % ________
for nx=1:NumNx % | |
for ny=1:NumNy % | |
ni=ni+1; % |________|
NodNdx(ni,:)=[nx,ny];
NodXY(ni,:)=[Xgrid(ny,nx),Ygrid(ny,nx)];
nf=nf+ndf; % 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(MaxDof,2); % (node, type) associated with dof
ndof=0; dd = 1:ndf;
for n=1:NumNod
for k=1:ndf
nf2nnt(ndof+k,:)=[n,k];
end
ndof=ndof+ndf;
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 = ndf*2*(NumNx+NumNy-2);
EBC.dof=zeros(MaxEBC,1); % Degree-of-freedom index
EBC.val=zeros(MaxEBC,1); % Dof value
X1=XNc(2); X2=XNc(NumNx-1);
nc=0;
for nf=1:MaxDof
ni=nf2nnt(nf,1);
x=NodXY(ni,1); y=NodXY(ni,2);
if(x==Xmin || x==Xmax || y==Ymin)
nt=nf2nnt(nf,2);
switch nt;
case {1, 2, 3}
nc=nc+1; EBC.dof(nc)=nf; EBC.val(nc)=0; % psi, u, v
end % switch (type)
elseif (y==Ymax)
nt=nf2nnt(nf,2);
switch nt;
case {1, 3}
nc=nc+1; EBC.dof(nc)=nf; EBC.val(nc)=0; % psi, v
case 2
nc=nc+1; EBC.dof(nc)=nf; EBC.val(nc)=Vlid; % u
end % switch (type)
end % if (boundary)
end % for nf
EBC.num=nc;
if (size(EBC.dof,1)>nc) % Truncate arrays if necessary
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(ndf,nnd);
Xe=zeros(2,nnd); % coordinates of element corners
ItType=0;
NLitr=0; ND=1:ndf;
while (NLitr<MaxNLit), NLitr=NLitr+1; % <<< BEGIN NONLINEAR ITERATION
% <<<<<<Newton?<<<<<<<<<<<<<<*******
if(ItType==0 && NLitr>1 && nv0(NLitr-1)<ItThreshold && ...
np0(NLitr-1)<ItThreshold)
ItType=1; disp(' >>> Begin Newton method >>>');
end
tclock=clock; % Start assembly time <<<<<<<<<
% Generate and assemble element matrices
Mat=spalloc(MaxDof,MaxDof,36*MaxDof);
RHS=spalloc(MaxDof,1,MaxDof);
% 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 diffusion mat
end
if (NLitr>1)
% Get stream function and velocities, loop over local element nodes
for n=1:nnd
Vdof(ND,n)=Q((nn2nft(Elcon(ne,n))-1)+ND,1);
end
% Fluid element convection matrix, first iteration uses Stokes equation
if ItType==0
% Convection term for simple iteration
[Emat,Rndx,Cndx] = CMatW(eu,Xe,Elcon(ne,:),nn2nft,Vdof);
else
% Convection term for Newton iteration
[Emat,Rndx,Cndx,Rcm,RcNdx] = CMatW(eu,Xe,Elcon(ne,:),nn2nft,Vdof);
RHS = RHS + sparse(RcNdx,1,Rcm,MaxDof,1);
end % if ItType...
% Global convection assembly
Mat=Mat+sparse(Rndx,Cndx,Emat,MaxDof,MaxDof);
end % if NLitr...
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);
[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+(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); 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;
disp(['(' num2str(NLitr) ') Solution time for linear system = '...
num2str(etime(clock,tclock)) ' sec,' ' dV = ' num2str(nv0(NLitr)) ...
', dP = ' num2str(np0(NLitr))]); % 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','Location','SouthWest');
% Plot nonlinear iteration correction contours
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 % if dP
% ********************** 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(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
% Cubic pressure
[P,Px,Py] = GetPres3W(eu,NodXY,Elcon,nn2nft,Q,EBCp,0);
% Copy pressure to structured mesh for plotting
P =reshape(P,NumNy,NumNx);
Px=reshape(Px,NumNy,NumNx);
Py=reshape(Py,NumNy,NumNx);
% ******************** 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');
pcolor(Xgrid,Ygrid,psi0);
shading interp %flat;
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(['Simple cubic pressure contours, Re=' num2str(Re)]);
% ********************* END FIGURE 1 *************************
disp(['Total elapsed time = '...
num2str(etime(clock,ETstart)/60) ' min']); % Elapsed time from start <<<
beep; pause(.25); beep; pause(.25); beep;
