%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright G. Allaire, Novembre 2000
%
% raffinement de maillage  
% avec l'equation de Burgers
% u,t + (u^2/2),x = 0
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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

plot(x3,uf3) ;
axis([0 lg -1.1 1.1]) ;
title('Lax-Friedrichs, cfl=0.9')

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

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

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

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

plot(x1,uf1,'o',x2,uf2,'*',x3,uf3,'d',x4,uf4) ;
axis([0 lg -1.1 1.1]) ;
title('comparaison de 4 maillages pour LF, cfl=0.9')
legend('50 mailles','100 mailles','200 mailles','500 mailles') ;
pause(0) ;

%
end


pause ;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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

plot(x3,uw3) ;
axis([0 lg -1.1 1.1]) ;
title('Lax-Wendroff, cfl=0.9')

pause(1) ;



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% marche en temps : Lax-Wendroff
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 

for n=1:nt1
%
% Lax-Wendroff
up1 = wshift2('1D',uw1,+1) ;
um1 = wshift2('1D',uw1,-1) ;
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 = wshift2('1D',uw2,+1) ;
um2 = wshift2('1D',uw2,-1) ;
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 = wshift2('1D',uw3,+1) ;
um3 = wshift2('1D',uw3,-1) ;
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 = wshift2('1D',uw4,+1) ;
um4 = wshift2('1D',uw4,-1) ;
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

plot(x1,uw1,'o',x2,uw2,'*',x3,uw3,'d',x4,uw4) ;
axis([0 lg -1.1 1.1]) ;
title('comparaison de 4 maillages pour LW, cfl=0.9')
legend('50 mailles','100 mailles','200 mailles','500 mailles') ;
pause(0) ;

%
end

pause ;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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

plot(x3,ug3) ;
axis([0 lg -1.1 1.1]) ;
title('Godunov, cfl=0.9')

pause(1) ;



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% marche en temps : Godunov
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 

for n=1:nt1
%
% Godunov
% evaluation du flux a l'interface
f1 = 0.5*ug1.^2 ;
up1 = wshift2('1D',ug1,+1) ;
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 = wshift2('1D',f1,-1) ;
v1 = ug1 - (dt1/dx1)*(f1-fm1) ;
ug1 = v1 ; 

for n2=1:(nx2/nx1)
%
% Godunov
f2 = 0.5*ug2.^2 ;
up2 = wshift2('1D',ug2,+1) ;
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 = wshift2('1D',f2,-1) ;
v2 = ug2 - (dt2/dx2)*(f2-fm2) ;
ug2 = v2 ;
end

for n3=1:(nx3/nx1)
%
% Godunov
f3 = 0.5*ug3.^2 ;
up3 = wshift2('1D',ug3,+1) ;
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 = wshift2('1D',f3,-1) ;
v3 = ug3 - (dt3/dx3)*(f3-fm3) ;
ug3 = v3 ;
end

for n4=1:(nx4/nx1)
%
% Godunov
f4 = 0.5*ug4.^2 ;
up4 = wshift2('1D',ug4,+1) ;
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 = wshift2('1D',f4,-1) ;
v4 = ug4 - (dt4/dx4)*(f4-fm4) ;
ug4 = v4 ;
end

plot(x1,ug1,'o',x2,ug2,'*',x3,ug3,'d',x4,ug4) ;
axis([0 lg -1.1 1.1]) ;
title('comparaison de 4 maillages pour Godunov, cfl=0.9')
legend('50 mailles','100 mailles','200 mailles','500 mailles') ;
pause(0) ;

end

pause(1) ;

plot(x1,ug1,'o',x4,uf4,x4,uw4) ;
axis([0 lg -1.1 1.1]) ;
title('comparaison de 3 schemas, cfl=0.9')
legend('Godunov 50 mailles','Lax-Friedrichs 500 mailles','Lax-Wendroff 500 mailles') ;

