PFV4 program script
From CFD-Wiki
Revision as of 19:03, 15 March 2013 by Jonas Holdeman (Talk | contribs)
%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;