Code covered by the BSD License  

Highlights from
Files Associated with a FREE Finite Volume Textbook

Files Associated with a FREE Finite Volume Textbook

by

 

28 Feb 2012 (Updated )

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

FV2DSWE
function FV2DSWE
% File name: FV2DSWE.m
% Author Clive Mingham 
% Date: 11 April 2011
% This m-file is associated with the free text book: 'Introductory Finite 
% Volume Methods for PDEs' which may be downloaded from: BookBooN
%                     not yet published
%
% Description:  Finite volume solvers for the 2D PDE, dU/dt + div(H) = 0.
% Solvers are all first order in time.
% This equation corresponds to the finite volume form:
% doubleintegral(dU/dt dA) + lineintegral(H.n ds) = 0.
% The equation is the (homogeneous) shallow water equation so:
%
% h=h(t,x,y)is water depth. 
% g is acceleration due to gravity.
% fi=h*g is the geopotential.
% v = v(t,x,y) = vx i + vy j = <vx,vy> is water velocity.
% U=[ fi        
%    fi*vx         
%    fi*vy]  is the matrix of conserved variables.        
%
% H=Fi+Gj is the flux density where i and j are the usual Cartesian 
% basis vectors and,
%
% F=[  fi*vx        ,   G=[   fi*vy
%    fi*vx^2+fi^2/2         fi*vx*vy
%     fi*vx*vy   ]         fi*vy^2+fi^2/2 ].
%
% The program is written in a structured way using subfunctions so that it
% can easily be modified for other equations, schemes and meshes.
% subfunctions: freadmesh, fcellarea, fsidevectors, finitialu,
%               fIflux, fghostcells, fcalcdt, fplotresults, fdisplay,
%               fdrawmesh,fgradients,fRiemann,flimiter
%
% This program was based on the following programs written by Clive Mingham: 
% 1) SWE.m, a finite difference 1D SWE solver and 
% 2) FV2Dlinearadvection.m, a structured finite volume solver for the 2D 
%    linear advection equation.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Verification: ???
% 
% Comments: This program is written for clarity using explicit for loops.  
% Compute speed could be greatly increased by using vector operations.
% 
%%
clc; clf; clear all;
% Set run parameters:
Ndim=1; % Number of dimensions: 1 or 2
ICname='Dambreak' % initial condition name: 'Dambreak', 'other'
scheme='FOU'; % scheme name: 'FOU', 'other'
runtime=1.0; % required runtime [s]
maxtimesteps=9999; % maximum number of time steps
stop=0; % stop flag = 1 to stop time marching at correct time.
t=0.0; % initial time [s]
tlevel=0; % time level
%
% Compute mesh data (includes ghost cells):
[x,y,A,xcen,ycen,solid,R,LI,Iu,LJ,Ju,LIl,LIr,LJl,LJr]=freadmesh(Ndim);
disp('mesh read in and mesh data computed')
%
S=fsidevectors(x,y); % compute and store cell side vectors in an
                     % (NI+2)x(NJ+2)x4x2 array which contains ghost cells.
disp('calculated cell side vectors')
U=finitialU(ICname,xcen,ycen); % put initial cell centre values of U in NIxNJx3 array
Uinitial=U          % store for plotting initial conditions
disp('inserted initial conditions')
%
% Append extra cells around arrays to store any ghost values.
U=finsertghostcells(U); % U is now (NI+2)x(NJ+2)x3
[m,n,d]=size(U);  % d is number of rows of U
NI=m-2; % number of computational cells in i direction.
NJ=n-2; % number of computational cells in j direction.
Unext=zeros(size(U))   % (NI+2)x(NJ+2)x3 array for U values at next time level
%%
disp('time marching starts')   
for timestep=1:maxtimesteps 
 U=fbcs(U);  % Implement boundary conditions using ghost cells. 
             % U is (NI+2)x(NJ+2)x3. Computational cells are 
             % indexed i=2 to NI+1, j=2 to NJ+1.
 [iGrad,jGrad,iGradL,iGradR,jGradL,jGradR]=fgradients(U,LI,Iu,LJ,Ju,LIl,LIr,LJl,LJr); % gets gradient vectors
 IH=fIflux(scheme,U,iGrad,jGrad,iGradL,iGradR,jGradL,jGradR,Iu,Ju,R,S); % gets all interface flux vectors
 dt=fcalcdt(Ndim,A,S,U); % Finds stable time step for each iteration.
 t=t+dt      % output current time 
 if (t>runtime)
     dt=runtime-(t-dt); % reduce last time step to stop at t = runtime
     t=runtime;
     stop=1; % flag to stop program at t = runtime
 end
 timestep  % output current time step
 for j=2:NJ+1 
   for i=2:NI+1 
       %
       s1=[S(i,j,1,1),S(i,j,1,2)]; % side vector for side 1
       s2=[S(i,j,2,1),S(i,j,2,2)]; % side vector for side 2
       s3=[S(i,j,3,1),S(i,j,3,2)]; % side vector for side 3
       s4=[S(i,j,4,1),S(i,j,4,2)]; % side vector for side 4
       %
       dtoverA=dt/A(i,j);
       %
     for k=1:d
       %IH is (NI+2)x(NJ+2)xdx4x2;
       IH1=[IH(i,j,d,1,1),IH(i,j,d,1,2)]; % interface flux vector for side 1
       IH2=[IH(i,j,d,2,1),IH(i,j,d,2,2)]; % interface flux vector for side 2
       IH3=[IH(i,j,d,3,1),IH(i,j,d,3,2)]; % interface flux vector for side 3
       IH4=[IH(i,j,d,4,1),IH(i,j,d,4,2)]; % interface flux vector for side 4      
       % FV scheme
       % compute total flux out of cell (i,j) for row k
       totalfluxout=dot(IH1,s1)+dot(IH2,s2)+dot(IH3,s3)+dot(IH4,s4);
       Unext(i,j,k)=U(i,j,k)-dtoverA*totalfluxout; 
     end % k loop
     %
   end % of i loop
 end % of j loop        
 U(2:NI+1,2:NJ+1,:)=Unext(2:NI+1,2:NJ+1,:);  % update U values
 U
 %
 if (stop==1)
     disp('time marching ends')
     t
     break % leave time marching loop
 end
end  % of time marching loop
%
%Uexact=fexactU(ICname,xcen,ycen,t); % exact solution at time t
% remove ghost cells
x=fremoveghostcells(x);
y=fremoveghostcells(y);
solid=fremoveghostcells(solid);
xcen=fremoveghostcells(xcen);
ycen=fremoveghostcells(ycen);
U=fremoveghostcells(U); % extract final computational values
%
%fdrawmesh(x,y,solid) % draws mesh if fplotresults is commented out
%fdisplay(U(:,:,1); % print results as a matrix for a small mesh if needed
fplotresults(xcen,ycen,U,Ndim) % plot results: comment in or out
%
disp('finished main program, see figures')
end % FV2Dlinearadvection
%%%%%%%%%%%%%%%%%%%%%%%%%%% subfunctions  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [x,y,A,xcen,ycen,solid,R,LI,Iu,LJ,Ju,LIl,LIr,LJl,LJr]=freadmesh(Ndim)
% The mesh is structured and has NI computational cells in the i direction 
% and NJ computational cells in the j direction (so each cell has 4 sides).
% 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)x(NJ+1) before 
% ghost cells re appended.  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)). A holds cell
% areas, xcen and ycen hold cell centres and solid holds solid flags - if cell (i,j) 
% is solid then solid(i,j)=1 otherwise solid(i,j)=0.
% R holds normal vectors from cell centres to each of the 4 cell faces.
% Ghost cells are appended to all arrays increasing each dimension index by 2.
%
% NI=51,NJ=51,dx=2,dy=2 makes mesh centre (51,51)
NI=10;
NJ=51;
dx=2;
dy=2;
if (Ndim==1)
    % 1D problem
    NJ=1;
    dy=1;
end
x=zeros(NI+1,NJ+1); % allocate correct sized array
y=zeros(NI+1,NJ+1); % allocate correct sized array
% create cell vertices
for i=1:NI+1
    for j=1:NJ+1
        x(i,j)=(i-1)*dx;
        y(i,j)=(j-1)*dy;
    end % of j loop
end % of i loop
nonCartesian=0; % 1 for non-Cartesian mesh
if (nonCartesian==1)
   disp('mesh is non-Cartesian')
   % Add 'random' numbers to mesh coordinates to produce a non-Cartesian mesh
   % ensure that the abs of these numbers are less than min(dx/2,dy/2).
   % Create non-Cartesian mesh
   for i=1:NI+1
     for j=1:NJ+1
         x(i,j)=x(i,j)+(dx/6)*cos(1.727*(i+j));
         if (Ndim==2)
         y(i,j)=y(i,j)+(dy/5)*sin(2.46723*i*j);
         end
     end % of j loop
   end % of i loop
else
   disp('mesh is Cartesian') 
end % end of if statement create mesh
%
% Create ghost cells around whole mesh (for gradient calcs)
% First embed x and y in larger arrays
tempx=zeros(NI+3,NJ+3);
tempx(2:NI+2,2:NJ+2)=x;
x=tempx;
clear tempx
tempy=zeros(NI+3,NJ+3);
tempy(2:NI+2,2:NJ+2)=y;
y=tempy;
clear tempy
% Now fill in ghost cell vertices by reflecting cell vertices
% i.e. if z is the vector from vertex (i,2) to vertex (i,3) then 
% ghost vertex (i,1) has position vector (i,2)-z = 2(i,2)-(i,3)
for i=2:NI+2
    x(i,1)=2*x(i,2)-x(i,3);
    y(i,1)=2*y(i,2)-y(i,3);
    x(i,NJ+3)=2*x(i,NJ+2)-x(i,NJ+1);
    y(i,NJ+3)=2*y(i,NJ+2)-y(i,NJ+1);
end
%
for j=2:NJ+2
    x(1,j)=2*x(2,j)-x(3,j);
    y(1,j)=2*y(2,j)-y(3,j);
    x(NI+3,j)=2*x(NI+2,j)-x(NI+1,j);
    y(NI+3,j)=2*y(NI+2,j)-y(NI+1,j);
end
%
% Find cell areas
A=zeros(NI+2,NJ+2);  % allocate correct size array
for i=1:NI+2
    for j=1:NJ+2
        d1=[x(i+1,j+1)-x(i,j),y(i+1,j+1)-y(i,j),0];
        d2=[x(i,j+1)-x(i+1,j),y(i,j+1)-y(i+1,j),0];
        d=cross(d1,d2);   
        A(i,j)=0.5*sqrt(dot(d,d));
    end % of j loop
end % of i loop
%
% Find cell centres by averaging coordinates of vertices
xcen=zeros(NI+2,NJ+2);  % allocate correct sized array
ycen=zeros(NI+2,NJ+2);  % allocate correct sized array
for i=1:NI+2
    for j=1:NJ+2
        xcen(i,j)=(x(i,j)+x(i+1,j)+x(i+1,j+1)+x(i,j+1))/4;
        ycen(i,j)=(y(i,j)+y(i+1,j)+y(i+1,j+1)+y(i,j+1))/4;
    end % of j loop
end % of i loop
% Input solid flags
solid=zeros(NI+2,NJ+2); % allocate correct sized array
% Cells are all fluid for now.
%
% Find normal vector from cell centre to each of 4 cell faces
% using a vector based formula R=p-q, where q is the cell centre position
% vector, p is the position vector of the point on the cell side,
% p = a + t (b-a) where a and b are position vectors of cell vertices
% defining the side and t = (q-a).(b-a)/(b-a).(b-a)
R=zeros(NI+2,NJ+2,4,2);  % create correct sized array
for i=1:NI+2
    for j=1:NJ+2
        q=[xcen(i,j),ycen(i,j)];
        % side 1
        a=[x(i,j),y(i,j)];
        b=[x(i+1,j),y(i+1,j)];
        bminusa=b-a;
        qminusa=q-a;
        t=dot(qminusa,bminusa)/dot(bminusa,bminusa);
        p=a+t*bminusa;
        answer=p-q;
        R(i,j,1,1)=answer(1);
        R(i,j,1,2)=answer(2);
        % side 2
        a=[x(i+1,j),y(i+1,j)];
        b=[x(i+1,j+1),y(i+1,j+1)];
        bminusa=b-a;
        qminusa=q-a;
        t=dot(qminusa,bminusa)/dot(bminusa,bminusa);
        p=a+t*bminusa;
        answer=p-q;
        R(i,j,2,1)=answer(1);
        R(i,j,2,2)=answer(2);
        % side 3
        a=[x(i+1,j+1),y(i+1,j+1)];
        b=[x(i,j+1),y(i,j+1)];
        bminusa=b-a;
        qminusa=q-a;
        t=dot(qminusa,bminusa)/dot(bminusa,bminusa);
        p=a+t*bminusa;
        answer=p-q;
        R(i,j,3,1)=answer(1);
        R(i,j,3,2)=answer(2);
        % side 4
        a=[x(i,j+1),y(i,j+1)];
        b=[x(i,j),y(i,j)];
        bminusa=b-a;
        qminusa=q-a;
        t=dot(qminusa,bminusa)/dot(bminusa,bminusa);
        p=a+t*bminusa;
        answer=p-q;
        R(i,j,4,1)=answer(1);
        R(i,j,4,2)=answer(2);
    end % of j loop
end % of i loop
% 
% Find:
% LI is the distance from centre of cell (i-1,j) to centre of 
% cell (i+1,j). Iu(i,j) holds the unit vector. Similarly for LJ and Ju. 
% These are used to give a limited gradient vectors.
LI=zeros(NI+2,NJ+2); % create correct sized array
Iu=zeros(NI+2,NJ+2,2); % create correct sized array
LJ=zeros(NI+2,NJ+2); % create correct sized array
Ju=zeros(NI+2,NJ+2,2); % create correct sized array
for i=2:NI+1 % indices for computational cells
    for j=2:NJ+1 % indices for computational cells
        % i direction
        a=[xcen(i-1,j),ycen(i-1,j)];
        b=[xcen(i+1,j),ycen(i+1,j)];
        n=b-a;  % vector joining cell centres
        LI(i,j)=sqrt(dot(n,n)); % length of vector n
        nunit=n/LI(i,j); % unit vector joining cell centres
        Iunit(i,j,1)=nunit(1); % store i component
        Iunit(i,j,2)=nunit(2); % store j component
        % j direction
        a=[xcen(i,j-1),ycen(i,j-1)];
        b=[xcen(i,j+1),ycen(i,j+1)];
        n=b-a;  % vector joining cell centres
        LJ(i,j)=sqrt(dot(n,n)); % length of vector n
        nunit=n/LJ(i,j); % unit vector joining cell centres
        Junit(i,j,1)=nunit(1); % store i component
        Junit(i,j,2)=nunit(2); % store j component       
    end % of j loop
end % of i loop
%
% Find:
% LIl is the distance between the centres of cells (i-1,j) and (i,j)
% LIr is the distance between the centres of cells (i,j) and (i+1,j)
% these are used for finding left and right I gradients in cell (i,j)
% LJl is the distance between the centres of cells (i,j-1) and (i,j)
% LJr is the distance between the centres of cells (i,j) and (i,j+1)
% these are used for finding left and right J gradients in cell (i,j)
LIl=zeros(NI+2,NJ+2); % create correct sized array
LIr=zeros(NI+2,NJ+2); % create correct sized array
LJl=zeros(NI+2,NJ+2); % create correct sized array
LJr=zeros(NI+2,NJ+2); % create correct sized array
for i=2:NI+1 % indices for computational cells
    for j=2:NJ+1 % indices for computational cells
        % i direction
        a=[xcen(i-1,j),ycen(i-1,j)];
        b=[xcen(i,j),ycen(i,j)];
        c=[xcen(i+1,j),ycen(i+1,j)]; 
        leftvector=b-a;
        rightvector=c-b;
        LIl(i,j)=sqrt(dot(leftvector,leftvector));
        LIr(i,j)=sqrt(dot(rightvector,rightvector));
        % j direction
        a=[xcen(i,j-1),ycen(i,j-1)];
        b=[xcen(i,j),ycen(i,j)];
        c=[xcen(i,j+1),ycen(i,j+1)]; 
        lowervector=b-a;
        uppervector=c-b;
        LJl(i,j)=sqrt(dot(lowervector,lowervector));
        LJr(i,j)=sqrt(dot(uppervector,uppervector));
    end % of j loop
end % of i loop
%       
end % freadmesh
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
function S=fsidevectors(x,y)
% Verification: correct for single diamond shaped cell.
%
% Finds the 4 side (S) vectors for each cell.  
% Cell sides are numbered 1 to 4 anti-clockwise.  For cell (i,j) side 1 
% is the vector s1 from (x(i,j), y(i,j)) to (x(i+1,j),y(i+1,j)).
% if s1=<a,b> then the corresponding side vector S1=<b,-a>.
% The side vector array contains ghost cells.
%
% Note that vectors are 1 by 2 arrays so care must be taken because side
% vectors are stored in higher dimensional arrays.
%
[m,n]=size(x);
NI=m-2; NJ=n-2;
S=zeros(NI+2,NJ+2,4,2);
for i=2:NI+1
    for j=2:NJ+1
        s1=[x(i+1,j)-x(i,j),y(i+1,j)-y(i,j)];
        S(i,j,1,1)=s1(2);
        S(i,j,1,2)=-s1(1);
        %
        s2=[x(i+1,j+1)-x(i+1,j),y(i+1,j+1)-y(i+1,j)];
        S(i,j,2,1)=s2(2);
        S(i,j,2,2)=-s2(1);
        %
        s3=[x(i,j+1)-x(i+1,j+1),y(i,j+1)-y(i+1,j+1)];
        S(i,j,3,1)=s3(2);
        S(i,j,3,2)=-s3(1);
        %
        s4=[x(i,j)-x(i,j+1),y(i,j)-y(i,j+1)];
        S(i,j,4,1)=s4(2);
        S(i,j,4,2)=-s4(1);    
    end % of j loop
end % of i loop
% [S(2,2,1,1),S(2,2,1,2)]
% [S(2,2,2,1),S(2,2,2,2)] they are correct for a diamond
% [S(2,2,3,1),S(2,2,3,2)]
% [S(2,2,4,1),S(2,2,4,2)]
end % fsidevectors
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
function U=finitialU(ICname,xcen,ycen)
% Verification:  correct.
% Inserts initial U values in computational cells (U is 3xNIxNJ)
[m,n]=size(xcen); % xcen contains ghost cells
NI=m-2; % number of computational cells in the i direction
NJ=n-2; % number of computational cells in the j direction
%
xcen=xcen(2:NI+1,2:NJ+1); % remove ghost cells from xcen
ycen=ycen(2:NI+1,2:NJ+1); % remove ghost cells from yecen
%
U=zeros(NI,NJ,3); % create correct sized array (without ghost values)
midI=floor(NI/2); % central i value
midIplus1=midI+1;
g=9.81;
%
switch ICname
 case 'Dambreak'
    % pseudo 1D Dam break.  Left and right states are either side 
    % of central i value in the computational domain.
    % left states
    hl=1.0/g; % left water depth
    vl=[0,0]; % left water velocity vector
    fil=hl*g; % left geopotential
    % right states
    hr=0.1/g; % right water depth
    vr=[0,0]; % right water velocity vector
    fir=hr*g; % right geopotential
    %
    % left state
    U(1:midI,:,1)=fil;
    U(1:midI,:,2)=fil*vl(1);
    U(1:midI,:,3)=fil*vl(2);
    % right state
    U(midIplus1:NI,:,1)=fir;
    U(midIplus1:NI,:,2)=fir*vr(1);
    U(midIplus1:NI,:,3)=fir*vr(2);
%% end of Dambreak   
 otherwise
   disp('Unknown initial conditions, program ends')
   exit
 end % of switch/case statements
%
end % finitialU
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function H=getHfromU(U)
% U = [   fi
%       fi*vx
%       fi*vy ]
%
% H=<     fi*vx    ,  fi*vy           
%    fi*vx^2+fi^2/2,  fi*vx*vy              
%         fi*vx*vy ,  fi*vy^2+fi^2/2 > 
[m,n,d]=size(U);
fi=zeros(m,n);
vx=zeros(m,n);
vy=zeros(m,n);
fivxvy=fi.*vx.*vy;
fi2over2=fi.*fi/2;
H=zeros(m,n,d,2); % flux density vectors
%
fi=U(:,:,1);
vx=U(:,:,2)./fi;
vy=U(:,:,3)./fi;
%
H(:,:,1,1)=U(:,:,2);
H(:,:,1,2)=U(:,:,3);
%
H(:,:,2,1)=fi.*vx.*vx+fi2over2;
H(:,:,2,2)=fivxvy;
%
H(:,:,3,1)=fivxvy;
H(:,:,3,2)=fi.*vy.*vy+fi2over2;
end % getHfromU
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function IH=fIflux(scheme,U,iGrad,jGrad,iGradL,iGradR,jGradL,jGradR,Iu,Ju,R,S)
% Calculates the 4 interface fluxes on the 4 sides of cell (i,j) 
% for each component of U.
% Cell sides are numbered 1 to 4 anti-clockwise.  For cell (i,j) side 1 
% is the vector s1 from (x(i,j), y(i,j)) and to (x(i+1,j),y(i+1,j)).
% This means that sides 2 and 4 are in the i direction and sides 1 and 3
% are in the j direction.
% iGrad and jGrad are gradient vectors in i and j directions.
% R holds normal vectors from cell centres to each of the 4 cell faces.
% S holds side vectors for each of the 4 cell faces.
% Interface fluxes are denoted IH1, IH2, IH3 and IH4 and are 2D vectors.
% Flux H depends on U, i.e. H = H(U). 
% For the 2D shallow water equations:
% H=<     fi*vx    ,  fi*vy
%    fi*vx^2+fi^2/2,  fi*vx*vy
%         fi*vx*vy ,  fi*vy^2+fi^2/2 >  
% 
% There is considerable choice for interface flux estimation: different 
% choices give different FV schemes.   Some schemes are coded below.
% 
[m,n,d]=size(U);  % d=number of rows of U
NI=m-2; NJ=n-2;  % number of computational cells in i and j directions
IH=zeros(NI+2,NJ+2,d,4,2); % array to store the 4 interface flux vectors 
%                            for each cell.
%
% 
switch scheme  % name of the scheme
 case 'FOU'
  %% FOU scheme
  disp('FOU scheme')
  % Interface flux is the flux at the neighbouring upwind cell centre
  H=getHfromU(U); % gets H at all cell centres from U values there
  for i=2:NI+1
    for j=2:NJ+1
        fi=U(i,j,1);
        vx=U(i,j,2)/fi;
        vy=U(i,j,3)/fi;
        v=[vx,vy];
        % i direction side 2
        S2=[S(i,j,2,1),S(i,j,2,2)];
        if (dot(v,S2)>0) % flow is direction of increasing i
           for k=1:d
             % loop over each row
             IH(i,j,k,2,1)=H(i,j,k,1);  % row k of IH2 in the x direction
             IH(i,j,k,2,2)=H(i,j,k,2);  % row k of IH2 in the y direction
           end % for loop                 
        else % flow is reversed
            for k=1:d
             % loop over each row
             IH(i,j,k,2,1)=H(i+1,j,k,1);  % row k of IH2 in the x direction
             IH(i,j,k,2,2)=H(i+1,j,k,2);  % row k of IH2 in the y direction
            end % for loop  
        end
        % i direction side 4
        S4=-[S(i,j,4,1),S(i,j,4,2)]; % reverse direction
        if (dot(v,S4)>0) % flow is in direction of increasing i
            for k=1:d
             % loop over each row
             IH(i,j,k,4,1)=H(i-1,j,k,1);  % row k of IH4 in the x direction
             IH(i,j,k,4,2)=H(i-1,j,k,2);  % row k of IH4 in the y direction
            end % for loop        
        else % flow is reversed
            for k=1:d
             % loop over each row
             IH(i,j,k,4,1)=H(i,j,k,1);  % row k of IH4 in the x direction
             IH(i,j,k,4,2)=H(i,j,k,2);  % row k of IH4 in the y direction
            end % for loop        
        end
        % j direction side 3
        S3=[S(i,j,3,1),S(i,j,3,2)];
        if (dot(v,S3)>0) % flow is in direction of increasing j
            for k=1:d
             % loop over each row
             IH(i,j,k,3,1)=H(i,j,k,1);  % row k of IH3 in the x direction
             IH(i,j,k,3,2)=H(i,j,k,2);  % row k of IH3 in the y direction
            end % for loop        
        else % flow is reversed
            for k=1:d
             % loop over each row
             IH(i,j,k,3,1)=H(i,j+1,k,1);  % row k of IH3 in the x direction
             IH(i,j,k,3,2)=H(i,j+1,k,2);  % row k of IH3 in the y direction
            end % for loop
        end
        % j direction side 1
        S1=-[S(i,j,1,1),S(i,j,1,2)]; % reverse direction
        if (dot(v,S1)>0) % flow is in direction of increasing j
            for k=1:d
             % loop over each row
             IH(i,j,k,1,1)=H(i,j-1,k,1);  % row k of IH1 in the x direction
             IH(i,j,k,1,2)=H(i,j-1,k,2);  % row k of IH1 in the y direction
            end % for loop
        else % flow is reversed
            for k=1:d
             % loop over each row
             IH(i,j,k,1,1)=H(i,j,k,1);  % row k of IH1 in the x direction
             IH(i,j,k,1,2)=H(i,j,k,2);  % row k of IH1 in the y direction
            end % for loop
        end 
    end % of j loop
   end % of i loop
%% end of FOU scheme
% case 'HighRes'
%   %% HighRes scheme
%   disp('HighRes scheme')
%     for i=2:NI+1
%      for j=2:NJ+1      
%         % cell side 2: interface (i+1/2,j) 
%           s=[S(i,j,2,1),S(i,j,2,2)];
%         % Left side of interface (i+1/2,j)      
%           r=[R(i,j,2,1),R(i,j,2,2)];
%           iG=[iGrad(i,j,1),iGrad(i,j,2)]; % u gradient vector in i direction
%           uL=u(i,j)+dot(r,iG); % piecewise linear reconstruction
%         % Right side of interface (i+1/2,j)              
%           r=[R(i+1,j,4,1),R(i+1,j,4,2)];
%           iG=[iGrad(i+1,j,1),iGrad(i+1,j,2)]; % u gradient vector in i direction
%           uR=u(i+1,j)+dot(r,iG); % piecewise linear reconstruction
%           %
%           Hstar=fRiemann(uL,uR,v,s); % Riemann interface flux
%           IH(i,j,2,1)=Hstar(1);  % component of IH2 in the x direction
%           IH(i,j,2,2)=Hstar(2);  % component of IH2 in the y direction
%         %  
%         % cell side 4: interface (i-1/2,j) 
%           s=-[S(i,j,4,1),S(i,j,4,2)]; % reverse direction 
%         % Left side of interface (i-1/2,j)         
%           r=[R(i-1,j,2,1),R(i-1,j,2,2)];
%           iG=[iGrad(i-1,j,1),iGrad(i-1,j,2)]; % u gradient vector in i direction
%           uL=u(i-1,j)+dot(r,iG); % piecewise linear reconstruction          
%         % Right side of interface (i-1/2,j)
%           r=[R(i,j,4,1),R(i,j,4,2)];
%           iG=[iGrad(i,j,1),iGrad(i,j,2)]; % u gradient vector in i direction
%           uR=u(i,j)+dot(r,iG); % piecewise linear reconstruction
%           %
%           Hstar=fRiemann(uL,uR,v,s); % Riemann interface flux
%           IH(i,j,4,1)=Hstar(1);  % component of IH4 in the x direction
%           IH(i,j,4,2)=Hstar(2);  % component of IH4 in the y direction
%         %
%         % cell side 3: interface (i,j+1/2)
%           s=[S(i,j,3,1),S(i,j,3,2)];
%           % Left side of interface (i,j+1/2)
%            r=[R(i,j,3,1),R(i,j,3,2)];
%            jG=[jGrad(i,j,1),jGrad(i,j,2)]; % u gradient vector in j direction
%            uL=u(i,j)+dot(r,jG); % piecewise linear reconstruction 
%           % Right side of interface (i,j+1/2)  
%            r=[R(i,j+1,1,1),R(i,j+1,1,2)];
%            jG=[jGrad(i,j+1,1),jGrad(i,j+1,2)]; % u gradient vector in j direction
%            uR=u(i,j+1)+dot(r,jG); % piecewise linear reconstruction
%            %
%            Hstar=fRiemann(uL,uR,v,s); % Riemann interface flux
%            IH(i,j,3,1)=Hstar(1);  % component of IH3 in the x direction
%            IH(i,j,3,2)=Hstar(2);  % component of IH3 in the x direction
%         %
%         % cell side 1
%           s=-[S(i,j,1,1),S(i,j,1,2)]; % reverse direction
%           % Left side of interface (i,j-1/2)
%            r=[R(i,j-1,3,1),R(i,j-1,3,2)];
%            jG=[jGrad(i,j-1,1),jGrad(i,j-1,2)]; % u gradient vector in j direction
%            uL=u(i,j-1)+dot(r,jG); % piecewise linear reconstruction
%            % Right side of interface (i,j-1/2)
%            r=[R(i,j,1,1),R(i,j,1,2)];
%            jG=[jGrad(i,j,1),jGrad(i,j,2)]; % u gradient vector in j direction
%            uR=u(i,j)+dot(r,jG); % piecewise linear reconstruction
%            %
%            Hstar=fRiemann(uL,uR,v,s); % Riemann interface flux
%            IH(i,j,1,1)=Hstar(1);  % component of IH1 in the x direction
%            IH(i,j,1,2)=Hstar(2);  % component of IH1 in the y direction
%            %
%      end % of j loop
%     end % of i loop  
%   %% end of HighRes scheme
%  case 'Halfhancock'
%      disp('Halfhancock scheme')
%     % Method is based on HighRes and includes limiters on gradients.
%     % First find limited gradient vectors iG, jG in each cell.
%     for i=2:NI+1
%      for j=2:NJ+1
%          % I direction
%          g=flimiter(iGradR(i,j),iGradL(i,j)); % limited u gradient
%          iG=g*[Iu(i,j,1),Iu(i,j,2)]; % u gradient vector in i direction
%          % J direction
%          g=flimiter(jGradR(i,j),jGradL(i,j)); % limited u gradient
%          jG=g*[Ju(i,j,1),Ju(i,j,2)]; % u gradient vector in j direction
%      end
%     end
%     % Now find left and right u values at interfaces and solve 
%     % the associated Riemann problem.
%     for i=2:NI+1
%      for j=2:NJ+1      
%         % cell side 2: interface (i+1/2,j) 
%           s=[S(i,j,2,1),S(i,j,2,2)];
%         % Left side of interface (i+1/2,j)      
%           r=[R(i,j,2,1),R(i,j,2,2)];
%           uL=u(i,j)+dot(r,iG); % piecewise linear reconstruction
%         % Right side of interface (i+1/2,j)              
%           r=[R(i+1,j,4,1),R(i+1,j,4,2)];
%           uR=u(i+1,j)+dot(r,iG); % piecewise linear reconstruction
%           %
%           Hstar=fRiemann(uL,uR,v,s); % Riemann interface flux
%           IH(i,j,2,1)=Hstar(1);  % component of IH2 in the x direction
%           IH(i,j,2,2)=Hstar(2);  % component of IH2 in the y direction
%         %  
%         % cell side 4: interface (i-1/2,j) 
%           s=-[S(i,j,4,1),S(i,j,4,2)]; % reverse direction 
%         % Left side of interface (i-1/2,j)         
%           r=[R(i-1,j,2,1),R(i-1,j,2,2)];
%           uL=u(i-1,j)+dot(r,iG); % piecewise linear reconstruction          
%         % Right side of interface (i-1/2,j)
%           r=[R(i,j,4,1),R(i,j,4,2)];
%           uR=u(i,j)+dot(r,iG); % piecewise linear reconstruction
%           %
%           Hstar=fRiemann(uL,uR,v,s); % Riemann interface flux
%           IH(i,j,4,1)=Hstar(1);  % component of IH4 in the x direction
%           IH(i,j,4,2)=Hstar(2);  % component of IH4 in the y direction
%         %
%         % cell side 3: interface (i,j+1/2)
%           s=[S(i,j,3,1),S(i,j,3,2)];
%           % Left side of interface (i,j+1/2)
%            r=[R(i,j,3,1),R(i,j,3,2)];
%            uL=u(i,j)+dot(r,jG); % piecewise linear reconstruction 
%           % Right side of interface (i,j+1/2)  
%            r=[R(i,j+1,1,1),R(i,j+1,1,2)];
%            uR=u(i,j+1)+dot(r,jG); % piecewise linear reconstruction
%            %
%            Hstar=fRiemann(uL,uR,v,s); % Riemann interface flux
%            IH(i,j,3,1)=Hstar(1);  % component of IH3 in the x direction
%            IH(i,j,3,2)=Hstar(2);  % component of IH3 in the x direction
%         %
%         % cell side 1
%           s=-[S(i,j,1,1),S(i,j,1,2)]; % reverse direction
%           % Left side of interface (i,j-1/2)
%            r=[R(i,j-1,3,1),R(i,j-1,3,2)];
%            uL=u(i,j-1)+dot(r,jG); % piecewise linear reconstruction
%            % Right side of interface (i,j-1/2)
%            r=[R(i,j,1,1),R(i,j,1,2)];
%            uR=u(i,j)+dot(r,jG); % piecewise linear reconstruction
%            %
%            Hstar=fRiemann(uL,uR,v,s); % Riemann interface flux
%            IH(i,j,1,1)=Hstar(1);  % component of IH1 in the x direction
%            IH(i,j,1,2)=Hstar(2);  % component of IH1 in the y direction
%            %
%     end % of j loop
%    end % of i loop    
% %% end of Halfhancock scheme    
 otherwise
   disp('Unknown scheme, program ends')
   exit
end % of switch/case statements
%
end % fIflux
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
function outarray=finsertghostcells(inarray)
% Assuming that ghost cells are needed around the whole domain an mxnxd
% array is embedded into (m+2)x(n+2)xd array of zeros.  Hence computational 
% indices go from i=2 to m+1 and j=2 to n+1.
% verification: Correct.
[m,n,d]=size(inarray); % d is number of PDEs
dummy=zeros(m+2,n+2);
for k=1:d
    dummy(2:m+1,2:n+1)=inarray(:,:,k);
    outarray(:,:,k)=dummy;
end % 
end % finsertghostcells
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
function outarray=fremoveghostcells(inarray)
% Removes ghost cells so that a mxnxd array becomes (m-2)x(n-2)xd
% verification: Correct.
[m,n,d]=size(inarray); % d is number of PDEs
outarray=zeros((m-2),(n-2),d);
for k=1:d
    outarray(:,:,k)=inarray(2:m-1,2:n-1,k);
end % 
end % fremoveghostcells
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function U=fbcs(U)
% Implements boundary conditions using ghost cells index by i=1, i=NI+2
% j=1, j=NJ+2.
[m,n,d]=size(U);
%
% Dirichlet: U=0 everywhere.
% Don't need to do anything as U values in ghost cells are already zero.
end % bcs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function dt=fcalcdt(Ndim,A,S,U)
% Verification:  needs checking, looks wrong
% Finds allowable time step dt using heuristic formula (which makes sense
% 0<F<1 is a tunable safety factor.
[M,N]=size(A);  % A and S now enlarged to include ghost cells 
dt=9999999;  % large value of dt to start
if (Ndim==2)
   F=0.4;
else
  % 1D
    F=0.9;
end
 for i=2:M-1
    for j=2:N-1
        % put side vectors into 1x2 arrays  
        s2=[S(i,j,2,1),S(i,j,2,2)];
        norms2=norm(s2); % length of side vector s2
        s3=[S(i,j,3,1),S(i,j,3,2)];
        norms3=norm(s3); % length of side vector s3
        %
        fi=U(i,j,1);
        cel=sqrt(fi);
        vx=U(i,j,2)/fi;
        vy=U(i,j,3)/fi;
        v=[vx,vy]; % velocity
        %
        v2=dot(v,s2); % velocity component through side 2 x side length
        v3=dot(v,s3); % velocity component through side 3 x side length
        speed2=max(abs(v2+cel*norms2),abs(v2-cel*norms2));
        speed3=max(abs(v3+cel*norms3),abs(v3-cel*norms3));
        dum=A(i,j)/max(speed2,speed3);
        % take smallest value of dt
        if (dum<dt)
           dt=dum;
        end
    end % j loop
 end % i loop
dt=F*dt
end  % fcalcdt
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function fdrawmesh(x,y,solid)
% 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.
% 
[m,n]=size(x);
NI=m-1;
NJ=n-1;
%
% Plot the mesh
hold on   % keeps all plots on the same axes
% draw lines in the i direction
for i=1:NI+1
    plot(x(i,:),y(i,:))
end
% draw lines in the j direction
for j=1:NJ+1
    plot(x(:,j),y(:,j))   
end
title('computational mesh')
xlabel('x')
ylabel('y')
% Fill in solid cells
for i=1:NI
    for j=1:NJ
        if (solid(i,j)==1)
            solidx=[x(i,j),x(i+1,j),x(i+1,j+1),x(i,j+1),x(i,j)];
            solidy=[y(i,j),y(i+1,j),y(i+1,j+1),y(i,j+1),y(i,j)];
            fill(solidx,solidy,'k')      
        end
    end % of j loop
end % of i loop
end % fdrawmesh
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function fplotresults(xcen,ycen,U,Ndim)
% Plots results on a structured mesh - even works for non-Cartesian!
%
g=9.81;  % acceleration due to gravity
if (Ndim==2)
  subplot(3,1,1), contour(xcen,ycen,U(:,:,1)/g)
  title('numerical results: contour plot of water depth')
  xlabel('x [m]')
  ylabel('y [m]')
  %
  subplot(3,1,2), contour(xcen,ycen,U(:,:,2)/U(:,:,1))
  title('numerical results: contour plot of water speed in x')
  xlabel('x [m]')
  ylabel('y [m]')
  %
  subplot(3,1,3), contour(xcen,ycen,U(:,:,3)/U(:,:,1))
  title('numerical results: contour plot of water speed in y')
  xlabel('x [m]')
  ylabel('y [m]')
else
  % Ndim =1
  subplot(3,1,1), plot(xcen(:,1),U(:,1,1)/g)
  title('numerical solution for water depth')
  xlabel('x [m]')
  ylabel('water depth [m]')
  %
  subplot(3,1,2), plot(xcen(:,1),U(:,1,2)./U(:,1,1))
  title('numerical solution for x component of water velocity')
  xlabel('x [m]')
  ylabel('speed [m/s]')
  %
  subplot(3,1,3), plot(xcen(:,1),U(:,1,3)./U(:,1,1))
  title('numerical solution for y component of water velocity')
  xlabel('x [m]')
  ylabel('speed [m/s]')
end % of if
%
end % fplotresults
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function fdisplay(in)
 % displays matrices 'correctly' - WORKS
 % Using index i for x which are columns and j for y which are rows in 
 % reverse order so this routine displays them in this fashion
 [NX,NY]=size(in);
 temp=transpose(in);
 out=temp(NY:-1:1,:);
 disp('printing as a matrix')
 out
end % fdisplay
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
function [iGrad,jGrad,iGradL,iGradR,jGradL,jGradR]=fgradients(U,LI,Iu,LJ,Ju,LIl,LIr,LJl,LJr)
% Verification: 
%
% Finds the gradients for each component of U in the i and j directions for each cell.
% iGrad(i,j,d) = [(U(i+1,j,d)-U(i-1,j,d))/LI(i,j)] Iu(i,j,:).
% iGradL(i,j,d) = [(U(i,j,d)-U(i-1,j,d))]/LIl(i,j).  etc.
% (ghost cell gradients are zero). 
%
[m,n,d]=size(U); % d is number of PDEs
NI=m-2; NJ=n-2;  % number of computational cells in i and j directions
iGrad=zeros(NI+2,NJ+2,3,2); % gradient vector
jGrad=zeros(NI+2,NJ+2,3,2); % gradient vector
iGradL=zeros(NI+2,NJ+2,3); % scalar gradient
jGradL=zeros(NI+2,NJ+2,3); % scalar gradient
iGradR=zeros(NI+2,NJ+2,3); % scalar gradient
jGradR=zeros(NI+2,NJ+2,3); % scalar gradient
Udiff=zeros(d);
Udiffoverd=zeros(d);
for i=2:NI+1
    for j=2:NJ+1
        % i direction 
        % gradient vectors
      for k=1:d
        Udiff(k)=U(i+1,j,k)-U(i-1,j,k);
        Udiffoverd(k)=Udiff(k)/LI(i,j);
        iGrad(i,j,k,1)=Udiffoverd(k)*Iu(i,j,1);
        iGrad(i,j,k,2)=Udiffoverd(k)*Iu(i,j,2);
        % left and right scalar gradients
        iGradL(i,j,k)=(U(i,j,k)-U(i-1,j,k))/LIl(i,j);
        iGradR(i,j,k)=(U(i+1,j,k)-U(i,j,k))/LIr(i,j);
        % j direction
        % gradient vectors
        Udiff(k)=U(i,j+1,k)-U(i,j-1,k);
        Udiffoverd(k)=Udiff(k)/LJ(i,j);
        jGrad(i,j,k,1)=Udiffoverd(k)*Ju(i,j,1);
        jGrad(i,j,k,2)=Udiffoverd(k)*Ju(i,j,2);
        % lower and upper scalar gradients
        jGradL(i,j,k)=(U(i,j,k)-U(i,j-1,k))/LJl(i,j);
        jGradR(i,j,k)=(U(i,j+1,k)-U(i,j,k))/LJr(i,j); 
      end % for loop
    end % of j loop
end % of i loop
%
end % fgradients
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Hstar=fRiemann(uL,uR,v,s)
% Verification: 
%
% Finds the Riemann solution, Hstar=v*ustar, for left and right states uL, uR 
% at a cell interface.
% v is the flow velocity at the cell interface
% s is the sidevector for cell interface (i+1/2,j) or (i,j+1/2) and
% -sidevector for cell interface (i-1/2,j) or (i,j-1/2).
%
% Solution for ustar is simply the upwind u value.
if (dot(v,s)>0) % flow is in the direction of increasing i (or j)
    ustar=uL;   % uL is upwind value
else
    ustar=uR;   % uR is upwind value
end 
Hstar=v*ustar;
%
end % fRiemann
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function g=flimiter(x,y)
% Verification: 
% Finds limited scalar gradient g.
%
k=2.0;  % k=2 is Superbee, k=1 is minmod
if (x>0) 
    g=max([0,min(k*x,y),min(x,k*y)]);
else
    g=min([0,max(k*x,y),max(x,k*y)]);
end 
%
end % flimiter
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Contact us