function x0new = Advection(InVec,CVec,p,dt,nsteps) % setup the mesh on which to solve % p=40; m=p^2; m=p^2; dx = 1/(p+1); x0 = InVec; % --------------------------------------------------------------------------------------------------------------- % compute the initial solution profile x0 (soln to heat equation) Ab=zeros(m,p+1); Ab(:,1)=4; Ab(:,2)=-1; Ab(p:p:m,2)=0; Ab(:,p+1)=-1; A=spdiags(Ab(:,[1,2,p+1]),-[0,1,p],m,m); A=A+tril(A,-1)'; % f=ones(m,1)/(p+1)^2; % x0=A\f; % --------------------------------------------------------------------------------------------------------------- % Calculate the first derivative matrices Ac_xdirn=zeros(m,p+1); Ac_xdirn(:,1)=1; Ac_xdirn(p+1:p:m+1,1)=0; Ac_xdirn(:,2)=-1; Ac_xdirn(p:p:m,2)=0; Ac_xdirn(:,3)=0; Ac_xdirn(1:p:m-p+1,3)=1; Ac_xdirn(:,4)=0; Ac_xdirn(1:p:m-p+1,4)=-1; Ac_ydirn=zeros(m,2); Ac_ydirn(:,1)=-1; Ac_ydirn(:,2)=1; A_yfirstderiv=spdiags(Ac_ydirn(:,[1,1,2,2]),[m-p,-p,p-m,p],m,m); A_xfirstderiv=spdiags(Ac_xdirn(:,[2,3]),[-1,-p+1],m,m); A_xfirstderiv = A_xfirstderiv - tril(A_xfirstderiv,-1)'; % disp('full(A_xfirstderiv)') % full(A_xfirstderiv) % disp('full(A_yfirstderiv)') % full(A_yfirstderiv) % --------------------------------------------------------------------------------------------------------------- % Calculate the second derivative matrices d2_xdirn=zeros(m,p+1); d2_xdirn(:,1)=-2; d2_xdirn(p+1:p:m+1,1)=0; d2_xdirn(:,2)=1; d2_xdirn(p:p:m,2)=0; d2_xdirn(:,3)=0; d2_xdirn(1:p:m-p+1,3)=1; d2_xdirn(:,4)=0; d2_xdirn(1:p:m-p+1,4)=-1; d2_xderiv=spdiags(d2_xdirn(:,[1,2,3]),[0,-1,-p+1],m,m); d2_xderiv = d2_xderiv - tril(d2_xderiv,-1)'; d2_ydirn=zeros(m,3); d2_ydirn(:,1)=1; d2_ydirn(:,2)=1; d2_ydirn(:,3)=-2; d2_yderiv=spdiags(d2_ydirn(:,[1,1,2,2,3]),[m-p,-p,p-m,p,0],m,m); % disp('full(d2_xderiv)') % full(d2_xderiv) % disp('full(d2_yderiv)') % full(d2_yderiv) % --------------------------------------------------------------------------------------------------------------- % Compute gradients of 'c' to obtain x and y advection velocities % do this tomorrow.... ConvectionVel_xdir = A_xfirstderiv*CVec; % norm(ConvectionVel_xdir) ConvectionVel_ydir = A_yfirstderiv*CVec; % norm(ConvectionVel_ydir) ConvectionVel_xdir = 15*ConvectionVel_xdir; ConvectionVel_ydir = 15*ConvectionVel_ydir; % ConvectionVel_xdir = ones(m,1); % ConvectionVel_ydir = ones(m,1); % ConvectionVel_xdir % ConvectionVel_ydir % alph = 0.01; % xVel = 0.01; % yVel = 0.0; % yVel = 0.01; % DiffCoeff1x = (xVel.*dt)/(2*dx); % DiffCoeff1y = (yVel.*dt)/(2*dx); % alpha = 0.01; alpha = 1; DiffCoeff1x = (alpha*dt)/(2*dx); DiffCoeff1y = (alpha*dt)/(2*dx); DiffCoeff2x = 2*(DiffCoeff1x.^2); DiffCoeff2y = 2*(DiffCoeff1y.^2); % DiffCoeff1 = Vel*dt/(2.*dx); % DiffCoeff2 = Vel*2*(DiffCoeff1.^2); % full(A_xfirstderiv) % full(A_yfirstderiv) % full(A_xsecondderiv) % full(A_ysecondderiv) % I = sparse(1:m,1:m,1); % ConvectionVel_xdir = 10; % ConvectionVel_ydir = 10; for i = 1:nsteps % disp('did an advection!') x0new = x0 - DiffCoeff1x*A_xfirstderiv*(ConvectionVel_xdir.*x0) - DiffCoeff1y*A_yfirstderiv*(ConvectionVel_ydir.*x0) - DiffCoeff2x*d2_xderiv*(ConvectionVel_xdir.*ConvectionVel_xdir.*x0) - DiffCoeff2y*d2_yderiv*(ConvectionVel_ydir.*ConvectionVel_ydir.*x0); % x0new = x0 - DiffCoeff1x*A_xfirstderiv*(ConvectionVel_xdir.*x0) - DiffCoeff1y*A_yfirstderiv*(ConvectionVel_ydir.*x0) - DiffCoeff2x*d2_xderiv*(ConvectionVel_xdir.*ConvectionVel_xdir.*x0) - DiffCoeff2y*d2_yderiv*(ConvectionVel_ydir.*ConvectionVel_ydir.*x0); % x0new = (I - DiffCoeff1x*A_xfirstderiv - DiffCoeff1y*A_yfirstderiv - DiffCoeff2x*d2_xderiv - DiffCoeff2y*d2_yderiv)*x0; % close; % if mod(i,100) == 0 % PlotData(p,x0new) % pause(0.1) % end x0 = x0new; end