/////////////////////////////////////////////////////////////////////////////
// Paola Goatin, Octobre 2006
//
// schema de Roe pour les equations d'Euler isentropiques en 1-D
// correction entropique ou non
// 
///////////////////////////////////////////////////////////////////////////
//
xset("font",2,3) ;
xset("thickness",2) ;
exec('shiftn.sci');
//
// 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)
sol = zeros(2,nx) ; 
// flux = flux par maille
flux = zeros(2,nx) ;
fluxp = zeros(2,nx) ;
// glux = flux numerique a l'interface
glux = zeros(2,nx) ;
gluxm = zeros(2,nx) ;
absroe = zeros(2,nx) ;
// vp = vecteurs propres de la jacobienne
vp1 = ones(2,nx) ;
vp2 = ones(2,nx) ;
// u = vitess, 
u = zeros(1,nx) ;
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// le tube a choc de Sod
// sans correctgion entropique
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// etats initiaux a gauche (l) et a droite (r)
// rho = densite, u = vitesse
rhol = 1. ;
ul = 0.75 ;
rhor = 0.125 ;
ur = 0. ;
// coeff de la pression p(rho) = a^2 rho, gama=1
a=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;
  else
     sol(1,i) = rhor ;     
     sol(2,i) = rhor*ur ;
  end
end
r = sol(1,:) ;
u = sol(2,:)./sol(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,1]);
plot2d(x,r,2,"000")
plotframe([0,liminf,lg,1.5],tics,[%f,%t],["Vitesse","",""],[0.5,0,0.5,1])
plot2d(x,u,3,"000")

   f=get("current_figure") ;//get the handle of the current figure : 
                                           //if none exists, create a figure and return the corresponding handle
  // f.figure_position
   f.figure_size= [596,240]; 
   f.figure_name="Euler 2x2 (Roe)";
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 a partir des variables conservatives
if sign(sol(1,:)) > 0
   u = sol(2,:)./sol(1,:) ;
else 
  'densite negative ou nulle'
  return
end
 
// moyennes de Roe
        
w = sqrt(sol(1,:));
wp = shiftn('+1',w) ;
up = shiftn('+1',u) ; 
solp = shiftn('+1',sol) ;
utild = (w.*u + wp.*up)./(w + wp);
rhotild = w.*wp ;

// c2son = vitesse du son au carre
c2son =  a^2 ;

// cson = vitesse du son
if sign(a) >= 0 
   cson = a ;
else 
   'vitesse du son negative'
   return
end

// vecteurs propres de la jacobienne

vp1(2,:) = utild - cson ;
vp2(2,:) = utild + cson ;

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

// // ce sont les coefficients alphai dans ur-ul = alpha1 vp1 + alpha2 vp2
// // alpha1 = 
// // alpha2 = 


// calcul du flux

flux(1,:) =     sol(2,:)      ;
flux(2,:) = u.* sol(2,:) + c2son.* sol(1,:) ;
fluxp = shiftn('+1',flux) ;

// terme de viscosite numerique dans le flux numerique

// // absroe =

// 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,:) ;
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,1]);
plot2d(x,r,2,"000")
plotframe([0,liminf,lg,1.5],tics,[%f,%t],["Vitesse","",""],[0.5,0,0.5,1])
plot2d(x,u,3,"000")
xpause(100000) ;

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

halt() ;


////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// tube a choc 
// 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
rhol = 1. ;
ul = 0.75 ;
rhor = 0.125 ;
ur = 0. ;
//
// 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;
  else
     sol(1,i) = rhor ;     
     sol(2,i) = rhor*ur ;
  end
end
r = sol(1,:) ;
u = sol(2,:)./sol(1,:) ;
    
// trace des conditions initiales
liminf = -0.1 ;
limmax = 1.1 ;
xinit

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

   f=get("current_figure") ;//get the handle of the current figure : 
                                           //if none exists, create a figure and return the corresponding handle
  // f.figure_position
   f.figure_size= [596,240]; 
   f.figure_name="Euler 2x2 (Roe avec correction entropique)";

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

// c2son = vitesse du son au carre
c2son = a^2 ;

// cson = vitesse du son
if sign(a) >= 0 
   cson = a ;
else 
   'vitesse du son negative'
   return
end

// vecteurs propres de la jacobienne

vp1(2,:) = utild - cson ;
vp2(2,:) = utild + cson ;

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

// // ce sont les coefficients alphai dans ur-ul = alpha1 vp1 + alpha2 vp2
// // alpha1 = 
// // alpha2 = 

// calcul du flux

flux(1,:) =     sol(2,:)      ;
flux(2,:) = u.* sol(2,:) + c2son.* sol(1,:) ;
fluxp = shiftn('+1',flux) ;

// um=solm état intermeriaire entre ul=sol et ur=solp

solm = zeros(2,nx) ; 
// // solm(1,:) = 
// // solm(2,:) = 

// valeurs propres de Df (matrice Jacobienne du problème exacte)

// // lambda1=lambda1(ul)
// // lambda1 = 
lambda1p = shiftn('+1',lambda1) ;  // lambda1(ur)
// // lambda1m=lambda1(um)
// // lambda1m = 
// // lambda1=lambda1(ul)
// // lambda2 = 
lambda2p = shiftn('+1',lambda2) ; // lambda1(ur)
// // lambda2m=lambda2(um)
// // lambda2m = 

// correction des vitesses

lambda1tildl = zeros(1,nx) ;
lambda1tildr = zeros(1,nx) ;
lambda2tildl = zeros(1,nx) ;
lambda2tildr = zeros(1,nx) ;

for i=1:nx
  if max(lambda1m(i),0) - min(lambda1(i),0) > 10^(-10)
    lambda1tildl(i) = min(lambda1(i),0).* (max(lambda1m(i),0) - (utild(i)-cson)) ./ (max(lambda1m(i),0) - min(lambda1(i),0)) ;
    lambda1tildr(i) = max(lambda1m(i),0).* ((utild(i)-cson) - min(lambda1(i),0)) ./ (max(lambda1m(i),0) - min(lambda1(i),0)) ;
  else
//    // cas d'un choc sonique, on garde les valeurs propres du problème linéarisé
// //     lambda1tildl(i) = 
// //     lambda1tildr(i) = 
  end

  if max(lambda2p(i),0) - min(lambda2m(i),0) > 10^(-10)
    lambda2tildl(i) = min(lambda2m(i),0).* (max(lambda2p(i),0) - (utild(i)+cson)) ./ (max(lambda2p(i),0) - min(lambda2m(i),0)) ;
    lambda2tildr(i) = max(lambda2p(i),0).* ((utild(i)+cson) - min(lambda2m(i),0)) ./ (max(lambda2p(i),0) - min(lambda2m(i),0)) ;
  else
//    // cas d'un choc sonique, on garde les valeurs propres du problème linéarisé  
// //    lambda2tildl(i) =
// //    lambda2tildr(i) = 
  end
end    

// terme de viscosite numerique dans le flux numerique

absroe(1,:) = (lambda1tildl-lambda1tildr).*alpha1.*vp1(1,:) ;
absroe(2,:) = (lambda1tildl-lambda1tildr).*alpha1.*vp1(2,:) ;

absroe(1,:) = absroe(1,:) + (lambda2tildl-lambda2tildr).*alpha2.*vp2(1,:) ;
absroe(2,:) = absroe(2,:) + (lambda2tildl-lambda2tildr).*alpha2.*vp2(2,:) ;

// 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,:) ;
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,1]);
plot2d(x,r,2,"000")
plotframe([0,liminf,lg,1.5],tics,[%f,%t],["Vitesse","",""],[0.5,0,0.5,1])
plot2d(x,u,3,"000")
xpause(100000) ;
end
//
if t>=tmax
 break
end
//
end