function temple8024_godunov_shallow_water
%TEMPLE8024_GONDUNOV_SHALLOW_WATER
%   Solves the 1d shallow water equations using Godunov's
%   method with Roe averaging. This version of the method
%   does not have an entropy fix in it, so spurious jumps
%   can arise at transsonic rarefactions.
%
%   02/2017 by Benjamin Seibold
%   http://www.math.temple.edu/~seibold
%   
%   This code may be used and modified for teaching,
%   learning, and research, as long as its orgin is credited.

%===========================================================
ax = [-1,1]; % domain
n =  1000; % number of grid cells
tf =  .5; % final time
cfl = .99; % CFL number
%===========================================================
h = diff(ax)/n; x = ax(1)+h/2:h:ax(2)-h/2; % cell centers
u = ic(x); u0 = u; % initial conditions
t = 0; % initial time
while t<tf % time loop
    % Godunov step
    uL = u(:,[1 1:end]); uR = u(:,[1:end end]); % extend at boundaries
    u_avg = roe_avg(uL,uR); % Roe-averaged stages on cell interfaces
    delta = uR-uL; % jump across cell interface
    [lambda1,lambda2,r1,l1] = spectral(u_avg); % eigs of Jacobian
    c = max(max(abs(lambda1),abs(lambda2))); % maximum wave speed
    k = min(cfl/c*h,tf-t); % time step
    alpha = sum(l1.*delta); % jumps in characteristic variables
    uM = uL+([1;1]*alpha).*r1; % middle states
    us = ([1;1]*(lambda1>=0)).*uL+... % Godunov state at x/t=0
        ([1;1]*(lambda1<0&lambda2>0)).*uM+...
        ([1;1]*(lambda2<=0)).*uR;
    fs = f(us); % nonlinear flux through cell interfaces
    u = u-k/h*(fs(:,2:end)-fs(:,1:end-1)); % update state
    t = t+k; % update time
    % plot
    clf
    subplot(2,1,1), plot(x,u0(1,:),'c.-',x,u(1,:),'r.-')
    title(sprintf('Water height at t = %0.2f',t))
    xlabel('x'), ylabel('h')
    subplot(2,1,2), plot(x,u0(2,:)./u0(1,:),'c.-',x,u(2,:)./u(1,:),'r.-')
    title('Water velocity'), xlabel('x'), ylabel('v')
    drawnow
end

%===========================================================
% Problem specific routines: shallow water equations
%===========================================================
function y = ic(x) % initial conditions
y = [.1+.8*(x<0);x*0];

function y = f(u) % flux function
y = [u(2,:);u(2,:).^2./u(1,:)+1/2*u(1,:).^2];

function u = roe_avg(uL,uR) % Roe-averaging
vL = uL(2,:)./uL(1,:); vR = uR(2,:)./uR(1,:); % velocities
wL = sqrt(uL(1,:)); wR = sqrt(uR(1,:)); % Roe weights
h_avg = (uL(1,:)+uR(1,:))/2; % averaged heights
v_avg = (wL.*vL+wR.*vR)./(wL+wR); % averaged velocities
u = [h_avg;v_avg.*h_avg]; % averaged state vector

function [lambda1,lambda2,r1,l1] = spectral(u) % properties of Jacobian
v = u(2,:)./u(1,:); % velocity
c = sqrt(u(1,:)); % speed of sound
e = c*0+1; % all ones
lambda1 = v-c; lambda2 = v+c; % eigenvalues
r1 = [e;lambda1]; l1 = [lambda2./(2*c);-e./(2*c)]; % first eigenvectors
