function temple8024_weno_claw
%TEMPLE8024_WENO_CLAW
%   Solves one-dimensional scalar conservation laws
%   u_t+f(u)_x = 0 by a semidiscrete approach, in which
%   5th order WENO is used for the reconstruction of the
%   Riemann states at cell boundaries, and the 3rd order
%   SSP Shu-Osher scheme is used for the time stepping.
%
%   10/2010 by Benjamin Seibold
%   http://www.math.temple.edu/~seibold
%   
%   This code may be used and modified for teaching,
%   learning, and research, given its orgin is credited.

%===========================================================
n =   100; % number of grid cells
tf =    5; % final time
cfl = .99; % CFL number
%===========================================================
h = 1/n; x = h/2:h:1-h/2; % cell centers
u = ic(x); % initial conditions
c = max(abs(f1(u))); % maximum characteristic velocity
nt = ceil(tf*c/h/cfl); k = tf/nt; la = k/h;
uax = [min(u),max(u)]+[-1,1]*(max(u)-min(u))*.1;
for j = 0:nt
    if j>0 % compute
        u = shu_osher_step(u,la);
    end
    % plot
    xp = reshape([x-h/2;x+h/2],1,[]);
    [vL,vR] = weno_twoside(u);
    up = reshape([vL;vR],1,[]);
    clf, plot(xp,up,'b-'), axis([0 1 uax])
    title(sprintf('t = %0.1f',j*k))
    xlabel('x'), ylabel('u')
    if j==0, pause(1.5), else drawnow, end
end

%===========================================================
% Problem specific routines
%===========================================================
function v = f(u)
% Flux function
% Here: Burger's equation
v = u.^2/2;

function u = ic(x)
% Initial conditions
u = .4*cos(2*pi*x)+.3*sin(6*pi*x)+.8*(.3<x&x<.7);

function v = shl(u)
% Shift solution left by one cell, considering
% boundary conditions
%v = u([1 1:end-1]); % Neumann b.c.
v = u([end 1:end-1]); % periodic b.c.

function v = shr(u)
% Shift solution right by one cell, considering
% boundary conditions
%v = u([2:end end]); % Neumann b.c.
v = u([2:end 1]); % periodic b.c.

%===========================================================
% Computational routines
%===========================================================
function unew = shu_osher_step(u,la)
% Shu-Osher SSP scheme with la = dt/h
u1 = euler_step(u,la);
u2 = 3/4*u+1/4*euler_step(u1,la);
unew = 1/3*u+2/3*euler_step(u2,la);

function unew = euler_step(u,la)
% Euler step with la = dt/h
unew = u+la*rhs(u);

function du = rhs(u)
% Right hand side, multiplied by h
[vL,vR] = weno_twoside(u);
du = -(f_num(vR,shr(vL))-f_num(shl(vR),vL));

function [vL,vR] = weno_twoside(u)
% WENO approximation to values at cell boundaries
uL = shl(u); uLL = shl(uL);
uR = shr(u); uRR = shr(uR);
vR = weno(uLL,uL,u,uR,uRR); % right cell boundary
vL = weno(uRR,uR,u,uL,uLL); % left cell boundary

function v = weno(uLL,uL,u,uR,uRR)
% 5th order WENO approximation to value at right cell
% boundary of central cell, computed from cell averages
% u, shifted by -2 to 2.
%
% smoothness indicators
betaL = 13/12*(uLL-2*uL+u  ).^2+1/4*(  uLL-4*uL+3*u  ).^2;
beta  = 13/12*(uL -2*u +uR ).^2+1/4*(  uL      -  uR ).^2;
betaR = 13/12*(u  -2*uR+uRR).^2+1/4*(3*u  -4*uR+  uRR).^2;
% regularization
epsilon = 1e-8*max(u,max(max(uL.^2,uLL.^2),...
               max(uR.^2,uRR.^2)))+1e-99;
% weights
wL = .1./(epsilon+betaL);
w  = .6./(epsilon+beta );
wR = .3./(epsilon+betaR);
ws = wL+w+wR;
wL = wL./ws; wR = wR./ws; w = 1-wL-wR;
% quadratic reconstructions
pL = (2*uLL-7*uL+11*u  )/6;
p  = (  -uL+5*u + 2*uR )/6;
pR = ( 2*u +5*uR-   uRR)/6;
% final estimate between u and uR
v = wL.*pL+w.*p+wR.*pR;

function z = f_num(a,b)
% Numerical flux function
% Here: local Lax-Friedrichs flux
c = max(abs(f1(a)),abs(f1(b))); % assuming f convex
z = (f(a)+f(b)-c.*(b-a))/2; % central plus diffusion

function z = f1(u)
% Derivative of flux function
e = 1e-7;
z = (f(u+e)-f(u-e))/(2*e);

