-
Notifications
You must be signed in to change notification settings - Fork 0
/
2D Laplace Equation.py
112 lines (71 loc) · 2.71 KB
/
2D Laplace Equation.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
import numpy
from matplotlib import pyplot, cm
from mpl_toolkits.mplot3d import Axes3D
lineSingle = '------------------------------------------------'
print("Solving Laplace Equation using Finite Difference Method\n")
#Function for plotting initial & steady state solution
def plot2D(x, y, p):
fig = pyplot.figure(figsize=(11,7),dpi=100)
ax = fig.gca(projection='3d')
X,Y=numpy.meshgrid(x,y) #Generating 2D Mesh
surf = ax.plot_surface(X,Y,p[:],rstride=1,cstride=1,cmap=cm.viridis,linewidth=0,antialiased=False)
ax.set_xlim(0, 2)
ax.set_ylim(0, 1)
ax.set_xlabel('X Spacing')
ax.set_ylabel('Y Spacing')
ax.set_zlabel('Velocity')
ax.view_init(30, 225)
#Function for solving Laplace Equation
def laplace2d(p, y, dx, dy, residual_target):
residual = 1 #innitial error
pm = numpy.empty_like(p)
iteration = 0
while residual > residual_target: #Convergence Criteria
pn = p.copy()
p[1:-1,1:-1] = ((dy**2*(pn[2:,1:-1]+pn[0:-2,1:-1]) +
dx**2*(pn[1:-1,2:]+pn[1:-1,0:-2]))/
(2*(dx**2 + dy**2)))
#Boundary Condition
p[:,0] = 0 # p = 0 @ x = 0
p[:,-1] = y # p = y @ x = 2
p[0,:] = p[1,:] # dp/dy = 0 @ y = 0
p[-1,:] = p[-2,:] # dp/dy = 0 @ y = 1
residual = (numpy.sum(numpy.abs(p[:]) - numpy.abs(pn[:]))/
numpy.sum(numpy.abs(pn[:])))
iteration += 1
print('number of iteration :', iteration)
return p
nx = 31 #Grid Point Along X
ny = 31 #Grid Point Along Y
#Grid Spacing
dx = 2 / (nx-1)
dy = 2 / (ny-1)
#initial condition
p = numpy.zeros((ny,nx)) # innitial guess, pressure everywhere is 0
#array
x = numpy.linspace(0,2,nx)
y = numpy.linspace(0,1,ny)
#Innitial Condition
p[:,0] = 0 # p = 0 @ x = 0
p[:,-1] = y # p = y @ x = 2
p[0,:] = p[1,:] # dp/dy = 0 @ y = 0
p[-1,:] = p[-2,:] # dp/dy = 0 @ y = 1
#Calling function to plot initial state solution
plot2D(x,y,p)
print(lineSingle)
print("Plotting Innitial Solution")
print(lineSingle)
pyplot.show()
print(lineSingle)
print("Calculating Numerical Solution......")
print(lineSingle)
print(lineSingle)
print("Solution Converged!")
print(lineSingle)
p = laplace2d(p, y, dx, dy, 1e-5) #Calling function to start the simulation
print(lineSingle)
print("Plotting Numerical Solution")
print(lineSingle)
#Calling function to plot stead state solution
plot2D(x,y,p)
pyplot.show()