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

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

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
        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
        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

  • $\begingroup$ Have you attempted to use a simpler time stepper? Like Heun's method which is a second order method and is basically two lines. physics.utah.edu/~detar/phys6720/handouts/ode/ode/node3.html $\endgroup$ Commented Apr 27, 2015 at 23:46
  • $\begingroup$ I have not. I was using the Euler step as my "simpler time stepper" as it's just the first step of the RK3 with the 1/3 removed. $\endgroup$
    – cbell14
    Commented Apr 28, 2015 at 2:50

1 Answer 1


I got the code working, the full solution can be found in my github repo https://github.com/cbell14/CFD_Class. The file is Taylor_Green_Vortices.ipynb.

There were 2 main issues in the code posted above. The first was in my implementation of the Dirichlet Boundary Condition (BC) in the A matrix. Implementing the BC in the A matrix does not work as it effects more than just the point you wish to make Dirichlet. Instead after solving for the pressure, I subtract the value at the middle from the entire field. This is necessary because all the boundary conditions are Neumann which yields an infinite number of solutions so the pressure must be anchored to a particular value.

The second issue was that I was incorporating the dx and dy terms from the denominator of the numerical scheme in the A matrix. This was proving to be difficult, so I incorporated them into the RHS of the linear system (my b matrix). This solved all my issues. Hope this helps anyone else who has a similar issue.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.