/////////////////////////////////////////////////////////////////
// Copyright G. Allaire, Octobre 2000, Janvier 2012
//
// schemas numeriques pour le probleme de Riemann 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 = 1. ;
// nx = nombre de mailles
nx = 100 ;
// dx = pas d'espace
dx = lg/nx ;
// etat gauche et droit
ul = 0.8 ;
ur = 0.2 ;
liminf = min(ul,ur) - 0.1 ;
limmax = max(ul,ur) + 0.1 ;
// 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=1./abs(ul+ur) ;
else 
   tt=0.5/max(abs(ul),abs(ur)) ;
end
tt=0.9*tt ;
nt = round(tt/dt) ;
//
// initialisation
//
x=zeros(1,nx) ;
u=zeros(1,nx) ;
v=zeros(1,nx) ;
w=zeros(1,nx) ;
up=zeros(1,nx) ;
um=zeros(1,nx) ;
// x = points du maillage
// u = donnee initiale
for i=1:nx
  x(i) = (i-1)*dx  ;
  if x(i) < lg/2
     u(i) = ul ;
  else
     u(i) = ur ;
  end
end
w = u ;

clf()
tics=[4,10,4,10];
plotframe([0,liminf,lg,limmax],tics);
plot2d(x,u,1,"000")
xtitle ('Lax-Friedrichs, cfl=0.9' ,' ',' ') ;

halt() ;
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// marche en temps : Lax-Friedrichs
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// 
for n=1:nt
//
up = shiftn('+1',u) ;
um = shiftn('-1',u) ;
v  = (up+um)/2 - (0.5*dt/dx)*(up.^2-um.^2)/2 ;
u = v ; 
if modulo(n,5) == 0
  w=solbu(w,nx,dx,dt,lg,ul,ur,n) ;
  clf() 
  plotframe([0,liminf,lg,limmax],tics);
  xset("thickness",2) ;
  plot2d(x,u,[1,1],"100","Lax-Friedrichs")
  xset("thickness",1) ;
  plot2d(x,w,[2,4],"100","exact")
  xtitle('Lax-Friedrichs, cfl=0.9',' ',' ');
  xpause(100000) ;
end
//
end

halt() ;

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// nouveau schema: Lax-Friedrichs avec cfl=0.3
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//
dt = 0.3*cfl ;
nt = round(tt/dt) ;
//
// initialisation
//
x=zeros(1,nx) ;
u=zeros(1,nx) ;
v=zeros(1,nx) ;
w=zeros(1,nx) ;
up=zeros(1,nx) ;
um=zeros(1,nx) ;
for i=1:nx
  x(i) = (i-1)*dx  ;
  if x(i) < lg/2
     u(i) = ul ;
  else
     u(i) = ur ;
  end
end
w = u ;

clf()
tics=[4,10,4,10];
plotframe([0,liminf,lg,limmax],tics);
plot2d(x,u,1,"000")
xtitle ('Lax-Friedrichs, cfl=0.3' ,' ',' ') ;

xpause(800000) ;
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// marche en temps : Lax-Friedrichs
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// 
for n=1:nt
//
up = shiftn('+1',u) ;
um = shiftn('-1',u) ;
v  = (up+um)/2 - (0.5*dt/dx)*(up.^2-um.^2)/2 ;
u = v ; 
if modulo(n,5) == 0
  w=solbu(w,nx,dx,dt,lg,ul,ur,n) ;
  clf() 
  plotframe([0,liminf,lg,limmax],tics);
  xset("thickness",2) ;
  plot2d(x,u,[1,1],"100","Lax-Friedrichs")
  xset("thickness",1) ;
  plot2d(x,w,[2,4],"100","exact")
  xtitle('Lax-Friedrichs, cfl=0.3',' ',' ');
  xpause(100000) ;
end
//
end

halt() ;

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// nouveau schema: centre explicite avec cfl=0.1 (explose)
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//
dt = 0.1*cfl ;
nt = round(tt/dt) ;
//
// initialisation
//
x=zeros(1,nx) ;
u=zeros(1,nx) ;
v=zeros(1,nx) ;
w=zeros(1,nx) ;
for i=1:nx
  x(i) = (i-1)*dx  ;
  if x(i) < lg/2
     u(i) = ul ;
  else
     u(i) = ur ;
  end
end
w = u ;

clf()
tics=[4,10,4,10];
plotframe([0,liminf-0.2,lg,limmax+0.2],tics);
plot2d(x,u,1,"000")
xtitle ('centre explicite, cfl=0.1' ,' ',' ') ;

xpause(800000) ;
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// marche en temps : centre explicite
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// 
for n=1:nt
//
up = shiftn('+1',u) ;
um = shiftn('-1',u) ;
v  = u - (0.5*dt/dx)*(up.^2-um.^2)/2 ;
u = v ; 
if modulo(n,5) == 0
  w=solbu(w,nx,dx,dt,lg,ul,ur,n) ;
  clf() 
  plotframe([0,liminf-0.2,lg,limmax+0.2],tics);
  xset("thickness",2) ;
  plot2d(x,u,[1,1],"100","centre explicite")
  xset("thickness",1) ;
  plot2d(x,w,[2,4],"100","exact")
  xtitle('centre explicite, cfl=0.1',' ',' ');
  xpause(100000) ;
end
//
end

halt() ;

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// nouveau schema: Lax-Wendroff avec cfl=0.9
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//
dt = 0.9*cfl ;
nt = round(tt/dt) ;
//
// initialisation
//
x=zeros(1,nx) ;
u=zeros(1,nx) ;
v=zeros(1,nx) ;
w=zeros(1,nx) ;
for i=1:nx
  x(i) = (i-1)*dx  ;
  if x(i) < lg/2
     u(i) = ul ;
  else
     u(i) = ur ;
  end
end
w = u ;

clf()
tics=[4,10,4,10];
plotframe([0,liminf,lg,limmax],tics);
plot2d(x,u,1,"000")
xtitle ('Lax-Wendroff, cfl=0.9' ,' ',' ') ;

xpause(800000) ;
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// marche en temps : Lax-Wendroff
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// 
for n=1:nt
//
up = shiftn('+1',u) ;
um = shiftn('-1',u) ;
v  = u - (0.5*dt/dx)*(up.^2-um.^2)/2 + ... 
(0.5*dt*dt/(dx*dx))*( (u+up)*0.5.*(up.^2-u.^2) ...
        - (u+um)*0.5.*(u.^2-um.^2) )/2 ;
u = v ; 
if modulo(n,5) == 0
  w=solbu(w,nx,dx,dt,lg,ul,ur,n) ;
  clf() 
  plotframe([0,liminf,lg,limmax],tics);
  xset("thickness",2) ;
  plot2d(x,u,[1,1],"100","Lax-Wendroff")
  xset("thickness",1) ;
  plot2d(x,w,[2,4],"100","exact")
  xtitle('Lax-Wendroff, cfl=0.9',' ',' ');
  xpause(100000) ;
end
//
end

halt() ;

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// nouveau schema: Lax-Wendroff avec cfl=0.3
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//
dt = 0.3*cfl ;
nt = round(tt/dt) ;
//
// initialisation
//
x=zeros(1,nx) ;
u=zeros(1,nx) ;
v=zeros(1,nx) ;
w=zeros(1,nx) ;
for i=1:nx
  x(i) = (i-1)*dx  ;
  if x(i) < lg/2
     u(i) = ul ;
  else
     u(i) = ur ;
  end
end
w = u ;

clf()
tics=[4,10,4,10];
plotframe([0,liminf,lg,limmax],tics);
plot2d(x,u,1,"000")
xtitle ('Lax-Wendroff, cfl=0.3' ,' ',' ') ;

xpause(800000) ;
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// marche en temps : Lax-Wendroff
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// 
for n=1:nt
//
up = shiftn('+1',u) ;
um = shiftn('-1',u) ;
v  = u - (0.5*dt/dx)*(up.^2-um.^2)/2 + ... 
(0.5*dt*dt/(dx*dx))*( 0.5*(u+up).*(up.^2-u.^2) ...
        - 0.5*(u+um).*(u.^2-um.^2) )/2 ;
u = v ; 
if modulo(n,5) == 0
  w=solbu(w,nx,dx,dt,lg,ul,ur,n) ;
  clf() 
  plotframe([0,liminf,lg,limmax],tics);
  xset("thickness",2) ;
  plot2d(x,u,[1,1],"100","Lax-Wendroff")
  xset("thickness",1) ;
  plot2d(x,w,[2,4],"100","exact")
  xtitle('Lax-Wendroff, cfl=0.3',' ',' ');
  xpause(100000) ;
end
//
end

halt() ;

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// nouveau schema: non conservatif avec cfl=0.9 (fausse solution)
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//
dt = 0.9*cfl ;
nt = round(tt/dt) ;
//
// initialisation
//
x=zeros(1,nx) ;
u=zeros(1,nx) ;
v=zeros(1,nx) ;
w=zeros(1,nx) ;
for i=1:nx
  x(i) = (i-1)*dx  ;
  if x(i) < lg/2
     u(i) = ul ;
  else
     u(i) = ur ;
  end
end
w = u ;

clf()
tics=[4,10,4,10];
plotframe([0,liminf,lg,limmax],tics);
plot2d(x,u,1,"000")
xtitle ('non conservatif, cfl=0.9' ,' ',' ') ;

xpause(800000) ;
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// marche en temps : non conservatif
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// 
for n=1:nt
//
up = shiftn('+1',u) ;
um = shiftn('-1',u) ;
v  = (up+um)/2 - (0.5*dt/dx)*um.*(up-um) ;
u = v ; 
if modulo(n,1) == 0
  w=solbu(w,nx,dx,dt,lg,ul,ur,n) ;
  clf() 
  plotframe([0,liminf,lg,limmax],tics);
  xset("thickness",2) ;
  plot2d(x,u,[1,1],"100","non conservatif")
  xset("thickness",1) ;
  plot2d(x,w,[2,4],"100","exact")
  xtitle('non conservatif, cfl=0.9',' ',' ');
  xpause(100000) ;
end
//
end

halt() ;

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// nouveau schema: Godunov avec cfl=0.9 
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//
dt = 0.9*cfl ;
nt = round(tt/dt) ;
//
// initialisation
//
x=zeros(1,nx) ;
u=zeros(1,nx) ;
v=zeros(1,nx) ;
w=zeros(1,nx) ;
f=zeros(1,nx) ;
sigma=zeros(1,nx) ;
for i=1:nx
  x(i) = (i-1)*dx  ;
  if x(i) < lg/2
     u(i) = ul ;
  else
     u(i) = ur ;
  end
end
w = u ;

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

xpause(800000) ;
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// marche en temps : Godunov
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// 
for n=1:nt
//
// evaluation du flux a l'interface
f = 0.5*u.^2 ;
up = shiftn('+1',u) ;
for i=1:nx
if u(i) > up(i)
// choc
     sigma = 0.5*(u(i)+up(i)) ;
     if sigma < 0.
        f(i) = 0.5*up(i)^2 ;
     end
elseif u(i) < up(i)
// detente
     if up(i) < 0.
        f(i) = 0.5*up(i)^2 ;
     elseif u(i) > 0.
        f(i) = 0.5*u(i)^2 ;
     else 
        f(i) = 0. ;
     end
end
end
fm = shiftn('-1',f) ;
v = u - (dt/dx)*(f-fm) ;
u = v ; 
if modulo(n,5) == 0
  w=solbu(w,nx,dx,dt,lg,ul,ur,n) ;
  clf() 
  plotframe([0,liminf,lg,limmax],tics);
  xset("thickness",2) ;
  plot2d(x,u,[1,1],"100","Godunov")
  xset("thickness",1) ;
  plot2d(x,w,[2,4],"100","exact")
  xtitle('Godunov, cfl=0.9',' ',' ');
  xpause(100000) ;
end
//
end