Code covered by the BSD License  

Highlights from
Files Associated with a FREE Finite Volume Textbook

Files Associated with a FREE Finite Volume Textbook



28 Feb 2012 (Updated )

A series of finite volume m-files for solving the linear advection and shallow water equations.

function FVlinearadvectionQUAD1D
% File name: FVlinearadvectionQUAD1D.m
% Author: Clive Mingham
% Date: 7 July 2011
% Description:  Solves the single 1D PDE, dU/dt + div(H) = 0,
% using the QUAD finite volume scheme on a uniform Cartesian mesh which is
% a single row of cells in the x direction with sides dx and dy.
% This equation corresponds to the finite volume form:
% doubleintegral(dU/dt dA) + lineintegral(H.n ds) = 0.
% In this case the equation is the linear advection equation so: 
% velocity: v = vx i + vy j where vx is constant and vy = 0.
% flux density: H = vU = vx U i
% There are NI cells in the i (x) direction and 1 cell in the j (y)
% direction. Since vy = 0 there is flow only through the left and right 
% sides of each cell. Cell areas and side vectors are the same for each cell.
% The program is written in a structured way using subfunctions so that it
% can be modified easily to the pseudo-1D and 2D cases on non-Cartesian
% structured meshes.
% Initial conditions:  Gaussian at cell centres.
% Boundary conditions: Dirichlet, u=0 in the 2 left and 1 right ghost cells. 
% subfunctions: freadmesh, fcellarea, fsidevectors, finitialu,
%               finterfaceflux, fghostcells, fcalcdt, fplotresults, 
%               fdrawmesh.
% Verified:  Not verified.
clc; clf; clear all;
runtime=130.0;  % runtime in seconds
t=0;          % current time
timelevel=0;  % current time level
vx=0.5;        % velocity component in x direction [m/s]
vy=0;          % zero velocity component in y direction 
v=[vx,vy];     % velocity vector = vx i + 0 j
[x,y,xcen,solid]=freadmesh; % gets coordinates of the 1D computational mesh 
                            % and solid flags (excluding ghost cells)
% Note that there are no solid interior cells in 1D.
disp('mesh read in')
NI=m-1; % number of computational cells in i direction.
A=fcellarea(x,y);   % computes constant cell area
disp('calculated cell area')
S=fsidevectors(x,y); % compute and store cell side vectors                                         
disp('calculated cell side vectors')
u0=finitialu(xcen,t,vx); % Put initial cell centre values of u in a NIx1 array
uinitial=u0; % store initial profile for plotting
disp('inserted initial conditions')
% Append extra cells around arrays to store any ghost values.
u0=finsertghostcells(u0);  % u0 is now (NI+2)x1
u1=zeros(size(u0));  % (NI+2)x1 array for u values at next time level
disp('time marching starts')   
 u0=fbcs(u0);  % Implement boundary conditions using ghost cells. 
               % u0 is (NI+2)x1 and each computational cell is 
               % indexed i=2 to NI+1.
 dt=fcalcdt(A,S,v); % Finds stable time step for each iteration.
 t=t+dt  % update time  
   for i=3:NI+2 
       IH=finterfaceflux(v,u0,i); % gets the left and right interface fluxes for cell i
       IHright=[IH(1,1),IH(1,2)]; % interface flux vector right side 
       IHleft=[IH(2,1),IH(2,2)];  % interface flux vector for left side
       sright=[S(1,1),S(1,2)]; % side vector for right side
       sleft=[S(2,1),S(2,2)];  % side vector for left side
       % FV scheme
       % compute total flux out of cell i
       % totalfluxout=vx*u0(i,1)*dy+vx*u0(i-1,1)*(-dy);
   end % of i loop
   u0=u1;  % update u values
end  % of while loop
disp('1D Quad scheme ended at time t')
%fdisplay(u0); % print results for a small mesh as a matrix
uexact=finitialu(xcen,t,vx); % exact solution
fplotresults(xcen,u0,uinitial,uexact) % plot results
disp('finished 1D QUAD scheme, see figures')
end % FVlinearadvectionQUAD1D
%%%%%%%%%%%%%%%%%%%%%%%%%%% subfunctions  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [x,y,xcen,solid]=freadmesh
% Verification:  not verified
% The mesh is structured and has NI cells in the i (x) direction and 
% 1 cell in the j (y) direction.
% The x and y coordinates of the lower left hand corner of cell (i,j) are
% held in arrays x and y respectively which are both (NI+1)x2.  In
% this way the 4 vertices of cell (i,j) are (x(i),y(j)), (x(i+1),y(j)),
% (x(i+1),y(j+1)) and (x(i),y(j+1)).  solid is an NI by 1 array which 
% flags solid cells.  If cell (i,j) is solid then solid(i,j)=1 otherwise
% solid(i,j)=0.
x=zeros(NI+1,2); % allocate correct sized array
y=zeros(NI+1,2); % allocate correct sized array
xcen=zeros(NI,1); % allocate correct sized array
solid=zeros(NI,1); % allocate correct sized array
dx=2; % cell length in x direction
dy=1; % cell length in y direction (arbitrary)
for i=1:NI+1
    for j=1:2
    end % of j loop
end % of i loop
% find cell centres by averaging coordinates of vertices 
% (this works for a general structured mesh)
for i=1:NI
end % of i loop
end % freadmesh
function A=fcellarea(x,y)  
% Verification:  verified.
% In a uniform Cartesian mesh cell area = dx dy 
end % fcellarea
function S=fsidevectors(x,y)
% Verification: not verified
% For each cell:
% right side vector = dy i + 0 j 
% left side vector = -dy i + 0 j
% Right side vectors
% Left side vectors
end % fsidevectors
function u=finitialu(xcen,t,vx) 
% Verification:  verified
% Inserts initial or exact u values at time t. 
% Initial 1D Gaussian function based on the centre, x=xc, of the 
% computational domain.
xmax=max(max(xcen));  % max x value in mesh
xmin=min(min(xcen));  % min x value in mesh
xc=(xmax+xmin)/2;  % approx x coord of mesh centre
cutoff=(xmax-xc)/2;  % cut off radius for Gaussian
for i=1:NI
        d=abs((x-xc)); % distance from centre
        if (d<cutoff)
           u(i,1)=exp(-0.01*d^2);  % Gaussian
end % of i loop
end % finitialu
function IH=finterfaceflux(v,u0,i)
% Verification:  not verified
% Calculates right and left interface fluxes for each cell.
% Flux H depends on U, i.e. H = H(U). 
% For the 1D linear advection equation, H(U) = vx U i + 0 j.
%% The following QUAD scheme 
%  obtains interface u values at i+1/2 by quadratic interpolation using
%  the 2 neighbouring upwind cell centre values and the neighbouring
%  downwind cell centre value.  For vx > 0 these are indexed by i-1, i and
%  i+1.
% it can be shown that the interpolated interface value of u at i+1/2 is
% (-u(i-1)+6*u(i)+3*u(i+1))/8;
%  so  IHright = v ((-u0(i-1,1)+6*u0(i,1)+3*u0(i+1,1))/8;
%  and IHleft  = v ((-u0(i-2,1)+6*u0(i-1,1)+3*u0(i,1))/8;
%  so 2 left hand and 1 right hand ghost values for u are needed
IH=zeros(2,2); % array to store left and right interface flux vectors
% Right
IH(1,1)=vx*((-u0(i-1,1)+6*u0(i,1)+3*u0(i+1,1))/8);  % i component of right interface flux
IH(1,2)=0;           % j component of right interface flux
% Left
IH(2,1)=vx*((-u0(i-2,1)+6*u0(i-1,1)+3*u0(i,1))/8);  % i component of left interface flux
IH(2,2)=0;             % j component of left interface flux
end % fIflux
function outarray=finsertghostcells(inarray)
% Verification:  not verified
% For the quad scheme 2 ghost cells are needed at left edge of domain, 
% and 1 at right.  array is embedded into (m+3)x1 array of zeros.  
% Hence computational indices go from i=3 to i=m+2.
end % finsertghostcells
function u=fbcs(u)
% Verification:  not verified
% Implements boundary conditions using ghost cells indexed by i=1,2,m+2.
% Dirichlet: u=0 at each end.
% Don't need to do anything as u values in ghost cells are already zero.
end % bcs
function dt=fcalcdt(A,S,v)
% Verification:  not verified
% Finds allowable time step dt using heuristic formula.
F=0.01; % tuning factor (F=1 gives exact results)
dt=A/abs(dot(v,rightside)); % heuristic formula
end  % fcalcdt
function fdrawmesh(x,y,solid)
% Verification:  not verified
% Description:  Draws a structured 2D finite volume mesh and fills in any
% solid cells. 
% Date structures:
% The mesh has NI cells in the i direction and NJ cells in the j direction.
% The x and y coordinates of the lower left hand corner of cell (i,j) are
% held in arrays x and y respectively which are both NI+1 by NJ+1.  In
% this way the 4 vertices of cell (i,j) are (x(i),y(j)), (x(i+1),y(j)),
% (x(i+1),y(j+1)) and (x(i),y(j+1)).  solid is an NI by NJ array which 
% flags solid cells.  If cell (i,j) is solid then solid(i,j)=1 otherwise
% solid(i,j)=0.
NI=m-1; % number of cells in i direction
NJ=1; % number of cells in j direction
% Plot the mesh
hold on   % keeps all plots on the same axes
% draw lines in the i direction
for i=1:NI+1
% draw lines in the j direction
for j=1:NJ+1
title('computational mesh')
% Fill in solid cells
for i=1:NI
    for j=1:NJ
        if (solid(i,j)==1)
    end % of j loop
end % of i loop
end % fdrawmesh
function fplotresults(xcen,u0,uinitial,uexact)
% Verification:  not verified
% Dispays results on regular Cartesian mesh
u0=u0(3:NI+2,1); % extract computational values
title('plots of initial profile (--), Quad numerical(.) and exact solutions (solid line)')
xlabel('x [m]')
ylabel('U [kg/m]')
end % fplotresults

Contact us