/////////////////////////////////////////////////////////////////////////////
// Copyright G. Allaire, Decembre 2002, Janvier 2012
//
// schema de Roe pour les equations d'Euler en 1-D
// correction entropique ou non
// 
///////////////////////////////////////////////////////////////////////////
//
xset("font",2,3) ;
xset("thickness",2) ;
//
// 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
gama = 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/(gama-1.) + 0.5*rhol*ul*ul ;
  else
     sol(1,i) = rhor ;     
     sol(2,i) = rhor*ur ;
     sol(3,i) = pr/(gama-1.) + 0.5*rhor*ur*ur ;
  end
end
r = sol(1,:) ;
u = sol(2,:)./sol(1,:) ;
p = (gama-1.)*(sol(3,:) - 0.5*sol(1,:).*u.*u) ;
e = (p./r)/(gama-1.) ;
    
// trace des conditions initiales
liminf = -0.1 ;
limmax = 1.1 ;
clf()

tics=[4,10,4,10];
rect = [0,liminf,lg,limmax];
plotframe([0,liminf,lg,limmax],tics,[%f,%t],...
["Roe: Densité","",""],[0,0,0.5,0.5]);
plot2d(x,r,2,"000")
plotframe([0,liminf,lg,limmax],tics,[%f,%t],["Vitesse","",""],[0.5,0,0.5,0.5])
plot2d(x,u,3,"000")
plotframe([0,1.5,lg,3],tics,[%f,%t],["Energie interne","",""],[0,0.5,0.5,0.5])
plot2d(x,e,4,"000")
plotframe(rect,tics,[%f,%t],["Pression","",""],[0.5,0.5,0.5,0.5])
plot2d(x,p,4,"000")

halt() ;

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// marche en temps : schema de Roe
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//
// si u est le vecteur des valeurs u(i) alors 
// shiftn('+1',u) est le vecteur des u(i+1) 
// shiftn('-1',u) est le vecteur des u(i-1)
//
for n=1:nt

// on recupere u et p a partir des variables conservatives
if sign(sol(1,:)) > 0
   u = sol(2,:)./sol(1,:) ;
   p = (gama-1.)*(sol(3,:) - 0.5*sol(1,:).*u.*u) ;
else 
  'densite negative ou nulle'
  return
end
 
// moyennes de Roe
        
w = sqrt(sol(1,:));
wp = shiftn('+1',w) ;
up = shiftn('+1',u) ; 
pp = shiftn('+1',p) ;
solp = shiftn('+1',sol) ;
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 = (gama - 1.)*(htild - 0.5 * utild.*utild) ;

// cson = vitesse du son
if sign(c2son) >= 0 
   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/(gama-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 = shiftn('+1',flux) ;

// 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 = shiftn('-1',glux) ;

// 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 modulo(n,8) == 0
r = sol(1,:) ;
u = sol(2,:)./sol(1,:) ;
p = (gama-1.)*(sol(3,:) - 0.5*sol(1,:).*u.*u) ;
e = (p./r)/(gama-1.) ;
clf()
tics=[4,10,4,10];
rect = [0,liminf,lg,limmax];
plotframe([0,liminf,lg,limmax],tics,[%f,%t],...
["Roe: Densité","",""],[0,0,0.5,0.5]);
plot2d(x,r,2,"000")
plotframe([0,liminf,lg,limmax],tics,[%f,%t],["Vitesse","",""],[0.5,0,0.5,0.5])
plot2d(x,u,3,"000")
plotframe([0,1.5,lg,3],tics,[%f,%t],["Energie interne","",""],[0,0.5,0.5,0.5])
plot2d(x,e,4,"000")
plotframe(rect,tics,[%f,%t],["Pression","",""],[0.5,0.5,0.5,0.5])
plot2d(x,p,4,"000")
xpause(100000) ;
//   title(['Roe: densite au temps =',num2str(t)])

end
//
if t>=tmax
 break
end
//
end

halt() ;

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// 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/(gama-1.) + 0.5*rhol*ul*ul ;
  else
     sol(1,i) = rhor ;     
     sol(2,i) = rhor*ur ;
     sol(3,i) = pr/(gama-1.) + 0.5*rhor*ur*ur ;
  end
end
r = sol(1,:) ;
u = sol(2,:)./sol(1,:) ;
p = (gama-1.)*(sol(3,:) - 0.5*sol(1,:).*u.*u) ;
e = (p./r)/(gama-1.) ;
    
// trace des conditions initiales
liminf = -0.1 ;
limmax = 1.1 ;
clf()

tics=[4,10,4,10];
rect = [0,liminf,lg,limmax];
plotframe([0,liminf,lg,limmax],tics,[%f,%t],...
["Roe: Densité","",""],[0,0,0.5,0.5]);
plot2d(x,r,2,"000")
plotframe([0,liminf,lg,1.5],tics,[%f,%t],["Vitesse","",""],[0.5,0,0.5,0.5])
plot2d(x,u,3,"000")
plotframe([0,1.9,lg,4],tics,[%f,%t],["Energie interne","",""],[0,0.5,0.5,0.5])
plot2d(x,e,4,"000")
plotframe(rect,tics,[%f,%t],["Pression","",""],[0.5,0.5,0.5,0.5])
plot2d(x,p,4,"000")

halt() ;

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// marche en temps : schema de Roe
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//
// si u est le vecteur des valeurs u(i) alors 
// shiftn('+1',u) est le vecteur des u(i+1) 
// shiftn('-1',u) est le vecteur des u(i-1)
//
for n=1:nt

// on recupere u et p a partir des variables conservatives
if sign(sol(1,:)) > 0
   u = sol(2,:)./sol(1,:) ;
   p = (gama-1.)*(sol(3,:) - 0.5*sol(1,:).*u.*u) ;
else 
  'densite negative ou nulle'
  return
end
 
// moyennes de Roe
        
w = sqrt(sol(1,:));
wp = shiftn('+1',w) ;
up = shiftn('+1',u) ; 
pp = shiftn('+1',p) ;
solp = shiftn('+1',sol) ;
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 = (gama - 1.)*(htild - 0.5 * utild.*utild) ;

// cson = vitesse du son
if sign(c2son) >= 0 
   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/(gama-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 = shiftn('+1',flux) ;

// 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 = shiftn('-1',glux) ;

// 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 modulo(n,8) == 0
r = sol(1,:) ;
u = sol(2,:)./sol(1,:) ;
p = (gama-1.)*(sol(3,:) - 0.5*sol(1,:).*u.*u) ;
e = (p./r)/(gama-1.) ;
clf()
tics=[4,10,4,10];
rect = [0,liminf,lg,limmax];
plotframe([0,liminf,lg,limmax],tics,[%f,%t],...
["Roe: Densité","",""],[0,0,0.5,0.5]);
plot2d(x,r,2,"000")
plotframe([0,liminf,lg,1.5],tics,[%f,%t],["Vitesse","",""],[0.5,0,0.5,0.5])
plot2d(x,u,3,"000")
plotframe([0,1.9,lg,4],tics,[%f,%t],["Energie interne","",""],[0,0.5,0.5,0.5])
plot2d(x,e,4,"000")
plotframe(rect,tics,[%f,%t],["Pression","",""],[0.5,0.5,0.5,0.5])
plot2d(x,p,4,"000")
xpause(100000) ;
end
//
if t>=tmax
 break
end
//
end


halt() ;

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// 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/(gama-1.) + 0.5*rhol*ul*ul ;
  else
     sol(1,i) = rhor ;     
     sol(2,i) = rhor*ur ;
     sol(3,i) = pr/(gama-1.) + 0.5*rhor*ur*ur ;
  end
end
r = sol(1,:) ;
u = sol(2,:)./sol(1,:) ;
p = (gama-1.)*(sol(3,:) - 0.5*sol(1,:).*u.*u) ;
e = (p./r)/(gama-1.) ;
    
// trace des conditions initiales
liminf = -0.1 ;
limmax = 1.1 ;
clf()

tics=[4,10,4,10];
rect = [0,liminf,lg,limmax];
plotframe([0,liminf,lg,limmax],tics,[%f,%t],...
["Roe: Densité","",""],[0,0,0.5,0.5]);
plot2d(x,r,2,"000")
plotframe([0,liminf,lg,1.5],tics,[%f,%t],["Vitesse","",""],[0.5,0,0.5,0.5])
plot2d(x,u,3,"000")
plotframe([0,1.9,lg,4],tics,[%f,%t],["Energie interne","",""],[0,0.5,0.5,0.5])
plot2d(x,e,4,"000")
plotframe(rect,tics,[%f,%t],["Pression","",""],[0.5,0.5,0.5,0.5])
plot2d(x,p,4,"000")

halt() ;


////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// marche en temps : schema de Roe
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//
// si u est le vecteur des valeurs u(i) alors 
// shiftn('+1',u) est le vecteur des u(i+1) 
// shiftn('-1',u) est le vecteur des u(i-1)
//
for n=1:nt

// on recupere u et p a partir des variables conservatives
if sign(sol(1,:)) > 0
   u = sol(2,:)./sol(1,:) ;
   p = (gama-1.)*(sol(3,:) - 0.5*sol(1,:).*u.*u) ;
else 
  'densite negative ou nulle'
  return
end
 
// moyennes de Roe
        
w = sqrt(sol(1,:));
wp = shiftn('+1',w) ;
up = shiftn('+1',u) ; 
pp = shiftn('+1',p) ;
solp = shiftn('+1',sol) ;
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 = (gama - 1.)*(htild - 0.5 * utild.*utild) ;

// cson = vitesse du son
if sign(c2son) >= 0 
   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/(gama-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 = shiftn('+1',flux) ;

// 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 = shiftn('-1',glux) ;

// 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 modulo(n,8) == 0
r = sol(1,:) ;
u = sol(2,:)./sol(1,:) ;
p = (gama-1.)*(sol(3,:) - 0.5*sol(1,:).*u.*u) ;
e = (p./r)/(gama-1.) ;
clf()
tics=[4,10,4,10];
rect = [0,liminf,lg,limmax];
plotframe([0,liminf,lg,limmax],tics,[%f,%t],...
["Roe: Densité","",""],[0,0,0.5,0.5]);
plot2d(x,r,2,"000")
plotframe([0,liminf,lg,1.5],tics,[%f,%t],["Vitesse","",""],[0.5,0,0.5,0.5])
plot2d(x,u,3,"000")
plotframe([0,1.9,lg,4],tics,[%f,%t],["Energie interne","",""],[0,0.5,0.5,0.5])
plot2d(x,e,4,"000")
plotframe(rect,tics,[%f,%t],["Pression","",""],[0.5,0.5,0.5,0.5])
plot2d(x,p,4,"000")
xpause(100000) ;
end
//
if t>=tmax
 break
end
//
end