% This Matlab script solves the one-dimensional convection
% equation using a finite volume algorithm. The
% discretization forward Euler in time.

clear all;
close all;
clc;
clf;

global a;


a=1;
Nx = 128;
L=2;
x = linspace(0,L,Nx);

CFL = 0.5;

s.dx = L/(Nx-1);
s.dt = CFL*s.dx/abs(a);

cr=1:(Nx-1);
cl=circshift(cr,1);
cll=circshift(cr,2);
crr=circshift(cr,-1);

s.faces={};
for i=1:Nx
   s.faces{i}.Ul=0;
   s.faces{i}.Ur=0;
   s.faces{i}.l=1;
   s.faces{i}.n=[1;0;0];
   s.faces{i}.nl=i-1;
   s.faces{i}.nr=i;
   if i==Nx
       s.faces{i}.cll=cll(1);
       s.faces{i}.cl=cl(1);
       s.faces{i}.cr=cr(1);
       s.faces{i}.crr=crr(1);   
   else
       s.faces{i}.cll=cll(i);
       s.faces{i}.cl=cl(i);
       s.faces{i}.cr=cr(i);
       s.faces{i}.crr=crr(i);
   end
   if i==Nx
      s.faces{i}.nr=-1; 
   end
   if i==1
      s.faces{i}.nl=-1; 
   end
end

s.cells={};
for i=1:(Nx-1)
    s.cells{i}.Q=0;
    s.cells{i}.V=s.dx;
    s.cells{i}.x=x(i)+s.dx/2;
    xc(i)=s.cells{i}.x;
end
    

Q=zeros(Nx-1,1);

Q(x>=0&x<=0.25)=1;

Qu=Q;
Qglf=Q;
Qllf=Q;
Qlw=Q;
Qroe=Q;
Qhoroe=Q;


figure(1);
hold on;
p1=plot(xc,Qglf,'-','DisplayName','Upwind');
p2=plot(xc,Qglf,'-','DisplayName','Global Lax-Friedrich (GLF)');
p3=plot(xc,Qllf,'-','DisplayName','Local Lax-Friedrich (LLF)');
p4=plot(xc,Qlw,'-','DisplayName','Lax-Wendroff (LW)');
p5=plot(xc,Qroe,'-','DisplayName','Roe');
p6=plot(xc,Qhoroe,'-','DisplayName','High-order Roe');
hold off;
ylim([-0.25,1.25]);
xlim([0,L]);
grid on;
legend show;
    
t = 0;
tend=0.7;
while(t<tend)
    
   
    
    s.method=0;
    s.high_order=0;
    R=residual(s,Qu);
    Qu=Qu+s.dt*R;
    
    s.method=1;
    s.high_order=0;
    R=residual(s,Qglf);
    Qglf=Qglf+s.dt*R;
    
    s.method=2;
    R=residual(s,Qllf);
    Qllf=Qllf+s.dt*R;
    
    s.method=3;
    s.high_order=0;
    R=residual(s,Qlw);
    Qlw=Qlw+s.dt*R;
    
    s.method=4;
    s.high_order=0;
    R=residual(s,Qroe);
    Qroe=Qroe+s.dt*R;
    
    s.method=4;
    s.high_order=1;
    s.limiter='minmod';
    R=residual(s,Qhoroe);
    Qhoroe=Qhoroe+s.dt*R;
    
    t=t+s.dt;
    
    set(p1,'Ydata',Qu);
    set(p2,'Ydata',Qglf);
    set(p3,'Ydata',Qllf);
    set(p4,'Ydata',Qlw);
    set(p5,'Ydata',Qroe);
    set(p6,'Ydata',Qhoroe);
    drawnow;
end



% Residual function
%%%%%%%%%%%%%%%%%%%%%%

function R=residual(s,Q)
    
    R=zeros(size(Q));
    
    for i=1:length(s.cells)
        s.cells{i}.Q=Q(i);
        s.cells{i}.R=0;
    end
    
    for i=1:length(s.faces)
        cl=s.faces{i}.cl;
        cr=s.faces{i}.cr;
        Ql=Q(cl);
        Qr=Q(cr);
        s.faces{i}.Ql=Ql;
        s.faces{i}.Qr=Qr;
    end
    
    %Reconstruction
    %%%%%%%%%%%%%%%%%%%%%%%%%%
    
    if s.high_order==1
        for i=1:length(s.faces)
            
            cll=s.faces{i}.cll;
            cl=s.faces{i}.cl;
            cr=s.faces{i}.cr;
            crr=s.faces{i}.cll;
            
            Qll=s.cells{cll}.Q; %i-2
            Ql=s.cells{cl}.Q; %i-1
            Qr=s.cells{cr}.Q; %i
            Qrr=s.cells{crr}.Q; %i+1
            
            rl=(Qr-Ql)/(Ql-Qll);
            rr=(Ql-Qr)/(Qr-Qrr);
            
            phil=limiter(rl,s.limiter);
            phir=limiter(rr,s.limiter);
            
            grad_Ql=(Ql-Qll)/s.dx;
            grad_Qr=(Qr-Qrr)/s.dx;
            
            s.faces{i}.Ql=Ql+phil*grad_Ql*s.dx/2;
            s.faces{i}.Qr=Qr-phir*grad_Qr*s.dx/2;
         
        end
    end
    
    %Fluxes
    %%%%%%%%%%%%%%%%%%%%%%%%%%
    for i=1:length(s.faces)
        
        nl=s.faces{i}.nl;
        nr=s.faces{i}.nr;
        
        Ql=s.faces{i}.Ql;
        Qr=s.faces{i}.Qr;
        
        n=s.faces{i}.n;
        l=s.faces{i}.l;
        
        F=flux(Ql,Qr,n,s.method,s.dx,s.dt);
        
        if nl~=-1
            s.cells{nl}.R=s.cells{nl}.R+F*l;
        end
        
        if nr~=-1
            s.cells{nr}.R=s.cells{nr}.R-F*l;
        end
        
    end
    
    for i=1:length(s.cells)
       R(i)=-s.cells{i}.R/s.cells{i}.V; 
    end
    
end

% Flux function
%%%%%%%%%%%%%%%%%%%%%%

function F=flux(Ql,Qr,n,method,dx,dt)
     
    global a;

    nx=n(1);
    
    delta=Qr-Ql; %Jump in conservarive variables
    
    ul=Ql; %Left u velocity
    ur=Qr; %Right u velocity
    
    Fxl=a*ul; %Left x-flux
    Fxr=a*ur; %Right x-flux
    
    Fl=Fxl*nx; %Left flux 
    Fr=Fxr*nx; %Right flux 
    
    %Upwind
    if method==0
       ubar=0.5*(ul+ur);
       if ubar>0
           F=Fl;
       else
           F=Fr;
       end
    end
    
    %Global Lax-Friedrichs 
    if method==1
       w=dx/dt*delta;
       F=0.5*(Fl+Fr)-0.5*w;
    end
    
    %Local Lax-Friedrichs
    if method==2
        %w=max(|F'(ul)|,|F'(ur)|)
        wl=abs(a*nx);
        wr=abs(a*nx);
        w=max(wl,wr)*delta;
        F=0.5*(Fl+Fr)-0.5*w;
    end
    
    %Lax-Wendroff 
    if method==3
        A=a;
        %w=dt/dx*A*(Fr-Fl)=dt/dx*A^2*delta
        w=dt/dx*A^2*delta;
        F=0.5*(Fl+Fr)-0.5*w;
    end
    
    %Roe
    if method==4
        A=a; %Jacobian A is a since f'(u)=a. Remember that f(u)=au
        w=abs(A)*delta;
        F=0.5*(Fl+Fr)-0.5*w;
    end
    
end

% Limiter function
%%%%%%%%%%%%%%%%%%%%%%

function phi=limiter(r,type)
    switch type
        case 'minmod'
            phi=max([0, min([r,1])]);
        case 'monotonizedcentral'
            phi=max([0, min([2*r,0.5*(r+1),2])]);
        case 'osher'
            phi=max([0,min(r,1.5)]);
        case 'sweby'
            phi=max([0, min([1.5*r,1]),min([r,1.5])]);
        case 'superbee'
            phi=max([0,min(2*r,1),min(r,2)]);
        case 'vanalbada1'
            phi=(r^2+r)/(r^2+1);
        case 'vanalbada2'
            phi=2*r/(r^2+1);
        case 'vanleer'
            phi=(r+abs(r))/(1+abs(r));
        case 'zero'
            phi=0;
        otherwise
            phi=1;
    end
end