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]
%vx=0;vy=10;
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')
[m,n]=size(x); 
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')   
while(t<runtime) 
 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 
 tlevel=tlevel+1;  
 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)
       totalfluxout=dot(IH1,s1)+dot(IH2,s2)+dot(IH3,s3)+dot(IH4,s4);
       unext(i,j)=u(i,j)-(dt/A(i,j))*totalfluxout; 
   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)
NI=61;
NJ=41;
dx=2;
dy=4;
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
        x(i,j)=(i-1)*dx;
        y(i,j)=(j-1)*dy;
    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
        x(i,j)=x(i,j)+(dx/6)*cos(1.727*(i+j));
        y(i,j)=y(i,j)+(dy/5)*sin(2.46723*i*j);
    end % of j loop
   end % of i loop
else
   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
        xcen(i,j)=(x(i,j)+x(i+1,j)+x(i+1,j+1)+x(i,j+1))/4;
        ycen(i,j)=(y(i,j)+y(i+1,j)+y(i+1,j+1)+y(i,j+1))/4;
    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
[m,n]=size(x);
NI=m-1; NJ=n-1;
A=zeros(NI,NJ);
for i=1:NI
    for j=1:NJ
        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,j)=0.5*sqrt(dot(d,d));
    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.
%
[m,n]=size(x);
NI=m-1; NJ=n-1;
S=zeros(NI+2,NJ+2,4,2);
for i=1:NI
    for j=1:NJ
        s1=[x(i+1,j)-x(i,j),y(i+1,j)-y(i,j)];
        S(i+1,j+1,1,1)=s1(2);
        S(i+1,j+1,1,2)=-s1(1);
        %
        s2=[x(i+1,j+1)-x(i+1,j),y(i+1,j+1)-y(i+1,j)];
        S(i+1,j+1,2,1)=s2(2);
        S(i+1,j+1,2,2)=-s2(1);
        %
        s3=[x(i,j+1)-x(i+1,j+1),y(i,j+1)-y(i+1,j+1)];
        S(i+1,j+1,3,1)=s3(2);
        S(i+1,j+1,3,2)=-s3(1);
        %
        s4=[x(i,j)-x(i,j+1),y(i,j)-y(i,j+1)];
        S(i+1,j+1,4,1)=s4(2);
        S(i+1,j+1,4,2)=-s4(1);    
    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) 
[NI,NJ]=size(xcen);
u=zeros(NI,NJ);  % create correct sized array (without ghost values)
vx=v(1);
vy=v(2);
%
% 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
        x=xcen(i,j)-vx*t; 
        y=ycen(i,j)-vy*t;
        d=sqrt((x-xc)^2+(y-yc)^2); % distance from centre
        if (d<cutoff)
           u(i,j)=exp(-0.01*d^2);
        else
           u(i,j)=0;
        end
    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
%
vx=v(1);
vy=v(2);
%
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.
[m,n]=size(inarray);
dummy=zeros(m+2,n+2);
dummy(2:m+1,2:n+1)=inarray;
outarray=dummy;
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.
[m,n]=size(u);
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
F=0.4;
 for i=2:NI+1
    for j=2:NJ+1
        % put side vectors into 1x2 arrays  
        s2=[S(i,j,2,1),S(i,j,2,2)];
        s3=[S(i,j,3,1),S(i,j,3,2)];
        %
        v2=abs(dot(v,s2)); % side 2
        v3=abs(dot(v,s3)); % side 3
        dum=A(i,j)/max(v2,v3);
        % take smallest value of dt
        if (dum<dt)
           dt=dum;
        end
    end % j loop
 end % i loop
dt=F*dt;
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.
% 
[m,n]=size(x);
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
    plot(x(i,:),y(i,:))
end
% draw lines in the j direction
for j=1:NJ+1
    plot(x(:,j),y(:,j))   
end
title('computational mesh')
xlabel('x')
ylabel('y')
% Fill in solid cells
for i=1:NI
    for j=1:NJ
        if (solid(i,j)==1)
            solidx=[x(i,j),x(i+1,j),x(i+1,j+1),x(i,j+1),x(i,j)];
            solidy=[y(i,j),y(i+1,j),y(i+1,j+1),y(i,j+1),y(i,j)];
            fill(solidx,solidy,'k')      
        end
    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
 [NX,NY]=size(in);
 temp=transpose(in);
 out=temp(NY:-1:1,:);
 disp('printing as a matrix')
 out
end % fdisplay
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%