/////////////////////////////////////////////////////////////
// Copyright G. Allaire, Janvier 2008, Janvier 2012
//
// schemas numeriques sur l'equation de transport lineaire 
// avec condition aux limites periodiques
// u,t + ax u,x + ay u,y = 0
// 
/////////////////////////////////////////////////////////////
//
xset("font",2,3) ; 
xset("thickness",2) ;
lgx = 1. ; // lgx = longueur du domaine
lgy = 1. ; // lgy = hauteur du domaine
dx = 0.01 ; // dx = pas d'espace
dy = 0.01 ; // dy = pas d'espace
nx = lgx/dx ; // nx = nombre de mailles en x
ny = lgy/dy ; // ny = nombre de mailles en y 
cfl = 1. ; // cfl
ax = 1. ; // ax = vitesse selon x
ay = 0.5 ; // ay = vitesse selon y
if (ax~=0.)
if (ay~=0.)
   dt = cfl*min(dx/ax,dy/ay) ; // dt = pas de temps
else
   dt = cfl*dx/ax ;
end
else
   dt = cfl*dx/ay ;
end
nt = 100 ; // nt = nombre de pas de temps effectues
//
// initialisation
//
x=zeros(1,nx) ;
y=zeros(1,ny) ;
u=zeros(nx,ny) ;
v=zeros(nx,ny) ;
w=zeros(nx,ny) ;
un=zeros(nx,ny) ;
us=zeros(nx,ny) ;
ue=zeros(nx,ny) ;
uo=zeros(nx,ny) ;
// (x,y) = points du maillage
// u = donnee initiale
for i=1:nx
  x(i) = (i-1)*dx  ;
end
for j=1:ny
  y(j) = (j-1)*dy  ;
end
for i=1:nx
for j=1:ny
  r = sqrt((x(i)-lgx/4)^2 + (y(j)-lgy/2)^2) ;
  if r<=0.1 
     u(i,j) = 1. ;
  end
end
end


clf()
leg='x@y@u' ;
flag=[0,1,4] ;
ebox=[0.,lgx,0.,lgy,0.,1.] ;
plot3d(x,y,u,35,45,leg,flag,ebox) ;
xtitle ('donnée initiale' ,' ',' ') ;

halt() ;

////////////////////////////////////////////////
// marche en temps : amont avec coin
////////////////////////////////////////////////
// 
for n=1:nt
//
// un = u decale vers le nord
un = shiftp2d('n',u) ;
// us = u decale vers le sud
us = shiftp2d('s',u) ;
// uo = u decale vers l'ouest
uo = shiftp2d('o',u) ;
// ue = u decale vers l'est
ue = shiftp2d('e',u) ;
// uc = u decale vers le sud-ouest
uc = shiftp2d('so',u) ;

ca = (dt/dx)*ax ;
cb = (dt/dy)*ay ;
v = (1-ca-cb+ca*cb)*u + cb*(1-ca)*us + ca*(1-cb)*uo... 
  + ca*cb*uc ;
u = v ; 

if modulo(n,5) == 0
        clf()
//        flag=[0,0,4] ;
        plot3d(x,y,u,35,45,leg,flag,ebox) ;
        titre =  msprintf('Schema amont avec coin, temps=%.2f, nombre de pas de temps=%i',n*dt,n) ; 
        xtitle(titre,' ',' ');
        xpause(100000) ;
end
//
end

halt() ;


//////////////////////////////////////////////////////////////
// nouveau schema: transport amont avec coins du 2eme ordre
// (on suppose a et b > 0.)
//////////////////////////////////////////////////////////////
//
//
// initialisation
//
x=zeros(1,nx) ;
y=zeros(1,ny) ;
u=zeros(nx,ny) ;
v=zeros(nx,ny) ;
w=zeros(nx,ny) ;
un=zeros(nx,ny) ;
us=zeros(nx,ny) ;
ue=zeros(nx,ny) ;
uo=zeros(nx,ny) ;
// (x,y) = points du maillage
// u = donnee initiale
for i=1:nx
  x(i) = (i-1)*dx  ;
end
for j=1:ny
  y(j) = (j-1)*dy  ;
end
for i=1:nx
for j=1:ny
  r = sqrt((x(i)-lgx/4)^2 + (y(j)-lgy/2)^2) ;
  if r<=0.1 
     u(i,j) = 1. ;
  end
end
end

// 
for n=1:nt
//
// un = u decale vers le nord
un = shiftp2d('n',u) ;
// us = u decale vers le sud
us = shiftp2d('s',u) ;
// uo = u decale vers l'ouest
uo = shiftp2d('o',u) ;
// ue = u decale vers l'est
ue = shiftp2d('e',u) ;
// uso = u decale vers le sud-ouest
uso = shiftp2d('so',u) ;
// use = u decale vers le sud-est
use = shiftp2d('se',u) ;
// uno = u decale vers le nord-ouest
uno = shiftp2d('no',u) ;

ca = (dt/dx)*ax ;
cb = (dt/dy)*ay ;
v = (1-0.5*ca*ca*(2-cb)-0.5*cb*cb*(2-ca))*u...
  - 0.25*ca*(1-ca)*(2-cb)*ue... 
  - 0.25*cb*(1-cb)*(2-ca)*un...
  - 0.25*ca*cb*(1-ca)*use... 
  - 0.25*ca*cb*(1-cb)*uno...
  + 0.25*(ca*(1+ca)*(2-cb)-2*ca*cb*cb)*uo...
  + 0.25*(cb*(1+cb)*(2-ca)-2*ca*ca*cb)*us...
  + 0.25*ca*cb*(2+ca+cb)*uso ;
u = v ; 

if modulo(n,5) == 0
        clf()
//        flag=[0,0,4] ;
        plot3d(x,y,u,35,45,leg,flag,ebox) ;
        titre =  msprintf('Schema amont avec coin 2eme ordre, temps=%.2f, nombre de pas de temps=%i',n*dt,n) ; 
        xtitle(titre,' ',' ');
        xpause(100000) ;
end
//
end

halt() ;


//////////////////////////////////////////////////////////////
// nouveau schema: Lax-Friedrichs par splitting 1-D
//////////////////////////////////////////////////////////////
//
// initialisation
//
x=zeros(1,nx) ;
y=zeros(1,ny) ;
u=zeros(nx,ny) ;
v=zeros(nx,ny) ;
w=zeros(nx,ny) ;
un=zeros(nx,ny) ;
us=zeros(nx,ny) ;
ue=zeros(nx,ny) ;
uo=zeros(nx,ny) ;
// (x,y) = points du maillage
// u = donnee initiale
for i=1:nx
  x(i) = (i-1)*dx  ;
end
for j=1:ny
  y(j) = (j-1)*dy  ;
end
for i=1:nx
for j=1:ny
  r = sqrt((x(i)-lgx/4)^2 + (y(j)-lgy/2)^2) ;
  if r<=0.1 
     u(i,j) = 1. ;
  end
end
end

// 
for n=1:nt
//
// uo = u decale vers l'ouest
uo = shiftp2d('o',u) ;
// ue = u decale vers l'est
ue = shiftp2d('e',u) ;

v = (uo+ue)/2 - (0.5*dt/dx)*ax*(ue-uo) ;
u = v ; 

// un = u decale vers le nord
un = shiftp2d('n',u) ;
// us = u decale vers le sud
us = shiftp2d('s',u) ;

v = (un+us)/2 - (0.5*dt/dy)*ay*(un-us) ;
u = v ; 

if modulo(n,5) == 0
        clf()
//        flag=[0,0,4] ;
        plot3d(x,y,u,35,45,leg,flag,ebox) ;
        titre =  msprintf('Schema Lax-Friedrichs splitting 1-D, temps=%.2f, nombre de pas de temps=%i',n*dt,n) ; 
        xtitle(titre,' ',' ');
        xpause(100000) ;
end
//
end

halt() ;

////////////////////
// nouveau schema: Lax-Wendroff par splitting 1-D
////////////////////
//
//
// initialisation
//
x=zeros(1,nx) ;
y=zeros(1,ny) ;
u=zeros(nx,ny) ;
v=zeros(nx,ny) ;
w=zeros(nx,ny) ;
un=zeros(nx,ny) ;
us=zeros(nx,ny) ;
ue=zeros(nx,ny) ;
uo=zeros(nx,ny) ;
// (x,y) = points du maillage
// u = donnee initiale
for i=1:nx
  x(i) = (i-1)*dx  ;
end
for j=1:ny
  y(j) = (j-1)*dy  ;
end
for i=1:nx
for j=1:ny
  r = sqrt((x(i)-lgx/4)^2 + (y(j)-lgy/2)^2) ;
  if r<=0.1 
     u(i,j) = 1. ;
  end
end
end

// 
for n=1:nt
//
// uo = u decale vers l'ouest
uo = shiftp2d('o',u) ;
// ue = u decale vers l'est
ue = shiftp2d('e',u) ;

v = u - (0.5*dt/dx)*ax*(ue-uo) + (0.5*dt*dt/(dx*dx))*ax*ax*(uo-2*u+ue) ;
u = v ; 

// un = u decale vers le nord
un = shiftp2d('n',u) ;
// us = u decale vers le sud
us = shiftp2d('s',u) ;

v = u - (0.5*dt/dy)*ay*(un-us) + (0.5*dt*dt/(dy*dy))*ay*ay*(un-2*u+us) ;
u = v ; 

if modulo(n,5) == 0
        clf()
//        flag=[0,0,4] ;
        plot3d(x,y,u,35,45,leg,flag,ebox) ;
        titre =  msprintf('Schema Lax-Wendroff splitting 1-D, temps=%.2f, nombre de pas de temps=%i',n*dt,n) ; 
        xtitle(titre,' ',' ');
        xpause(100000) ;
end
//
end

halt() ;


////////////////////////
// nouveau schema: TVD superbee par splitting 1-D
//////////////////////
//
//
// initialisation
//
x=zeros(1,nx) ;
y=zeros(1,ny) ;
u=zeros(nx,ny) ;
v=zeros(nx,ny) ;
w=zeros(nx,ny) ;
un=zeros(nx,ny) ;
us=zeros(nx,ny) ;
ue=zeros(nx,ny) ;
uo=zeros(nx,ny) ;
du=zeros(nx,ny) ;
dum=zeros(nx,ny) ;
phi=zeros(nx,ny) ;
phim=zeros(nx,ny) ;
// (x,y) = points du maillage
// u = donnee initiale
for i=1:nx
  x(i) = (i-1)*dx  ;
end
for j=1:ny
  y(j) = (j-1)*dy  ;
end
for i=1:nx
for j=1:ny
  r = sqrt((x(i)-lgx/4)^2 + (y(j)-lgy/2)^2) ;
  if r<=0.1 
     u(i,j) = 1. ;
  end
end
end

// 
for n=1:nt
//
// uo = u decale vers l'ouest
uo = shiftp2d('o',u) ;
// ue = u decale vers l'est
ue = shiftp2d('e',u) ;
du = u-uo ;
du = du + 1.e-6*(1+sign(du)).*(1-sign(du)) ;
dum = shiftp2d('o',du) ;
phi = dum./du ;
phi = max(min(2*phi,1.),min(2.,phi)) ;
phi = max(phi,0.) ;
phim = shiftp2d('o',phi) ;

v = u - (dt/dx)*ax*(u-uo) + 0.5*(dt/dx)*ax*...
    ((dt/dx)*ax-1.)*(phi.*(ue-u)-phim.*(u-uo)) ;
u = v ; 

// un = u decale vers le nord
un = shiftp2d('n',u) ;
// us = u decale vers le sud
us = shiftp2d('s',u) ;
du = u-us ;
du = du + 1.e-6*(1+sign(du)).*(1-sign(du)) ;
dum = shiftp2d('n',du) ;
phi = dum./du ;
phi = max(min(2*phi,1.),min(2.,phi)) ;
phi = max(phi,0.) ;
phim = shiftp2d('n',phi) ;

v = u - (dt/dy)*ay*(u-us) + 0.5*(dt/dy)*ay*...
    ((dt/dy)*ay-1.)*(phi.*(un-u)-phim.*(u-us)) ;
u = v ; 

if modulo(n,5) == 0
        clf()
//        flag=[0,0,4] ;
        plot3d(x,y,u,35,45,leg,flag,ebox) ;
        titre =  msprintf('Schema TVD superbee splitting 1-D, temps=%.2f, nombre de pas de temps=%i',n*dt,n) ; 
        xtitle(titre,' ',' ');
        xpause(100000) ;
end
//
end