% This Matlab script solves the two-dimensional convection
% equation using a finite volume algorithm.  
%

clear all;

% Set flag for plotting
doplot = 1;

if (doplot == 1),
  figure(1); clf;
  figure(2); clf;
  figure(3); clf;
end

xymult = 2;

% Specify x range and number of points
x0 = -2;
x1 =  2;
Nx = xymult*20;

% Specify y range and number of points
y0 = -2;
y1 =  2;
Ny = xymult*20;

% Construct mesh
x       = linspace(x0,x1,Nx+1);
y       = linspace(y0,y1,Ny+1);
[xg,yg] = ndgrid(x,y);
  
% Construct mesh needed for plotting
xp = zeros(4,Nx*Ny);
yp = zeros(4,Nx*Ny);
n = 0;
for j = 1:Ny,
  for i = 1:Nx,
    
    n = n + 1;
    xp(1,n) = x(i);
    yp(1,n) = y(j);

    xp(2,n) = x(i+1);
    yp(2,n) = y(j);
    
    xp(3,n) = x(i+1);
    yp(3,n) = y(j+1);

    xp(4,n) = x(i);
    yp(4,n) = y(j+1);

  end
end

% Calculate midpoint values in each control volume
xmid = 0.5*(x(1:Nx) + x(2:Nx+1));
ymid = 0.5*(y(1:Ny) + y(2:Ny+1));

% Calculate cell size in control volumes (assumed equal)
dx = x(2) - x(1);
dy = y(2) - y(1);
A  = dx*dy;

% Set velocity 
u =  1;
v = -1;

% Set tolerances on iteration and residual tolerance
ntol = 3000;
Rtol = -13;

% Set timestep
CFL = 1.0;
dt = CFL/(abs(u)/dx + abs(v)/dy);

% Set initial conditions
U = zeros(Nx,Ny);
n = 0;
Rmax = 0;

% Loop until one of the termination criterion are met
while ((n < ntol) & (Rmax > Rtol)),
  
  % The following implement the bc's by creating a larger array
  % for U and putting the appropriate values in the first and last
  % columns or rows to set the correct bc's
  Ubc(2:Nx+1,2:Ny+1) = U; 
  Ubc(   1,2:Ny+1) = exp(-10*(ymid-1).^2); %
  Ubc(Nx+2,2:Ny+1) = zeros(1, Ny); %
  Ubc(2:Nx+1,   1) = zeros(Nx, 1); %
  Ubc(2:Nx+1,Ny+2) = exp(-10*(y1-1)^2); %
  
  % Calculate the flux at each interface

  % First the i interfaces
  F =   0.5*    u *( Ubc(2:Nx+2,2:Ny+1) + Ubc(1:Nx+1,2:Ny+1)) ...
      - 0.5*abs(u)*( Ubc(2:Nx+2,2:Ny+1) - Ubc(1:Nx+1,2:Ny+1));
  
  % Now the j interfaces
  G =   0.5*    v *( Ubc(2:Nx+1,2:Ny+2) + Ubc(2:Nx+1,1:Ny+1)) ...
      - 0.5*abs(v)*( Ubc(2:Nx+1,2:Ny+2) - Ubc(2:Nx+1,1:Ny+1));
  
  % Add contributions to residuals from fluxes
  R = (F(2:Nx+1,:) - F(1:Nx,:))*dy + (G(:,2:Ny+1) - G(:,1:Ny))*dx;
  
  % Calculate maximum residual
  Rmax = log10(max(max(abs(R))));
  
  % Forward Euler step
  U = U - (dt/A)*R;

  % Increment iteration
  n = n + 1;
  fprintf('Iter: %i, Log10 Max Res = %f\n',n,Rmax);

  if (doplot == 1),
    
    % Plot current solution
    Up = reshape(U,1,Nx*Ny);
    figure(1);
    clf;
    [Hp] = patch(xp,yp,Up);
    set(Hp,'EdgeAlpha',0);
    axis('equal');
    title(sprintf('Log10 Max Res = %f\n',Rmax));
    xlabel('x'); ylabel('y');
    caxis([0,1]);
    colorbar;
    
    % Plot max residual
    figure(2);
    hold on; plot(n,Rmax,'*'); hold off;
    xlabel('Iteration');
    ylabel('Log10 Rmax');
    
    % Plot solution at bottom of domain
    figure(3);
    stairs(x,[U(:,1)', U(Nx,1)]);
    xlabel('x');
    ylabel('U');
    axis([x0,x1,0,1]);
    title(sprintf('U at y = %f\n',y0));
  
    drawnow;
    
  end

end