%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright G. Allaire, Novembre 2000
%
% schemas numeriques pour l'equation de transport en 2-D 
% u,t + a u,x + b u,y = 0
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% lgx, lgy = longueurs du domaine
lgx = 1. ;
lgy = 1. ;
% nx, ny = nombre de mailles
nx = 100 ;
ny = 100. ;
% dx, dy = pas d'espace
dx = lgx/nx ;
dy = lgy/ny ;
% a, b = vitesse de transport
a = 0. ;
b = 1.5 ;
% condition cfl
cfl = min(dx/a,dy/b) ;
% dt = pas de temps
dt = 0.5*cfl ;
nt = 200 ;
%
% 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

surf(x,y,u) ;
mesh(x,y,u) ;
axis([0 lgx 0 lgy -0.1 1.1]) ;
%axis equal ;
view(-37.5,30) ;
colormap('autumn') ;
shading('interp') ;
title('transport 2-D') ;

pause(1) ;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% dt = pas de temps
dt = 1.0*cfl ;
nt = 100 ;
%
% 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

surf(x,y,u) ;
mesh(x,y,u) ;
axis([0 lgx 0 lgy -0.1 1.1]) ;
%axis equal ;
view(-37.5,30) ;
colormap('autumn') ;
shading('interp') ;
title('transport 2-D') ;

pause(1) ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% marche en temps : amont avec coin
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
for n=1:nt
%
% un = u decale vers le nord
un = wshift('2D',u,[0 +1]) ;
% us = u decale vers le sud
us = wshift('2D',u,[0 -1]) ;
% uo = u decale vers l'ouest
uo = wshift('2D',u,[-1 0]) ;
% ue = u decale vers l'est
ue = wshift('2D',u,[+1 0]) ;
% uc = u decale vers le sud-ouest
uc = wshift('2D',u,[-1 -1]) ;

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

if mod(n,5) == 0
surf(x,y,u) ;
mesh(x,y,u) ;
axis([0 lgx 0 lgy -0.1 1.1]) ;
%axis equal ;
view(-37.5,30) ;
colormap('autumn') ;
shading('interp') ;
title('schema amont avec coin, cfl=1.0')
pause(0) ;
end
%
end

pause ;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% nouveau schema: transport amont avec coins du 2eme ordre
% (on suppose a et b > 0.)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% dt = pas de temps
dt = 1.0*cfl ;
nt = 100 ;
%
% 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

surf(x,y,u) ;
mesh(x,y,u) ;
axis([0 lgx 0 lgy -0.1 1.1]) ;
%axis equal ;
view(-37.5,30) ;
colormap('autumn') ;
shading('interp') ;
title('transport 2-D') ;

pause(1) ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% marche en temps : amont avec coin du 2eme ordre
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
for n=1:nt
%
% un = u decale vers le nord
un = wshift('2D',u,[0 +1]) ;
% us = u decale vers le sud
us = wshift('2D',u,[0 -1]) ;
% uo = u decale vers l'ouest
uo = wshift('2D',u,[-1 0]) ;
% ue = u decale vers l'est
ue = wshift('2D',u,[+1 0]) ;
% uso = u decale vers le sud-ouest
uso = wshift('2D',u,[-1 -1]) ;
% use = u decale vers le sud-est
use = wshift('2D',u,[+1 -1]) ;
% uno = u decale vers le nord-ouest
uno = wshift('2D',u,[-1 +1]) ;

ca = (dt/dx)*a ;
cb = (dt/dy)*b ;
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 mod(n,5) == 0
surf(x,y,u) ;
mesh(x,y,u) ;
axis([0 lgx 0 lgy -0.1 1.1]) ;
%axis equal ;
view(-37.5,30) ;
colormap('autumn') ;
shading('interp') ;
title('schema amont avec coin 2eme ordre, cfl=1.0')
pause(0) ;
end
%
end

pause ;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% nouveau schema: Lax-Friedrichs par splitting 1-D
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% dt = pas de temps
dt = 1.0*cfl ;
nt = 200 ;
%
% 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

surf(x,y,u) ;
mesh(x,y,u) ;
axis([0 lgx 0 lgy -0.1 1.1]) ;
%axis equal ;
view(-37.5,30) ;
colormap('autumn') ;
shading('interp') ;
title('transport 2-D') ;

pause(1) ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% marche en temps : Lax-Friedrichs 1-D
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
for n=1:nt
%
% uo = u decale vers l'ouest
uo = wshift('2D',u,[-1 0]) ;
% ue = u decale vers l'est
ue = wshift('2D',u,[+1 0]) ;

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

% un = u decale vers le nord
un = wshift('2D',u,[0 +1]) ;
% us = u decale vers le sud
us = wshift('2D',u,[0 -1]) ;

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

if mod(n,5) == 0
surf(x,y,u) ;
mesh(x,y,u) ;
axis([0 lgx 0 lgy -0.1 1.1]) ;
%axis equal ;
view(-37.5,30) ;
colormap('autumn') ;
shading('interp') ;
title('Lax-Friedrichs splitting 1-D, cfl=1.0')
pause(0) ;
end
%
end

pause ;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% nouveau schema: Lax-Wendroff par splitting 1-D
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% dt = pas de temps
dt = 1.0*cfl ;
nt = 200 ;
%
% 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

surf(x,y,u) ;
mesh(x,y,u) ;
axis([0 lgx 0 lgy -0.1 1.1]) ;
%axis equal ;
view(-37.5,30) ;
colormap('autumn') ;
shading('interp') ;
title('transport 2-D') ;

pause(1) ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% marche en temps : Lax-Wendroff 1-D
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
for n=1:nt
%
% uo = u decale vers l'ouest
uo = wshift('2D',u,[-1 0]) ;
% ue = u decale vers l'est
ue = wshift('2D',u,[+1 0]) ;

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

% un = u decale vers le nord
un = wshift('2D',u,[0 +1]) ;
% us = u decale vers le sud
us = wshift('2D',u,[0 -1]) ;

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

if mod(n,5) == 0
surf(x,y,u) ;
mesh(x,y,u) ;
axis([0 lgx 0 lgy -0.1 1.1]) ;
%axis equal ;
view(-37.5,30) ;
colormap('autumn') ;
shading('interp') ;
title('Lax-Wendroff splitting 1-D, cfl=1.0')
pause(0) ;
end
%
end

pause ;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% nouveau schema: TVD superbee par splitting 1-D
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% dt = pas de temps
dt = 1.0*cfl ;
nt = 200 ;
%
% 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

surf(x,y,u) ;
mesh(x,y,u) ;
axis([0 lgx 0 lgy -0.1 1.1]) ;
%axis equal ;
view(-37.5,30) ;
colormap('autumn') ;
shading('interp') ;
title('transport 2-D') ;

pause(1) ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% marche en temps : TVD superbee 1-D
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
for n=1:nt
%
% uo = u decale vers l'ouest
uo = wshift('2D',u,[-1 0]) ;
% ue = u decale vers l'est
ue = wshift('2D',u,[+1 0]) ;
du = u-uo ;
du = du + 1.e-6*(1+sign(du)).*(1-sign(du)) ;
dum = wshift('2D',du,[-1 0]) ;
phi = dum./du ;
phi = max(min(2*phi,1.),min(2.,phi)) ;
phi = max(phi,0.) ;
phim = wshift('2D',phi,[-1 0]) ;

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

% un = u decale vers le nord
un = wshift('2D',u,[0 +1]) ;
% us = u decale vers le sud
us = wshift('2D',u,[0 -1]) ;
du = u-us ;
du = du + 1.e-6*(1+sign(du)).*(1-sign(du)) ;
dum = wshift('2D',du,[-1 0]) ;
phi = dum./du ;
phi = max(min(2*phi,1.),min(2.,phi)) ;
phi = max(phi,0.) ;
phim = wshift('2D',phi,[-1 0]) ;

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

if mod(n,5) == 0
surf(x,y,u) ;
mesh(x,y,u) ;
axis([0 lgx 0 lgy -0.1 1.1]) ;
%axis equal ;
view(-37.5,30) ;
colormap('autumn') ;
shading('interp') ;
title('TVD superbee splitting 1-D, cfl=1.0')
pause(0) ;
end
%
end

