%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright G. Allaire, Novembre 2000
%
% schemas de Godunov et de Van Leer  
% pour l'equation de Burgers
% u,t + (u^2/2),x = 0
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% lg = longueur du domaine
lg = 5. ;
% nx = nombre de mailles
nx = 50 ;
% dx = pas d'espace
dx = lg/nx ;
% etat gauche et droit
ul = 0.8 ;
ur = 0.  ;
liminf = min(ul,ur) - 0.2 ;
limmax = max(ul,ur) + 0.2 ;
% condition cfl
cfl = dx/max(abs(ul),abs(ur)) ;
% dt = pas de temps
dt = 0.9*cfl ;
% tt = temps de simulation
if ul > ur 
   tt=0.8*lg/abs(ul+ur) ;
else 
   tt=0.4*lg/max(abs(ul),abs(ur)) ;
end
nt = round(tt/dt) ;
%
% initialisation
%
% x = points du maillage
x=zeros(1,nx+1) ;
% ug = solution de Godunov
ug=zeros(1,nx+1) ;
% uv = solution de Van Leer
uv=zeros(1,nx+1) ;
% ue = solution exacte
ue=zeros(1,nx+1) ;
% f = flux
f=zeros(1,nx+1) ;
% p = pente
p=zeros(1,nx+1) ;
% du = increment u(i+1)-u(i)
du=zeros(1,nx+1) ;
% tableaux auxiliaires 
v=zeros(1,nx+1) ;
w=zeros(1,nx+1) ;
up=zeros(1,nx+1) ;
um=zeros(1,nx+1) ;
uvp=zeros(1,nx+1) ;
uvm=zeros(1,nx+1) ;
uvmp=zeros(1,nx+1) ;
fm=zeros(1,nx+1) ;
dum=zeros(1,nx+1) ;
sdu=zeros(1,nx+1) ;
sdum=zeros(1,nx+1) ;
% donnee initiale
for i=1:nx+1
  x(i) = (i-1)*dx  ;
  if x(i) < lg/2
     ue(i) = ul ;
  else
     ue(i) = ur ;
  end
end
ug = ue ;
uv = ue ;

plot(x,ue) ;
axis([0 lg liminf limmax]) ;
title('comparaison Godunov - Van Leer, cfl=0.9')

pause ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% marche en temps : 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
for n=1:nt
%
% Godunov
% evaluation du flux a l'interface
f = 0.5*ug.^2 ;
for i=1:nx
  if ug(i) > ug(i+1)
% choc
     sigma = 0.5*(ug(i)+ug(i+1)) ;
     if sigma < 0.
        f(i) = 0.5*ug(i+1)^2 ;
     end
  elseif ug(i) < ug(i+1)
% detente
     if ug(i+1) < 0.
        f(i) = 0.5*ug(i+1)^2 ;
     elseif ug(i) > 0.
     else 
        f(i) = 0. ;
     end
  end
end
fm(2:nx+1) = f(1:nx) ;
v(2:nx) = ug(2:nx) - (dt/dx)*(f(2:nx)-fm(2:nx)) ;
v(1) = ug(1) ;
v(nx+1) = ug(nx+1) ;
ug = v ; 


% Van Leer
up(1:nx) = uv(2:nx+1) ;
um(2:nx+1) = uv(1:nx) ;
du = up-uv ;
dum(2:nx+1) = du(1:nx) ;
sdu = sign(du) ;
sdum = sign(dum) ;
% evaluation de la pente par la fonction minmod
p = max(0,sdu.*sdum).*sdu.*min(abs(du),abs(dum)) ;
% uvm = solution en (i+1/2)-
% uvp = solution en (i-1/2)+
uvm = uv + 0.5*p ;
uvp = uv - 0.5*p ;
% evaluation de la solution au demi pas de temps
v(2:nx) = uvm(2:nx) - (0.5*dt/dx)*(uvm(2:nx).^2-uvp(2:nx).^2)/2 ;
w(2:nx) = uvp(2:nx) - (0.5*dt/dx)*(uvm(2:nx).^2-uvp(2:nx).^2)/2 ;
v(1) = uv(1) ;
v(nx+1) = uv(nx+1) ;
w(1) = uv(1) ;
w(nx+1) = uv(nx+1) ;
uvm = v ;
uvp = w ;
uvpp = wshift3('1D',uvp,+1) ;
uvpp(1:nx) = uvp(2:nx+1) ;
% resolution du probleme de Riemann entre uvm et uvpp
% evaluation du flux a l'interface
f = 0.5*uvm.^2 ;
for i=1:nx
  if uvm(i) > uvpp(i)
% choc
     sigma = 0.5*(uvm(i)+uvpp(i)) ;
     if sigma < 0.
        f(i) = 0.5*uvpp(i)^2 ;
     end
  elseif uvm(i) < uvpp(i)
% detente
     if uvpp(i) <= 0.
        f(i) = 0.5*uvpp(i)^2 ;
     elseif uvm(i) >= 0.
     else 
        f(i) = 0. ;
     end
  end
end
fm(2:nx+1) = f(1:nx) ;
v(2:nx) = uv(2:nx) - (dt/dx)*(f(2:nx)-fm(2:nx)) ;
v(1) = uv(1) ;
v(nx+1) = uv(nx+1) ;
uv = v ; 


if mod(n,1) == 0
ue=solbu(ue,nx,dx,dt,lg,ul,ur,n) ;
plot(x,ug,'o',x,uv,'*',x,ue) ;
axis([0 lg liminf limmax]) ;
title('comparaison Godunov - Van Leer, cfl=0.9')
legend('Godunov','Van Leer','exact') ;
pause(0) ;
end
%
end

pause ;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% conditions aux limites periodiques
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
lg = 1. ;
% nx = nombre de mailles
nx = 50 ;
% dx = pas d'espace
dx = lg/nx ;
% dt = pas de temps
% (on sait que la vitesse maximale sera 2.)
cfl = dx/2 ;
dt = 0.9*cfl ;
% nt = nombre de pas de temps effectues
nt = 300 ;
%
% initialisation
%
% x = points du maillage
x=zeros(1,nx) ;
% ug = solution de Godunov
ug=zeros(1,nx) ;
% uv = solution de Van Leer
uv=zeros(1,nx) ;
% f = flux
f=zeros(1,nx) ;
% p = pente
p=zeros(1,nx) ;
% du = increment u(i+1)-u(i)
du=zeros(1,nx) ;
% tableaux auxiliaires 
v=zeros(1,nx) ;
w=zeros(1,nx) ;
up=zeros(1,nx) ;
um=zeros(1,nx) ;
uvp=zeros(1,nx) ;
uvm=zeros(1,nx) ;
uvpp=zeros(1,nx) ;
fm=zeros(1,nx) ;
dum=zeros(1,nx) ;
sdu=zeros(1,nx) ;
sdum=zeros(1,nx) ;
% donnee initiale
for i=1:nx
  x(i) = (i-1)*dx  ;
  ug(i) = 1. + sin(x(i)*2*pi) ;
end
uv = ug ;

% affichage de la donnee initiale
plot(x,ug,'*') ;
axis([0 lg -0.1 2.1]) ;
title('comparaison Godunov - Van Leer, cfl=0.9')

pause(1) ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% marche en temps : 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
for n=1:nt
%
% Godunov
% evaluation du flux a l'interface
f = 0.5*ug.^2 ;
up = wshift2('1D',ug,+1) ;
for i=1:nx
  if ug(i) > up(i)
% choc
     sigma = 0.5*(ug(i)+up(i)) ;
     if sigma < 0.
        f(i) = 0.5*up(i)^2 ;
     end
  elseif ug(i) < up(i)
% detente
     if up(i) <= 0.
        f(i) = 0.5*up(i)^2 ;
     elseif ug(i) >= 0.
     else 
        f(i) = 0. ;
     end
  end
end
fm = wshift2('1D',f,-1) ;
v = ug - (dt/dx)*(f-fm) ;
ug = v ; 

% Van Leer
um = wshift2('1D',uv,-1) ;
up = wshift2('1D',uv,+1) ;
du = up-uv ;
dum = wshift2('1D',du,-1) ;
sdu = sign(du) ;
sdum = sign(dum) ;
% evaluation de la pente par la fonction minmod
p = max(0,sdu.*sdum).*sdu.*min(abs(du),abs(dum)) ;
% uvm = solution en (i+1/2)-
% uvp = solution en (i-1/2)+
uvm = uv + 0.5*p ;
uvp = uv - 0.5*p ;
% evaluation de la solution au demi pas de temps
v = uvm - (0.5*dt/dx)*(uvm.^2-uvp.^2)/2 ;
w = uvp - (0.5*dt/dx)*(uvm.^2-uvp.^2)/2 ;
uvm = v ;
uvp = w ;
uvpp = wshift2('1D',uvp,+1) ;
% resolution du probleme de Riemann entre uvm et uvpp
% evaluation du flux a l'interface
f = 0.5*uvm.^2 ;
for i=1:nx
  if uvm(i) > uvpp(i)
% choc
     sigma = 0.5*(uvm(i)+uvpp(i)) ;
     if sigma < 0.
        f(i) = 0.5*uvpp(i)^2 ;
     end
  elseif uvm(i) < uvpp(i)
% detente
     if uvpp(i) <= 0.
        f(i) = 0.5*uvpp(i)^2 ;
     elseif uvm(i) >= 0.
     else 
        f(i) = 0. ;
     end
  end
end
fm = wshift2('1D',f,-1) ;
v = uv - (dt/dx)*(f-fm) ;
uv = v ; 

if mod(n,1) == 0
plot(x,ug,'o',x,uv,'*') ;
axis([0 lg -0.1 2.1]) ;
title('comparaison Godunov - Van Leer, cfl=0.9')
legend('Godunov','Van Leer') ;
pause(0) ;
end
%
end


