%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright G. Allaire, Decembre 2000
%
% schema de Roe pour les equations d'Euler en 1-D
% correction entropique ou non
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% lg = longueur du domaine
lg = 1. ;
% nx = nombre de mailles
nx = 1000 ;
dx = lg/nx ;
% t = temps, tmax = temps maximal, nt = nombre max de pas de temps
t = 0. ;
nt = 400 ;
tmax = 0.15 ;
% loi d'etat du gaz parfait
gamma = 1.4 ;
% sol = variables conservatives = (rho, rho*u, E)
sol = zeros(3,nx) ; 
% flux = flux par maille
flux = zeros(3,nx) ;
fluxp = zeros(3,nx) ;
% glux = flux numerique a l'interface
glux = zeros(3,nx) ;
gluxm = zeros(3,nx) ;
absroe = zeros(3,nx) ;
% vp = vecteurs propres de la jacobienne
vp1 = ones(3,nx) ;
vp2 = ones(3,nx) ;
vp3 = ones(3,nx) ;
% u = vitess, p = pression
u = zeros(1,nx) ;
p = zeros(1,nx) ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% le tube a choc de Sod
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% etats initiaux a gauche (l) et a droite (r)
% rho = densite, u = vitesse, p = pression
rhol = 1. ;
ul = 0. ;
pl = 1. ;
rhor = 0.125 ;
ur = 0. ;
pr = 0.1 ;
%
% initialisation
%
x=zeros(1,nx) ;

for i=1:nx
  x(i) = (i-1)*dx  ;
  if x(i) < lg/2
     sol(1,i) = rhol ;     
     sol(2,i) = rhol*ul;
     sol(3,i) = pl/(gamma-1.) + 0.5*rhol*ul*ul ;
  else
     sol(1,i) = rhor ;     
     sol(2,i) = rhor*ur ;
     sol(3,i) = pr/(gamma-1.) + 0.5*rhor*ur*ur ;
  end
end
u = sol(2,:)./sol(1,:) ;
p = (gamma-1.)*(sol(3,:) - 0.5*sol(1,:).*u.*u) ;
    
% trace des conditions initiales
liminf = -0.1 ;
limmax = 1.1 ;
subplot(3,1,1), plot(x,sol(1,:)) ;
axis([0 lg 0. limmax]) ;
title('Roe: densite')
subplot(3,1,2), plot(x,u) ;
axis([0 lg liminf limmax]) ;
title('Roe: vitesse')
subplot(3,1,3), plot(x,p) ;
axis([0 lg 0. limmax]) ;
title('Roe: pression')
pause ;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% marche en temps : schema de Roe
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% 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

% on recupere u et p a partir des variables conservatives
if all(sign(sol(1,:))+1)
   u = sol(2,:)./sol(1,:) ;
   p = (gamma-1.)*(sol(3,:) - 0.5*sol(1,:).*u.*u) ;
else 
  'densite negative ou nulle'
  return
end
 
% moyennes de Roe
        
w = sqrt(sol(1,:));
wp = wshift3('1D',w,+1) ;
up = wshift3('1D',u,+1) ; 
pp = wshift3('1D',p,+1) ;
solp = wshift3('2D',sol,[0 +1]) ;
utild = (w.*u + wp.*up)./(w + wp);
htild = (w.*(sol(3,:) + p) + wp.*(solp(3,:) + pp))./...
        (w.* sol(1,:)      + wp.* solp(1,:)) ;
rhotild = w.*wp ;

% c2son = vitesse du son au carre
c2son = (gamma - 1.)*(htild - 0.5 * utild.*utild) ;

% cson = vitesse du son
if all(sign(c2son)+1.) 
   cson = sqrt(c2son) ;
else 
   'perte d hyperbolicite numerique'
   'carre de la vitesse du son negative'
   'dans les mailles'
   for i=1:nx
      if c2son(i) <= 0. 
         i
      end
   end
   return
end

% vecteurs propres de la jacobienne

vp1(2,:) = utild - cson ;
vp1(3,:) = htild - utild.*cson ;
vp2(2,:) = utild ;
vp2(3,:) = htild - c2son/(gamma-1.) ;
vp3(2,:) = utild + cson ;
vp3(3,:) = htild + utild.*cson ;

% coefficients de l'increment dans la base 
% des vecteurs propres

alpha1 = ( (pp-p) - cson.*rhotild.*(up-u) )./(2*c2son) ;
alpha2 = (solp(1,:)-sol(1,:)) - (pp-p)./c2son ;
alpha3 = ( (pp-p) + cson.*rhotild.*(up-u) )./(2*c2son) ;

% calcul du flux

flux(1,:) =     sol(2,:)      ;
flux(2,:) = u.* sol(2,:) + p  ;
flux(3,:) = u.*(sol(3,:) + p) ;
fluxp = wshift3('2D',flux,[0 +1]) ;

% terme de viscosite numerique dans le flux numerique

absroe(1,:) = abs(utild-cson).*alpha1.*vp1(1,:) ;
absroe(2,:) = abs(utild-cson).*alpha1.*vp1(2,:) ;
absroe(3,:) = abs(utild-cson).*alpha1.*vp1(3,:) ;

absroe(1,:) = absroe(1,:) + abs(utild).*alpha2.*vp2(1,:) ;
absroe(2,:) = absroe(2,:) + abs(utild).*alpha2.*vp2(2,:) ;
absroe(3,:) = absroe(3,:) + abs(utild).*alpha2.*vp2(3,:) ;

absroe(1,:) = absroe(1,:) + abs(utild+cson).*alpha3.*vp3(1,:) ;
absroe(2,:) = absroe(2,:) + abs(utild+cson).*alpha3.*vp3(2,:) ;
absroe(3,:) = absroe(3,:) + abs(utild+cson).*alpha3.*vp3(3,:) ;

% flux numerique a l'interface i+1/2

glux = 0.5*( flux + fluxp - absroe ) ;
gluxm = wshift3('2D',glux,[0 -1]) ;

% calcul de la condition CFL
cfl = max(max(abs(utild+cson)),max(abs(utild-cson))) ;

% mise a jour

if cfl <= 0. 
   'cfl nulle'
   return
else
   dt = 0.9*dx/cfl ;
end

t = t + dt ;

aux1 = sol(:,1) ;
aux2 = sol(:,nx) ;

sol = sol - (dt/dx)*(glux-gluxm) ;

% correction des effets de bord
sol(:,1) = aux1 ;
sol(:,nx) = aux2 ;

% trace de la solution
if mod(n,5) == 0
   u = sol(2,:)./sol(1,:) ;
   p = (gamma-1.)*(sol(3,:) - 0.5*sol(1,:).*u.*u) ;
                
   subplot(3,1,1), plot(x,sol(1,:)) ;
   axis([0 lg 0. limmax]) ;
   title(['Roe: densite au temps =',num2str(t)])
   subplot(3,1,2), plot(x,u) ;
   axis([0 lg liminf limmax]) ;
   title(['Roe: vitesse au temps =',num2str(t)])
   subplot(3,1,3), plot(x,p) ;
   axis([0 lg 0. limmax]) ;
   title(['Roe: pression au temps =',num2str(t)])
   pause(0) ;
end
%
if t>=tmax
 break
end
%
end

pause ;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% un autre tube a choc avec point sonique
% schema de Roe sans correction entropique
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
t = 0. ;
% etats initiaux a gauche (l) et a droite (r)
% rho = densite, u = vitesse, p = pression
rhol = 1. ;
ul = 0.75 ;
pl = 1.0 ;
rhor = 0.125 ;
ur = 0. ;
pr = 0.1 ;
%
% initialisation
%
x=zeros(1,nx) ;

for i=1:nx
  x(i) = (i-1)*dx  ;
  if x(i) < lg/2
     sol(1,i) = rhol ;     
     sol(2,i) = rhol*ul;
     sol(3,i) = pl/(gamma-1.) + 0.5*rhol*ul*ul ;
  else
     sol(1,i) = rhor ;     
     sol(2,i) = rhor*ur ;
     sol(3,i) = pr/(gamma-1.) + 0.5*rhor*ur*ur ;
  end
end
u = sol(2,:)./sol(1,:) ;
p = (gamma-1.)*(sol(3,:) - 0.5*sol(1,:).*u.*u) ;
    
% trace des conditions initiales
liminf = -0.1 ;
limmax = 1.1 ;
subplot(3,1,1), plot(x,sol(1,:)) ;
axis([0 lg 0. 1.1]) ;
title('Roe: densite')
subplot(3,1,2), plot(x,u) ;
axis([0 lg liminf 1.6]) ;
title('Roe: vitesse')
subplot(3,1,3), plot(x,p) ;
axis([0 lg 0. 1.1]) ;
title('Roe: pression')
pause ;



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% marche en temps : schema de Roe
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% 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

% on recupere u et p a partir des variables conservatives
if all(sign(sol(1,:))+1)
   u = sol(2,:)./sol(1,:) ;
   p = (gamma-1.)*(sol(3,:) - 0.5*sol(1,:).*u.*u) ;
else 
  'densite negative ou nulle'
  return
end
 
% moyennes de Roe
        
w = sqrt(sol(1,:));
wp = wshift3('1D',w,+1) ;
up = wshift3('1D',u,+1) ; 
pp = wshift3('1D',p,+1) ;
solp = wshift3('2D',sol,[0 +1]) ;
utild = (w.*u + wp.*up)./(w + wp);
htild = (w.*(sol(3,:) + p) + wp.*(solp(3,:) + pp))./...
        (w.* sol(1,:)      + wp.* solp(1,:)) ;
rhotild = w.*wp ;

% c2son = vitesse du son au carre
c2son = (gamma - 1.)*(htild - 0.5 * utild.*utild) ;

% cson = vitesse du son
if all(sign(c2son)+1.) 
   cson = sqrt(c2son) ;
else 
   'perte d hyperbolicite numerique'
   'carre de la vitesse du son negative'
   'dans les mailles'
   for i=1:nx
      if c2son(i) <= 0. 
         i
      end
   end
   return
end

% vecteurs propres de la jacobienne

vp1(2,:) = utild - cson ;
vp1(3,:) = htild - utild.*cson ;
vp2(2,:) = utild ;
vp2(3,:) = htild - c2son/(gamma-1.) ;
vp3(2,:) = utild + cson ;
vp3(3,:) = htild + utild.*cson ;

% coefficients de l'increment dans la base 
% des vecteurs propres

alpha1 = ( (pp-p) - cson.*rhotild.*(up-u) )./(2*c2son) ;
alpha2 = (solp(1,:)-sol(1,:)) - (pp-p)./c2son ;
alpha3 = ( (pp-p) + cson.*rhotild.*(up-u) )./(2*c2son) ;

% calcul du flux

flux(1,:) =     sol(2,:)      ;
flux(2,:) = u.* sol(2,:) + p  ;
flux(3,:) = u.*(sol(3,:) + p) ;
fluxp = wshift3('2D',flux,[0 +1]) ;

% terme de viscosite numerique dans le flux numerique

absroe(1,:) = abs(utild-cson).*alpha1.*vp1(1,:) ;
absroe(2,:) = abs(utild-cson).*alpha1.*vp1(2,:) ;
absroe(3,:) = abs(utild-cson).*alpha1.*vp1(3,:) ;

absroe(1,:) = absroe(1,:) + abs(utild).*alpha2.*vp2(1,:) ;
absroe(2,:) = absroe(2,:) + abs(utild).*alpha2.*vp2(2,:) ;
absroe(3,:) = absroe(3,:) + abs(utild).*alpha2.*vp2(3,:) ;

absroe(1,:) = absroe(1,:) + abs(utild+cson).*alpha3.*vp3(1,:) ;
absroe(2,:) = absroe(2,:) + abs(utild+cson).*alpha3.*vp3(2,:) ;
absroe(3,:) = absroe(3,:) + abs(utild+cson).*alpha3.*vp3(3,:) ;

% flux numerique a l'interface i+1/2

glux = 0.5*( flux + fluxp - absroe ) ;
gluxm = wshift3('2D',glux,[0 -1]) ;

% calcul de la condition CFL
cfl = max(max(abs(utild+cson)),max(abs(utild-cson))) ;

% mise a jour

if cfl <= 0. 
   'cfl nulle'
   return
else
   dt = 0.9*dx/cfl ;
end

t = t + dt ;

aux1 = sol(:,1) ;
aux2 = sol(:,nx) ;

sol = sol - (dt/dx)*(glux-gluxm) ;

% correction des effets de bord
sol(:,1) = aux1 ;
sol(:,nx) = aux2 ;

% trace de la solution
if mod(n,5) == 0
   u = sol(2,:)./sol(1,:) ;
   p = (gamma-1.)*(sol(3,:) - 0.5*sol(1,:).*u.*u) ;
                
   subplot(3,1,1), plot(x,sol(1,:)) ;
   axis([0 lg 0. limmax]) ;
   title(['Roe: densite au temps =',num2str(t)])
   subplot(3,1,2), plot(x,u) ;
   axis([0 lg liminf 1.6]) ;
   title(['Roe: vitesse au temps =',num2str(t)])
   subplot(3,1,3), plot(x,p) ;
   axis([0 lg 0. limmax]) ;
   title(['Roe: pression au temps =',num2str(t)])
   pause(0) ;
end
%
if t>=tmax
 break
end
%
end




pause ;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% un autre tube a choc avec point sonique
% schema de Roe avec correction entropique
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
t = 0. ;
eps = 5.e-2 ;
% etats initiaux a gauche (l) et a droite (r)
% rho = densite, u = vitesse, p = pression
rhol = 1. ;
ul = 0.75 ;
pl = 1.0 ;
rhor = 0.125 ;
ur = 0. ;
pr = 0.1 ;
%
% initialisation
%
x=zeros(1,nx) ;

for i=1:nx
  x(i) = (i-1)*dx  ;
  if x(i) < lg/2
     sol(1,i) = rhol ;     
     sol(2,i) = rhol*ul;
     sol(3,i) = pl/(gamma-1.) + 0.5*rhol*ul*ul ;
  else
     sol(1,i) = rhor ;     
     sol(2,i) = rhor*ur ;
     sol(3,i) = pr/(gamma-1.) + 0.5*rhor*ur*ur ;
  end
end
u = sol(2,:)./sol(1,:) ;
p = (gamma-1.)*(sol(3,:) - 0.5*sol(1,:).*u.*u) ;
    
% trace des conditions initiales
liminf = -0.1 ;
limmax = 1.1 ;
subplot(3,1,1), plot(x,sol(1,:)) ;
axis([0 lg 0. 1.1]) ;
title('Roe: densite')
subplot(3,1,2), plot(x,u) ;
axis([0 lg liminf 1.6]) ;
title('Roe: vitesse')
subplot(3,1,3), plot(x,p) ;
axis([0 lg 0. 1.1]) ;
title('Roe: pression')
pause ;



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% marche en temps : schema de Roe
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% 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

% on recupere u et p a partir des variables conservatives
if all(sign(sol(1,:))+1)
   u = sol(2,:)./sol(1,:) ;
   p = (gamma-1.)*(sol(3,:) - 0.5*sol(1,:).*u.*u) ;
else 
  'densite negative ou nulle'
  return
end
 
% moyennes de Roe
        
w = sqrt(sol(1,:));
wp = wshift3('1D',w,+1) ;
up = wshift3('1D',u,+1) ; 
pp = wshift3('1D',p,+1) ;
solp = wshift3('2D',sol,[0 +1]) ;
utild = (w.*u + wp.*up)./(w + wp);
htild = (w.*(sol(3,:) + p) + wp.*(solp(3,:) + pp))./...
        (w.* sol(1,:)      + wp.* solp(1,:)) ;
rhotild = w.*wp ;

% c2son = vitesse du son au carre
c2son = (gamma - 1.)*(htild - 0.5 * utild.*utild) ;

% cson = vitesse du son
if all(sign(c2son)+1.) 
   cson = sqrt(c2son) ;
else 
   'perte d hyperbolicite numerique'
   'carre de la vitesse du son negative'
   'dans les mailles'
   for i=1:nx
      if c2son(i) <= 0. 
         i
      end
   end
   return
end

% vecteurs propres de la jacobienne

vp1(2,:) = utild - cson ;
vp1(3,:) = htild - utild.*cson ;
vp2(2,:) = utild ;
vp2(3,:) = htild - c2son/(gamma-1.) ;
vp3(2,:) = utild + cson ;
vp3(3,:) = htild + utild.*cson ;

% coefficients de l'increment dans la base 
% des vecteurs propres

alpha1 = ( (pp-p) - cson.*rhotild.*(up-u) )./(2*c2son) ;
alpha2 = (solp(1,:)-sol(1,:)) - (pp-p)./c2son ;
alpha3 = ( (pp-p) + cson.*rhotild.*(up-u) )./(2*c2son) ;

% calcul du flux

flux(1,:) =     sol(2,:)      ;
flux(2,:) = u.* sol(2,:) + p  ;
flux(3,:) = u.*(sol(3,:) + p) ;
fluxp = wshift3('2D',flux,[0 +1]) ;

% terme de viscosite numerique dans le flux numerique

absroe(1,:) = sqrt((utild-cson).^2+eps^2).*alpha1.*vp1(1,:) ;
absroe(2,:) = sqrt((utild-cson).^2+eps^2).*alpha1.*vp1(2,:) ;
absroe(3,:) = sqrt((utild-cson).^2+eps^2).*alpha1.*vp1(3,:) ;

absroe(1,:) = absroe(1,:) + sqrt(utild.^2+eps^2).*alpha2.*vp2(1,:) ;
absroe(2,:) = absroe(2,:) + sqrt(utild.^2+eps^2).*alpha2.*vp2(2,:) ;
absroe(3,:) = absroe(3,:) + sqrt(utild.^2+eps^2).*alpha2.*vp2(3,:) ;

absroe(1,:) = absroe(1,:) + sqrt((utild+cson).^2+eps^2).*alpha3.*vp3(1,:) ;
absroe(2,:) = absroe(2,:) + sqrt((utild+cson).^2+eps^2).*alpha3.*vp3(2,:) ;
absroe(3,:) = absroe(3,:) + sqrt((utild+cson).^2+eps^2).*alpha3.*vp3(3,:) ;

% flux numerique a l'interface i+1/2

glux = 0.5*( flux + fluxp - absroe ) ;
gluxm = wshift3('2D',glux,[0 -1]) ;

% calcul de la condition CFL
cfl = max(max(abs(utild+cson)),max(abs(utild-cson))) ;

% mise a jour

if cfl <= 0. 
   'cfl nulle'
   return
else
   dt = 0.9*dx/cfl ;
end

t = t + dt ;

aux1 = sol(:,1) ;
aux2 = sol(:,nx) ;

sol = sol - (dt/dx)*(glux-gluxm) ;

% correction des effets de bord
sol(:,1) = aux1 ;
sol(:,nx) = aux2 ;

% trace de la solution
if mod(n,5) == 0
   u = sol(2,:)./sol(1,:) ;
   p = (gamma-1.)*(sol(3,:) - 0.5*sol(1,:).*u.*u) ;
                
   subplot(3,1,1), plot(x,sol(1,:)) ;
   axis([0 lg 0. limmax]) ;
   title(['Roe: densite au temps =',num2str(t)])
   subplot(3,1,2), plot(x,u) ;
   axis([0 lg liminf 1.6]) ;
   title(['Roe: vitesse au temps =',num2str(t)])
   subplot(3,1,3), plot(x,p) ;
   axis([0 lg 0. limmax]) ;
   title(['Roe: pression au temps =',num2str(t)])
   pause(0) ;
end
%
if t>=tmax
 break
end
%
end