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