///////////////////////////////////////////////////////////
// Copyright G. Allaire, Octobre 2002, Janvier 2012
//
// comparaison de schemas numeriques  
// avec l'equation de Burgers
// u,t + (u^2/2),x = 0
//
///////////////////////////////////////////////////////////
//
xset("font",2,3) ; 
xset("thickness",2) ;
//
pi = %pi ;
//
// lg = longueur du domaine
lg = 5. ;
// nx = nombre de mailles
nx = 200 ;
// dx = pas d'espace
dx = lg/nx ;
// etat gauche et droit
ul = 0.1 ;
ur = 0.5 ;
liminf = min(ul,ur) - 0.2 ;
limmax = max(ul,ur) + 0.2 ;
// condition cfl
cfl = dx/max(abs(ul),abs(ur)) ;
// dt = pas de temps
dt = 0.9*cfl ;
// tt = temps de simulation
if ul > ur 
   tt=0.8*lg/abs(ul+ur) ;
else 
   tt=0.4*lg/max(abs(ul),abs(ur)) ;
end
nt = round(tt/dt) ;
//
// initialisation
//
x=zeros(1,nx) ;
v=zeros(1,nx) ;
f=zeros(1,nx) ;
// uf = solution de Lax-Friedrichs
uf=zeros(1,nx) ;
// uw = solution de Lax-Wendroff
uw=zeros(1,nx) ;
// ug = solution de Godunov
ug=zeros(1,nx) ;
// ue = solution exacte
ue=zeros(1,nx) ;
up=zeros(1,nx) ;
um=zeros(1,nx) ;
// x = points du maillage
// ue = donnee initiale
for i=1:nx
  x(i) = (i-1)*dx  ;
  if x(i) < lg/2
     ue(i) = ul ;
  else
     ue(i) = ur ;
  end
end
uf = ue ;
uw = ue ;
ug = ue ;

clf()
tics=[4,10,4,10];
plotframe([0,liminf,lg,limmax],tics);
plot2d(x,ue,1,"000")
xtitle ('LF, LW, Godunov, cfl=0.9' ,' ',' ') ;

halt() ;
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// marche en temps : 
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// 
for n=1:nt
//
// Lax-Friedrichs
up = shiftn('+1',uf) ;
um = shiftn('-1',uf) ;
v  = (up+um)/2 - (0.5*dt/dx)*(up.^2-um.^2)/2 ;
uf = v ; 
// Lax-Wendroff
up = shiftn('+1',uw) ;
um = shiftn('-1',uw) ;
v  = uw - (0.5*dt/dx)*(up.^2-um.^2)/2 + ... 
(0.5*dt*dt/(dx*dx))*( (uw+up)*0.5.*(up.^2-uw.^2) ...
        - (uw+um)*0.5.*(uw.^2-um.^2) )/2 ;
uw = v ; 
// Godunov
// evaluation du flux a l'interface
f = 0.5*ug.^2 ;
up = shiftn('+1',ug) ;
for i=1:nx
if ug(i) > up(i)
// choc
     sigma = 0.5*(ug(i)+up(i)) ;
     if sigma < 0.
        f(i) = 0.5*up(i)^2 ;
     end
elseif ug(i) < up(i)
// detente
     if up(i) < 0.
        f(i) = 0.5*up(i)^2 ;
     elseif ug(i) > 0.
        f(i) = 0.5*ug(i)^2 ;
     else 
        f(i) = 0. ;
     end
end
end
fm = shiftn('-1',f) ;
v = ug - (dt/dx)*(f-fm) ;
ug = v ; 
if modulo(n,5) == 0
  ue=solbu(ue,nx,dx,dt,lg,ul,ur,n) ;
  clf() 
  plotframe([0,liminf,lg,limmax],tics);
  xset("thickness",2) ;
  plot2d(x,uf,[1,1],"100","Lax-Friedrichs")
  xset("thickness",2) ;
  plot2d(x,uw,[2,3],"100","Lax-Wendroff")
  xset("thickness",2) ;
  plot2d(x,ug,[3,4],"100","Godunov")
  xset("thickness",1) ;
  plot2d(x,ue,[4,6],"100","exacte")
  xtitle('comparaison de 3 schemas, cfl=0.9',' ',' ');
  xpause(100000) ;
end
//
end

halt() ;

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// conditions aux limites periodiques
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// 
lg = 1. ;
// nx = nombre de mailles
nx = 100 ;
// dx = pas d'espace
dx = lg/nx ;
cfl = dx ;
// dt = pas de temps
// (on sait que la vitesse maximale sera 1.)
dt = 0.9*cfl ;
// nt = nombre de pas de temps effectues
nt = 200 ;
//
// initialisation
//
x=zeros(1,nx) ;
v=zeros(1,nx) ;
f=zeros(1,nx) ;
// uf = solution de Lax-Friedrichs
uf=zeros(1,nx) ;
// uw = solution de Lax-Wendroff
uw=zeros(1,nx) ;
// ug = solution de Godunov
ug=zeros(1,nx) ;
up=zeros(1,nx) ;
um=zeros(1,nx) ;
u=zeros(1,nx) ;
v=zeros(1,nx) ;
up=zeros(1,nx) ;
um=zeros(1,nx) ;
sigma=zeros(1,nx) ;
// x = points du maillage
// u = donnee initiale
for i=1:nx
  x(i) = (i-1)*dx  ;
  uf(i) = sin(x(i)*2*pi) ;
end
uw = uf ;
ug = uf ;

// affichage de la donnee initiale
liminf = -1.1 ;
limmax = 1.1 ;
clf()
tics=[4,10,4,10];
plotframe([0,liminf,lg,limmax],tics);
plot2d(x,uf,1,"000")
xtitle ('comparaison LF, LW, Godunov, cfl=0.9' ,' ',' ') ;

halt() ;
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// marche en temps : 
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// 
for n=1:nt
//
// Lax-Friedrichs
up = shiftp('+1',uf) ;
um = shiftp('-1',uf) ;
v = (up+um)/2 - (0.5*dt/dx)*(up.^2-um.^2)/2 ;
uf = v ; 
// Lax-Wendroff
up = shiftp('+1',uw) ;
um = shiftp('-1',uw) ;
v = uw - (0.5*dt/dx)*(up.^2-um.^2)/2 + ... 
(0.5*dt*dt/(dx*dx))*( (uw+up)*0.5.*(up.^2-uw.^2) ...
        - (uw+um)*0.5.*(uw.^2-um.^2) )/2 ;
uw = v ; 
// Godunov
// evaluation du flux a l'interface (i+1/2)
// ug = etat gauche, up = etat droit
f = 0.5*ug.^2 ;
up = shiftp('+1',ug) ;
for i=1:nx
if ug(i) > up(i)
// choc
     sigma = 0.5*(ug(i)+up(i)) ;
     if sigma < 0.
        f(i) = 0.5*up(i)^2 ;
     end
elseif ug(i) < up(i)
// detente
     if up(i) < 0.
        f(i) = 0.5*up(i)^2 ;
     elseif ug(i) > 0.
        f(i) = 0.5*ug(i)^2 ;
     else 
        f(i) = 0. ;
     end
end
end
fm = shiftp('-1',f) ;
v = ug - (dt/dx)*(f-fm) ;
ug = v ; 
if modulo(n,5) == 0
  clf() 
  plotframe([0,liminf,lg,limmax],tics);
  xset("thickness",2) ;
  plot2d(x,uf,[1,1],"100","Lax-Friedrichs")
  xset("thickness",2) ;
  plot2d(x,uw,[2,3],"100","Lax-Wendroff")
  xset("thickness",2) ;
  plot2d(x,ug,[3,4],"100","Godunov")
  xtitle('comparaison de 3 schemas, cfl=0.9',' ',' ');
  xpause(100000) ;
end
//
end