% godunov.m: Solve the Burgers' equation using the first-order
%            Godunov method.
%
% Date:   March 7, 2011
% Author: John Stockie
%         Department of Mathematics
%         Simon Fraser University
%

function godunov( nx, dt, ictype )

if nargin < 1, nx    = 100; end;
if nargin < 2, dt    = 0.1; end;
if nargin < 3, ictype= 1;   end;  % 1 = shock; 2 = expansion; 
                                  % 3 = sonic expansion; 4 = cosine.

tend = 4;     % end time
xmax = 8;     % domain length [0,xmax]
dx = xmax/nx; % mesh spacing (constant)
x  = [0 : dx : xmax];
nt = floor(tend/dt);
dt = tend / nt;
ntprint = 50; % for printing

% Set up the initial solution values.
u0   = uinit(x,ictype);
u    = u0;
unew = 0*u;
us   = unew(1:end-1);

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

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

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

ip = 1;
figure(1)
for i = 1 : nt,
  t = i*dt;

  us = rpbu(u(1:end-1), u(2:end));
  unew(2:end-1) = u(2:end-1) + dt/dx * (f(us(1:end-1)) - f(us(2:end)));
  unew(1)   = u(1);
  unew(end) = u(end);
    
  % Plot the solution profiles.
  if t >= ip*dtprint,
    plot(x, unew)
    xlabel('x'), ylabel('u')
    title( ['Solution at time t=', num2str(t,'%9.4f')] )
    grid on, shg
    pause(0.1)
    ip = ip + 1;
    uall(ip,:) = unew;
  end
  u = unew;
end

uall(ip+1,:) = unew;  % Be sure to plot the final solution.

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Burgers flux.
function ret = f( u )
ret = 0.5 * u.^2;

% Derivative of Burgers flux.
function ret = df( u )
ret = u;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Initial conditions.
function ret = uinit( x, ictype )

xshift = 1.0;
switch ictype
 case 1    % shock (uL > uR)
  uL = 1.0;
  uR = 0.5;
  ret = uR + (uL-uR) * ((x-xshift) <= 0.0);
 case 2    % expansion (uL < uR)
  uL = 0.5;
  uR = 1.0;
  ret = uR + (uL-uR) * ((x-xshift) <= 0.0);
 case 3    % sonic expansion (uL < 0 < uR)
  uL =-0.25; 
  uR = 0.5;
  xshift = 4.0;
  ret = uR + (uL-uR) * ((x-xshift) <= 0.0);
 otherwise % cosine (forms shock later)
  ret = 0.5*(1 + cos(pi*(x-xshift))) .* (abs(x-xshift) <= 1);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compact version of the Rieman solver for Burgers' equation.
% Returns the exact solution along x/t = 0.
function ret = rpbu( uL, uR )
usonic = 0.0;        % sonic point
s = 0.5 * (uL + uR); % shock speed from Rankine-Hugoniot
issonic = (uL < 0) & (0 < uR);
ret = issonic * usonic + (~issonic) .* (uL .* (s >= 0) + uR .* (s < 0));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Long version of the Rieman solver for Burgers' equation.
% Returns the exact solution along x/t = 0.
function ret = rpbu2( uL, uR )
if uL > 0,
  if uR > 0,
    ret = uL;
  else
    s = 0.5 * (uL + uR);
    if s > 0,
      ret = uL;
    else
      ret = uR;
    end
  end
else
  if uR > 0,
    ret = 0;
  else
    ret = uR;
  end
end
