% upwind.m: Solve the linear advection equation with
%           the upwind scheme (first order in space and time)
%
% Date:   January 19, 2011
% Author: John Stockie
%         Department of Mathematics
%         Simon Fraser University
%

function upwind( CFL )

if nargin == 0,
  disp( 'Usage:  upwind( CFL )' )
  return;
end

a    = 1.0;   % advection speed
tend = 5;     % end time
xmax = 10;    % domain length [0,xmax]
nx   = 100;   % number of grid points

dx = xmax/nx;   % mesh spacing (constant)
x  = [0 : dx : xmax];
ntprint = 50;   % for plotting
xex= [0 : xmax/1000 : xmax];

% Calculate the largest time step (dt) allowed 
% so that the CFL condition is still satisfied.
% A constant dt is used (the largest that still 
% satisfies the CFL condition) so that the actual 
% CFL number may be different than what is input.
maxdt= CFL * dx / a;
nt   = ceil( tend / maxdt );   
dt   = tend / nt;      % constant time step
sigma= a * dt / dx;    % actual CFL number used

% Set up the initial solution values.
u0   = uexact( x, 0, a );
u    = u0;
unew = 0*u;

disp( ['    a = ', num2str(a)] );
disp( ['   dx = ', num2str(dx)] );
disp( ['   dt = ', num2str(dt)] );
disp( ['  CFL = ', num2str(sigma)] );

ntprint = min(nt, ntprint);
dtprint = tend / ntprint;

% Step the solution in time.
uall = zeros(ntprint+1,nx+1);
uall(1,:) = u0;
uerr = zeros(1,nt);
ip = 1;

figure(1)

for i = 1 : nt,

  t = i*dt;

  % Apply the upwind method.
  unew(2:end) = u(2:end) - sigma * (u(2:end) - u(1:end-1));
  unew(1) = 0.0;  % boundary condition
  
  % Calculate the error in the L1-norm.
  uerr(i) = sum( dx*abs(unew - uexact(x,t,a)) );
  
  % Plot the solution profiles.
  if t >= ip*dtprint,
    uex = uexact( xex, t, a );
    plot(x, unew)
    xlabel('x'), ylabel('u')
    title( ['Solution at time t=', num2str(t,'%9.4f')] )
    hold on, plot(xex,uex,'r:'), hold off
    grid on, shg
    pause(0.1)
    ip = ip + 1;
    uall(ip,:) = unew;
  end
  
  u = unew;

end

% Plot a sequence of profiles throughout the run.
figure(2)
nskip = 3;
plot(x,uall(1:nskip:end,:));
xlabel('x'), ylabel('u')
grid on

% Plot the L1-error.
figure(3)
plot(dt*[1:nt], uerr);
xlabel('t'), ylabel('L1-error')
% Uncomment the next line to get a log-log plot of error
% (or replace 'plot' with 'loglog').
%%set(gca,'Xscale','log','Yscale','log') 
grid on

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% uexact: This is the exact solution (also used for IC's)
function uex = uexact( x, t, a )
u0shift = 2;
uex = 0.5 * (1 + sin(2*pi*(x-1/4-u0shift-a*t))) .* ...
      (abs(x-0.5-u0shift-a*t) < 0.5);
