I am trying to write a 2D Navier Stokes solver using an RK3 for time advancement in python. For debugging, I have converted the RK3 to an Euler step for simplicity. Checking my divergence for my predictor velocity I find that the divergence is acceptably close to zero. I then go to my pressure solver and the solution is incorrect. It does not even follow the same pattern as the exact solution. Can anyone see where my error is?
My code isn't the most efficient, but I'm not worried about that currently just trying to get a physically correct answer. The analytical solutions for the u, v, and p are given by $$u = -e^{-2t}\cos x \sin y$$ $$v = -e^{-2t}\sin x \cos y$$ $$p = -\frac{e^{-4t}}{4} \left(\cos 2x + \sin 2y \right)$$
def Poisson_Solver_Neumann(u, v, Nx, Ny, dt, dx, dy, T, xp, yp):
"""Solves the 2D Poisson equation implicitly on a staggered grid
using Neumann Boundary Conditions
Params:
------
u, v 2D array of float, x and y velocities
Nx, Ny float, Number of segments in x and y
dt float, time step size
T float, current time
X, Y 2D array of float, meshgrid
Returns:
-------
ANeum 2D array of float, A matrix with Neumann conditions
f_RHSn 1D array of float, f(x,y) for Neumann conditions
"""
#Building A
ANeum = np.zeros((Nx*Ny,Nx*Ny),dtype=float)
a = Nx * Ny
b = Nx * Ny
c_int = -4 * (Nx * Ny)
c_edge = -3 * (Nx * Ny)
c_corner = -2 * (Nx * Ny)
d = Nx * Ny
e = Nx * Ny
#Set corner points
ANeum[0,0] = c_corner
ANeum[-1,-1] = c_corner
ANeum[Nx*Ny-Ny,Nx*Ny-Ny] = c_corner
ANeum[Ny-1,Ny-1] = c_corner
#Set edges in first block
for j in range(1,Ny-1):
ANeum[j,j] = c_edge
j +=j
#Set edges in last block
for j in range((Nx*Ny)-Ny,(Nx*Ny)-2):
ANeum[j+1,j+1] = c_edge
j +=j
#Set edges along main diagonal except for first block
for j in range(Nx+1,Ny*Nx):
if j % Nx ==0:
ANeum[j-1,j-1] = c_edge
j +=j
#Set edges on main diagonal except for last block
for j in range(Nx,(Ny*Nx)-Nx):
if j % Nx ==0:
ANeum[j,j] = c_edge
j +=j
#Second diagonal above and below diagonal
for j in range(Ny,Ny*Nx):
ANeum[j,j-Ny] = a
ANeum[j-Nx,j] = e
j +=j
#first diagonal below main diagonal
for j in range(1,Ny*Nx):
if j % Ny ==0:
ANeum[j,j-1] = 0
else:
ANeum[j,j-1] = b
j +=j
#first diagonal above main diagonal
for j in range(0,Ny*Nx):
if j % Nx ==0:
ANeum[j-1,j] = 0
else:
ANeum[j-1,j] = d
j +=j
#Main Diagonal
for j in range(0,Ny*Nx):
if ANeum[j,j] ==0:
ANeum[j,j] = c_int
#Setting random point to Dirichlet
ANeum[0,:] = 0
ANeum[0,0] = 1
return ANeum
My Euler step is here:
for t in range(0,nt):
un = np.empty_like(u)
vn = np.empty_like(v)
#copy info from last loop, needed for current loop into variable_n
un = u.copy()
vn = v.copy()
#Pn = p.copy()
T = T_0 + t*dt
####STEP 1
### t -----> t + (1/3)*dt
un, vn = BCs(un, vn)
G1 = Build_G1(un, vn, dx, dy, nu)
G2 = Build_G2(un, vn, dx, dy, nu)
u1 = un + (dt*G1) #put 1/3 back in
v1 = vn + (dt*G2)
u1, v1 = BCs(u1, v1)
u_div = (u1[1:-1,2:-1] - u1[1:-1,1:-2])/dx
v_div = (v1[2:-1,1:-1] - v1[1:-2,1:-1])/dy
#Solve for Pressure Field
b1 = np.zeros((Ny+2,Nx+2), dtype=float)
b1[1:-1,1:-1] = (u1[1:-1,2:-1] - u1[1:-1,1:-2])/dx +\
(v1[2:-1,1:-1] - v1[1:-2,1:-1])/dy
b1 = np.reshape(b1[1:-1,1:-1], Ny*Nx)
A1 = Poisson_Solver_Neumann(u1, v1, Nx, Ny, dt, dx, dy, T, xp, yp)
b1 = b1[:]*(1.0/dt)
temp1 = np.linalg.solve(A1,b1)
p_star = np.reshape(temp1, (Ny,Nx))
F1_Pstar, F2_Pstar = Pressure_Terms(p_star, dx, dy)
#calc predictor velocities
u_star = u1.copy()
v_star = v1.copy()
u_star[1:-1,2:-2] = u1[1:-1,2:-2] - (dt*F1_Pstar)
v_star[2:-2,1:-1] = v1[2:-2,1:-1] - (dt*F2_Pstar)
u = u_star.copy()
v = v_star.copy()
Update: My code is now working for the first time step but subsequently blows up.
EDIT: Some code changes
