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

PFVT thermal program script

From CFD-Wiki

Revision as of 20:48, 15 March 2013 by Jonas Holdeman (Talk | contribs)
(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)
Jump to: navigation, search
%TC44SW            THERMAL-DRIVEN CAVITY 
%
% Finite element solution of the 2D Navier-Stokes equation for the 
%  bouyancy-driven thermal cavity problem using quartic Hermite elements 
%   and SEGMENTED SOLUTIONS (T-V split).
% 
% This could be characterized as a VELOCITY-STREAM FUNCTION or 
%   STREAM FUNCTION-VELOCITY method.
%
% Reference:  "Computation of incompressible thermal flows using Hermite finite 
%    elements", Comput. Methods Appl. Mech. Engrg, 199, P3297-3304 (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. 
% Segregated version 
%
%This script direcly calls the user-defined functions:
%  ELS4424r      - class with definitions of velocity shape functions
%  ELG3412r      - class with definitions of temperature & pressure fns
%  regrade       - to regrade the mesh (optional)
%  DMatW         - to evaluate element velocity diffusion matrix 
%  CMatW         - to evaluate element velocity convection matrix
%  TDMatSW       - to evaluate element thermal diffusion matrix 
%  TCMatSW       - to evaluate element thermal convection matrix
%  BMatSW        - to evaluate element bouyancy matrix
%  GetPres3W     - to evaluate the pressure 
%
% Uses:
%   ilu          - incomplete LU preconditioner for solver
%   gmres        - sparse linear solver to solve the sparse system 
% Indirectly uses:
%   Gquad2       - Gauss integraion rules
%   ELG3412r     - class of pressure basis functions
%
% Jonas Holdeman,   December 2008   revised March 2013

clear all;
eu = ELS4424r;    % class of functions for velocity
et = ELG3412r;    % class fo functions for emperaure

disp('Bouyancy-driven thermal cavity, T-V split');
disp([' Four-node, 24 dof ' eu.name 'basis, 12 dof cubic thermal basis.']);

% ------------------ Fixed constants -------------------------- 
nnd = eu.nnodes;      % nnodes = number of element nodes
nV = eu.nndofs;       % nndofs = number of velocity dofs per node, (6); 
nT = et.nndofs;       % nndofs = number of temperature dofs per node, (3); 
nn = eu.nn;           % defines local nodal order, [-1 -1; 1 -1; 1 1; -1 1]
nV2 = nV*nV;    % Number velocity DOFs per node squared. 
nT2 = nT*nT;    % Number temperature DOFs per node squared.
nD = nV+nT; nD2=nD*nD; % Total number of DOFs per node
% -------------------------------------------------------------
% Tolerance parameters for the GMRES iterative sparse solver
GM.VTol = 2.e-12;    % convergence tolerance for velociy solutions
GM.TTol = 1.e-13;    % convegence tolerance for temperature solutions
GM.MIter = 30;       % number of iterations allowed
GM.MRstrt = 6;       % number of restarts
% parameters for ilu preconditioner
% Decrease su.droptol if ilu preconditioner fails
su.type='ilutp';
su.droptol=1.e-5;
% -------------------------------------------------------------

% Fixed parameters: 
% Temperature on left, right side
TL=1;    TR=0;
% Prandtl number
Pr=.71;
% Dimensionless equation form to be used 
% Use EquationForm=1 for large Ra, EquationForm=0 for Ra->0 
EquationForm=1;   % (Choose 0 or 1)

ETstart=clock;

% Define the problem geometry:
Xmin = 0.0; Xmax = 1.0;    Ymin = 0.0; Ymax = 1.0;  

% Mesh grading parameters
xgrd = 1.0; ygrd=1.0; %
%xgrd = .75; ygrd=.75; %

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

%  DEFINE THE MESH: 
% Number of elements in the x-direction
NumEx = 18; 
% Number of elements in the y-direction
NumEy = 18; 

NumEl = NumEx*NumEy;

% Number of nodes
NumNx = NumEx+1;  NumNy = NumEy+1;
NumNod = NumNx*NumNy;

% --------------------------------------
% Suggested relaxation factors for steady flow
% Ra         10^3    10^4    10^5    10^6
% RelxFac   .94     .67     .35     .04

% Rayleigh number Ra
  Ra = 10^4;
% Relaxation factor ( <= 1 )
  RelxFac=.67; 

% Number of nonlinear iterations
MaxNLit=20; 
  
if (xgrd~=1 || ygrd~=1), meshsp='graded'; else meshsp='uniform'; end
disp(['Rayleigh number = ' num2str(Ra) ',  ' num2str(NumEx) 'x' ...
       num2str(NumEy) ' element ' meshsp ' mesh' ]);
disp(['Maximum number of nonlinear iterations = ' num2str(MaxNLit)]);
disp(' ');
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 towards the ends
%   if agrd<1 then refine 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),.5*(YNc(1)+YNc(2)),YNc(2:end-1),(YNc(end-1)+YNc(end))/2.,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.
 
% 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:NumNy-1)=.5*(YNc(1:NumNy-2)-YNc(3:NumNy));
wy(NumNy)=.5*(YNc(NumNy-1)-YNc(NumNy));
Wa=wy'*wx;

% Initial stream function and velocity, temperature & gradient 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);
T0=zeros(NumNy,NumNx);
Tx0=zeros(NumNy,NumNx);
Ty0=zeros(NumNy,NumNx);

%------------------Begin grid plot--------------------
% Plot the grid 
figure(1);
clf;
orient portrait;  orient tall;
subplot(2,2,2);
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 thermal cavity']);

%-------------------End grid plot---------------------

MaxVdof=nV*NumNod;
MaxTdof=nT*NumNod;
MaxDof=MaxVdof+MaxTdof;        % maximum number of degrees of freedom
VBC.Mxdof=MaxVdof;             % maximum number V of degrees of freedom
TBC.Mxdof=MaxTdof;             % maximum number of T degrees of freedom

NodXY =zeros(NumNod,2);
nn2vft=zeros(NumNod,2); % node number -> nfv & nt
nn2tft=zeros(NumNod,2); % node number -> nft & nt
NodNdx=zeros(NumNod,2);
% Generate lists of active nodal indices, freedom number & type 
ni=0;   nt=1; 
nfv=-nV+1;  nft=-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)];
      nfv=nfv+nV;               % all V nodes have 6 dofs 
      nn2vft(ni,:)=[nfv,nt];    % dof number & type (all nodes type 1)
      nft=nft+nT;               % all T nodes have 3 dofs 
      nn2tft(ni,:)=[nft,nt];    % dof number & type (all nodes type 1)
   end;
end;
%NumNod=ni;     % total number of nodes

nvf2nnt=zeros(MaxVdof,2);    % (node, type) associated with dof
ntf2nnt=zeros(MaxTdof,2);    % (node, type) associated with dof
nfv=0; nft=0;
for n=1:NumNod
  for k=1:nV
    nvf2nnt(nfv+k,:)=[n,k];
  end
  nfv=nfv+nV;
  for k=1:nT
    ntf2nnt(nft+k,:)=[n,k];
  end
  nft=nft+nT;
end

% Generate element connectivity, from upper left to lower right. 
Elcon=zeros(NumEl,4);
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 (Dirichlet) V boundary conditions
MaxVBC=2*(5*NumNy+5*(NumNx-2));  % Allocate space for VBC
VBC.dof=zeros(MaxVBC,1);         % Degree-of-freedom index  
VBC.val=zeros(MaxVBC,1);         % Dof value 
nc=0;
for nf=1:MaxVdof
  nt=nvf2nnt(nf,2);  ni=nvf2nnt(nf,1);
  nx=NodNdx(ni,1);   ny=NodNdx(ni,2);
  x=XNc(nx);         y=YNc(ny); 
  if((x==Xmin || x==Xmax) && (y==Ymax || y==Ymin))    % corner 
      nc=nc+1; VBC.dof(nc)=nf; VBC.val(nc)=0;   % psi, u, v, pxx, pxy, pyy
  elseif(x==Xmin || x==Xmax)        % left or right boundary 
    switch nt;
    case {1, 2, 3, 5, 6}            % psi, u, v, pxy, pyy (pxx=-v_x free)
      nc=nc+1; VBC.dof(nc)=nf; VBC.val(nc)=0; 
    end  % switch (type)
  elseif (y==Ymax || y==Ymin)       % top, bottom
      switch nt
      case {1, 2, 3, 4, 5}          % psi, u, v, pxx, pxy (pyy=u_y free)
        nc=nc+1; VBC.dof(nc)=nf; VBC.val(nc)=0;  
      end  % switch (type)
  end  % if (boundary)
end  % for nf 
VBC.num=nc; 
if (size(VBC.dof,1)>nc)   % Truncate V arrays if necessary 
   VBC.dof=VBC.dof(1:nc);
   VBC.val=VBC.val(1:nc);
end   % truncate V
% End ESSENTIAL (Dirichlet) V boundary conditions

% Begin ESSENTIAL (Dirichlet) T boundary conditions
MaxTBC=2*(2*NumNy+(NumNx-2));    % Allocate space for TBC
TBC.dof=zeros(MaxTBC,1);         % Degree-of-freedom index  
TBC.val=zeros(MaxTBC,1);         % Dof value 
nc=0;
for nf=1:MaxTdof
  nt=ntf2nnt(nf,2);   ni=ntf2nnt(nf,1);
  nx=NodNdx(ni,1);    ny=NodNdx(ni,2);
  x=XNc(nx);          y=YNc(ny); 
  if(x==Xmin)                     % left boundary 
    switch nt;
    case 1
      nc=nc+1; TBC.dof(nc)=nf; TBC.val(nc)=TL;     % T 
    case 3
      nc=nc+1; TBC.dof(nc)=nf; TBC.val(nc)=0;      % Ty 
    end  % switch (type)
  elseif (x==Xmax)                % right boundary 
    switch nt;
    case 1
      nc=nc+1; TBC.dof(nc)=nf; TBC.val(nc)=TR;      % T 
    case 3
      nc=nc+1; TBC.dof(nc)=nf; TBC.val(nc)=0;       % Ty 
    end  % switch (type)
  elseif (y==Ymax || y==Ymin)       % top, bottom
    switch nt
    case 3
      nc=nc+1; TBC.dof(nc)=nf; TBC.val(nc)=0;       % Ty 
    end  % switch (type)
  end  % if (boundary)
end  % for nf 
TBC.num=nc; 
if (size(TBC.dof,1)>nc)   % Truncate T arrays if necessary 
   TBC.dof=TBC.dof(1:nc);
   TBC.val=TBC.val(1:nc);
end   % truncate T  
% End ESSENTIAL (Dirichlet) T boundary conditions

% Number active DOFs
ADOFs=MaxDof-VBC.num-TBC.num;
disp(['Number of active DOFs = ' num2str(ADOFs)]);

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

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

Qv=zeros(MaxVdof,1); % Allocate space for velocity solution (dof) vector
Qt=zeros(MaxTdof,1); % Allocate space for temperature solution (dof) vector
   
  for k = 1:MaxTdof
    nn=ntf2nnt(k,1);
    nx=NodNdx(nn,1);   ny=NodNdx(nn,2);
    switch ntf2nnt(k,2) % switch on DOF type 
     case 1
        Qt(k)=T0(ny,nx);
      case 2
        Qt(k)=Tx0(ny,nx);
      case 3
        Qt(k)=Ty0(ny,nx);
    end;  % switch (type)
  end   % loop over TDOFs

% Initialize fields to boundary conditions
KB=1:VBC.num;
Qv(VBC.dof(KB))=VBC.val(KB);
KB=1:TBC.num;
Qt(TBC.dof(KB))=TBC.val(KB);

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

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

if EquationForm==1 
   Ctd=1/sqrt(Ra*Pr); Cvd=Ctd*Pr; Cbf=1; % Coefficients for large Ra
else 
   Ctd=1; Cvd=Pr; Cbf=Pr*Ra;             % Coefficients for small Ra
end

DMat = spalloc(MaxVdof,MaxVdof,54*MaxVdof);   % for the diffusion matrix
TDMat = spalloc(MaxTdof,MaxTdof,30*MaxTdof);  % for thermal diffusion mat.
BMat = spalloc(MaxVdof,MaxTdof,18*MaxVdof);   % for the bouyancy matrix
MatV = [];   % Velocity submatrix
MatT = [];   % Temperature submatrix

Vdof=zeros(6,nnd);
Xe=zeros(2,nnd);     % coordinates of element corners 
NLitr=0;
NV=1:nV;
ItType=0;          % Initialy simple iteration

while NLitr<MaxNLit, NLitr=NLitr+1;   % <<< BEGIN NONLINEAR ITERATION 
      
% Generate and assemble element matrices
CMat = spalloc(MaxVdof,MaxVdof,54*MaxVdof);   % for fluid convection matrix
TCMat = spalloc(MaxTdof,MaxTdof,30*MaxTdof);  % for thermal convection mat.
RHS=spalloc(MaxVdof,1,MaxVdof);

% Copy fields to dof vector
Qv0=Qv;  Qt0=Qt;
tclock=clock;   % Start assembly time <<<<<<<<<

for ne=1:NumEl   % BEGIN GLOBAL MATRIX ASSEMBLY
  
  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,:),nn2vft);
     DMat = DMat + sparse(Rndx,Cndx,Emat,MaxVdof,MaxVdof); 
     
    % Element thermal diffusion matrix
     [Emat,Rndx,Cndx] = TDMatSW(et,Xe,Elcon(ne,:),nn2tft); 
     TDMat = TDMat + sparse(Rndx,Cndx,Emat,MaxTdof,MaxTdof); 
    
     % Global fluid bouyancy matrix
     [Emat,Rndx,Cndx] = BMatSW(eu,Xe,Elcon(ne,:),nn2vft,nn2tft);
     BMat = BMat + sparse(Rndx,Cndx,Emat,MaxVdof,MaxTdof); 
   end 
  
   if (NLitr>1) 
%   First iteration uses Stokes equation.     
% Get stream function and velocities for linearized Navier-Stokes
   for n=1:nnd                  % Loop over local nodes of element
      Vdof(NV,n)=Qv((nn2vft(Elcon(ne,n),1)-1)+NV); 
   end
   
%     Fluid element convection matrix   
       [Emat,Rndx,Cndx] = CMatW(eu,Xe,Elcon(ne,:),nn2vft,Vdof); 
      CMat = CMat + sparse(Rndx,Cndx,Emat,MaxVdof,MaxVdof);  
      
%     Thermal convection matrix       
      [Emat,Rndx,Cndx] = TCMatSW(eu,Xe,Elcon(ne,:),nn2tft,Vdof); 
      TCMat = TCMat + sparse(Rndx,Cndx,Emat,MaxTdof,MaxTdof); 
   end;  % NLitr>1 
   
end;   % END GLOBAL MATRIX ASSEMBLY

MatT = -TCMat-Ctd*TDMat;
MatV = -CMat -Cvd*DMat;

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

% Solve system
tclock=clock;    % Start solution time  <<<<<<<<<<<<<<

RHSt=sparse(MaxTdof,1);
RHSt=RHSt(TBC.ndro)-MatT(TBC.ndro,TBC.dof)*TBC.val;
MatT=MatT(TBC.ndro,TBC.ndro);
Qs=Qt(TBC.ndro);

%  Solve for temperatures
[Lm,Um] = ilu(MatT,su);      % incomplete LU
Qr = gmres(MatT,RHSt,GM.MIter,GM.TTol,GM.MRstrt,Lm,Um,Qs);	% GMRES

Qt=[Qr;TBC.val];       % Augment active dofs with esential (Dirichlet) dofs
Qt=Qt(TBC.p_vec_undo);     % Restore natural order

RHS=RHS-Cbf*BMat*Qt;
RHS=RHS(VBC.ndro)-MatV(VBC.ndro,VBC.dof)*VBC.val;
MatV=MatV(VBC.ndro,VBC.ndro);
Qs=Qv(VBC.ndro);

%  Solve for velocities
[Lm,Um] = ilu(MatV,su);      % incomplete LU
Qr = gmres(MatV,RHS,GM.MIter,GM.VTol,GM.MRstrt,Lm,Um,Qs);	% GMRES

Qv=[Qr;VBC.val];       % Augment active dofs with esential (Dirichlet) dofs
Qv=Qv(VBC.p_vec_undo);       % Restore natural order

% --------------------------------------------
   Qv=RelxFac*Qv+(1-RelxFac)*Qv0;
% --------------------------------------------
% Compute change and copy dofs to field arrays
dsqp=0; dsqv=0; dsqt=0;
for k=1:MaxVdof
  ni=nvf2nnt(k,1);
  nx=NodNdx(ni,1);  ny=NodNdx(ni,2);
  switch nvf2nnt(k,2) % switch on dof type 
    case 1
      dsqp=dsqp+Wa(ny,nx)*(Qv(k)-Qv0(k))^2;  psi0(ny,nx)=Qv(k);
      errpsi(ny,nx)=Qv0(k)-Qv(k);
    case 2
      dsqv=dsqv+Wa(ny,nx)*(Qv(k)-Qv0(k))^2;   u0(ny,nx)=Qv(k);
    case 3
      dsqv=dsqv+Wa(ny,nx)*(Qv(k)-Qv0(k))^2;   v0(ny,nx)=Qv(k);
    case 4
      pxx0(ny,nx)=Qv(k);
    case 5
      pxy0(ny,nx)=Qv(k);
    case 6
      pyy0(ny,nx)=Qv(k);
  end  % switch on dof type 
end  % for 
for k=1:MaxTdof
  ni=ntf2nnt(k,1);
  nx=NodNdx(ni,1);  ny=NodNdx(ni,2);
  switch ntf2nnt(k,2) % switch on dof type 
    case 1
      dsqt=dsqt+Wa(ny,nx)*(Qt(k)-Qt0(k))^2;   T0(ny,nx)=Qt(k);
    case 2
      Tx0(ny,nx)=Qt(k);
    case 3
      Ty0(ny,nx)=Qt(k);
  end  % switch on dof type 
end  % for 
np0(NLitr)=sqrt(dsqp); 
nv0(NLitr)=sqrt(dsqv); 
nT0(NLitr)=sqrt(dsqt); 

disp(['Solution time for linear system = '...
     num2str(etime(clock,tclock)) ' sec,   nv0=' num2str(nv0(NLitr)) ...
     ',  nT0=' num2str(nT0(NLitr))]); 

if (np0(NLitr)<=1e-15||nv0(NLitr)<=1e-15||nT0(NLitr)<=1e-15) 
  MaxNLit=NLitr; np0=np0(1:MaxNLit); nv0=nv0(1:MaxNLit); nT0=nT0(1:MaxNLit);  end;

%---------- Begin plot of intermediate results ----------
   
figure(2);
orient portrait;
orient tall;

% Stream function (intermediate) 
subplot(2,2,1);
contour(Xgrid,Ygrid,psi0,10,'k');  % Plot contours (trajectories)
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(['Thermal cavity stream lines;  Ra=' num2str(Ra)]);

% Temperature (intermediate) 
subplot(2,2,3);
contour(Xgrid,Ygrid,T0,10,'k');
axis([Xmin,Xmax,Ymin,Ymax]);
hold on;
plot([Xmin,Xmin,Xmax,Xmax,Xmin],[Ymax,Ymin,Ymin,Ymax,Ymax],'k');
hold off;
axis equal;  axis image;
title(['Thermal cavity isotherms, Ra=' num2str(Ra)]);

% Convergence 
subplot(2,2,2);
semilogy(1:MaxNLit,nT0,'k-x',1:MaxNLit,nv0,'k-+',1:MaxNLit,np0,'k-o');
xlabel('Nonlinear iteration number');
ylabel('Nonlinear correction');
axis square;
title(['Iteration conv.,  Ra=' num2str(Ra)]);
legend('T','U','Psi');

% Error
subplot(2,2,4);
axis([Xmin,Xmax,Ymin,Ymax]);
contour(Xgrid,Ygrid,errpsi,8,'k');  % Plot contours (errors)
hold on;
pcolor(Xgrid,Ygrid,errpsi);
shading interp  %flat;
plot([Xmin,Xmin,Xmax,Xmax,Xmin],[Ymax,Ymin,Ymin,Ymax,Ymax],'k');
hold off;
axis([Xmin,Xmax,Ymin,Ymax]);
axis image;
title('Iteration correction, psi');
pause(1);

%----------  End plot of intermediate results  ---------
 
end;   % <<< END NONLINEAR ITERATION

format short g;
disp('Convergence results by iteration: temperature, velocity, stream function');
disp(['nT0:  ' num2str(nT0)]); 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,nn2vft,Qv,EBCp,Cvd);
P=reshape(P,NumNy,NumNx); 
Px = reshape(Px,NumNy,NumNx); 
Py = reshape(Py,NumNy,NumNx); 
  
%----------Plot of final result-------------

figure(1);

% Stream function    (final)
subplot(2,2,1);
contour(Xgrid,Ygrid,psi0,10,'k');  % Plot contours (trajectories)
hold on;
pcolor(Xgrid,Ygrid,psi0);
shading interp  %flat;
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, Ra=' num2str(Ra)]);

% Temperature        (final)
subplot(2,2,3);
contour(Xgrid,Ygrid,T0,10,'k');
hold on;
pcolor(Xgrid,Ygrid,T0);
shading interp  %flat;
plot([Xmin,Xmin,Xmax,Xmax,Xmin],[Ymax,Ymin,Ymin,Ymax,Ymax],'k');
hold off;
axis([Xmin,Xmax,Ymin,Ymax]);
axis equal;  axis image;
title(['Isotherms, Ra=' num2str(Ra)]);...
  
% Plot pressure contours (final)
subplot(2,2,4);
contour(Xgrid,Ygrid,P,10,'k');  % Plot pressure contours
hold on;
pcolor(Xgrid,Ygrid,P);
shading interp  %flat;
plot([Xmin,Xmin,Xmax,Xmax,Xmin],[Ymax,Ymin,Ymin,Ymax,Ymax],'k');
hold off;
axis([Xmin,Xmax,Ymin,Ymax]);
axis equal;   axis image;
title(['Pressure contours, Ra=' num2str(Ra)]);

% ************* FIGURE 2 ***********************

figure (2);
clf;
orient tall;

% U-velocity    (final)
subplot(2,2,1);
contour(Xgrid,Ygrid,u0,10,'k');     % Plot vector field 
hold on;
pcolor(Xgrid,Ygrid,u0);
shading interp  %flat;
plot([Xmin,Xmin,Xmax,Xmax,Xmin],[Ymax,Ymin,Ymin,Ymax,Ymax],'k');
hold off;
axis([Xmin,Xmax,Ymin,Ymax]);
axis equal;  axis image;
title('U-velocity field');...

% V-velocity    (final)
subplot(2,2,3);
contour(Xgrid,Ygrid,v0,10,'k');       % Plot vector field 
hold on;
pcolor(Xgrid,Ygrid,v0);
shading interp  %flat;
plot([Xmin,Xmin,Xmax,Xmax,Xmin],[Ymax,Ymin,Ymin,Ymax,Ymax],'k');
hold off;
axis([Xmin,Xmax,Ymin,Ymax]);
axis equal;  axis image;
title('V-velocity field');...
  
% Convergence 
subplot(2,2,2);
semilogy(1:MaxNLit,nT0,'k-x',1:MaxNLit,nv0,'k-+',1:MaxNLit,np0,'k-o');
xlabel('Nonlinear iteration number');
ylabel('Nonlinear correction');
axis square;
title(['Iteration conv.,  Ra=' num2str(Ra)]);
legend('T','U','Psi');

subplot(2,2,4);
nc=fix((NumNy+1)/2);
plot(Xgrid(nc,:),T0(nc,:),'k-x');
title('Centerline temperature');...

% ---------- End final plot ------------

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