function FVlinearadvectionFOU2D
% File name: FVlinearadvectionFOU2D.m
% Author: Clive Mingham
% Date: 24 Feb 2011
% Description:  Solves the 2D PDE, dU/dt + div(H) = 0
% using the FOU finite volume scheme.
% This equation corresponds to the finite volume form:
% doubleintegral(dU/dt dA) + lineintegral(H.n ds) = 0.
% In this case the equation is the linear advection equation so:
% flux density: H = vU 
% velocity: v = vx i + vy j = <vx,vy> where vx and vy are constant.
% 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,
%               fdrawmesh
% Verification:  (on uniform Cartesian mesh)
% for NI=61,NJ=81,dx=2,dy=4, initial profile centre is at (61,82) 
% contour plots of numerical agree with exact for: 
% vx=10,vy=0, t=1.04, new centre (71, 82) is correct
% vx=0, vy=10, t=1.12, new centre (61, 94) is correct
% vx=10,vy=10, t=1.04, new centre (71, 94) is correct
% Comments:  Although program appears to work the numerical scheme is highly
%            dissipative.

clc; clf; clear all;
runtime=1.0;     % target runtime [s]
t=0.0;           % time [s]
tlevel=0;        % time level
vx=10;   % velocity component in x direction [m/s]
vy=10;   % velocity component in y direction [m/s]
v=[vx,vy];       % velocity vector <vx,vy>
[x,y,xcen,ycen,solid]=freadmesh; % gets coordinates of the 2D computational mesh 
                     % cell centres and solid flags (excluding ghost cells)
disp('read in mesh')
NI=m-1; % number of computational cells in i direction.
NJ=n-1; % number of computational cells in j direction.
A=fcellarea(x,y);    % compute and store cell areas in an NIxNJ array
disp('calculated cell areas')
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(xcen,ycen,0,v); % put initial cell centre values of u in NIxNJ 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)
unext=zeros(size(u));   % (NI+2)x(NJ+2) array for u values at next time level
A=finsertghostcells(A); % A is now (NI+2)x(NJ+2)
disp('time marching starts')   
 u=fbcs(u);  % Implement boundary conditions using ghost cells. 
             % u is (NI+2)x(NJ+2) and each computational cell is 
             % indexed i=2 to NI+1, j=2 to NJ+1.
 dt=fcalcdt(A,S,v); % Finds stable time step for each iteration.
 t=t+dt  % update time 
 for j=2:NJ+1 
   for i=2:NI+1 
       IH=fIflux(v,u,i,j); % gets interface fluxes for cell (i,j)
       IH1=[IH(1,1),IH(1,2)]; % interface flux vector for side 1
       IH2=[IH(2,1),IH(2,2)]; % interface flux vector for side 2
       IH3=[IH(3,1),IH(3,2)]; % interface flux vector for side 3
       IH4=[IH(4,1),IH(4,2)]; % interface flux vector for side 4
       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
       % FV scheme
       % compute total flux out of cell (i,j)
   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
end  % of while loop
disp('time marching ends') 
uexact=finitialu(xcen,ycen,t,v); % exact solution at time t
u=u(2:NI+1,2:NJ+1); % extract final computational values
fdrawmesh(x,y,solid) % draws mesh if fplotresults is commented out
%fdisplay(u); % print results as a matrix for a small mesh if needed
fplotresults(xcen,ycen,uinitial,u,uexact) % plot results: comment in or out
disp('finished main program, see figures')
end % FVlinearadvectionFOU2D
%%%%%%%%%%%%%%%%%%%%%%%%%%% subfunctions  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [x,y,xcen,ycen,solid]=freadmesh
% The mesh is structured and has NI cells in the i direction and 
% NJ 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).  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)).  xcen and ycen are NI by NJ arrays
% which hold cell centres. 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.
% % test using diamond
% x(1,1)=0; y(1,1)=-1;
% x(2,1)=1; y(2,1)=0;
% x(2,2)=0; y(2,2)=1;
% x(1,2)=-1; y(1,2)=0;
% Test using regular Cartesian mesh starting at (0,0).  i direction is
% x direction.
% NI=61,NJ=41,dx=2,dy=4 makes mesh centre (61,82)
x=zeros(NI+1,NJ+1); % allocate correct sized array
y=zeros(NI+1,NJ+1); % allocate correct sized array
xcen=zeros(NI,NJ);  % allocate correct sized array
ycen=zeros(NI,NJ);  % allocate correct sized array
solid=zeros(NI,NJ); % allocate correct sized array
% create cell vertices
for i=1:NI+1
    for j=1:NJ+1
    end % of j loop
end % of i loop
nonCartesian=1; % 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
    end % of j loop
   end % of i loop
   disp('mesh is Cartesian') 
end % end of if statement create mesh
% find cell centres by averaging coordinates of vertices
for i=1:NI
    for j=1:NJ
    end % of j loop
end % of i loop
end % freadmesh
function A=fcellarea(x,y)  
% Verification:  correct for single diamond
% cell area = half the modulus of cross product of cell diagonal vectors.
% These vectors are 2D so append 0 k component to make 3D for cross product
NI=m-1; NJ=n-1;
for i=1:NI
    for j=1:NJ
    end % of j loop
end % of i loop
end % fcellarea
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 has ghost cells around it so that side vector 
% indices (i+1,j+1) refer to cell (i,j).
% Note that vectors are 1 by 2 arrays so care must be taken because side
% vectors are stored in higher dimensional arrays.
NI=m-1; NJ=n-1;
for i=1:NI
    for j=1:NJ
    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(xcen,ycen,t,v)   
% Inserts exact u values (and initial conditions for t=0) 
u=zeros(NI,NJ);  % create correct sized array (without ghost values)
% Initial 2D Gaussian function based on the centre (xc,yc) of the 
% computational domain.
xmax=max(max(xcen));  % max x value in mesh
ymax=max(max(ycen));  % max y value in mesh
xmin=min(min(xcen));  % min x value in mesh
ymin=min(min(ycen));  % min y value in mesh
xc=(xmax+xmin)/2;  % approx x coord of mesh centre
yc=(ymax+ymin)/2;  % approx y coord of mesh centre
cutoff=min(xmax-xc,ymax-yc)/2;  % cut off radius for Gaussian
for i=1:NI
    for j=1:NJ
        d=sqrt((x-xc)^2+(y-yc)^2); % distance from centre
        if (d<cutoff)
    end % of j loop
end % of i loop
end % finitialu
function IH=fIflux(v,u,i,j)
% Calculates the 4 interface fluxes on the 4 sides of cell (i,j). 
% 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)).
% 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 linear advection equation, H(U) = <vx U, vy U>.
% There is considerable choice for interface flux estimation: different 
% choices give different FV schemes.
%% The following scheme simply takes:
%  IH1=H(u(i,j-1)) = v u(i,j-1) 
%  IH2=H(u(i,j)) = v u(i,j)
%  IH3=H(u(i,j)) = v u(i,j)
%  IH4=H(u(i-1,j)) = v u(i-1,j)
% This is an upwind scheme for vx > 0 and vy >0 and corresponds to the 
% FOU finite difference scheme on a Cartesian mesh.
IH=zeros(4,2); % array to store all 4 interface flux vectors for a cell
IH(1,1)=vx*u(i,j-1);  % component of IH1 in the x direction
IH(1,2)=vy*u(i,j-1);  % component of IH1 in the y direction 
IH(2,1)=vx*u(i,j);    % component of IH2 in the x direction
IH(2,2)=vy*u(i,j);    % component of IH2 in the y direction
IH(3,1)=vx*u(i,j);    % component of IH3 in the x direction
IH(3,2)=vy*u(i,j);    % component of IH3 in the y direction 
IH(4,1)=vx*u(i-1,j);  % component of IH4 in the x direction
IH(4,2)=vy*u(i-1,j);  % component of IH4 in the y direction
end % fIflux
function outarray=finsertghostcells(inarray)
% Assuming that ghost cells are needed around the whole domain an mxn
% array is embedded into (m+2)x(n+2) array of zeros.  Hence computational 
% indices go from i=2 to m+1 and j=2 to n+1.
end % finsertghostcells
% function outarray=fremoveghostcells(inarray)
% % Removes ghost cells so that a mxn array becomes (m-2)x(n-2)
% [m,n]=size(inarray);
% outarray=inarray(2:m-1,2:n-1);
% end % fremoveghostcells
function u=fbcs(u)
% Implements boundary conditions using ghost cells index by i=1, i=NI+2
% j=1, j=NJ+2.
NI=m-2; NJ=n-2;
% Dirichlet: u=0 everywhere.
% Don't need to do anything as u values in ghost cells are already zero.
end % bcs
function dt=fcalcdt(A,S,v)
% Verification:
% 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
NI=m-2; NJ=n-2; 
dt=9999999;  % large value of dt to start
 for i=2:NI+1
    for j=2:NJ+1
        % put side vectors into 1x2 arrays  
        v2=abs(dot(v,s2)); % side 2
        v3=abs(dot(v,s3)); % side 3
        % take smallest value of dt
        if (dum<dt)
    end % j loop
 end % i loop
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.
NI=m-1; % number of cells in i direction
NJ=n-1; % number of cells in j direction
% Plot the mesh
hold on   % keeps all plots on the same axes
% draw lines in the i direction
for i=1:NI+1
% draw lines in the j direction
for j=1:NJ+1
title('computational mesh')
% Fill in solid cells
for i=1:NI
    for j=1:NJ
        if (solid(i,j)==1)
    end % of j loop
end % of i loop
end % fdrawmesh
function fplotresults(x,y,uinitial,u,uexact)
% Plots results on a structured mesh - even works for non-Cartesian!
subplot(3,1,1), contour(x,y,uinitial)
title('Contour plot of initial conditions')
xlabel('x [m]')
ylabel('y [m]')
subplot(3,1,2), contour(x,y,u)
title('Contour plot of numerical solution')
xlabel('x [m]')
ylabel('y [m]')
subplot(3,1,3), contour(x,y,uexact)
title('Contour plot of exact solution')
xlabel('x [m]')
ylabel('y [m]')
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
 disp('printing as a matrix')
end % fdisplay