%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright G. Allaire, Octobre 2000
%
% schemas numeriques pour le probleme de Riemann avec l'equation de Burgers
% u,t + (u^2/2),x = 0
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% lg = longueur du domaine
lg = 1. ;
% nx = nombre de mailles
nx = 100 ;
% dx = pas d'espace
dx = lg/nx ;
% etat gauche et droit
ul = 0.2 ;
ur = 0.8 ;
liminf = min(ul,ur) - 0.1 ;
limmax = max(ul,ur) + 0.1 ;
% 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=1./abs(ul+ur) ;
else 
   tt=0.5/max(abs(ul),abs(ur)) ;
end
tt=0.9*tt ;
nt = round(tt/dt) ;
%
% initialisation
%
x=zeros(1,nx) ;
u=zeros(1,nx) ;
v=zeros(1,nx) ;
w=zeros(1,nx) ;
up=zeros(1,nx) ;
um=zeros(1,nx) ;
% x = points du maillage
% u = donnee initiale
for i=1:nx
  x(i) = (i-1)*dx  ;
  if x(i) < lg/2
     u(i) = ul ;
  else
     u(i) = ur ;
  end
end
w = u ;

plot(x,u) ;
axis([0 lg liminf limmax]) ;
title('Lax-Friedrichs, cfl=0.9')

pause ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% marche en temps : Lax-Friedrichs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
for n=1:nt
%
up = wshift3('1D',u,+1) ;
um = wshift3('1D',u,-1) ;
v  = (up+um)/2 - (0.5*dt/dx)*(up.^2-um.^2)/2 ;
u = v ; 
if mod(n,5) == 0
w=solbu(w,nx,dx,dt,lg,ul,ur,n) ;
plot(x,u,'*',x,w) ;
axis([0 lg liminf limmax]) ;
title('Lax-Friedrichs, cfl=0.9')
pause(0) ;
end
%
end
pause ;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% nouveau schema: Lax-Friedrichs avec cfl=0.3
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
dt = 0.3*cfl ;
nt = round(tt/dt) ;
%
% initialisation
%
x=zeros(1,nx) ;
u=zeros(1,nx) ;
v=zeros(1,nx) ;
w=zeros(1,nx) ;
up=zeros(1,nx) ;
um=zeros(1,nx) ;
for i=1:nx
  x(i) = (i-1)*dx  ;
  if x(i) < lg/2
     u(i) = ul ;
  else
     u(i) = ur ;
  end
end
w = u ;

plot(x,u) ;
axis([0 lg liminf limmax]) ;
title('Lax-Friedrichs, cfl=0.3')

pause(1) ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% marche en temps : Lax-Friedrichs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
for n=1:nt
%
up = wshift3('1D',u,+1) ;
um = wshift3('1D',u,-1) ;
v  = (up+um)/2 - (0.5*dt/dx)*(up.^2-um.^2)/2 ;
u = v ; 
if mod(n,5) == 0
w=solbu(w,nx,dx,dt,lg,ul,ur,n) ;
plot(x,u,'*',x,w) ;
axis([0 lg liminf limmax]) ;
title('Lax-Friedrichs, cfl=0.3')
pause(0) ;
end
%
end
pause ;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% nouveau schema: centre explicite avec cfl=0.1 (explose)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
dt = 0.1*cfl ;
nt = round(tt/dt) ;
%
% initialisation
%
x=zeros(1,nx) ;
u=zeros(1,nx) ;
v=zeros(1,nx) ;
w=zeros(1,nx) ;
for i=1:nx
  x(i) = (i-1)*dx  ;
  if x(i) < lg/2
     u(i) = ul ;
  else
     u(i) = ur ;
  end
end
w = u ;

plot(x,u) ;
axis([0 lg liminf limmax+0.2]) ;
title('centre explicite, cfl=0.1')

pause(1) ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% marche en temps : centre explicite
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
for n=1:nt
%
up = wshift3('1D',u,+1) ;
um = wshift3('1D',u,-1) ;
v  = u - (0.5*dt/dx)*(up.^2-um.^2)/2 ;
u = v ; 
if mod(n,5) == 0
w=solbu(w,nx,dx,dt,lg,ul,ur,n) ;
plot(x,u,'*',x,w,x,u,':') ;
axis([0 lg liminf limmax+0.2]) ;
title('centre explicite, cfl=0.1')
pause(0) ;
end
%
end
pause ;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% nouveau schema: Lax-Wendroff avec cfl=0.9
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
dt = 0.9*cfl ;
nt = round(tt/dt) ;
%
% initialisation
%
x=zeros(1,nx) ;
u=zeros(1,nx) ;
v=zeros(1,nx) ;
w=zeros(1,nx) ;
for i=1:nx
  x(i) = (i-1)*dx  ;
  if x(i) < lg/2
     u(i) = ul ;
  else
     u(i) = ur ;
  end
end
w = u ;

plot(x,u) ;
axis([0 lg liminf limmax+0.2]) ;
title('Lax-Wendroff, cfl=0.9')

pause(1) ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% marche en temps : Lax-Wendroff
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
for n=1:nt
%
up = wshift3('1D',u,+1) ;
um = wshift3('1D',u,-1) ;
v  = u - (0.5*dt/dx)*(up.^2-um.^2)/2 + ... 
(0.5*dt*dt/(dx*dx))*( (u+up)*0.5.*(up.^2-u.^2) ...
        - (u+um)*0.5.*(u.^2-um.^2) )/2 ;
u = v ; 
if mod(n,5) == 0
w=solbu(w,nx,dx,dt,lg,ul,ur,n) ;
plot(x,u,'*',x,w,x,u,':') ;
axis([0 lg liminf limmax+0.2]) ;
title('Lax-Wendroff, cfl=0.9')
pause(0) ;
end
%
end
pause ;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% nouveau schema: Lax-Wendroff avec cfl=0.3
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
dt = 0.3*cfl ;
nt = round(tt/dt) ;
%
% initialisation
%
x=zeros(1,nx) ;
u=zeros(1,nx) ;
v=zeros(1,nx) ;
w=zeros(1,nx) ;
for i=1:nx
  x(i) = (i-1)*dx  ;
  if x(i) < lg/2
     u(i) = ul ;
  else
     u(i) = ur ;
  end
end
w = u ;

plot(x,u) ;
axis([0 lg liminf limmax+0.2]) ;
title('Lax-Wendroff, cfl=0.3')

pause(1) ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% marche en temps : Lax-Wendroff
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
for n=1:nt
%
up = wshift3('1D',u,+1) ;
um = wshift3('1D',u,-1) ;
v  = u - (0.5*dt/dx)*(up.^2-um.^2)/2 + ... 
(0.5*dt*dt/(dx*dx))*( 0.5*(u+up).*(up.^2-u.^2) ...
        - 0.5*(u+um).*(u.^2-um.^2) )/2 ;
u = v ; 
if mod(n,5) == 0
w=solbu(w,nx,dx,dt,lg,ul,ur,n) ;
plot(x,u,'*',x,w,x,u,':') ;
axis([0 lg liminf limmax+0.2]) ;
title('Lax-Wendroff, cfl=0.3')
pause(0) ;
end
%
end
pause ;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% nouveau schema: non conservatif avec cfl=0.9 (fausse solution)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
dt = 0.9*cfl ;
nt = round(tt/dt) ;
%
% initialisation
%
x=zeros(1,nx) ;
u=zeros(1,nx) ;
v=zeros(1,nx) ;
w=zeros(1,nx) ;
for i=1:nx
  x(i) = (i-1)*dx  ;
  if x(i) < lg/2
     u(i) = ul ;
  else
     u(i) = ur ;
  end
end
w = u ;

plot(x,u) ;
axis([0 lg liminf limmax+0.2]) ;
title('non conservatif, cfl=0.9')

pause(1) ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% marche en temps : non conservatif
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
for n=1:nt
%
up = wshift3('1D',u,+1) ;
um = wshift3('1D',u,-1) ;
v  = (up+um)/2 - (0.5*dt/dx)*um.*(up-um) ;
u = v ; 
if mod(n,1) == 0
w=solbu(w,nx,dx,dt,lg,ul,ur,n) ;
plot(x,u,'*',x,w,x,u,':') ;
axis([0 lg liminf limmax]) ;
title('non conservatif, cfl=0.9')
pause(0) ;
end
%
end

pause ;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% nouveau schema: Godunov avec cfl=0.9 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
dt = 0.9*cfl ;
nt = round(tt/dt) ;
%
% initialisation
%
x=zeros(1,nx) ;
u=zeros(1,nx) ;
v=zeros(1,nx) ;
w=zeros(1,nx) ;
f=zeros(1,nx) ;
sigma=zeros(1,nx) ;
for i=1:nx
  x(i) = (i-1)*dx  ;
  if x(i) < lg/2
     u(i) = ul ;
  else
     u(i) = ur ;
  end
end
w = u ;

plot(x,u) ;
axis([0 lg liminf limmax+0.2]) ;
title('Godunov, cfl=0.9')

pause(1) ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% marche en temps : Godunov
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
for n=1:nt
%
% evaluation du flux a l'interface
f = 0.5*u.^2 ;
up = wshift3('1D',u,+1) ;
for i=1:nx
if u(i) > up(i)
% choc
     sigma = 0.5*(u(i)+up(i)) ;
     if sigma < 0.
        f(i) = 0.5*up(i)^2 ;
     end
elseif u(i) < up(i)
% detente
     if up(i) < 0.
        f(i) = 0.5*up(i)^2 ;
     elseif u(i) > 0.
        f(i) = 0.5*u(i)^2 ;
     else 
        f(i) = 0. ;
     end
end
end
fm = wshift3('1D',f,-1) ;
v = u - (dt/dx)*(f-fm) ;
u = v ; 
if mod(n,5) == 0
w=solbu(w,nx,dx,dt,lg,ul,ur,n) ;
plot(x,u,'*',x,w,x,u,':') ;
axis([0 lg liminf limmax]) ;
title('Godunov, cfl=0.9')
pause(0) ;
end
%
end









