Code covered by the BSD License

# Files Associated with a FREE Finite Volume Textbook

### Clive Mingham (view profile)

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
%                     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,
%
% 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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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):
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.
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')
%%%%%%%%%%%%%%%%%%%%%%%%%%% subfunctions  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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;
% 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;
% 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;
% 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;
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
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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.
% 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)];
%           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 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
%          iG=g*[Iu(i,j,1),Iu(i,j,2)]; % u gradient vector in i direction
%          % J direction
%          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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Verification:
%
% Finds the gradients for each component of U in the i and j directions for each cell.
% (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
Udiff=zeros(d);
Udiffoverd=zeros(d);
for i=2:NI+1
for j=2:NJ+1
% i direction
for k=1:d
Udiff(k)=U(i+1,j,k)-U(i-1,j,k);
Udiffoverd(k)=Udiff(k)/LI(i,j);
% left and right scalar gradients
% j direction
Udiff(k)=U(i,j+1,k)-U(i,j-1,k);
Udiffoverd(k)=Udiff(k)/LJ(i,j);
% lower and upper scalar gradients
end % for loop
end % of j loop
end % of i loop
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%```