clear

for N = 1:4
    for k = 1:7
        K = 2.^k;
        L2err(k) = advec1D(N,K);
        h(k) = 2/K;
    end
    loglog(h,L2err,'o--','linewidth',2)
    hold on
    loglog(h,(L2err(1)/h(1)^(N+1))*h.^(N+1),'k--') % line with slope (N+1)
end
xlabel('h','fontsize',15)
ylabel('L2 error','fontsize',15)


function L2err = advec1D(Nin,K)

% clear

% demo for 1D advection
global Dr LIFT N

if nargin==0
    N = 4; % order of polynomial
    K = 16; % # elements
else
    N = Nin;
end
FinalTime = 1.1;

% set wavespeed
global a
a = 1;

% make grid
VX = linspace(-1,1,K+1);

% make nodal operators
Np = N+1;
r = JacobiGL(0,0,N); % interpolation points
V = Vandermonde1D(N,r);
Vr = GradVandermonde1D(N,r);
Dr = Vr/V;
M = inv(V')*inv(V); % inv(V*V');
Mf = zeros(N+1,2); Mf(1,1) = 1; Mf(N+1,2) = 1;
LIFT = M\Mf;

% make nodes: x(r) = xk + .5*(1+r)*h;
x = zeros(Np,K);
for e = 1:K
    he = VX(e+1)-VX(e);
    x(:,e) = VX(e) + .5*he*(1+r);
end

% make node maps
global mapM mapP
Nfp = 1; Nfaces = 2;
mapM = reshape(1:Nfp*Nfaces*K,Nfp*Nfaces,K);
mapP = mapM;
for e = 2:K-1
    mapP(1,e) = mapM(2,e-1);
    mapP(2,e) = mapM(1,e+1);
end
mapP(2,1) = mapM(1,2);
mapP(1,K) = mapM(2,K-1);

% periodic mapping
mapP(1) = mapM(end);
mapP(end) = mapM(1);

% compute geometric factors: drdx
% x(r) = xk + .5*(1+r)*h, dx/dr = .5*h, dr/dx*dx/dr = 1, dr/dx = 2/h
global rx J nx Fscale
xr = Dr*x; J = xr; rx = 1./xr;
nx = zeros(2,K);
nx(1,:) = -1; nx(2,:) = 1;
Fscale = 1./J([1 N+1],:);

% set initial condition
u = zeros(Np,K);
x0 = -.5;
u0 = @(x) exp(-100*(x-x0).^2);
u0 = @(x) sin(pi*x);
u = u0(x);
% u = abs(x+.25) < .5;
% plot(x,u,'o')
% return

%% timestepping

% Low storage Runge-Kutta coefficients
rk4a = [            0.0 ...
    -567301805773.0/1357537059087.0 ...
    -2404267990393.0/2016746695238.0 ...
    -3550918686646.0/2091501179385.0  ...
    -1275806237668.0/842570457699.0];
rk4b = [ 1432997174477.0/9575080441755.0 ...
    5161836677717.0/13612068292357.0 ...
    1720146321549.0/2090206949498.0  ...
    3134564353537.0/4481467310338.0  ...
    2277821191437.0/14882151754819.0];
rk4c = [             0.0  ...
    1432997174477.0/9575080441755.0 ...
    2526269341429.0/6820363962896.0 ...
    2006345519317.0/3224310063776.0 ...
    2802321613138.0/2924317926251.0];

% timestep size
dxmin = min(abs(x(1,:)-x(2,:)));
dt = .125*dxmin / a;
Nsteps = ceil(FinalTime/dt);
dt = FinalTime/Nsteps;

resu = zeros(Np,K);
% figure(1)
for i = 1:Nsteps
    
    time = i*dt;
    for INTRK = 1:5
        timelocal = time + rk4c(INTRK)*dt;
        [rhsu] = advecRHS(u,timelocal);
        resu = rk4a(INTRK)*resu + dt*rhsu;
        u = u + rk4b(INTRK)*resu;
    end
    
    % sum_k (J*u'*M*u) = ||u||^2_[-1,1]
    unorm(i) = sum(sum((u.*J).*(M*u)));
    
    if nargin==0 && mod(i,10)==0
        plot(x,u,'o--')
        title(sprintf('time = %f\n',i*dt))
        axis([-1,1, -2 2])
        drawnow
    end
end

% interpolate to quadrature
[rq wq] = JacobiGQ(0,0,N+1);
Vq = Vandermonde1D(N,rq)/V;
xq = Vq*x;
wJq = diag(wq)*(Vq*J);

uex = u0(xq-FinalTime);

L2err = sqrt(sum(sum(wJq.*(Vq*u - uex).^2)));

end

function rhs = advecRHS(u,time)

global Dr LIFT N
global rx Fscale nx
global mapM mapP
global a

uf = u([1 N+1],:);
uP = uf(mapP);
% uP(1) = 0*exp(-100*(time-.25).^2); % no boundary condition for periodic

% compute flux term: (avg(au)*n) - jump(u)
du = .5*(uP+uf(mapM)) - uf(mapM);
ujump = uP-uf(mapM);
flux = a*du.*nx - .5*abs(a)*ujump; % multiply flux by nx

rhs = -(a*rx.*(Dr*u) + LIFT*(Fscale.*flux));

end