function FVSWE1Dv1
% File name: FVSWE1Dv1.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.
%
% Note: This prog is version one of a sequence of programs which gradually
% build up to coding a FV solver on a general 2D structured mesh.
% A simple approach is used in this code which can be generalized later.
%
% 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: finitialU, finterfaceflux, finsertghostcells, fcalcdt,
% fsolver, fremoveghostcells, fplotresults, fgetHifromUi.
%
% Verification: Results match those from the equivalent finite difference
% (Lax-Friedrichs) code which ought to be the case as the
% code is actually the same on a uniform Cartesian mesh.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
clc; clf; clear all; % clear everything
runtime=1.2; % runtime in seconds
t=0; % current time
timelevel=0; % current time level
p=0; % start of computational domain [m]
q=100; % end of computational domain [m]
NI=100; % Number of uniform rectangular dx by dy cells in the x direction
dx=(q-p)/NI; % cell length in x direction
dy=3; % cell length in y direction (arbitrary as it cancels out)
xcen=dx/2:dx:(NI*dx-dx/2); % cell centres
A=dx*dy; % constant cell areas [m^2]
SRight=[dy,0]; % constant right cell side vectors
SLeft=[-dy,0]; % constant left cell side vectors
%
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)
%%
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,SRight,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,SRight,SLeft); % 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 % FVSWE1Dv1
%%%%%%%%%%%%%%%%%%%%%%%%%%% subfunctions %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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,SRight,SLeft)
% 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
[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)*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,SRight,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
a=dot(v,SRight);
b=cel*norm(SRight);
%
speed=max(abs(a+b),abs(a-b)); % max abs wave speed
localdt=A/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('FVSWE1Dv1.m: numerical solution for water depth')
xlabel('x [m]')
ylabel('water depth [m]')
%
subplot(2,1,2), plot(xcen,waterspeed)
title('FVSWE1Dv1.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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%