function FVSWE1Dv2
% File name: FVSWE1Dv2.m
% Author: Clive Mingham
% Date: 7 Oct 2011
% This m-file is associated with the free text book: 'Introductory Finite
% Volume Methods for PDEs' which may be downloaded from bookboon.com
% when published.
%
% Notes:
% 1. This prog is version 2 of a sequence of programs which gradually
% build up to coding a FV solver on a general 2D structured mesh.
%
% 2. This code is based on FVSWE1Dv1.m by Clive Mingham. The changes are
% that mesh cells are given by x and y coordinates of their vertices
% and cell areas and side vectors are stored in arrays and calculated
% for general quadrilateral cells although only left and right side
% vectors are found in this 1D case.
%
% Description: Solves the system of 1D PDEs, dU/dt + div(H) = 0,
% using a finite volume scheme on a single row of uniform Cartesian cells
% in the x direction.
%
% This equation corresponds to the finite volume form:
%
% doubleintegral(dU/dt dA) + lineintegral(H.n ds) = 0.
%
% In this case the equations are the 1D shallow water equations (SWE) so:
% Water depth: h=h(t,x) [m].
% Water velocity: v = vx i + 0 j where vx=vx(t,x) [m/s].
% Acceleration due to gravity: g=9.81 [m/s^2].
% Geopotential: phi = gh.
% Matrix of conserved variables: U = [ phi
% phi vx].
% Flux density: H = [ phi v i + [0 j
% phi vx^2 + 0.5*phi^2 ] 0]
%
% subfunctions: freadmesh, fcellareas, fsidevectors, finitialU,
% finterfaceflux, finsertghostcells, fcalcdt,
% fsolver, fremoveghostcells, fplotresults, fgetHifromUi.
%
% Verification: results look exactly the same as for FVSWE1Dv1.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
clc; clf; clear all; % clear everything
runtime=1.2; % runtime [s]
t=0; % current time [s]
timelevel=0; % current time level
[x,y,xcen,ycen]=freadmesh; % gets coordinates of the computational mesh
% and cell centres
A=fcellareas(x,y); % compute and store areas for each cell. Includes
% left and right ghost cells.
S=fsidevectors(x,y); % compute and store cell side vectors in an
% (NI+2)x2x2 array which contains left and
% right ghost cells.
%
U=finitialU(xcen); % Put initial cell centre values of U in a 2xNI array
%
Uinitial=U; % Store initial profile for plotting if required
U=finsertghostcells(U); % Append ghost cells to each end of the domain.
% This means that indices of computational cell
% start at i=2 and end at i = NI+1.
Unext=zeros(size(U)); % Creates an array for the updated U values
[d,N]=size(U); % d is the number of PDEs (=rows of U)
NI=N-2; % number of computational cells in the x direction
%%
disp('time marching starts')
while(t<runtime)
timelevel=timelevel+1
U=fbcs(U); % Insert boundary conditions into ghost cells at each end
% of the domain. U becomes 2x(NI+2) so each computational
% cell is indexed i=2 to NI+1.
dt=fcalcdt(A,S,U); % Finds stable time step before each iteration.
t=t+dt; % updates time
if (t>runtime) % ensures program stops at exactly t=runtime [s]
dt=dt-(t-runtime);
t=runtime;
end
t % Outputs current time to the Command Window
dt
Unext=fsolver(U,A,dt,S); % implements a time step of the solver
%
U(:,2:NI+1)=Unext(:,2:NI+1); % update computational U values prior to
% the next time iteration
end % of while loop
disp('1D finite volume scheme ended at time t')
t
U=fremoveghostcells(U);
fplotresults(xcen,U) % plot results
disp('finished main program, see figures')
end % FVSWE1Dv2
%%%%%%%%%%%%%%%%%%%%%%%%%%% subfunctions %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [x,y,xcen,ycen]=freadmesh
% Verification: verified.
% The mesh has NI uniform quadrilateral cells in the x direction and
% NJ=1 cell in the j direction (for a 1D problem).
% 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). 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)).
% The mesh starts at the origin so x(1,1)=0=y(1,1).
%
NI=100; % number of cells in the i direction
NJ=1; % number of cells in the j direction (=1 for a 1D problem)
x=zeros(NI+1,NJ+1); % allocate correct sized array for cell vertices
y=zeros(NI+1,NJ+1); % allocate correct sized array for cell vertices
xcen=zeros(NI,NJ); % allocate correct sized array for cell centres
ycen=zeros(NI,NJ); % allocate correct sized array for cell centres
%
dx=1.0; % cell length in x direction
dy=2.0; % cell length in y direction (arbitrary)
%
x(1,1)=0; % first row starting point for x
y(1,1)=0; % first row starting point for x
x(1,2)=0; % second row starting point for x
y(1,2)=dy; % second row starting point for y
% create uniform mesh
for i=2:NI+1
x(i,1)=x(1,1)+(i-1)*dx;
y(i,1)=y(1,1);
x(i,2)=x(1,2)+(i-1)*dx;
y(i,2)=y(1,2);
end % of i loop
% find cell centres by averaging coordinates of vertices
% (this works for a general structured mesh)
for i=1:NI
xcen(i,1)=(x(i,1)+x(i+1,1)+x(i+1,2)+x(i,2))/4;
ycen(i,1)=(y(i,1)+y(i+1,1)+y(i+1,2)+y(i,2))/4;
end % of i loop
end % freadmesh
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function A=fcellareas(x,y)
% General way to find quadrilateral cell areas.
% Defines A to include left and right ghost cells.
[m,n]=size(x);
NI=m-1; % number of computational cells in the x direction
NJ=n-1; % number of computational cells in the y direction
A=zeros(NI+2); % allocate correct size array including ghost cells
% inclusion of a left ghost cell means that indices of computational cells
% go from i=2 to i=NI+1
for i=1:NI
for j=1:NJ % j goes from 1 to 1
% write the diagonal vectors as 3D vectorrs for the cross product
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+1)=0.5*norm(d); % note shift of i index by 1 in A
end % of j loop
end % of i loop
end % of fcellareas
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function S=fsidevectors(x,y)
% Verification:
%
% Finds the left and right side vectors for each cell and stores in array
% S. The first component of S holds the cell index. The second component
% of S indicates the right (entry = 1) or left (entry = 2) side vector.
% The third component of S indicates the i (entry = 1) or j (entry = 2)
% component of the side vector.
%
% If a i + b j is the vector representing the cell side then the
% corresponding side vector is b i - a j.
% Notes:
% 1. S contains left and right ghost cells.
%
% 2. 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-1; % number of computational cells in the x direction
NJ=n-1; % number of computational cells in the y direction
S=zeros(NI+2,2,2); % includes ghost cells
% ghost cells shift the i index by 1 so computational indices
% go from 2 to NI+1
for i=1:NI
for j=1:NJ
rightcellside=[x(i+1,j+1)-x(i+1,j),y(i+1,j+1)-y(i+1,j)];
S(i+1,1,1)=rightcellside(2); % note shift of i index
S(i+1,1,2)=-rightcellside(1);
%
leftcellside=[x(i,j)-x(i,j+1),y(i,j)-y(i,j+1)];
S(i+1,2,1)=leftcellside(2);
S(i+1,2,2)=-leftcellside(1);
end % of j loop
end % of i loop
end % fsidevectors
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function U=finitialU(xcen)
% This sets up a dam break problem.
% Verification: correct.
% Inserts initial U values in computational cells (U is 2xNI)
NI=length(xcen); % number of computational cells in x direction
%
d=2; % number of PDEs
U=zeros(d,NI); % create correct sized array (without ghost cells)
dx=xcen(2)-xcen(1);
lengthofdomain=NI*dx;
xmidpoint=lengthofdomain/2;
g=9.81;
%
% 1D Dam break. Left and right states are either side
% of the centre of computational domain.
% Left states
hleft=10.0; % left water depth [m]
vleft=[0,0]; % left water velocity vector [m/s]
phileft=hleft*g; % left geopotential
% Right states
hright=5.0; % right water depth [m]
vright=[0,0]; % right water velocity vector [m/s]
phiright=hright*g; % right geopotential
%
% Fill in U values for each computational cell
for i=1:NI
if (xcen(i)<=xmidpoint)
U(1,i)=phileft;
U(2,i)=phileft*vleft(1);
else
U(1,i)=phiright;
U(2,i)=phiright*vright(1);
end % of if
end % of i loop
%
end % finitialU
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function U=finsertghostcells(U)
% Verification: verified.
% Inserts ghost cells at left and right ends of the computational domain.
% Hence computational indices go from i=2 to NI+1.
% Assumes that U has 2 rows (i.e. 2 PDEs) which is correct for the 1D SWE.
%
[d,NI]=size(U); % d is the number of equations (= rows of U).
leftghost=[0
0];
%
rightghost=[0
0];
%
U=[leftghost U rightghost];
end % finsertghostcells
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function U=fbcs(U)
% Verification: verified.
% Inserts boundary conditions in to the ghost cells at left and right ends
% of the domain. Note that ghost cells have already been appended to U and
% U is initially zero in each ghost cell.
%
% The following applies zero gradient boundary conditions everywhere.
[d,N]=size(U); % d is the number of equations (= rows of U).
% Left
U(:,1)=U(:,2); % copies the neighbouring interior values to the ghost cells
% Right
U(:,N)=U(:,N-1); % copies the neighbouring interior values to the ghost cells
end % fbcs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Unext=fsolver(U,A,dt,S)
% Implements the finite volume scheme for 1 time step
[d,N]=size(U); % includes ghost cells
NI=N-2; % number of computational cells
for i=2:NI+1 % loop over computational cells
SRight=[S(i,1,1),S(i,1,2)]; % right side vector
SLeft=[S(i,2,1),S(i,2,2)];
[IHRight,IHLeft]=finterfaceflux(U,i); % finds the left and right
% interface fluxes for cell i
% Update U to Unext
for k=1:d
% compute total flux out of cell i for equation k
IHR=[IHRight(k,1),IHRight(k,2)]; % right flux vector
IHL=[IHLeft(k,1),IHLeft(k,2)]; % left flux vector
totalfluxout=dot(IHR,SRight)+dot(IHL,SLeft);
% Uses Lax-Friedrichs
Unext(k,i)=0.5*(U(k,i-1)+U(k,i+1))-(dt/A(i))*totalfluxout;
end % k loop
end % of i loop
end % fsolver
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [IHRight,IHLeft]=finterfaceflux(U,i)
% Interface H values are based on linear interpolation of neighbouring cell
% centre data which is the same as averaging on a uniform mesh.
% These may be found either by finding the average value of U at the cell
% at the interface and then constructing H or by direct averaging of H
% components. Since the PDEs are non-linear these approaches are not the same.
% The latter appoach is analogous to the Lax-Friedrichs finite
% difference method and is the better method as the former method produces
% unwanted oscillations.
%
% Interface fluxes are denoted IHRight, IHLeft.
% For the 1D shallow water equations H is:
% [ phi*vx i + [ 0 j
% phi*vx^2+phi^2/2 ] 0 ]
[d,m]=size(U); % d=number of rows of U (= number of PDEs = 2)
NI=m-2; % number of computational cells
IHRight=zeros(d,2); % array to store the right interface flux for a cell.
% IHRight(,1) is the i component of H and
% IHRight(,2) is the j component of H which is zero.
IHLeft=zeros(d,2);
%
% Find interface values of H directly by averaging neighbouring components.
% Cell i
Hi=fgetHifromUi(U,i);
% Cell i-1
Himinus1=fgetHifromUi(U,i-1);
% Cell i plus 1
Hiplus1=fgetHifromUi(U,i+1);
%
% Right interface is the average of H components in cells i and i+1.
IHRight=(Hiplus1+Hi)/2;
% Left interface is the average of H components in cells i-1 and i.
IHLeft=(Himinus1+Hi)/2;
end % finterfaceflux
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function outarray=fremoveghostcells(inarray)
% Verification: Verified
% Removes left and right ghost cells.
[d,m]=size(inarray);
outarray=inarray(:,2:m-1);
end % fremoveghostcells
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function dt=fcalcdt(A,S,U)
% Verification: not verified.
% Finds allowable time step dt using heuristic formula.
% 0<F<1 is a tunable safety factor.
[nrows,M]=size(U); % U now enlarged to include ghost cells
dt=9999999; % large value of dt to start
F=0.9; % safety factor
for i=2:M-1 % indices go over computational cells
phi=U(1,i);
cel=sqrt(phi); % wave celerity
vx=U(2,i)/phi; % flow speed in x direction
v=[vx,0]; % flow velocity
SRight=[S(i,1,1),S(i,1,2)];
a=dot(v,SRight);
b=cel*norm(SRight);
%
speed=max(abs(a+b),abs(a-b)); % max abs wave speed
localdt=A(i)/speed;
% take smallest value of dt
if (localdt<dt)
dt=localdt;
end
end % i loop
dt=F*dt;
end % fcalcdt
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function fplotresults(xcen,U)
% Plots results on a structured mesh - even works for non-Cartesian!
%
g=9.81; % acceleration due to gravity
waterdepth=U(1,:)/g;
waterspeed=U(2,:)./U(1,:);
subplot(2,1,1), plot(xcen,waterdepth)
title('FVSWE1Dv2.m: numerical solution for water depth')
xlabel('x [m]')
ylabel('water depth [m]')
%
subplot(2,1,2), plot(xcen,waterspeed)
title('FVSWE1Dv2.m: numerical solution for x component of water velocity')
xlabel('x [m]')
ylabel('speed [m/s]')
%
end % fplotresults
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function H=fgetHifromUi(U,i)
% Constructs flux vector H at centre of cell i from U(i)
% U = [ phi
% fi*vx ]
%
% H=[ phi*vx , 0 Each row of H is a 2D vector
% phi*vx^2+fi^2/2, 0 ]
%
[d,n]=size(U); % d is the number of PDEs (= rows of U)
H=zeros(d,2); % flux density vectors
%
phi=U(1,i); % geopotential
vx=U(2,i)/phi; % velocity component in x direction
%
H(1,1)=U(2,i); % i component of first row of H
H(1,2)=0; % j component of first row of H
%
H(2,1)=phi*vx*vx+0.5*phi*phi; % i component of second row of H
H(2,2)=0; % j component of second row of H
%
end % fgetHifromUi
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%