//////////////////////////////////////////////////////////////////////////
// Copyright G. Allaire, Decembre 2002, Janvier 2012
//
// tube a choc de Sod pour les equations d'Euler
// etude de divers schemas simples
//
///////////////////////////////////////////////////////////////////////////
//
xset("font",2,3) ;
xset("thickness",2) ;
//
// 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.) ;
// limitation du pas de temps par
// la condition CFL
cfl = 0.9/(1. + sqrt(gam)) ;
dt = dx*cfl ;
// nt = nombre de pas de temps
nt = 400 ;
//
// initialisation
//
x=zeros(1,nx+1) ;
r=zeros(1,nx+1) ;
u=zeros(1,nx+1) ;
e=zeros(1,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
// p = pression
p = (gam-1.)*r.*e ;

// 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],["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 Lax-Friedrichs
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//
// 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
//
rp = shiftn('+1',r) ;
rm = shiftn('-1',r) ;
up = shiftn('+1',u) ;
um = shiftn('-1',u) ;
ep = shiftn('+1',e) ;
em = shiftn('-1',e) ;
pp = (gam-1.)*rp.*ep ;
pm = (gam-1.)*rm.*em ;
ip = rp.*(up.^2/2.+ep) ;
im = rm.*(um.^2/2.+em) ;
rn = (rp+rm)/2 - (0.5*dt/dx)*(rp.*up-rm.*um) ;
un = (rp.*up+rm.*um)/2 - (0.5*dt/dx)*...
     (rp.*up.^2+pp-rm.*um.^2-pm) ;
un = un./rn ;
en = (ip+im)/2 - (0.5*dt/dx)*... 
     ((ip+pp).*up-(im+pm).*um) ;
en = en./rn - un.^2/2. ;
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) ;
r = rn ;
u = un ; 
e = en ;
p = (gam-1.)*r.*e ;

if modulo(n,5) == 0
clf()
tics=[4,10,4,10];
rect = [0,liminf,lg,limmax];
plotframe([0,liminf,lg,limmax],tics,[%f,%t],...
["Lax-Friedrichs: 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() ;

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


////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// marche en temps : Van Leer
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//
for n=1:(5*nt/4)
//
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,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() ;

//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// schema de Richtmyer (Lax-Wendroff en deux etapes)
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//
// initialisation
//
cfl = 0.9/(1. + sqrt(gam)) ;
dt = dx*cfl ;
x=zeros(1,nx+1) ;
r=zeros(1,nx+1) ;
u=zeros(1,nx+1) ;
e=zeros(1,nx+1) ;
// f = flux en i, fp = flux en i+1, fm = flux en i-1
f=zeros(3,nx+1) ;
fp=zeros(3,nx+1) ;
fm=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


////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// marche en temps : Richtmyer
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//
for n=1:nt
//
rp = shiftn('+1',r) ;
up = shiftn('+1',u) ;
ep = shiftn('+1',e) ;
f = flux(gam,r,u,e,f) ;
fp = shiftn('+1',f) ;
// etat intermediaire (n+1/2)
rs = 0.5*(rp+r) - (0.5*dt/dx)*(fp(1,:)-f(1,:)) ;

us = 0.5*(rp.*up+r.*u) - (0.5*dt/dx)*(fp(2,:)-f(2,:)) ;
us = us./rs ;

es = 0.5*(rp.*(up.^2/2.+ep)+r.*(u.^2/2.+e)) - ...
     (0.5*dt/dx)*(fp(3,:)-f(3,:)) ;
es = es./rs - us.^2/2. ;
// etat final (n+1)
f = flux(gam,rs,us,es,f) ;
fm = shiftn('-1',f) ;
rn = r - (dt/dx)*(f(1,:)-fm(1,:)) ;

un = r.*u - (dt/dx)*(f(2,:)-fm(2,:)) ;
un = un./rn ;

en = r.*(u.^2/2.+e) - (dt/dx)*(f(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,limmax],tics,[%f,%t],...
["Richtmyer: 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() ;


//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// schema de Lax-Wendroff
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//
// initialisation
//
cc =1. ;
cfl = 0.9/(1. + sqrt(gam)) ;
dt = dx*cfl ;
x=zeros(1,nx+1) ;
r=zeros(1,nx+1) ;
u=zeros(1,nx+1) ;
e=zeros(1,nx+1) ;
// f = flux en i, fp = flux en i+1, fm = flux en i-1
f=zeros(3,nx+1) ;
fp=zeros(3,nx+1) ;
fm=zeros(3,nx+1) ;
gp=zeros(3,nx+1) ;
gm=zeros(3,nx+1) ;
// a = matrice jacobienne
a2=zeros(3,nx+1) ;
a3=zeros(3,nx+1) ;
for i=1:nx+1
  x(i) = (i-1)*dx  ;
  if x(i) < lg/2.01
     r(i) = rl ;
     u(i) = ul ;
     e(i) = el ;
  elseif x(i) > lg/1.99
     r(i) = rr ;
     u(i) = ur ;
     e(i) = er ;
  else
     r(i) = (rl+rr)/2 + 0.5*(rl-rr)*cos((x(i)-lg/2.01)*%pi/(lg/1.99-lg/2.01)) ;
     u(i) = (ul+ur)/2 + 0.5*(ul-ur)*cos((x(i)-lg/2.01)*%pi/(lg/1.99-lg/2.01)) ;
     e(i) = (el+er)/2 + 0.5*(el-er)*cos((x(i)-lg/2.01)*%pi/(lg/1.99-lg/2.01)) ;
  end
end
// p = pression
p = (gam-1.)*r.*e ;

// trace des conditions initiales
// il faut les regulariser sinon ca explose...
clf()
tics=[4,10,4,10];
rect = [0,liminf,lg,limmax];
plotframe([0,liminf,lg,limmax],tics,[%f,%t],...
["Lax-Wendroff: 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) ;

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// marche en temps : Lax-Wendroff
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//
for n=1:nt
//
f = flux(gam,r,u,e,f) ;
rp = shiftn('+1',r) ;
up = shiftn('+1',u) ;
ep = shiftn('+1',e) ;
fp = shiftn('+1',f) ;
fm = shiftn('-1',f) ;
// etat moyen (i+1/2)
rmoy = (r+rp)/2. ;
umoy = (r.*u+rp.*up)./(r+rp) ;
//umoy = (u+up)/2. ;
emoy = ( r.*(0.5*u.^2+e) + rp.*(0.5*up.^2+ep) )./...
       (r+rp) - 0.5*umoy.^2 ;
//emoy = (e+ep)/2. ;
// jacobien a l'etat moyen
a2(1,:) = 0.5*(gam-3.)*umoy.^2 ;
a2(2,:) = -(gam-3.)*umoy ;
a3(1,:) = (gam-1.)*umoy.^3 - gam*umoy.*(emoy+0.5*umoy.^2) ;
a3(2,:) = gam*(emoy+0.5*umoy.^2) - 1.5*(gam-1.)*umoy.^2 ;
a3(3,:) = gam*umoy ;
df = fp-f ;

gp(1,:) = df(2,:) ;
gp(2,:) = a2(1,:).*df(1,:) + a2(2,:).*df(2,:) + (gam-1)*df(3,:) ;
gp(3,:) = a3(1,:).*df(1,:) + a3(2,:).*df(2,:) + a3(3,:).*df(3,:) ;
gm = shiftn('-1',gp) ;

rn = r - (0.5*dt/dx)*(fp(1,:)-fm(1,:))...
       + cc*(0.5*dt*dt/(dx*dx))*(gp(1,:)-gm(1,:)) ;

un = r.*u - (0.5*dt/dx)*(fp(2,:)-fm(2,:))...
          + cc*(0.5*dt*dt/(dx*dx))*(gp(2,:)-gm(2,:)) ;
un = un./rn ;

en = r.*(u.^2/2.+e) - (0.5*dt/dx)*(fp(3,:)-fm(3,:))...
     + cc*(0.5*dt*dt/(dx*dx))*(gp(3,:)-gm(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,limmax],tics,[%f,%t],...
["Lax-Wendroff: 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