%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright G. Allaire, Octobre 2000
%
% probleme de Riemann pour les equations d'Euler
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% lg = longueur du domaine
lg = 1. ;
% nx = nombre de mailles
nx = 1000 ;
% dx = pas d'espace
dx = lg/nx ;
% gam = gamma de la loi d'etat du gaz parfait
gam = 1.4 ;
% etats initiaux a gauche (l) et a droite (r)
% r = densite, u = vitesse, e = energie interne
rl = 1. ;
ul = 0. ;
el = 1/(gam-1.) ;
rr = 0.125 ;
ur = 0. ;
er = 0.8/(gam-1.) ;
% limitation du pas de temps par 
% la condition CFL
cfl = 0.9/(1. + sqrt(gam)) ;
dt = dx*cfl ;
% nt = nombre de pas de temps
nt = 400 ;
%
% initialisation
%
x=zeros(1,nx+1) ;
r=zeros(1,nx+1) ;
u=zeros(1,nx+1) ;
e=zeros(1,nx+1) ;
for i=1:nx+1
  x(i) = (i-1)*dx  ;
  if x(i) < lg/2
     r(i) = rl ;
     u(i) = ul ;
     e(i) = el ;
  else
     r(i) = rr ;
     u(i) = ur ;
     e(i) = er ;
  end
end
% p = pression
p = (gam-1.)*r.*e ;

% trace des conditions initiales
liminf = -0.1 ;
limmax = 1.1 ;
subplot(3,1,1), plot(x,r) ;
axis([0 lg 0. limmax]) ;
title('Lax-Friedrichs: densite')
subplot(3,1,2), plot(x,u) ;
axis([0 lg liminf limmax]) ;
title('Lax-Friedrichs: vitesse')
subplot(3,1,3), plot(x,p) ;
axis([0 lg 0. limmax]) ;
title('Lax-Friedrichs: pression')
pause ;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% marche en temps : schema de Lax-Friedrichs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% si u est le vecteur des valeurs u(i) alors 
% wshift3('1D',u,+1) est le vecteur des u(i+1) 
% wshift3('1D',u,-1) est le vecteur des u(i-1)
%
for n=1:nt
%
rp = wshift3('1D',r,+1) ;
rm = wshift3('1D',r,-1) ;
up = wshift3('1D',u,+1) ;
um = wshift3('1D',u,-1) ;
ep = wshift3('1D',e,+1) ;
em = wshift3('1D',e,-1) ;
pp = (gam-1.)*rp.*ep ;
pm = (gam-1.)*rm.*em ;
ip = rp.*(up.^2/2.+ep) ;
im = rm.*(um.^2/2.+em) ;
rn = (rp+rm)/2 - (0.5*dt/dx)*(rp.*up-rm.*um) ;
un = (rp.*up+rm.*um)/2 - (0.5*dt/dx)*... 
     (rp.*up.^2+pp-rm.*um.^2-pm) ;
un = un./rn ;
en = (ip+im)/2 - (0.5*dt/dx)*... 
     ((ip+pp).*up-(im+pm).*um) ;
en = en./rn - un.^2/2. ;
rn(1) = r(1) ;
rn(nx+1) = r(nx+1) ;
un(1) = u(1) ;
un(nx+1) = u(nx+1) ;
en(1) = e(1) ;
en(nx+1) = e(nx+1) ;
r = rn ;
u = un ; 
e = en ;
p = (gam-1.)*r.*e ;

if mod(n,5) == 0
   subplot(3,1,1), plot(x,r) ;
   axis([0 lg 0. limmax]) ;
   title('Lax-Friedrichs: densite')
   subplot(3,1,2), plot(x,u) ;
   axis([0 lg liminf limmax]) ;
   title('Lax-Friedrichs: vitesse')
   subplot(3,1,3), plot(x,p) ;
   axis([0 lg 0. limmax]) ;
   title('Lax-Friedrichs: pression')
   pause(0) ;
end
%
end


pause ;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% schema de Van Leer (flux splitting)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% initialisation
%
cfl = 0.6 ;
x=zeros(1,nx+1) ;
r=zeros(1,nx+1) ;
u=zeros(1,nx+1) ;
e=zeros(1,nx+1) ;
% a = vitesse du son
a=zeros(1,nx+1) ;
% (fp,fm)   = decomposition du flux en i
% (fpp,fmp) = decomposition du flux en i+1
% (fpm,fmm) = decomposition du flux en i-1
fp=zeros(3,nx+1) ;
fm=zeros(3,nx+1) ;
fpp=zeros(3,nx+1) ;
fmp=zeros(3,nx+1) ;
fpm=zeros(3,nx+1) ;
fmm=zeros(3,nx+1) ;
for i=1:nx+1
  x(i) = (i-1)*dx  ;
  if x(i) < lg/2
     r(i) = rl ;
     u(i) = ul ;
     e(i) = el ;
  else
     r(i) = rr ;
     u(i) = ur ;
     e(i) = er ;
  end
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% marche en temps : Van Leer
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
for n=1:(5*nt/4)
%
umax = max(abs(u)) ;
a = sqrt(gam*(gam-1.)*e) ; 
amax = max(abs(a)) ;
dt = dx*cfl/(umax+amax) ;
[fp,fm] = fluvl(gam,r,u,a) ;
fpp = wshift3('2D',fp,[0 +1]) ;
fpm = wshift3('2D',fp,[0 -1]) ;
fmp = wshift3('2D',fm,[0 +1]) ;
fmm = wshift3('2D',fm,[0 -1]) ;

rn = r - (dt/dx)*( (fp(1,:)+fmp(1,:))...
       -           (fpm(1,:)+fm(1,:)) )  ;

un = r.*u - (dt/dx)*( (fp(2,:)+fmp(2,:))...
          -           (fpm(2,:)+fm(2,:)) )  ;
un = un./rn ;

en = r.*(u.^2/2.+e) - (dt/dx)*( (fp(3,:)+fmp(3,:))...
                    -           (fpm(3,:)+fm(3,:)) )  ;
en = en./rn - un.^2/2. ;

% corrections des effets de bord
rn(1) = r(1) ;
rn(nx+1) = r(nx+1) ;
un(1) = u(1) ;
un(nx+1) = u(nx+1) ;
en(1) = e(1) ;
en(nx+1) = e(nx+1) ;
% mise a jour des variables
r = rn ;
u = un ; 
e = en ;
p = (gam-1.)*r.*e ;
% trace de la solution
if mod(n,5) == 0
   subplot(3,1,1), plot(x,r) ;
   axis([0 lg 0. limmax]) ;
   title('Van Leer: densite')
   subplot(3,1,2), plot(x,u) ;
   axis([0 lg liminf limmax]) ;
   title('Van Leer: vitesse')
   subplot(3,1,3), plot(x,p) ;
   axis([0 lg 0. limmax]) ;
   title('Van Leer: pression')
   pause(0) ;
end
%
end


pause ;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% schema de Richtmyer (Lax-Wendroff en deux etapes)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% initialisation
%
cfl = 0.9/(1. + sqrt(gam)) ;
dt = dx*cfl ;
x=zeros(1,nx+1) ;
r=zeros(1,nx+1) ;
u=zeros(1,nx+1) ;
e=zeros(1,nx+1) ;
% f = flux en i, fp = flux en i+1, fm = flux en i-1
f=zeros(3,nx+1) ;
fp=zeros(3,nx+1) ;
fm=zeros(3,nx+1) ;
for i=1:nx+1
  x(i) = (i-1)*dx  ;
  if x(i) < lg/2
     r(i) = rl ;
     u(i) = ul ;
     e(i) = el ;
  else
     r(i) = rr ;
     u(i) = ur ;
     e(i) = er ;
  end
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% marche en temps : Richtmyer
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
for n=1:nt
%
rp = wshift3('1D',r,+1) ;
up = wshift3('1D',u,+1) ;
ep = wshift3('1D',e,+1) ;
f = flux(gam,r,u,e) ;
fp = wshift3('2D',f,[0 +1]) ;
% etat intermediaire (n+1/2)
rs = 0.5*(rp+r) - (0.5*dt/dx)*(fp(1,:)-f(1,:)) ;

us = 0.5*(rp.*up+r.*u) - (0.5*dt/dx)*(fp(2,:)-f(2,:)) ;
us = us./rs ;

es = 0.5*(rp.*(up.^2/2.+ep)+r.*(u.^2/2.+e)) - ... 
     (0.5*dt/dx)*(fp(3,:)-f(3,:)) ;
es = es./rs - us.^2/2. ;
% etat final (n+1)
f = flux(gam,rs,us,es) ;
fm = wshift3('2D',f,[0 -1]) ;
rn = r - (dt/dx)*(f(1,:)-fm(1,:)) ;

un = r.*u - (dt/dx)*(f(2,:)-fm(2,:)) ;
un = un./rn ;

en = r.*(u.^2/2.+e) - (dt/dx)*(f(3,:)-fm(3,:)) ;
en = en./rn - un.^2/2. ;
% corrections des effets de bord
rn(1) = r(1) ;
rn(nx+1) = r(nx+1) ;
un(1) = u(1) ;
un(nx+1) = u(nx+1) ;
en(1) = e(1) ;
en(nx+1) = e(nx+1) ;
% mise a jour des variables
r = rn ;
u = un ; 
e = en ;
p = (gam-1.)*r.*e ;
% trace de la solution
if mod(n,5) == 0
   subplot(3,1,1), plot(x,r) ;
   axis([0 lg 0. limmax]) ;
   title('Richtmyer: densite')
   subplot(3,1,2), plot(x,u) ;
   axis([0 lg liminf limmax]) ;
   title('Richtmyer: vitesse')
   subplot(3,1,3), plot(x,p) ;
   axis([0 lg 0. limmax]) ;
   title('Richtmyer: pression')
   pause(0) ;
end
%
end


pause ;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% schema de Lax-Wendroff
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% initialisation
%
cc =1. ;
cfl = 0.9/(1. + sqrt(gam)) ;
dt = dx*cfl ;
x=zeros(1,nx+1) ;
r=zeros(1,nx+1) ;
u=zeros(1,nx+1) ;
e=zeros(1,nx+1) ;
% f = flux en i, fp = flux en i+1, fm = flux en i-1
f=zeros(3,nx+1) ;
fp=zeros(3,nx+1) ;
fm=zeros(3,nx+1) ;
% a = matrice jacobienne
a2=zeros(3,nx+1) ;
a3=zeros(3,nx+1) ;
for i=1:nx+1
  x(i) = (i-1)*dx  ;
  if x(i) < lg/2.01
     r(i) = rl ;
     u(i) = ul ;
     e(i) = el ;
  elseif x(i) > lg/1.99
     r(i) = rr ;
     u(i) = ur ;
     e(i) = er ;
  else
     r(i) = (rl+rr)/2 + 0.5*(rl-rr)*cos((x(i)-lg/2.01)*pi/(lg/1.99-lg/2.01)) ;
     u(i) = (ul+ur)/2 + 0.5*(ul-ur)*cos((x(i)-lg/2.01)*pi/(lg/1.99-lg/2.01)) ;
     e(i) = (el+er)/2 + 0.5*(el-er)*cos((x(i)-lg/2.01)*pi/(lg/1.99-lg/2.01)) ;
  end
end
% p = pression
p = (gam-1.)*r.*e ;

% trace des conditions initiales 
% il faut les regulariser sinon ca explose...
subplot(3,1,1), plot(x,r) ;
axis([0 lg 0. limmax]) ;
title('Lax-Wendroff: densite')
subplot(3,1,2), plot(x,u) ;
axis([0 lg liminf limmax]) ;
title('Lax-Wendroff: vitesse')
subplot(3,1,3), plot(x,p) ;
axis([0 lg 0. limmax]) ;
title('Lax-Wendroff: pression')


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% marche en temps : Lax-Wendroff
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
for n=1:nt
%
f = flux(gam,r,u,e) ;
rp = wshift3('1D',r,+1) ;
up = wshift3('1D',u,+1) ;
ep = wshift3('1D',e,+1) ;
fp = wshift3('2D',f,[0 +1]) ;
fm = wshift3('2D',f,[0 -1]) ;
% etat moyen (i+1/2)
rmoy = (r+rp)/2. ;
umoy = (r.*u+rp.*up)./(r+rp) ;
%umoy = (u+up)/2. ;
emoy = ( r.*(0.5*u.^2+e) + rp.*(0.5*up.^2+ep) )./...
       (r+rp) - 0.5*umoy.^2 ;
%emoy = (e+ep)/2. ;
% jacobien a l'etat moyen
a2(1,:) = 0.5*(gam-3.)*umoy.^2 ;
a2(2,:) = -(gam-3.)*umoy ;
a3(1,:) = (gam-1.)*umoy.^3 - gam*umoy.*(emoy+0.5*umoy.^2) ;
a3(2,:) = gam*(emoy+0.5*umoy.^2) - 1.5*(gam-1.)*umoy.^2 ;
a3(3,:) = gam*umoy ;
df = fp-f ;

gp(1,:) = df(2,:) ;
gp(2,:) = a2(1,:).*df(1,:) + a2(2,:).*df(2,:) + (gam-1)*df(3,:) ;
gp(3,:) = a3(1,:).*df(1,:) + a3(2,:).*df(2,:) + a3(3,:).*df(3,:) ;
gm = wshift3('2D',gp,[0 -1]) ;

rn = r - (0.5*dt/dx)*(fp(1,:)-fm(1,:))...
       + cc*(0.5*dt*dt/(dx*dx))*(gp(1,:)-gm(1,:)) ;

un = r.*u - (0.5*dt/dx)*(fp(2,:)-fm(2,:))... 
          + cc*(0.5*dt*dt/(dx*dx))*(gp(2,:)-gm(2,:)) ;
un = un./rn ;

en = r.*(u.^2/2.+e) - (0.5*dt/dx)*(fp(3,:)-fm(3,:))...
     + cc*(0.5*dt*dt/(dx*dx))*(gp(3,:)-gm(3,:)) ;
en = en./rn - un.^2/2. ;

% corrections des effets de bord
rn(1) = r(1) ;
rn(nx+1) = r(nx+1) ;
un(1) = u(1) ;
un(nx+1) = u(nx+1) ;
en(1) = e(1) ;
en(nx+1) = e(nx+1) ;
% mise a jour des variables
r = rn ;
u = un ; 
e = en ;
p = (gam-1.)*r.*e ;
% trace de la solution
if mod(n,5) == 0
   subplot(3,1,1), plot(x,r) ;
   axis([0 lg 0. limmax]) ;
   title('Lax-Wendroff: densite')
   subplot(3,1,2), plot(x,u) ;
   axis([0 lg liminf limmax]) ;
   title('Lax-Wendroff: vitesse')
   subplot(3,1,3), plot(x,p) ;
   axis([0 lg 0. limmax]) ;
   title('Lax-Wendroff: pression')
   pause(0) ;
end
%
end








