Advection_1D.m

Tiedoston sisältö

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  Integrate the one-dimensional
%  advection equation
%
%   du/dt + c * du/dx = 0
%
%  using the leapfrog method and
%  assuming periodic boundary conditions
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear;   %  Clear the memory
clf;     %  Clear the graphics window

L = 10000;  %   Domain size (km)
M = 100;    %   Number of grid steps 
DX = L/M;   %   Grid length (km)

T = 5;      %   Time length of integration (days)
N = 100;    %   Number of time steps
DT = T/N;   %   Time step (days)

dx = 1000 * DX;        %   Grid length (m)
dt = (24*60*60) * DT;  %   Time step (s)

c = 10 ;   %    Constant phase speed (m/s)

lambda = c*dt/dx;  %  Courant number
fprintf(' CFL number: %g \n', lambda)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Arrays for the solution

un   = zeros(1,M); %  Solution at time     n*dt
unm1 = zeros(1,M); %  Solution at time (n-1)*dt
unp1 = zeros(1,M); %  Solution at time (n+1)*dt

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Indices for periodic domain

M0 = [1:M    ];  %   Indices for centre points
MM = [M 1:M-1];  %   Indices for points to left
MP = [2:M 1  ];  %   Indices for points to right

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Define the initial conditions

k = 2;    %   Number of wavelengths

un   = 10*sin(2*pi*k*M0/M);

plot(M0,un,'b',M0,un,'rx');
axis([0 M -12 +12]);
hold on; drawnow; pause

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  MAIN LOOP

% First time step is forward
  dudx = un(MP)-un(MM) ;
  unp1 = un - 0.5*(lambda)*dudx;

  % shuffle the solutions
  unm1 = un;
  un   = unp1;
  plot(M0,un); drawnow; pause

for n=2:N

  dudx = un(MP)-un(MM) ;
  unp1 = unm1 - lambda*dudx;

  % shuffle the solutions
  unm1 = un;
  un   = unp1;
  plot(M0,un); drawnow;  pause(0.1)

end

hold off;   %   Release hold on graphics window
return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%