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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%