%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright G. Allaire, Novembre 2000
%
% comparaison de schemas numeriques  
% avec l'equation de Burgers
% u,t + (u^2/2),x = 0
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% lg = longueur du domaine
lg = 5. ;
% nx = nombre de mailles
nx = 200 ;
% dx = pas d'espace
dx = lg/nx ;
% etat gauche et droit
ul = 0.1 ;
ur = 0.5 ;
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=zeros(1,nx) ;
v=zeros(1,nx) ;
f=zeros(1,nx) ;
% uf = solution de Lax-Friedrichs
uf=zeros(1,nx) ;
% uw = solution de Lax-Wendroff
uw=zeros(1,nx) ;
% ug = solution de Godunov
ug=zeros(1,nx) ;
% ue = solution exacte
ue=zeros(1,nx) ;
up=zeros(1,nx) ;
um=zeros(1,nx) ;
% x = points du maillage
% ue = donnee initiale
for i=1:nx
  x(i) = (i-1)*dx  ;
  if x(i) < lg/2
     ue(i) = ul ;
  else
     ue(i) = ur ;
  end
end
uf = ue ;
uw = ue ;
ug = ue ;

plot(x,ue) ;
axis([0 lg liminf limmax]) ;
title('LF, LW, Godunov, cfl=0.9')

pause(1) ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% marche en temps : 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
for n=1:nt
%
% Lax-Friedrichs
up = wshift3('1D',uf,+1) ;
um = wshift3('1D',uf,-1) ;
v  = (up+um)/2 - (0.5*dt/dx)*(up.^2-um.^2)/2 ;
uf = v ; 
% Lax-Wendroff
up = wshift3('1D',uw,+1) ;
um = wshift3('1D',uw,-1) ;
v  = uw - (0.5*dt/dx)*(up.^2-um.^2)/2 + ... 
(0.5*dt*dt/(dx*dx))*( (uw+up)*0.5.*(up.^2-uw.^2) ...
        - (uw+um)*0.5.*(uw.^2-um.^2) )/2 ;
uw = v ; 
% Godunov
% evaluation du flux a l'interface
f = 0.5*ug.^2 ;
up = wshift3('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.
        f(i) = 0.5*ug(i)^2 ;
     else 
        f(i) = 0. ;
     end
end
end
fm = wshift3('1D',f,-1) ;
v = ug - (dt/dx)*(f-fm) ;
ug = v ; 
if mod(n,10) == 0
ue=solbu(ue,nx,dx,dt,lg,ul,ur,n) ;
plot(x,uf,'o',x,uw,'*',x,ug,'d',x,ue) ;
axis([0 lg liminf limmax]) ;
title('comparaison de 3 schemas, cfl=0.9')
legend('Lax-Friedrichs','Lax-Wendroff','Godunov','exact') ;
pause(0) ;
end
%
end

pause ;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% conditions aux limites periodiques
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
lg = 1. ;
% nx = nombre de mailles
nx = 100 ;
% dx = pas d'espace
dx = lg/nx ;
cfl = dx ;
% dt = pas de temps
% (on sait que la vitesse maximale sera 1.)
dt = 0.9*cfl ;
% nt = nombre de pas de temps effectues
nt = 200 ;
%
% initialisation
%
x=zeros(1,nx) ;
v=zeros(1,nx) ;
f=zeros(1,nx) ;
% uf = solution de Lax-Friedrichs
uf=zeros(1,nx) ;
% uw = solution de Lax-Wendroff
uw=zeros(1,nx) ;
% ug = solution de Godunov
ug=zeros(1,nx) ;
up=zeros(1,nx) ;
um=zeros(1,nx) ;
u=zeros(1,nx) ;
v=zeros(1,nx) ;
up=zeros(1,nx) ;
um=zeros(1,nx) ;
sigma=zeros(1,nx) ;
% x = points du maillage
% u = donnee initiale
for i=1:nx
  x(i) = (i-1)*dx  ;
  uf(i) = sin(x(i)*2*pi) ;
end
uw = uf ;
ug = uf ;

% affichage de la donnee initiale
plot(x,uf,'*') ;
axis([0 lg -1.1 1.1]) ;
title('comparaison LF, LW, Godunov, cfl=0.9')
'old'
pause(1) ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% marche en temps : 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
for n=1:nt
%
% Lax-Friedrichs
up = wshift2('1D',uf,+1) ;
um = wshift2('1D',uf,-1) ;
v = (up+um)/2 - (0.5*dt/dx)*(up.^2-um.^2)/2 ;
uf = v ; 
% Lax-Wendroff
up = wshift2('1D',uw,+1) ;
um = wshift2('1D',uw,-1) ;
v = uw - (0.5*dt/dx)*(up.^2-um.^2)/2 + ... 
(0.5*dt*dt/(dx*dx))*( (uw+up)*0.5.*(up.^2-uw.^2) ...
        - (uw+um)*0.5.*(uw.^2-um.^2) )/2 ;
uw = v ; 
% Godunov
% evaluation du flux a l'interface (i+1/2)
% ug = etat gauche, up = etat droit
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.
        f(i) = 0.5*ug(i)^2 ;
     else 
        f(i) = 0. ;
     end
end
end
fm = wshift2('1D',f,-1) ;
v = ug - (dt/dx)*(f-fm) ;
ug = v ; 
if mod(n,5) == 0
plot(x,uf,'o',x,uw,'*',x,ug,'d') ;
axis([0 lg -1.1 1.1]) ;
title('comparaison de 3 schemas, cfl=0.9')
legend('Lax-Friedrichs','Lax-Wendroff','Godunov') ;
pause(0) ;
end
%
end


