////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// Copyright G. Allaire, Octobre 2002, Janvier 2012
//
// raffinement de maillage  
// 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
nx1 = 50 ;
nx2 = 100 ; 
nx3 = 200 ; 
nx4 = 500 ;
// dx = pas d'espace
dx1 = lg/nx1 ;
dx2 = lg/nx2 ;
dx3 = lg/nx3 ;
dx4 = lg/nx4 ;
// dt = pas de temps
// (on sait que la vitesse maximale sera 1.)
dt1 = 0.9*dx1 ;
dt2 = 0.9*dx2 ;
dt3 = 0.9*dx3 ;
dt4 = 0.9*dx4 ;
// nombre de pas de temps
nt1 = 100 ;
nt2 = nt1*nx2/nx1 ;
nt3 = nt1*nx3/nx1 ;
nt4 = nt1*nx4/nx1 ;
//
// initialisation
//
x1=zeros(1,nx1) ;
x2=zeros(1,nx2) ;
x3=zeros(1,nx3) ;
x4=zeros(1,nx4) ;
v1=zeros(1,nx1) ;
v2=zeros(1,nx2) ;
v3=zeros(1,nx3) ;
v4=zeros(1,nx4) ;
// uf = solution de Lax-Friedrichs
uf1=zeros(1,nx1) ;
uf2=zeros(1,nx2) ;
uf3=zeros(1,nx3) ;
uf4=zeros(1,nx4) ;
// uw = solution de Lax-Wendroff
uw1=zeros(1,nx1) ;
uw2=zeros(1,nx2) ;
uw3=zeros(1,nx3) ;
uw4=zeros(1,nx4) ;
// ug = solution de Lax-Friedrichs
ug1=zeros(1,nx1) ;
ug2=zeros(1,nx2) ;
ug3=zeros(1,nx3) ;
ug4=zeros(1,nx4) ;
// f = flux
f1=zeros(1,nx1) ;
f2=zeros(1,nx2) ;
f3=zeros(1,nx3) ;
f4=zeros(1,nx4) ;
// up, um = solution uf decalee
up1=zeros(1,nx1) ;
up2=zeros(1,nx2) ;
up3=zeros(1,nx3) ;
up4=zeros(1,nx4) ;
um1=zeros(1,nx1) ;
um2=zeros(1,nx2) ;
um3=zeros(1,nx3) ;
um4=zeros(1,nx4) ;
// x = points du maillage
// ue = donnee initiale
for i=1:nx1
  x1(i) = (i-1)*dx1  ;
  uf1(i) = sin(x1(i)*2*pi) ;
end
for i=1:nx2
  x2(i) = (i-1)*dx2  ;
  uf2(i) = sin(x2(i)*2*pi) ;
end
for i=1:nx3
  x3(i) = (i-1)*dx3  ;
  uf3(i) = sin(x3(i)*2*pi) ;
end
for i=1:nx4
  x4(i) = (i-1)*dx4  ;
  uf4(i) = sin(x4(i)*2*pi) ;
end

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

halt() ;
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// marche en temps : Lax-Friedrichs
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// 
for n=1:nt1
//
// Lax-Friedrichs
up1 = shiftp('+1',uf1) ;
um1 = shiftp('-1',uf1) ;
v1 = (up1+um1)/2 - (0.5*dt1/dx1)*(up1.^2-um1.^2)/2 ;
uf1 = v1 ; 

for n2=1:(nx2/nx1)
//
// Lax-Friedrichs
up2 = shiftp('+1',uf2) ;
um2 = shiftp('-1',uf2) ;
v2 = (up2+um2)/2 - (0.5*dt2/dx2)*(up2.^2-um2.^2)/2 ;
uf2 = v2 ; 
end

for n3=1:(nx3/nx1)
//
// Lax-Friedrichs
up3 = shiftp('+1',uf3) ;
um3 = shiftp('-1',uf3) ;
v3 = (up3+um3)/2 - (0.5*dt3/dx3)*(up3.^2-um3.^2)/2 ;
uf3 = v3 ; 
end

for n4=1:(nx4/nx1)
//
// Lax-Friedrichs
up4 = shiftp('+1',uf4) ;
um4 = shiftp('-1',uf4) ;
v4 = (up4+um4)/2 - (0.5*dt4/dx4)*(up4.^2-um4.^2)/2 ;
uf4 = v4 ; 
end

  clf() 
  plotframe([0,liminf,lg,limmax],tics);
  xset("thickness",2) ;
  plot2d(x1,uf1,[1,1],"100","50 mailles")
  xset("thickness",2) ;
  plot2d(x2,uf2,[2,3],"100","100 mailles")
  xset("thickness",2) ;
  plot2d(x3,uf3,[3,4],"100","200 mailles")
  xset("thickness",2) ;
  plot2d(x4,uf4,[4,6],"100","500 mailles")
  xtitle('comparaison de 4 maillages pour LF, cfl=0.9',' ',' ');
  xpause(100000) ;

//
end

halt() ;

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// Lax Wendroff
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// 
// x = points du maillage
// ue = donnee initiale
for i=1:nx1
  x1(i) = (i-1)*dx1  ;
  uw1(i) = sin(x1(i)*2*pi) ;
end
for i=1:nx2
  x2(i) = (i-1)*dx2  ;
  uw2(i) = sin(x2(i)*2*pi) ;
end
for i=1:nx3
  x3(i) = (i-1)*dx3  ;
  uw3(i) = sin(x3(i)*2*pi) ;
end
for i=1:nx4
  x4(i) = (i-1)*dx4  ;
  uw4(i) = sin(x4(i)*2*pi) ;
end

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

halt() ;

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// marche en temps : Lax-Wendroff
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// 

for n=1:nt1
//
// Lax-Wendroff
up1 = shiftp('+1',uw1) ;
um1 = shiftp('-1',uw1) ;
v1 = uw1 - (0.5*dt1/dx1)*(up1.^2-um1.^2)/2 + ... 
(0.5*dt1*dt1/(dx1*dx1))*( (uw1+up1)*0.5.*(up1.^2-uw1.^2) ...
        - (uw1+um1)*0.5.*(uw1.^2-um1.^2) )/2 ;
uw1 = v1 ; 

for n2=1:(nx2/nx1)
//
// Lax-Wendroff
up2 = shiftp('+1',uw2) ;
um2 = shiftp('-1',uw2) ;
v2 = uw2 - (0.5*dt2/dx2)*(up2.^2-um2.^2)/2 + ... 
(0.5*dt2*dt2/(dx2*dx2))*( (uw2+up2)*0.5.*(up2.^2-uw2.^2) ...
        - (uw2+um2)*0.5.*(uw2.^2-um2.^2) )/2 ;
uw2 = v2 ; 
end

for n3=1:(nx3/nx1)
//
// Lax-Wendroff
up3 = shiftp('+1',uw3) ;
um3 = shiftp('-1',uw3) ;
v3 = uw3 - (0.5*dt3/dx3)*(up3.^2-um3.^2)/2 + ... 
(0.5*dt3*dt3/(dx3*dx3))*( (uw3+up3)*0.5.*(up3.^2-uw3.^2) ...
        - (uw3+um3)*0.5.*(uw3.^2-um3.^2) )/2 ;
uw3 = v3 ; 
end

for n4=1:(nx4/nx1)
//
// Lax-Wendroff
up4 = shiftp('+1',uw4) ;
um4 = shiftp('-1',uw4) ;
v4 = uw4 - (0.5*dt4/dx4)*(up4.^2-um4.^2)/2 + ... 
(0.5*dt4*dt4/(dx4*dx4))*( (uw4+up4)*0.5.*(up4.^2-uw4.^2) ...
        - (uw4+um4)*0.5.*(uw4.^2-um4.^2) )/2 ;
uw4 = v4 ; 
end

  clf() 
  plotframe([0,liminf,lg,limmax],tics);
  xset("thickness",2) ;
  plot2d(x1,uw1,[1,1],"100","50 mailles")
  xset("thickness",1) ;
  plot2d(x2,uw2,[2,3],"100","100 mailles")
  xset("thickness",1) ;
  plot2d(x3,uw3,[3,4],"100","200 mailles")
  xset("thickness",1) ;
  plot2d(x4,uw4,[4,6],"100","500 mailles")
  xtitle('comparaison de 4 maillages pour LW, cfl=0.9',' ',' ');
  xpause(100000) ;

//
end

halt() ;

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// Godunov
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// 
// 
// x = points du maillage
// ue = donnee initiale
for i=1:nx1
  x1(i) = (i-1)*dx1  ;
  ug1(i) = sin(x1(i)*2*pi) ;
end
for i=1:nx2
  x2(i) = (i-1)*dx2  ;
  ug2(i) = sin(x2(i)*2*pi) ;
end
for i=1:nx3
  x3(i) = (i-1)*dx3  ;
  ug3(i) = sin(x3(i)*2*pi) ;
end
for i=1:nx4
  x4(i) = (i-1)*dx4  ;
  ug4(i) = sin(x4(i)*2*pi) ;
end

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

halt() ;

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// marche en temps : Godunov
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// 

for n=1:nt1
//
// Godunov
// evaluation du flux a l'interface
f1 = 0.5*ug1.^2 ;
up1 = shiftp('+1',ug1) ;
for i=1:nx1
  if ug1(i) > up1(i)
// choc
     sigma = 0.5*(ug1(i)+up1(i)) ;
     if sigma < 0.
        f1(i) = 0.5*up1(i)^2 ;
     end
  elseif ug1(i) < up1(i)
// detente
     if up1(i) <= 0.
        f1(i) = 0.5*up1(i)^2 ;
     elseif up1(i) >= 0.
        f1(i) = 0.5*ug1(i)^2 ;
     else 
        f1(i) = 0. ;
     end
  end
end
fm1 = shiftp('-1',f1) ;
v1 = ug1 - (dt1/dx1)*(f1-fm1) ;
ug1 = v1 ; 

for n2=1:(nx2/nx1)
//
// Godunov
f2 = 0.5*ug2.^2 ;
up2 = shiftp('+1',ug2) ;
for i=1:nx2
  if ug2(i) > up2(i)
// choc
     sigma = 0.5*(ug2(i)+up2(i)) ;
     if sigma < 0.
        f2(i) = 0.5*up2(i)^2 ;
     end
  elseif ug2(i) < up2(i)
// detente
     if up2(i) <= 0.
        f2(i) = 0.5*up2(i)^2 ;
     elseif up2(i) >= 0.
        f2(i) = 0.5*ug2(i)^2 ;
     else 
        f2(i) = 0. ;
     end
  end
end
fm2 = shiftp('-1',f2) ;
v2 = ug2 - (dt2/dx2)*(f2-fm2) ;
ug2 = v2 ;
end

for n3=1:(nx3/nx1)
//
// Godunov
f3 = 0.5*ug3.^2 ;
up3 = shiftp('+1',ug3) ;
for i=1:nx3
  if ug3(i) > up3(i)
// choc
     sigma = 0.5*(ug3(i)+up3(i)) ;
     if sigma < 0.
        f3(i) = 0.5*up3(i)^2 ;
     end
  elseif ug3(i) < up3(i)
// detente
     if up3(i) <= 0.
        f3(i) = 0.5*up3(i)^2 ;
     elseif up3(i) >= 0.
        f3(i) = 0.5*ug3(i)^2 ;
     else 
        f3(i) = 0. ;
     end
  end
end
fm3 = shiftp('-1',f3) ;
v3 = ug3 - (dt3/dx3)*(f3-fm3) ;
ug3 = v3 ;
end

for n4=1:(nx4/nx1)
//
// Godunov
f4 = 0.5*ug4.^2 ;
up4 = shiftp('+1',ug4) ;
for i=1:nx4
  if ug4(i) > up4(i)
// choc
     sigma = 0.5*(ug4(i)+up4(i)) ;
     if sigma < 0.
        f4(i) = 0.5*up4(i)^2 ;
     end
  elseif ug4(i) < up4(i)
// detente
     if up4(i) <= 0.
        f4(i) = 0.5*up4(i)^2 ;
     elseif up4(i) >= 0.
        f4(i) = 0.5*ug4(i)^2 ;
     else 
        f4(i) = 0. ;
     end
  end
end
fm4 = shiftp('-1',f4) ;
v4 = ug4 - (dt4/dx4)*(f4-fm4) ;
ug4 = v4 ;
end

  clf() 
  plotframe([0,liminf,lg,limmax],tics);
  xset("thickness",2) ;
  plot2d(x1,ug1,[1,1],"100","50 mailles")
  xset("thickness",2) ;
  plot2d(x2,ug2,[2,3],"100","100 mailles")
  xset("thickness",2) ;
  plot2d(x3,ug3,[3,4],"100","200 mailles")
  xset("thickness",2) ;
  plot2d(x4,ug4,[4,6],"100","500 mailles")
  xtitle('comparaison de 4 maillages pour Godunov, cfl=0.9',' ',' ');
  xpause(100000) ;

end

halt() ;

  clf() 
  plotframe([0,liminf,lg,limmax],tics);
  xset("thickness",2) ;
  plot2d(x1,ug1,[1,1],"100","Godunov 50 mailles")
  xset("thickness",2) ;
  plot2d(x4,uf4,[2,3],"100","Lax-Friedrichs 500 mailles")
  xset("thickness",2) ;
  plot2d(x4,uw4,[3,4],"100","Lax-Wendroff 500 mailles")
  xtitle('comparaison de 3 schémas à maillages différents',' ',' ');
  xpause(100000) ;