%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright G. Allaire, Decembre 2000
%
% tuyere pour les equations d'Euler isentropiques
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% lg = longueur du domaine
lg = 1. ;
% nx = nombre de mailles
nx = 1000 ;
% dx = pas d'espace
dx = lg/nx ;
% r = densite, u = vitesse
% gam = gamma de la loi d'etat du gaz parfait
% p = r^gam = pression
gam = 1.4 ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 1er calcul: entree supersonique
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% condition aux limites
% entree
r_entree = 2. ;
u_entree = 5. ;
% sortie
r_sortie = 1. ;
u_sortie = 0. ;
% etat initial
ri = r_sortie ;
ui = u_sortie ;
% limitation du pas de temps par 
% la condition CFL
cfl = 0.9 ;
t = 0. ;
tmax = 0.5 ;
% nt = nombre de pas de temps
nt = 4000 ;
%
% initialisation
%
x=zeros(1,nx+1) ;
s=zeros(1,nx+1) ;
ds=zeros(1,nx+1) ;
r=zeros(1,nx+1) ;
u=zeros(1,nx+1) ;
rs=zeros(1,nx+1) ;
rsu=zeros(1,nx+1) ;
for i=1:nx+1
  x(i) = (i-1)*dx  ;
  s(i) = cosh(2*(x(i)-lg/2)) - 0.5 ;
  ds(i) = sinh(x(i)-lg/2) ;
  r(i) = ri ;
  u(i) = ui ;
end
rs = r.*s ;
rsu = r.*s.*u ;
ps = rs.^gam./s.^(gam-1) ;

% 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,s) ;
axis([0 lg 0. limmax]) ;
title('section de la tuyere')
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:6000
%
umax = max(abs(rsu./rs)) ;
cmax = max(sqrt(gam*rs.^(gam-1)./s.^(gam-1))) ;
dt = dx*cfl/(umax+cmax) ;
t = t + dt ;
ps = rs.^gam./s.^(gam-1) ;
rsp = wshift3('1D',rs,+1) ;
rsm = wshift3('1D',rs,-1) ;
rsup = wshift3('1D',rsu,+1) ;
rsum = wshift3('1D',rsu,-1) ;
psp = wshift3('1D',ps,+1) ;
psm = wshift3('1D',ps,-1) ;
rsn = (rsp+rsm)/2 - (0.5*dt/dx)*(rsup-rsum) ;
rsun = (rsup+rsum)/2 - (0.5*dt/dx)*... 
     ((rsup.^2)./rsp+psp-(rsum.^2)./rsm-psm)...
     + dt*ps.*ds./s ;

rsn(1) = r_entree*s(1) ;
rsn(nx+1) = r_sortie*s(nx+1) ;
rsun(1) = u_entree*rsn(1) ;
rsun(nx+1) = rsu(nx)/rs(nx) ;
rs = rsn ;
rsu = rsun ; 
r = rs./s ;
p = r.^gam ;
u = rsu./rs ;

if mod(n,10) == 0
   subplot(3,1,1), plot(x,r) ;
   axis([0 lg 0. 10.]) ;
   title('Lax-Friedrichs: densite')
   subplot(3,1,2), plot(x,u) ;
   axis([0 lg 0. 7.]) ;
   title('Lax-Friedrichs: vitesse')
   subplot(3,1,3), plot(x,p) ;
   axis([0 lg 0. 23.]) ;
   title(['Lax-Friedrichs: pression au temps =',num2str(t)])
   pause(0) ;
end
%
if t>=tmax
 break
end
%
end

   subplot(3,1,1), plot(x,r) ;
   axis([0 lg 1. 5.]) ;
   title(['Lax-Friedrichs: densite au temps =',num2str(t)])
   subplot(3,1,2), plot(x,u) ;
   axis([0 lg 4.6 5.1]) ;
   title(['Lax-Friedrichs: vitesse au temps =',num2str(t)])
   subplot(3,1,3), plot(x,rsu) ;
   title(['Lax-Friedrichs: debit au temps =',num2str(t)])


pause ;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 2eme calcul: entree subsonique
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% condition aux limites
% entree
r_entree = 1. ;
u_entree = 0.3 ;
% sortie
r_sortie = 1. ;
u_sortie = 0. ;
% etat initial
ri = r_sortie ;
ui = u_sortie ;
% limitation du pas de temps par 
% la condition CFL
cfl = 0.9 ;
t = 0. ;
tmax = 5. ;
% nt = nombre de pas de temps
nt = 12000 ;
%
% initialisation
%
x=zeros(1,nx+1) ;
s=zeros(1,nx+1) ;
ds=zeros(1,nx+1) ;
r=zeros(1,nx+1) ;
u=zeros(1,nx+1) ;
rs=zeros(1,nx+1) ;
rsu=zeros(1,nx+1) ;
for i=1:nx+1
  x(i) = (i-1)*dx  ;
  s(i) = cosh(2*(x(i)-lg/2)) - 0.5 ;
  ds(i) = sinh(x(i)-lg/2) ;
  r(i) = ri ;
  u(i) = ui ;
end
rs = r.*s ;
rsu = r.*s.*u ;
% p = r.^gam = pression
ps = rs.^gam./s.^(gam-1) ;

% 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,s) ;
axis([0 lg 0. limmax]) ;
title('section de la tuyere')
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
%
umax = max(abs(rsu./rs)) ;
cmax = max(sqrt(gam*rs.^(gam-1)./s.^(gam-1))) ;
dt = dx*cfl/(umax+cmax) ;
t = t + dt ;
ps = rs.^gam./s.^(gam-1) ;
rsp = wshift3('1D',rs,+1) ;
rsm = wshift3('1D',rs,-1) ;
rsup = wshift3('1D',rsu,+1) ;
rsum = wshift3('1D',rsu,-1) ;
psp = wshift3('1D',ps,+1) ;
psm = wshift3('1D',ps,-1) ;
rsn = (rsp+rsm)/2 - (0.5*dt/dx)*(rsup-rsum) ;
rsun = (rsup+rsum)/2 - (0.5*dt/dx)*... 
     ((rsup.^2)./rsp+psp-(rsum.^2)./rsm-psm)...
     + dt*ps.*ds./s ;

rsn(1) = r_entree*s(1) ;
rsn(nx+1) = r_sortie*s(nx+1) ;
rsun(1) = u_entree*rsn(1) ;
rsun(nx+1) = rsu(nx)/rs(nx) ;
rs = rsn ;
rsu = rsun ; 
r = rs./s ;
p = r.^gam ;
u = rsu./rs ;

nnn = 10 ;
if n >= 3000
   nnn = 100 ;
end
if mod(n,nnn) == 0
   subplot(3,1,1), plot(x,r) ;
   axis([0 lg 0. 2.]) ;
   title('Lax-Friedrichs: densite')
   subplot(3,1,2), plot(x,u) ;
   axis([0 lg -0.3 0.6]) ;
   title('Lax-Friedrichs: vitesse')
   subplot(3,1,3), plot(x,p) ;
   axis([0 lg 0. 3.]) ;
   title(['Lax-Friedrichs: pression au temps =',num2str(t)])
   pause(0) ;
end
%
if t>=tmax
 break
end
%
end

   subplot(3,1,1), plot(x,r) ;
   axis([0 lg 0.9 1.6]) ;
   title(['Lax-Friedrichs: densite au temps =',num2str(t)])
   subplot(3,1,2), plot(x,u) ;
   axis([0 lg 0.2 0.5]) ;
   title(['Lax-Friedrichs: vitesse au temps =',num2str(t)])
   subplot(3,1,3), plot(x,rsu) ;
   title(['Lax-Friedrichs: debit au temps =',num2str(t)])

pause ;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 3eme calcul: choc au col
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% condition aux limites
% entree
r_entree = 1. ;
u_entree = 0.7 ;
% sortie
r_sortie = 1. ;
u_sortie = 0. ;
% etat initial
ri = r_sortie ;
ui = u_sortie ;
% limitation du pas de temps par 
% la condition CFL
cfl = 0.9 ;
t = 0. ;
tmax = 5. ;
% nt = nombre de pas de temps
nt = 15000 ;
%
% initialisation
%
x=zeros(1,nx+1) ;
s=zeros(1,nx+1) ;
ds=zeros(1,nx+1) ;
r=zeros(1,nx+1) ;
u=zeros(1,nx+1) ;
rs=zeros(1,nx+1) ;
rsu=zeros(1,nx+1) ;
for i=1:nx+1
  x(i) = (i-1)*dx  ;
  s(i) = cosh(2*(x(i)-lg/2)) - 0.5 ;
  ds(i) = sinh(x(i)-lg/2) ;
  r(i) = ri ;
  u(i) = ui ;
end
rs = r.*s ;
rsu = r.*s.*u ;
% p = r.^gam = pression
ps = rs.^gam./s.^(gam-1) ;

% 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,s) ;
axis([0 lg 0. limmax]) ;
title('section de la tuyere')
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
%
umax = max(abs(rsu./rs)) ;
cmax = max(sqrt(gam*rs.^(gam-1)./s.^(gam-1))) ;
dt = dx*cfl/(umax+cmax) ;
t = t + dt ;
ps = rs.^gam./s.^(gam-1) ;
rsp = wshift3('1D',rs,+1) ;
rsm = wshift3('1D',rs,-1) ;
rsup = wshift3('1D',rsu,+1) ;
rsum = wshift3('1D',rsu,-1) ;
psp = wshift3('1D',ps,+1) ;
psm = wshift3('1D',ps,-1) ;
rsn = (rsp+rsm)/2 - (0.5*dt/dx)*(rsup-rsum) ;
rsun = (rsup+rsum)/2 - (0.5*dt/dx)*... 
     ((rsup.^2)./rsp+psp-(rsum.^2)./rsm-psm)...
     + dt*ps.*ds./s ;

rsn(1) = r_entree*s(1) ;
rsn(nx+1) = r_sortie*s(nx+1) ;
rsun(1) = u_entree*rsn(1) ;
rsun(nx+1) = rsu(nx)/rs(nx) ;
rs = rsn ;
rsu = rsun ; 
r = rs./s ;
p = r.^gam ;
u = rsu./rs ;

nnn = 10 ;
if n >= 2000
   nnn = 50 ;
end
if mod(n,nnn) == 0
   subplot(3,1,1), plot(x,r) ;
   axis([0 lg 0. 3.]) ;
   title('Lax-Friedrichs: densite')
   subplot(3,1,2), plot(x,u) ;
   axis([0 lg -0.2 1.8]) ;
   title('Lax-Friedrichs: vitesse')
   subplot(3,1,3), plot(x,p) ;
   axis([0 lg 0. 4.]) ;
   title(['Lax-Friedrichs: pression au temps =',num2str(t)])
   pause(0) ;
end
%
if t>=tmax
 break
end
%
end

   subplot(3,1,1), plot(x,r) ;
   axis([0 lg 0.5 1.5]) ;
   title(['Lax-Friedrichs: densite au temps =',num2str(t)])
   subplot(3,1,2), plot(x,u) ;
   axis([0 lg 0.5 1.8]) ;
   title(['Lax-Friedrichs: vitesse au temps =',num2str(t)])
   subplot(3,1,3), plot(x,rsu) ;
   title(['Lax-Friedrichs: debit au temps =',num2str(t)])


