//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
// Paola Goatin, Octobre 2006
// Translated from a Matlab code by G. Allaire
// 
// schemas de Godunov et de Van Leer  
// pour l'equation de Burgers
// u,t + (u^2/2),x = 0
// 
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
//
exec('solbu.sci'); 
exec('wshift2.sci'); 
exec('wshift3.sci'); 
// 
// lg = longueur du domaine
lg = 5;
// nx = nombre de mailles
nx = 50;
// dx = pas d'espace
dx = lg/nx;
// etat gauche et droit
ul = 0.8;
ur = 0.;
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 = points du maillage
x = zeros(1,nx+1);
// ug = solution de Godunov
ug = zeros(1,nx+1);
// uv = solution de Van Leer
uv = zeros(1,nx+1);
// ue = solution exacte
ue = zeros(1,nx+1);
// f = flux
f = zeros(1,nx+1);
// p = pente
p = zeros(1,nx+1);
// du = increment u(i+1)-u(i)
du = zeros(1,nx+1);
// tableaux auxiliaires 
v = zeros(1,nx+1);
w = zeros(1,nx+1);
up = zeros(1,nx+1);
um = zeros(1,nx+1);
uvp = zeros(1,nx+1);
uvm = zeros(1,nx+1);
uvmp = zeros(1,nx+1);
fm = zeros(1,nx+1);
dum = zeros(1,nx+1);
sdu = zeros(1,nx+1);
sdum = zeros(1,nx+1);
// donnee initiale
for i = 1:nx+1;
  x(i) = (i-1)*dx;
  if x(i) < (lg/2) 
    ue(i) = ul;
  else
    ue(i) = ur;
  end;
end;
ug = ue;
uv = ue;
  
  clf()
  tics=[4,10,4,10];
  plotframe([0,liminf,lg,limmax],tics);
  plot2d(x,ue,1,"000")
  xtitle('comparaison Godunov - Van Leer, cfl=0.9');

halt();
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
// marche en temps : 
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
// 
for n = 1:nt;
  // 
  // Godunov
  // evaluation du flux a l'interface
  f = 0.5*(ug.^2);
  for i = 1:nx
    if ug(i) > ug(i+1) 
      // choc
//      // definir la vitesse du choc 'sigma'
//      // definir la valeur du flux 'f'
    elseif ug(i) < ug(i+1) 
      // detente
      if ug(i+1) < 0.
        f(i) = 0.5*ug(i+1)^2;
      elseif ug(i) > 0
        f(i) = 0.5*ug(i)^2;
      else
        f(i) = 0.;
      end
    end
  end
  fm(2:nx+1) = f(1:nx);
  v(2:nx) = ug(2:nx)-(dt/dx)*(f(2:nx)-fm(2:nx));
  v(1) = ug(1);
  v(nx+1) = ug(nx+1);
  ug = v;


  // Van Leer
  up(1:nx) = uv(2:nx+1);
  um(2:nx+1) = uv(1:nx);
  du = up-uv;
  dum(2:nx+1) = du(1:nx);
  sdu = sign(du);
  sdum = sign(dum);
  // evaluation de la pente par la fonction minmod
  p = max(0,sdu.*sdum) .*sdu .*min(abs(du),abs(dum));
  // uvm = solution en (i+1/2)-
  // uvp = solution en (i-1/2)+
  uvm = uv+0.5*p;
  uvp = uv-0.5*p;
  // evaluation de la solution au demi pas de temps
//  // mise a jour de 'uvm(2:nx)' avec 'v(2:nx)' 
//  // mise a jour de 'uvp(2:nx)' avec 'w(2:nx)'
  v(1) = uv(1);
  v(nx+1) = uv(nx+1);
  w(1) = uv(1);
  w(nx+1) = uv(nx+1);
  uvm = v;
  uvp = w;
  uvpp = wshift3('1D',uvp,1);
  uvpp(1:nx) = uvp(2:nx+1);
  // resolution du probleme de Riemann entre uvm et uvpp
  // evaluation du flux a l'interface
  f = 0.5*uvm.^2;
  for i = 1:nx
    if uvm(i) > uvpp(i)
      // choc
//      // definir la vitesse du choc 'sigma'
//      // definir la valeur du flux 'f' 
    elseif uvm(i) < uvpp(i)
      // detente
//      // definir la valeur du flux 'f' 
    end
  end
  fm(2:nx+1) = f(1:nx);
  v(2:nx) = uv(2:nx) - (dt/dx)*(f(2:nx)-fm(2:nx));
  v(1) = uv(1);
  v(nx+1) = uv(nx+1);
  uv = v;


  if pmodulo(n,1)==0 
    ue = solbu(ue,nx,dx,dt,lg,ul,ur,n);
    clf();
    plotframe([0,liminf,lg,limmax],tics);
    plot2d(x',[ug',uv',ue'],[-1,-3,3])
    xtitle('comparaison Godunov - Van Leer, cfl=0.9');
    legends(['Godunov','Van Leer','exact'],[-1,-3,3],opt=1);
    xpause(100000);
  end
  // 
end

halt();

//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
// conditions aux limites periodiques
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
// 
lg = 1;
// nx = nombre de mailles
nx = 50;
// dx = pas d'espace
dx = lg/nx;
// dt = pas de temps
// (on sait que la vitesse maximale sera 2.)
cfl = dx/2;
dt = 0.9*cfl;
// nt = nombre de pas de temps effectues
nt = 300;
// 
// initialisation
// 
// x = points du maillage
x = zeros(1,nx);
// ug = solution de Godunov
ug = zeros(1,nx);
// uv = solution de Van Leer
uv = zeros(1,nx);
// f = flux
f = zeros(1,nx);
// p = pente
p = zeros(1,nx);
// du = increment u(i+1)-u(i)
du = zeros(1,nx);
// tableaux auxiliaires 
v = zeros(1,nx);
w = zeros(1,nx);
up = zeros(1,nx);
um = zeros(1,nx);
uvp = zeros(1,nx);
uvm = zeros(1,nx);
uvpp = zeros(1,nx);
fm = zeros(1,nx);
dum = zeros(1,nx);
sdu = zeros(1,nx);
sdum = zeros(1,nx);
// donnee initiale
for i = 1:nx;
  x(i) = (i-1)*dx;
  ug(i) = 1. + sin(x(i)*2*%pi);
end;
uv = ug;

// affichage de la donnee initiale
  clf()
  tics=[4,10,4,10];
  plotframe([0,-0.1,lg,2.1],tics);
  plot2d(x,ug,1,"000")
  xtitle('comparaison Godunov - Van Leer, cfl=0.9');

halt();
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
// marche en temps : 
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
// 
for n = 1:nt;
  // 
  // Godunov
  // evaluation du flux a l'interface
  f = 0.5*ug.^2;
  up = wshift2('1D',ug,1);
  for i = 1:nx;
    if ug(i) > up(i)
      // choc
//  // completer
    elseif ug(i) < up(i) then
      // detente
//  // completer
    end
  end
  fm = wshift2('1D',f,-1);
  v = ug-(dt/dx)*(f-fm);
  ug = v;

  // Van Leer
  um = wshift2('1D',uv,-1);
  up = wshift2('1D',uv,1);
  du = up-uv;
  dum = wshift2('1D',du,-1);
  sdu = sign(du);
  sdum = sign(dum);
  // evaluation de la pente par la fonction minmod
  p = max(0,sdu.*sdum).*sdu.*min(abs(du),abs(dum));
  // uvm = solution en (i+1/2)-
  // uvp = solution en (i-1/2)+
  uvm = uv+0.5*p;
  uvp = uv-0.5*p;
  // evaluation de la solution au demi pas de temps
//  v = uvm - // completer
//  w = uvp - // completer
  uvm = v;
  uvp = w;
  uvpp = wshift2('1D',uvp,1);
  // resolution du probleme de Riemann entre uvm et uvpp
  // evaluation du flux a l'interface
  f = 0.5*uvm.^2;
  for i = 1:nx;
    if uvm(i) > uvpp(i) 
      // choc
 // completer
    elseif uvm(i) < uvpp(i)
      // detente
// completer
    end
  end
  fm = wshift2('1D',f,-1);
  v = uv-(dt/dx)*(f-fm);
  uv = v;

  if pmodulo(n,1) == 0 
    clf();
    plotframe([0,-0.1,lg,2.1],tics);
    plot2d(x',[ug',uv'],[-1,-3])
    xtitle('comparaison Godunov - Van Leer, cfl=0.9');
    legends(['Godunov','Van Leer'],[-1,-3],opt=1);
    xpause(100000);
  end
  // 
end