%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright G. Allaire, Decembre 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.) ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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
p = (gam-1.)*r.*e ;

liminf = -0.1 ;
limmaxd = 1.1 ;
limmaxv = 1.1 ;
limmaxp = 1.1 ;
subplot(3,1,1), plot(x,r) ;
axis([0 lg liminf limmaxd]) ;
title('Van Leer: densite')
subplot(3,1,2), plot(x,u) ;
axis([0 lg liminf limmaxv]) ;
title('Van Leer: vitesse')
subplot(3,1,3), plot(x,p) ;
axis([0 lg liminf limmaxp]) ;
title('Van Leer: pression')
pause ;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% marche en temps : Van Leer
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
for n=1:300
%
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. limmaxd]) ;
   title('Van Leer: densite')
   subplot(3,1,2), plot(x,u) ;
   axis([0 lg liminf limmaxv]) ;
   title('Van Leer: vitesse')
   subplot(3,1,3), plot(x,p) ;
   axis([0 lg 0. limmaxp]) ;
   title('Van Leer: pression')
   pause(0) ;
end
%
end


pause ;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% etats initiaux a gauche (l) et a droite (r)
% r = densite, u = vitesse, e = energie interne
vitesse = 1.2 ;
rl = 1. ;
ul = -vitesse ;
el = 0.4/(gam-1.) ;
rr = 1. ;
ur = vitesse ;
er = 0.4/(gam-1.) ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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
p = (gam-1.)*r.*e ;

liminf = -0.1 ;
limmax = 1.1 ;
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 -vitesse-0.1 vitesse+0.1]) ;
title('Van Leer: vitesse')
subplot(3,1,3), plot(x,p) ;
axis([0 lg 0. 0.5]) ;
title('Van Leer: pression')
pause ;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% marche en temps : Van Leer
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
for n=1:300
%
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 -vitesse-0.1 vitesse+0.1]) ;
   title('Van Leer: vitesse')
   subplot(3,1,3), plot(x,p) ;
   axis([0 lg 0. 0.5]) ;
   title('Van Leer: pression')
   pause(0) ;
end
%
end


pause ;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% etats initiaux a gauche (l) et a droite (r)
% r = densite, u = vitesse, e = energie interne
rl = 1. ;
ul = 0. ;
el = 0.01/(gam-1.) ;
rr = 1. ;
ur = 0. ;
er = 100./(gam-1.) ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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
p = (gam-1.)*r.*e ;


subplot(3,1,1), plot(x,r) ;
axis([0 lg 0. 6.]) ;
title('Van Leer: densite')
subplot(3,1,2), plot(x,u) ;
axis([0 lg -7. 0.5]) ;
title('Van Leer: vitesse')
subplot(3,1,3), plot(x,p) ;
axis([0 lg -1. 110.]) ;
title('Van Leer: pression')
pause ;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% marche en temps : Van Leer
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
for n=1:500
%
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. 6]) ;
   title('Van Leer: densite')
   subplot(3,1,2), plot(x,u) ;
   axis([0 lg -7. 0.5]) ;
   title('Van Leer: vitesse')
   subplot(3,1,3), plot(x,p) ;
   axis([0 lg -1. 110.]) ;
   title('Van Leer: pression')
   pause(0) ;
end
%
end


