2D Poisson Equation

Problem Description

$$\frac{\partial^2{p}}{\partial^2{x}}+\frac{\partial^2{p}}{\partial^2{y}} = c$$

Problem Formulation

Input Data:

  • niter = 51 - maximum number of iterations
  • nx = 21 - number of x spatial points
  • ny = 11 - number of y spatial points
  • xmax = 2
  • ymax = 1
  • infinity = 100 - a large number

ICs

  • $\forall{(n,x,y)}$: $n=0 \rightarrow p =0$

BCs

  • Dirichlet:

    • $\forall{(n,y)}$: $x=0, x=2, y=0, y=1\rightarrow p=0$

  • Neumann:

    • $\forall{(n,x)}$ if $y=0 \rightarrow \frac{\partial{p}}{\partial{y}}=0$
    • $\forall{(n,x)}$ if $y=1 \rightarrow \frac{\partial{p}}{\partial{y}}=0$

Output Data:

  • $\forall{(x, y, n)}p = ?$

Algorithm Design

Space-time Discretization

  • $i \rightarrow$ index of lattice in x
  • $j \rightarrow$ index of lattice in y
  • $n \rightarrow$ index of iterations
  • $m \rightarrow$ index of oddd numbers in analytical solution

Numerical Scheme

  • 1st order FD in time
  • 2nd order CD in space

Algorithm Implementation

In [1]:
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.cm as cm
In [2]:
def laplace(error_target, niter, nx, xmax, ymax):
    h = dy = dx = xmax/(nx-1)
    
    # Compute ny from dy and ymax - ny should be an integer
    ny = int(ymax/dy) + 1
    
    # Initialise data structures:
    p = b = np.zeros(((nx,ny,niter)))
    x = np.zeros(nx)
    y = np.zeros(ny)
    jpos = np.zeros(ny)
    jneg = np.zeros(ny)

    # X Loop
    for i in range(0,nx):
        x[i] = i*dx

    # Y Loop
    for j in range(0,ny):
        y[j] = j*dy

    # Initial conditions
    p[:,:,0] = 0

    # Dirichlet Boundary Conditions:

    # p boundary, left:
    p[0,:,:] = 0

    # p boundary, right:
    for j in range(0,ny):
        p[nx-1,j,:] = y[j] 

    # von Neumann Boundary Conditions:

    # Values for correction at boundaries:
    for j in range(0,ny):
        jpos[j] = j + 1
        jneg[j] = j - 1

    # Set Reflection: 
    jpos[ny-1] = ny-2 
    jneg[0] = 1 
    
    while True: 
        for n in range(0,niter-1):
            for i in range(1,nx-1):
                for j in range(0,ny):  
                    p[i,j,n+1] = 0.25*( p[i+1,j,n]+p[i-1,j,n]+p[i,int(jpos[j]),n]+p[i,int(jneg[j]),n])
            
        error = (np.sum(np.abs(p[i,j,n+1])-np.abs(p[i,j,n]) ) / 
                    np.sum(np.abs(p[i,j,n+1]) ))
        
        if(error < error_target):
            print ("Completed! with  n = " + str(n))
            break
        break
            
    return p, x, y
In [4]:
p1, x1, y1 = laplace(1.0e-6, 5000, 51, 2.0, 1.0)

Completed! with  n = 4998
In [5]:
def plot_3D(p,x,y,time,title,label):
    """
    Plots the 2D velocity field
    """

    import matplotlib.pyplot as plt
    from mpl_toolkits.mplot3d import Axes3D
    fig=plt.figure(figsize=(11,7),dpi=100)
    ax=fig.gca(projection='3d')
    ax.set_xlabel('x (m)')
    ax.set_ylabel('y (m)')
    ax.set_zlabel(label)
    ax.set_xlim(0,2)
    ax.set_ylim(0,1)
    ax.view_init(30,225)
    Y,X=np.meshgrid(y,x) #note meshgrid uses y,x not x,y!!!
    surf=ax.plot_surface(X,Y,p[:,:,time], rstride=1, cstride=1, cmap=cm.viridis,
            linewidth=0, antialiased=False)
    plt.title(title)
    plt.show()
In [6]:
plot_3D(p1,x1,y1,0,'Figure 1: Initial Conditions: n=0, nx=51','p (Pa)')
plot_3D(p1,x1,y1,4998,'Figure 2: Final Conditions: n=4998, nx=51','p (Pa)')
In [9]:
def poisson(error_target, niter, nx, xmax, ymax):
    """
    Returns the velocity field and distance for 2D Poisson Equation
    """
    # Increments in x and y are the same - h:
    h = dy = dx = xmax/(nx-1)
    
    # Compute ny from dy and ymax - ny should be an integer
    ny = int(ymax/dy) + 1
    
    # Initialise data structures:
    p = np.zeros(((nx,ny,niter)))
    b = np.zeros((nx,ny))
    x = np.zeros(nx)
    y = np.zeros(ny)

    b[int(nx/4),int(ny/4)] = 100
    b[int(3*nx/4),int(3*ny/4)] = -100 
    
    # X Loop
    for i in range(0,nx):
        x[i] = i*dx

    # Y Loop
    for j in range(0,ny):
        y[j] = j*dy

    # Initial conditions
    p[:,:,0] = 0

    # Dirichlet Boundary Conditions:
    p[0,:,:] = p[nx-1,:,:] = p[:,0,:] = p[:,ny-1,:] = 0

    #Loop
    while True: 
        for n in range(0,niter-1):
            for i in range(1,nx-1):
                for j in range(1,ny-1): 
                    p[i,j,n+1] = 0.25*( p[i+1,j,n]+p[i-1,j,n]+p[i,j+1,n]+p[i,j-1,n] - b[i,j]*h**2)
            
        error = (np.sum(np.abs(p[i,j,n+1])-np.abs(p[i,j,n]) ) / 
                    np.sum(np.abs(p[i,j,n+1]) ))
        
        if(error < error_target):
            print ("Completed! with  n = " + str(n))
            break
        break
            
    return p, x, y
In [10]:
p2, x2, y2 = poisson(1.0e-6, 3000, 51, 2.0, 1.0)

Completed! with  n = 2998
In [11]:
plot_3D(p2,x2,y2,0,'Figure 1: Initial Conditions: n=0, nx=51','p (Pa)')
plot_3D(p2,x2,y2,17,'Figure 2: Final Conditions: n=17, nx=51','p (Pa)')
In [ ]: