-
Notifications
You must be signed in to change notification settings - Fork 0
/
2D Linear Convection Equation .py
142 lines (99 loc) · 3.97 KB
/
2D Linear Convection 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
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
from mpl_toolkits.mplot3d import Axes3D #To plot a projected 3D result, make sure that you have added the Axes3D library
import numpy
from matplotlib import pyplot, cm
import time
lineSingle = '------------------------------------------------'
print("Solving 2D Linear Convection Equation using Finite Difference Method\n")
print(lineSingle)
print("METHOD - I")
print("Using Nested FOR Loop")
print(lineSingle)
#meshing
nx = 81 #grid points in x-Direction
ny = 81 #grid points in y-Direction
nt = 80 #number of time step
c = 1 #wave speed kept constant
#grid spacing
dx = 2 / (nx - 1)
dy = 2 / (ny - 1)
CFL = 0.2
dt = CFL*dx #timestep size
x = numpy.linspace(0, 2, nx) #array along x
y = numpy.linspace(0, 2, ny) #array along y
u = numpy.ones((ny, nx))
un = numpy.ones((ny, nx)) #2d temporaray array where we copy our velocity field
#innitial condition
u[int(0.5/dy):int(1/dy+1),int(0.5/dx):int(1/dx+1)] = 2 #Cuboidic Wave Profile
#plotting innitial condition
print(lineSingle)
print("Plotting Innitial Solution: Cuboidic Wave Profile")
print(lineSingle)
fig = pyplot.figure(figsize=(11, 7), dpi=100) #innitilize plot window
ax = fig.gca(projection = '3d') #defining axis is 3d
X,Y = numpy.meshgrid(x, y) #Generating 2D Mesh
#assign plot window the axes label ax, specifies its 3d projection plot
#plot_surface is regular plot command but it takes a grid of x and y for data point position
surf = ax.plot_surface(X, Y, u[:], cmap=cm.viridis)
ax.set_title('Innitial Velocity Field')
ax.set_xlabel('X Spacing')
ax.set_ylabel('Y Spacing')
ax.set_zlabel('Velocity')
pyplot.show()
start1 = time.time() #calculating solving time
start1 = numpy.around(start1, decimals = 2)
print(lineSingle)
print("Calculating Numerical Solution......")
print(lineSingle)
for n in range(nt + 1): #time marching
un = u.copy()
row, col = u.shape
#space marching
for j in range(1,row):
for i in range(1,col):
#Backward Difference Scheme
u[j,i] = (un[j,i] - (c * dt/dx * (un[j,i] - un[j,i-1])) - (c * dt/dx * (un[j,i] - un[j-1,i])))
#Boundary Conditions, U = 1 for x = 0,2 & Y = 0,2
u[0, :] = 1
u[-1, :] = 1
u[:, 0] = 1
u[:, -1] = 1
print("Solving time with Nested FOR loop: %.4s seconds"% (time.time() - start1))
print(lineSingle)
print("Plotting Numerical Solution")
print(lineSingle)
fig = pyplot.figure(figsize=(11,7), dpi=100)
ax = fig.gca(projection = '3d')
surf2 = ax.plot_surface(X, Y, u[:], cmap=cm.viridis)
ax.set_title('Method - I:Using Nested FOR Loop')
ax.set_xlabel('X Spacing')
ax.set_ylabel('Y Spacing')
ax.set_zlabel('Velocity')
pyplot.show()
#using array operation
print(lineSingle)
print("METHOD - II")
print("Using ARRAYS Operation")
print(lineSingle)
print(lineSingle)
print("Calculating Numerical Solution......")
print(lineSingle)
start2 = time.time()
for n in range(nt + 1):
un = u.copy() #itterative process
u[1:, 1:] = (un[1:,1:] - (c * dt / dx * (un[1:,1:] - un[1:,:-1])) - (c * dt / dy * (un[1:,1:] - un[:-1,1:]))) #Backward Difference Scheme
u[0, :] = 1
u[-1, :] = 1
u[:, 0] = 1
u[:, -1] = 1
print("Solving time with Array Operation: %.4s seconds"% (time.time() - start2))
print(lineSingle)
print("Plotting Numerical Solution")
print(lineSingle)
fig = pyplot.figure(figsize=(11,7), dpi = 100)
ax = fig.gca(projection = '3d')
surf3 = ax.plot_surface(X, Y, u[:], cmap = cm.viridis)
ax.set_title('Method - II: Using ARRAYS Operation')
ax.set_xlabel('X Spacing')
ax.set_ylabel('Y Spacing')
ax.set_zlabel('Velocity')
pyplot.show()