%===============================================================
% CodeC1-advChi.m
% A very simple Navier-Stokes solver for a drop falling in a
% rectangular box, using a conservative form of the equations. 
% A first-order explicit projection method and centered in space 
% discretizationa are used. The marker function is advected by 
% an advection diffusion equation. Last edited 7/6/2016
%===============================================================

Lx=1.0;Ly=1.0;gx=0.0;gy=-100.0;  % Domain size and
rho1=1.0; rho2=2.0; m0=0.01;     % physical variables
unorth=0;usouth=0;veast=0;vwest=0;time=0.0; 
rad=0.15;xc=0.5;yc=0.7; % Initial drop size and location

%-------------------- Numerical variables ----------------------
nx=32;ny=32;dt=0.00125;nstep=400; maxit=200;maxError=0.001;beta=1.5;

%-------------------- Zero various arrys -----------------------
u=zeros(nx+1,ny+2);  v=zeros(nx+2,ny+1);  p=zeros(nx+2,ny+2);
ut=zeros(nx+1,ny+2); vt=zeros(nx+2,ny+1); tmp1=zeros(nx+2,ny+2); 
uu=zeros(nx+1,ny+1); vv=zeros(nx+1,ny+1); tmp2=zeros(nx+2,ny+2);
r=zeros(nx+2,ny+2); chi=zeros(nx+2,ny+2);

dx=Lx/nx;dy=Ly/ny;                          % Set the grid 
for i=1:nx+2; x(i)=dx*(i-1.5);end; for j=1:ny+2; y(j)=dy*(j-1.5);end;

%-------------------- Initial Conditions -----------------------
r=zeros(nx+2,ny+2)+rho1;                       % Set density
for i=2:nx+1,for j=2:ny+1;                     % for the domain and the drop
  if((x(i)-xc)^2+(y(j)-yc)^2 <rad^2),r(i,j)=rho2;chi(i,j)=1.0;end, 
end,end
                                          
contour(x,y,chi'),axis equal,axis([0 Lx 0 Ly]);pause(0.01)               

%---------------------- START TIME LOOP ------------------------
for is=1:nstep,is

%--- ADVECT marker using centered difference plus diffusion ----
    chio=chi;
    for i=2:nx+1,for j=2:ny+1
          chi(i,j)=chio(i,j)-(0.5*dt/dx)*(u(i,j)*(chio(i+1,j)...
            +chio(i,j))-u(i-1,j)*(chio(i-1,j)+chio(i,j)))...
            -(0.5* dt/dy)*(v(i,j)*(chio(i,j+1)...
            +chio(i,j))-v(i,j-1)*(chio(i,j-1)+chio(i,j))  )...
            +(m0*dt/dx/dx)*(chio(i+1,j)-2.0*chio(i,j)+chio(i-1,j))...
            +(m0*dt/dy/dy)*(chio(i,j+1)-2.0*chio(i,j)+chio(i,j-1));
   end,end
       
%-------------------- Update the density -----------------------
    ro=r;
    for i=1:nx+2,for j=1:ny+2
      r(i,j)=rho1+(rho2-rho1)*chi(i,j);
    end,end 

%------------- Set tangential velocity at boundaries -----------	     
    u(1:nx+1,1)=2*usouth-u(1:nx+1,2);u(1:nx+1,ny+2)=2*unorth-u(1:nx+1,ny+1);
    v(1,1:ny+1)=2*vwest-v(2,1:ny+1);v(nx+2,1:ny+1)=2*veast-v(nx+1,1:ny+1);

%-------------- Find the predicted velocities ------------------	     
    for i=2:nx,for j=2:ny+1      % Temporary u-velocity
      ut(i,j)=(2.0/(r(i+1,j)+r(i,j)))*(0.5*(ro(i+1,j)+ro(i,j))*u(i,j)+ dt* (...
      -(0.25/dx)*(ro(i+1,j)*(u(i+1,j)+u(i,j))^2-ro(i,j)*(u(i,j)+u(i-1,j))^2)...
      -(0.0625/dy)*( (ro(i,j)+ro(i+1,j)+ro(i,j+1)+ro(i+1,j+1))*...
                                       (u(i,j+1)+u(i,j))*(v(i+1,j)+v(i,j)) ...
      -(ro(i,j)+ro(i+1,j)+ro(i+1,j-1)+ro(i,j-1))*(u(i,j)...
                                       +u(i,j-1))*(v(i+1,j-1)+v(i,j-1)))...
      +m0*((u(i+1,j)-2*u(i,j)+u(i-1,j))/dx^2+ (u(i,j+1)-2*u(i,j)+u(i,j-1))/dy^2)...
                               + 0.5*(ro(i+1,j)+ro(i,j))*gx ) );
    end,end

    for i=2:nx+1,for j=2:ny       % Temporary v-velocity 
      vt(i,j)=(2.0/(r(i,j+1)+r(i,j)))*(0.5*(ro(i,j+1)+ro(i,j))*v(i,j)+ dt* (...     
      -(0.0625/dx)*( (ro(i,j)+ro(i+1,j)+ro(i+1,j+1)+ro(i,j+1))*...
                                        (u(i,j)+u(i,j+1))*(v(i,j)+v(i+1,j)) ...
                  - (ro(i,j)+ro(i,j+1)+ro(i-1,j+1)+ro(i-1,j))*...
                                    (u(i-1,j+1)+u(i-1,j))*(v(i,j)+v(i-1,j)) )...                                 
      -(0.25/dy)*(ro(i,j+1)*(v(i,j+1)+v(i,j))^2-ro(i,j)*(v(i,j)+v(i,j-1))^2 )...
      +m0*((v(i+1,j)-2*v(i,j)+v(i-1,j))/dx^2+(v(i,j+1)-2*v(i,j)+v(i,j-1))/dy^2)...
                               + 0.5*(ro(i,j+1)+ro(i,j))*gy  ) );    
    end,end   

%------------------ Solve the Pressure Equation ----------------    
    rt=r; lrg=1000;   % Compute source term and the coefficient for p(i,j)
    rt(1:nx+2,1)=lrg;rt(1:nx+2,ny+2)=lrg;
    rt(1,1:ny+2)=lrg;rt(nx+2,1:ny+2)=lrg;

    for i=2:nx+1,for j=2:ny+1
      tmp1(i,j)= (0.5/dt)*( (ut(i,j)-ut(i-1,j))/dx+(vt(i,j)-vt(i,j-1))/dy );
      tmp2(i,j)=1.0/( (1./dx)*(1./(dx*(rt(i+1,j)+rt(i,j)))+   ...
                               1./(dx*(rt(i-1,j)+rt(i,j))) )+ ...
                      (1./dy)*(1./(dy*(rt(i,j+1)+rt(i,j)))+   ...
                               1./(dy*(rt(i,j-1)+rt(i,j))) )   );
    end,end

    for it=1:maxit	               % Solve for pressure by SOR
      oldArray=p;
      for i=2:nx+1,for j=2:ny+1
        p(i,j)=(1.0-beta)*p(i,j)+beta* tmp2(i,j)*(        ...
        (1./dx)*( p(i+1,j)/(dx*(rt(i+1,j)+rt(i,j)))+      ...
                  p(i-1,j)/(dx*(rt(i-1,j)+rt(i,j))) )+ ...
        (1./dy)*( p(i,j+1)/(dy*(rt(i,j+1)+rt(i,j)))+      ...
                  p(i,j-1)/(dy*(rt(i,j-1)+rt(i,j))) ) - tmp1(i,j));
      end,end
      if max(max(abs(oldArray-p))) <maxError, break, end
    end
                                      
    for i=2:nx,for j=2:ny+1   % Correct the u-velocity 
      u(i,j)=ut(i,j)-dt*(2.0/dx)*(p(i+1,j)-p(i,j))/(r(i+1,j)+r(i,j));
    end,end
      
    for i=2:nx+1,for j=2:ny   % Correct the v-velocity
      v(i,j)=vt(i,j)-dt*(2.0/dy)*(p(i,j+1)-p(i,j))/(r(i,j+1)+r(i,j));
    end,end

%------------------ Plot the results ---------------------------
  time=time+dt                   % plot the results
  uu(1:nx+1,1:ny+1)=0.5*(u(1:nx+1,2:ny+2)+u(1:nx+1,1:ny+1));
  vv(1:nx+1,1:ny+1)=0.5*(v(2:nx+2,1:ny+1)+v(1:nx+1,1:ny+1));
  for i=1:nx+1,xh(i)=dx*(i-1);end;     for j=1:ny+1,yh(j)=dy*(j-1);end
  hold off,contour(x,y,r'),axis equal,axis([0 Lx 0 Ly]);
  hold on;quiver(xh,yh,uu',vv','r');
  pause(0.01)

end                  % End of time step
