% Navier Stokes - du/dt + (u.D)u = -(1/rho)Dp + visLPu
% D = divergence LP = Laplacian  u = {u, v}
% in continuous casep satisifes Poisson equation -
%     -1/rho*(d2p/dx2 + d2p/dy2) = (du/dx)^2 + 2du/dydv/dx + (dv/dy)^2
% discretize in time
% u/n+1 = n/n + dt[-u/n.Du/n - (1/rho)Dp/n + visLPu/n]
% set Du/n+1 = 0 to get Poisson eq for p
% D.p/n+1 = D.u/n + dt[-D.u/n.Du/n - (1/rho)D^2p/n + visLP(Du/n)]
% set D.u/n+1 = 0 to get Poisson eq for p
% D^2p/n+1 = rhoD.u/n/dt + [-rhoD.u/n.Du/n  + rho*visLP(Du/n)]
% discretize Poisson for p
%(p/i+1,j - 2p/i,j + p/i-1/j)/dx^2 + (p/i,j+1 - 2p/i,j + p/i/j-1)/dy^2
%    = (u/i+1,j - u/i-1,j)/2dx)^2 
%     + (u/i+1,j - u/i-1,j))/2dy*(v/i+1,j - v/i-1,j))/2dx 
%     + (v/i+1,j - v/i-1,j)/2dy)^2 
%     + (rho/dt)((u/i+1,j - u/i-1,j)/2dx) + (v/i+1,j - v/i-1,j)/2dy)                                  
% ME 702 Video Lesson 10-11
% IC - u,v,p = 0 everywhere
% BC u=1 @ y=2
% u,v = 0  x=0,2 & y=0
% p = 0 @ y = 2, dp/dy = 0 @ y = 0 & dp/dx = 0 @ x = 0,2



% from step 8?
% un+1/i.j = un/i,j + dt*vis*(un/i+1,j - 2un/i,j + un/i-1,j)/dx^2 + dt*vis*(un/i,j+1 - 2un/i,j + un/i,j-1)/dy^2  

%
nx = 10; ny = 10;
nt = 100; nit=10; 
rho = 1;
nu = .1;
dt = .001;

dx = 2/(nx-1); dy = 2/(ny-1);

x=0:dx:2; y=0:dy:2;
[Y,X] = meshgrid(y,x);


u=zeros(nx,ny);
v=zeros(nx,ny);
p=zeros(nx,ny);
b=zeros(nx,ny);

for n=1:nt   % outer time step loop
    
    for i=2:nx-1    % fill b array (see step 10 for Poission eq)
        for j=2:ny-1
             b(i,j) = rho*(1/dt*((u(i+1,j)-u(i-1,j))/(2*dx)+ (v(i,j+1)-v(i,j-1))/(2*dy)) ...
             -  ((u(i+1,j)-u(i-1,j))/(2*dx))^2 ...
              - 2*((u(i,j+1)-u(i,j-1))/(2*dy)*(v(i+1,j)-v(i-1,j))/(2*dx)) ...
              -  ((v(i,j+1)-v(i,j-1))/(2*dy))^2);
              
       
        end
    end
    

    for iit=1:nit   % solve poission eq for p

       for i = 2:nx-1
           for j = 2:ny-1
               p(i,j) = [(p(i+1,j) + p(i-1,j))*dy^2 + (p(i,j+1) + p(i,j-1))*dx^2] ...
                          /  (2*(dx^2+dy^2)) ...
               - dx^2*dy^2/(2*(dx^2+dy^2))*b(i,j);
           end
       end
       p(nx,:) = p(nx-1,:);            % dp/dx = 0 at x = 2       
       p(1,:) = p(2,:);                % dp/dx = 0 at x = 0     
       p(:,1) = p(:,2);                % dp/dy = 0 at y = 0
       p(:,ny) = 0;                    % p = 0 at y = 2
             
    end

    un=u;
    vn=v;
    for i = 2:nx-1
        for j = 2:ny-1
          u(i,j) = un(i,j) ...
                  - un(i,j)*dt/dx*(un(i,j) - un(i-1,j)) ...
                  - vn(i,j)*dt/dx*(un(i,j) - un(i,j-1)) ...
                  - dt/(2*rho*dx)*(p(i+1,j) - p(i-1,j)) ...
                  + nu*(dt/dx^2*(un(i+1,j) - 2*un(i,j) + un(i-1,j)) ...
                  +    dt/dy^2*(un(i,j+1) - 2*un(i,j) + un(i,j-1)));
                      
          v(i,j) = vn(i,j) ...
                  - un(i,j)*dt/dx*(vn(i,j) - vn(i-1,j)) ...
                  - vn(i,j)*dt/dy*(vn(i,j) - vn(i,j-1)) ...
                  - dt/(2*rho*dy)*(p(i,j+1) - p(i,j-1)) ...
                  + nu*(dt/dx^2*(vn(i+1,j) - 2*vn(i,j) + vn(i-1,j)) ...
                  +    dt/dy^2*(vn(i,j+1) - 2*vn(i,j) + vn(i,j-1)));                       
        end
    end
   
    
    u(1,:)=0;     % u = 0 @ x = 0
    u(:,1) = 0;   % u = 0 @ y = 0 
    u(:,ny) = 1;  % u = 1 @ y = 2
    
    v(1,:)=0;     % v = 0 @ x = 0 
    v(nx,:)=0;    % v = 0 @ x = 2
    v(:,1) = 0;   % v = 0 @ y = 0 
    v(:,ny) = 0;  % v u = 0 @ y = 2
    
    u(nx,:)=0;    % u = 0 @ x = 2 
end

contourf(X,Y,p)
hold on
colorbar()
contour(X,Y,p)
quiver(X(1:2:nx,1:2:ny),Y(1:2:nx,1:2:ny),u(1:2:nx,1:2:ny),v(1:2:nx,1:2:ny))


hold off