from IPython.display import HTML
s="""<h1>2D Laplace Equation</h1></br><div id="toc"></div>
<script src="https://kmahelona.github.io/ipython_notebook_goodies/ipython_notebook_toc.js" defer></script>""";
h=HTML(s); h
What is the 2D pressure field for the Laplace equation when the initial conditions are zero everywhere and the x-boundary conditions are 0 and 1 and the y-boundary conditions are von Neumann (zero gradient)?
The Laplace Equation for pressure is described as follows:
This problem has no temporal component, so we use a number of iterations.
niter
= 51 (maximum number of iterations)nx
= 21 (number of x spatial points)ny
= 11 (number of y spatial points)xmax
= 2ymax
= 1infinity
= 100 (a large number)To use this as an iterative formula, $p_{i,j}^{n} = p_{i,j}^{n+1}$ This looks a bit dodgy but is the only thing that makes sense, unless we equate the RHS to forward differencing
niter = 51
nx = 21
xmax = 2
ny = 11
ymax = 1
dx = xmax/(nx-1)
dy = ymax/(ny-1)
infinity = 100
# Dirichlet Boundary Conditions
# p boundary, left:
p(0,:,:) = 0
# p boundary, right:
p(nx-1,:,:) = ymax
# von Neumann Boundary Conditions
# Values for correction at boundaries:
for j between 0 and ny-1
jpos(j) = j + 1
jneg(j) = j - 1
# Set Reflection:
# j: -1 0, 1,... ny-2, ny-1, ny
# jpos: start => => end
# jneg: start => => end
jpos(ny-1) = ny-2 # i.e. ny = ny-2
jneg(0) = 1 # i.e. -1 = 1
# Initial Conditions
p(:,:,0) = 0
# Numerical Computation
for n between 0 and niter-1
for i between 1 and nx-1
for j between 0 and ny
p(i,j,n+1) = { dy^2 * [ p(i+1,j,n)+p(i-1,j,n) ]
+ dx^2 * [ p(i,jpos(j),n)+p(i,jneg(j),n) ] }
/ {2(dx^2 + dy^2)}
# Analytical Solution
p_coefficient(i,j,1) = { 4 * [sinh(1*pi*x(i)]*[cos(1*pi*y(j)] }
/ { [(1*pi)**2]*sinh(2*pi*1) }
for m between 3 and infinity in odd steps
for i between 0 and nx-1
for j between 0 and ny-1
p_coefficient(i,j) += {[sinh(m*pi*x(i)]*[cos(m*pi*y(j)] }
/ { [(m*pi)**2]*sinh(2*pi*m) }
for i between 0 and nx-1
p_analytical(i,j)= (x(i) / 4) - 4*p_coefficient(i,j)
def laplace_equation_numerical(niter, nx, ny, xmax, ymax):
"""
Returns the velocity field and distance for 2D linear convection
"""
# Increments:
dx = xmax/(nx-1)
dy = ymax/(ny-1)
# Initialise data structures:
import numpy as np
p = 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-1):
jpos[j] = j + 1
jneg[j] = j - 1
# Set Reflection:
jpos[ny-1] = ny-2
jneg[0] = 1
# Loop
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] = (( dy**2 * ( p[i+1,j,n]+p[i-1,j,n] )
+ dx**2 * ( p[i,jpos[j],n]+p[i,jneg[j],n] ) )
/ (2*(dx**2 + dy**2)))
return p, x, y
p1, x1, y1 = laplace_equation_numerical(100, 51, 51, 2.0, 1.0)
def plot_2D(p,x,nt,title):
"""
Plots the 1D velocity field
"""
import matplotlib.pyplot as plt
import matplotlib.cm as cm
plt.figure()
ax=plt.subplot(111)
colour=iter(cm.rainbow(np.linspace(0,20,nt)))
for n in range(0,nt,20):
c=next(colour)
ax.plot(x,p[:,-1,n],linestyle='-',c=c,label='n='+str(n)+' numerical')
box=ax.get_position()
ax.set_position([box.x0, box.y0, box.width*1.5,box.height*1.5])
ax.legend( bbox_to_anchor=(1.02,1), loc=2)
plt.xlabel('x (m)')
plt.ylabel('p (Pa)')
plt.ylim([0,1.0])
plt.xlim([0,2.0])
plt.title(title)
plt.show()
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)
plt.title(title)
plt.show()
plot_3D(p1,x1,y1,0,'Figure 1: Initial Conditions: n=0, niter=100, nx=51, ny=51','p (Pa)')
plot_3D(p1,x1,y1,99,'Figure 2: Final Conditions: n=99, niter=100, nx=51, ny=51','p (Pa)')
plot_2D(p1,x1,99,'Figure 3: At Far Neumann Boundary: n=99, niter=100, nx=51, ny=51')
So now we need a target Relative Error, so choose $1 \times 10^{-2}$
This changes the pseudo-code a bit.
For safety we can still include the maximum number of iterations, but maybe increase it so that the while loop determines what happens?
We also need a new plotting function, as we need to step through less iterations in 2D.
error_target=1*10^-2
niter = 500
# Numerical Computation
while error is greater than error_target
for n between 0 and niter-1
for i between 1 and nx-1
for j between 1 and ny-1
p(i,j,n+1) = { dy^2 * [ p(i+1,j,n)+p(i-1,j,n) ] +
dx^2 * [ p(i,jpos(j),n)+p(i,jneg(j),n) ] }
/ {2(dx^2 + dy^2)}
error = [ sum(abs(p(i,j,n+1)-abs(p(i,j,n)) ] / sum(abs(p(i,j,n+1)))
def laplace_equation_numerical_2(error_target, niter, nx, ny, xmax, ymax):
"""
Returns the velocity field and distance for 2D linear convection
"""
# Increments:
dx = xmax/(nx-1)
dy = ymax/(ny-1)
# Initialise data structures:
import numpy as np
p = 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] = (( dy**2 * ( p[i+1,j,n]+p[i-1,j,n] )
+ dx**2 * ( p[i,jpos[j],n]+p[i,jneg[j],n] ) )
/ (2*(dx**2 + dy**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]) ))
print "n = " + str(n) + " completed"
if(error < error_target):
break
break
return p, x, y
p2, x2, y2 = laplace_equation_numerical_2(1.0e-2, 500, 51, 26, 2.0, 1.0)
n = 0 completed n = 1 completed n = 2 completed n = 3 completed n = 4 completed n = 5 completed n = 6 completed n = 7 completed n = 8 completed n = 9 completed n = 10 completed n = 11 completed n = 12 completed n = 13 completed n = 14 completed n = 15 completed n = 16 completed
def plot_2D_2(p,x,nt,title):
"""
Plots the 1D velocity field
"""
import matplotlib.pyplot as plt
import matplotlib.cm as cm
plt.figure()
ax=plt.subplot(111)
colour=iter(cm.rainbow(np.linspace(0,1,nt+1)))
for n in range(0,nt+1,1):
c=next(colour)
ax.plot(x,p[:,-1,n],linestyle='-',c=c,label='n='+str(n)+' numerical')
box=ax.get_position()
ax.set_position([box.x0, box.y0, box.width*1.5,box.height*1.5])
ax.legend( bbox_to_anchor=(1.02,1), loc=2)
plt.xlabel('x (m)')
plt.ylabel('p (Pa)')
plt.ylim([0,1.0])
plt.xlim([0,2.0])
plt.title(title)
plt.show()
plot_3D(p2,x2,y2,0,'Figure 1: Initial Conditions: n=0, niter=25, n=51, nx=51, ny=51','p (Pa)')
plot_3D(p2,x2,y2,16,'Figure 2: Final Conditions: n=16, niter=25, nx=51, ny=51','p (Pa)')
plot_2D_2(p1,x1,16,'Figure 3: At Far Neumann Boundary: n=16, niter=25, nx=51, ny=51')
def laplace_equation_analytical(INFINITY, NX, NY, XMAX, YMAX):
from math import pi as PI
import math as ma
# Initialise data structures
p_coefficient = np.zeros((NX, NY))
p_analytical = np.zeros((NX, NY))
x = np.zeros(NX)
y = np.zeros(NY)
# Constants
DX = XMAX/(NX-1)
DY = YMAX/(NY-1)
# X Loop
for i in range(0,NX):
x[i] = i*DX
# Y Loop
for j in range(0,NY):
y[j] = j*DY
# Analytical solution
for i in range(0, NX):
for j in range(0, NY):
p_coefficient[i,j] = ( (ma.sinh(1.0*PI*x[i]) * ma.cos(1.0*PI*y[j]) )
/ ( ((1.0*PI)**2.0) * ma.sinh(2.0*PI*1.0) ) )
for m in range(3, INFINITY, 2):
for i in range(0, NX):
for j in range(0, NY):
p_coefficient[i,j] += ( (ma.sinh(m*PI*x[i]) * ma.cos(m*PI*y[j]) )
/ ( ((m*PI)**2.0) * ma.sinh(2.0*PI*m) ) )
for i in range (0, NX):
for j in range(0, NY):
p_analytical[i,j]= (x[i] / 4.0) - (4.0 * p_coefficient[i,j])
return p_analytical, x, y
p3, x3, y3 = laplace_equation_analytical(100, 51, 26, 2.0, 1.0)
def plot_3D_2(p,x,y,title):
"""
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('p (Pa)')
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[:,:], rstride=1, cstride=1)
plt.title(title)
plt.show()
plot_3D_2(p3,x3,y3,'Figure 1: Analytical Solution: nx=51, ny=26')
where: $r = {{\Delta t} \over {h^2}}$
So the pressure is just the average of it's neighbours:
$$ 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})$$def laplace_equation_numerical_3(error_target, niter, nx, ny, xmax, ymax):
"""
Returns the velocity field and distance for 2D linear convection
"""
# Increments:
dx = xmax/(nx-1)
dy = ymax/(ny-1)
# Initialise data structures:
import numpy as np
p = 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,jpos[j],n]+p[i,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 "n = " + str(n) + " completed"
break
break
return p, x, y
p4, x4, y4 = laplace_equation_numerical_3(1.0e-6, 5000, 51, 26, 2.0, 1.0)
n = 3129 completed
plot_3D(p4,x4,y4,3129,'Figure 2: Final Conditions: n=3129, niter=3130, nx=51, ny=26','p (Pa)')
def plot_diffusion(p4,p3,x4,NY,TIME,TITLE):
"""
Plots the 1D velocity field
"""
import matplotlib.pyplot as plt
import matplotlib.cm as cm
plt.figure()
ax=plt.subplot(111)
colour=iter(cm.rainbow(np.linspace(0,5,NY)))
for j in range(0,NY,5):
c=next(colour)
ax.plot(x4,p4[:,j,TIME],'ko', markerfacecolor='none', alpha=0.5, label='j='+str(j)+' numerical')
ax.plot(x4,p3[:,j],linestyle='-',c=c,label='j='+str(j)+' analytical')
box=ax.get_position()
ax.set_position([box.x0, box.y0, box.width*1.5,box.height*1.5])
ax.legend( bbox_to_anchor=(1.02,1), loc=2)
plt.xlabel('x (m)')
plt.ylabel('p (Pa)')
#plt.ylim([0,8.0])
#plt.xlim([0,2.0*PI])
plt.title(TITLE)
plt.show()
plot_diffusion(p4,p3,x4,26,3129,'Analytical vs Numerical')
Laplace operator is typical for diffusion, diffusion is an isotropic phenomenon.
It has to be discretised with a central difference - numerics must be consistent with the physics that it represents.
Convection is anisotropic so numerical scheme must be directionally biased.
Second order CD in both x and y is the most widely used numerical scheme for $ \nabla^2 $
Also known as "five point difference operator"
This scheme is an iterative method for a steady state (artificial time variable)
Let $ \Delta x = \Delta y = h$ then
Linear system of equations in pentadiagonal coefficient matrix
Equivalent to Point Jacobi method