clear

% demo for 1D advection
global Dr LIFT N

N = 4; % order of polynomial
K = 16; % # elements 

FinalTime = 5;

% 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);
p = sin(pi*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;
Nsteps = ceil(FinalTime/dt);
dt = FinalTime/Nsteps;

resp = zeros(Np,K);
resu = zeros(Np,K);
figure(1)
for i = 1:Nsteps
    
    time = i*dt;
    for INTRK = 1:5
        timelocal = time + rk4c(INTRK)*dt;
        [rhsp rhsu] = waveRHS(p,u,timelocal);
        resp = rk4a(INTRK)*resp + dt*rhsp;
        resu = rk4a(INTRK)*resu + dt*rhsu;
        p = p + rk4b(INTRK)*resp;
        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 mod(i,10)==0        
        plot(x,p,'o--')
        title(sprintf('time = %f\n',i*dt))
        axis([-1,1, -2 2])
        drawnow
    end
end

% clf
% uexact = u0(x-a*FinalTime);
% plot(x,u,'o-')
% hold on
% plot(x,uexact,'--')
% return

% figure(2)
% plot(dt*(1:Nsteps),unorm,'--')
% hold on
% ylim([mean(unorm)-2 mean(unorm)+2])
% title('Norm of u over time','fontsize',15)

function [rhsp rhsu] = waveRHS(p,u,time)

global Dr LIFT N
global rx Fscale nx
global mapM mapP

tau = 1;

pf = p([1 N+1],:);
uf = u([1 N+1],:);

pP = pf(mapP);
uP = uf(mapP);
 
% compute flux term: (avg(au)*n) - jump(u)
ujump = uP-uf(mapM);
pjump = pP-pf(mapM);
pflux = .5*(ujump).*nx - .5*tau*pjump;
uflux = .5*(pjump).*nx - .5*tau*ujump;

rhsp = -(rx.*(Dr*u) + LIFT*(Fscale.*pflux));
rhsu = -(rx.*(Dr*p) + LIFT*(Fscale.*uflux));

end