//////////////////////////////////////////////////////////////////////////
// Copyright G. Allaire, Octobre 2002, Janvier 2012
//
// probleme de Riemann pour les equations d'Euler
// 
//////////////////////////////////////////////////////////////////////////
//
xset("font",2,3) ; 
xset("thickness",2) ;
//
pi = %pi ;
//
// 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

liminf = -0.1 ;
limmax = 1.1 ;

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// 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,fp,fm) ;
fpp = shiftn('+1',fp) ;
fpm = shiftn('-1',fp) ;
fmp = shiftn('+1',fm) ;
fmm = shiftn('-1',fm) ;

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 modulo(n,8) == 0
clf()
tics=[4,10,4,10];
rect = [0,liminf,lg,limmax];
plotframe([0,liminf,lg,limmax],tics,[%f,%t],...
["Van Leer: 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,5,"000")
xpause(100000) ;
end
//
end


halt() ;


// 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

liminf = 0. ;
limmax = 1. ;

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// 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,fp,fm) ;
fpp = shiftn('+1',fp) ;
fpm = shiftn('-1',fp) ;
fmp = shiftn('+1',fm) ;
fmm = shiftn('-1',fm) ;

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 modulo(n,5) == 0
clf()
tics=[4,10,4,10];
rect = [0,liminf,lg,limmax];
plotframe([0,liminf,lg,1.1],tics,[%f,%t],...
["Van Leer: Densité","",""],[0,0,0.5,0.5]);
plot2d(x,r,2,"000")
plotframe([0,-vitesse-0.1,lg,vitesse+0.1],tics,[%f,%t],["Vitesse","",""],[0.5,0,0.5,0.5])
plot2d(x,u,3,"000")
plotframe([0,0.4,lg,1.1],tics,[%f,%t],... 
["Energie interne","",""],[0,0.5,0.5,0.5])
plot2d(x,e,4,"000")
plotframe([0,0.,lg,0.5],tics,[%f,%t],["Pression","",""],[0.5,0.5,0.5,0.5])
plot2d(x,p,5,"000")
xpause(100000) ;
end
//
end


halt() ;

// 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

liminf = -0.1 ;
limmax = 1.1 ;

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// 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,fp,fm) ;
fpp = shiftn('+1',fp) ;
fpm = shiftn('-1',fp) ;
fmp = shiftn('+1',fm) ;
fmm = shiftn('-1',fm) ;

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 modulo(n,5) == 0
clf()
tics=[4,10,4,10];
rect = [0,liminf,lg,limmax];
plotframe([0,liminf,lg,6.],tics,[%f,%t],...
["Van Leer: Densité","",""],[0,0,0.5,0.5]);
plot2d(x,r,2,"000")
plotframe([0,-7.,lg,0.5],tics,[%f,%t],["Vitesse","",""],[0.5,0,0.5,0.5])
plot2d(x,u,3,"000")
plotframe([0,-1.,lg,270],tics,[%f,%t],...
["Energie interne","",""],[0,0.5,0.5,0.5])
plot2d(x,e,4,"000")
plotframe([0,-1.,lg,110.],tics,[%f,%t],["Pression","",""],[0.5,0.5,0.5,0.5])
plot2d(x,p,5,"000")
xpause(100000) ;
end
//
end