sandbox/easystab/advection_2D.m

    Advection 2D : f_t+Uf_x+Vf_y=0

    In this code, we do the marching in time of the advection equation \displaystyle f_{t}+Uf_{x}+Vf_{y}=0 with the boundary condition at the left by changing the code given vibrating_string.m.

    clear all; clf

% parameters
Nx=30; % number of gridpoints in x
Ny=30; % number of gridpoints in y 
Lx=10; % domain length x
Ly=10; % domain length y
U=1;% x velocity
V=0;% y velocity
dt=0.01; % time step
tmax=10; % final time

x0=Lx/4; % x-coordinate at initial time
y0=0;% y-coordinate at initial time
l0=1; % length width

% differentiation matrices
scale=-2/Lx;
[x,DM] = chebdif(Nx,2); 
dx=DM(:,:,1)*scale;        
x=(x-1)/scale; 

scale=-2/Ly;
[y,DM] = chebdif(Ny,1); 
dy=DM(:,:,1)*scale;    
y=(y-1)/scale; 

[X,Y]=meshgrid(x,y);
% 2D differentiation matrices
Dx=kron(dx,eye(Ny));
Dy=kron(eye(Nx),dy);
[X,Y]=meshgrid(x,y);


I=eye(Nx*Ny);

    System matrices

    Here we build the matrices that give the discretization of the equations. The equations is \displaystyle f_{t}+Uf_{x}+Vf_{y}=0 We can also write this equation in the following way \displaystyle f_{t}=-UD_xf-VD_yf

    We can put an identity matrix on the left of the equation, in this way, we can impose the boundary condition more easliy.

    So the equation becomes : \displaystyle If_{t}=-UD_xf-VD_yf

    Now we have a linear system with the form \displaystyle Ef_{t}=Af 

    % system matrices
E=I;
A=-U*Dx-V*Dy;

    Boundary conditions

    For 2D problems, by extracting the boundary cell locations (methode used in poisson2D.m), we can set our boundary conditions more efficiantly. Know that we have just one derivite on x in our equation, the only condition we have to set is the left side of the array.

    % boundary conditions

% locations at the boundary
dom=reshape(1:Nx*Ny,Ny,Nx);
top=dom(1,2:end-1); top=top(:); 
bot=dom(end,2:end-1); bot=bot(:); 
left=dom(1:end,1); left=left(:); 
right=dom(1:end,end); right=right(:); 
II=eye(Nx*Ny); ZZ=zeros(Nx*Ny,Nx*Ny);
loc=[left];


A(loc,:)=II(loc,:);
E(loc,:)=0;

    March in time

    Here we use the same methode with vibrating_string.m

    % march in time matrix 
M=(E-A*dt/2)\(E+A*dt/2);

    Initial condition

    Here we need to describe the initial state of the system. We say here that the surface is initially deformed as a bell curve, and that the velocity is zero. From that state we then perform a loop to repeatedly advance the state of a short forward step in time. We use drawnow to show the evolution of the simulation as a movie when running the code. We store for validation the string position at the midle of the domain, to do this without worrying about the the grid points are, we interpolate f with the function interp1.

    % initial condition
q=Y*0+exp(-((X-x0)/l0).^2-((Y-y0)/l0).^2); 
q=q(:);

% marching loop
tvec=dt:dt:tmax; 
Nt=length(tvec);
e=zeros(length(y),Nt);

for ind=1:Nt    
    q=M*q; % one step forward
    q_reshape=reshape(q,Ny,Nx);
    e(:,ind)=interp2(X,Y,q_reshape,Lx/2,y); % store center point
    
    % plotting
    figure (1)
    mesh(X,Y,q_reshape); 
   
    axis([0,Lx,0,Ly,0,1]);
    drawnow
end
legend('position'); title('Vibrating string')
xlabel('x'); ylabel('y');zlabel('f');

    Validation

    Now we show that the code does what we think it does, and also give ourselves the means to tell how precise and robust it is.In order to verify the solutions, we draw the position of the point in the middle of L at different times.

    % time evolution of central point
figure(2);
[TVEC,YT]=meshgrid(tvec,y);
mesh(TVEC,YT,e); 
axis([0,Lx,0,Ly,0,1]);

hold on

tt=linspace(0,tmax,500);
[TT,YTT]=meshgrid(tt,y);
etheo=YTT*0+exp(-((Lx/2-U*TT-x0)/l0).^2-((YTT-V*TT-y0)/l0).^2); 
surf(TT,YTT,etheo); 
axis([0,Lx,0,Ly,0,1]);
xlabel('temp/s');
zlabel('amplitude');
title('Valitation');
    validation validation

    Form the figure of validation, we find a superposition of the theory results and the numeral results

    ZHAO’s contribution

     See also advection_1D.m

     Or link to my contribution page zhao.m

    Notes

    De la part de Jérôme

    domaine	        valeur	note
connectivité 	2	2
recyclage	2	2
graphiques	2	1
théories	4	0
Originalité	4	0
note /14	14	5

    Ce que tu as fait est bien mais ça correspond seulement aux deux premières séances de TP. Dommage que tu as laissé tombé.