Our solver will solve the Navier Stokes equations in two dimensions to calculate the velocities at each grid point and then move to the next time step. The continuity equation will be used as a governing equation for pressure by some manipulation. Hence the we will solve the following equations in theory:
$$
\begin{equation}
    \cfrac{\partial u}{\partial x} +\cfrac{\partial v}{\partial y}=0
\end{equation}
\\
\begin{equation}
    \frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y} = -\frac{1}{\rho}\frac{\partial p}{\partial x}+\nu \left(\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2} \right)
\end{equation}\\
\begin{equation}
    \frac{\partial v}{\partial t}+u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial y} = -\frac{1}{\rho}\frac{\partial p}{\partial y}+\nu\left(\frac{\partial^2 v}{\partial x^2}+\frac{\partial^2 v}{\partial y^2}\right) 
\end{equation}
$$
But since we do not have an explicit governing equation for pressure, we will take the divergence of the momentum equations, add them and use the continuity equation to eliminate the $\cfrac{\partial u}{\partial x} +\cfrac{\partial v}{\partial y}$ terms in the equations, to get the pressure poisson equation, given by,
$$
\begin{equation}
    \frac{\partial^2 p}{\partial x^2}+\frac{\partial^2 p}{\partial y^2} = -\rho\left(\frac{\partial u}{\partial x}\frac{\partial u}{\partial x}+2\frac{\partial u}{\partial y}\frac{\partial v}{\partial x}+\frac{\partial v}{\partial y}\frac{\partial v}{\partial y} \right)
\end{equation}
$$
This equation essentially ensures that continuity is satisfied at the current time step and also couples the pressure and velocities together. 


## Constants and initial conditions


```python
import numpy as np
from matplotlib import pyplot as plt, cm
from solver import *

lx = 5      #length in the x direction
ly = 5      #length in the y direction
nx = 51     #grid points in x direction
ny = 51     #grid points in y direction
t = 10       #final time to calculate to
#nit = 10    #pseudo-time steps in each timestep
dx = lx / (nx - 1)
dy = ly / (ny - 1)
x = np.linspace(0, lx, nx)
y = np.linspace(0, ly, ny)
X, Y = np.meshgrid(x, y)

rho = 1
nu = 0.5
dt = .001 #time step size

nt=int(t/dt)
ds=int(100)

u = np.zeros((ny, nx)) #initial conditions 
v = np.zeros((ny, nx))
p = np.zeros((ny, nx)) 
#b = np.zeros((ny, nx))
```

## Function for calculations of pressure poisson equation


```python
u = np.zeros((ny, nx)) #comment out these lines if you want to run from last time step
v = np.zeros((ny, nx))
p = np.zeros((ny, nx))

vel_mag=(u**2+v**2)**0.5

np.savetxt('p0.csv', p, delimiter=',')
np.savetxt('u0.csv', u, delimiter=',')
np.savetxt('v0.csv', v, delimiter=',')
np.savetxt('Umag0.csv', vel_mag, delimiter=',')

#PRESSURE
pLeft=['N',0]
pRight=['N',0]
pTop=['D',0]
pBottom=['N',0]

pBCs=[pLeft,pRight,pTop,pBottom]
#pressureBC(pLeft,pRight,pTop,pBottom,p,dx,dy)

#VELOCITY

#x-component
uLeft=['D',0]
uRight=['D',0]
uTop=['D',10]
uBottom=['D',0]

uBCs=[uLeft,uRight,uTop,uBottom]
#XVelBC(uLeft,uRight,uTop,uBottom,u,dx,dy)

#y-component
vLeft=['D',0]
vRight=['D',0]
vTop=['D',0]
vBottom=['D',0]

#YVelBC(vLeft,vRight,vTop,vBottom,v,dx,dy)
vBCs=[vLeft,vRight,vTop,vBottom]

# solving
u, v, p = flow_solver(nt, u, v, dt, ds, dx, dy, p, rho, nu, pBCs,uBCs,vBCs)

```


```python
clearResults()
```

## Plotting the results


```python
fig = plt.figure(figsize=(11,7),dpi=100)
# plotting the pressure field as a contour
plt.contourf(X, Y, u, alpha=0.8, cmap=cm.viridis)  
plt.colorbar()
## plotting the pressure field outlines
#plt.contour(X, Y, p, cmap=cm.viridis)  
## plotting velocity field
plt.quiver(X[::lx, ::ly], Y[::lx, ::ly], u[::lx, ::ly], v[::lx, ::ly]) 
plt.xlabel('X')
plt.ylabel('Y')
ax = plt.gca()
ax.set_aspect(1)
plt.show()
```


```python
fig = plt.figure(figsize=(11,7),dpi=100)
# plotting the pressure field as a contour
#plt.contourf(X, Y, u, alpha=0.8, cmap=cm.viridis)
## plotting the pressure field outlines
plt.pcolormesh(X, Y, u, alpha=1, cmap=cm.viridis)
plt.colorbar()
## plotting velocity field
#plt.streamplot(X[::lx, ::2*ly], Y[::lx, ::2*ly], u[::lx, ::2*ly], v[::lx, ::2*ly])
plt.xlabel('X')
plt.ylabel('Y')
ax = plt.gca()
ax.set_aspect(1)
plt.show()
```




```python
import pandas as pd
import matplotlib.animation as animation

fig, ax = plt.subplots(dpi=100)
#levels=np.linspace(0,1.5,10)
def animate(i):
    tmpName='Umag'+str(i*ds)+'.csv'
    Tfile=pd.read_csv(r"Results/"+str(tmpName),header=None)
    Tarr=Tfile.to_numpy()
    # plotting the temperature field
    re = plt.contourf(X, Y, Tarr, 20, alpha=0.8, cmap=cm.viridis)
    #plt.xlabel('X')
    #plt.ylabel('Y')
    ax = plt.gca()
    ax.set_aspect(1)
    return re

anim = animation.FuncAnimation(fig, animate, frames=nt//ds, repeat=False)
anim.save(r'Gifs/cavity_Umag.gif', writer='imagemagick', fps=10)

```
