-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy path2DTPFA.py
123 lines (91 loc) · 3.44 KB
/
2DTPFA.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
import numpy as np
import matplotlib.pyplot as plt
import sympy as sym
from mpl_toolkits import mplot3d
import math
bottom_left = (0,0)
top_right = (1,1)
num_nodes_x = 12
num_nodes_y = 6
num_nodes = num_nodes_x*num_nodes_y
u_h = np.zeros(num_nodes)
nodes_x, nodes_y = np.meshgrid(np.linspace(bottom_left[0],top_right[0],num=num_nodes_x),np.linspace(bottom_left[1],top_right[1],num=num_nodes_y))
nodes = np.stack([nodes_x,nodes_y],axis=2)
x = sym.symbols('x')
y = sym.symbols('y')
u_fabric = x*y*(1-x)*(1-y)
f = 2*x*(1-x)+2*y*(1-y)
u_fabric_vec = np.zeros(num_nodes)
u_fabric_nodes = np.zeros((num_nodes_y,num_nodes_x))
f_vec = np.zeros(num_nodes)
A = np.zeros((num_nodes,num_nodes))
def meshToVec(j,i)->int:
return i*num_nodes_y + j
def vecToMesh(h)->(int,int):
return (h % num_nodes_y, math.floor(h/num_nodes_y))
def computeFlux(j,i):
#computes the fluxes(-grad u) corresponding to cell j,i with forward euler.
vec_i = meshToVec(j,i)
north_i = meshToVec(j+1,i)
south_i = meshToVec(j-1,i)
east_i = meshToVec(j,i+1)
west_i = meshToVec(j,i-1)
north_face = (nodes[j,i+1,0]-nodes[j,i-1,0])/2
south_face = north_face
east_face = (nodes[j+1,i,1]-nodes[j-1,i,1])/2
west_face = east_face
dy = nodes[j+1,i,1]-nodes[j,i,1]
A[vec_i,north_i] += -1*north_face*1/dy
A[vec_i,vec_i] += 1*north_face*1/dy
dx = nodes[j,i+1,0]-nodes[j,i,0]
A[vec_i,east_i] += -1*east_face*1/dx
A[vec_i,vec_i] += 1*east_face*1/dx
dy = nodes[j,i,1]-nodes[j-1,i,1]
A[vec_i,south_i] += -1*south_face*1/dy
A[vec_i,vec_i] += 1*south_face*1/dy
dx = nodes[j,i,0]-nodes[j,i-1,0]
A[vec_i,west_i] += -1*west_face*1/dx
A[vec_i,vec_i] += 1*west_face*1/dx
return
def computeSource(j,i,x_d,y_d):
#Computes the source term by simple midpoint quadrature
north_face = (nodes[j,i+1,0]-nodes[j,i-1,0])/2
south_face = north_face
east_face = (nodes[j+1,i,1]-nodes[j-1,i,1])/2
west_face = east_face
f_vec[meshToVec(j,i)] = f.subs([(x,x_d),(y,y_d)])*north_face*east_face
#Loop trough the mesh and assemble the matrix.
for i in range(num_nodes_x):
for j in range(num_nodes_y):
vec_i = meshToVec(j,i)
#Check if we are looking at a boundary cell
if (i==0) or (i==num_nodes_x-1) or (j==0) or (j==num_nodes_y-1):
A[vec_i,vec_i] = 1
f_vec[vec_i] = 0
continue
computeFlux(j,i)
x_d = nodes[j,i,0]
y_d = nodes[j,i,1]
computeSource(j,i,x_d,y_d)
u_fabric_vec[vec_i] = u_fabric.subs([(x,x_d),(y,y_d)])
u_fabric_nodes[j,i] = u_fabric.subs([(x,x_d),(y,y_d)])
print(A)
u_vec = np.linalg.solve(A,f_vec)
u_nodes = u_fabric_nodes.copy()
err_nodes = np.zeros(num_nodes)
err_nodes_max = 0
for i in range(num_nodes):
u_nodes[vecToMesh(i)] = u_vec[i]
err_nodes[i] = abs(u_fabric_vec[i]-u_vec[i])
err_nodes_max = np.max(err_nodes)
print('max error at nodes',err_nodes_max)
fig = plt.figure(figsize=plt.figaspect(0.5))
ax = fig.add_subplot(1,2,1,projection='3d')
ax.plot_surface(nodes[:,:,0],nodes[:,:,1],u_nodes,cmap='viridis', edgecolor='none')
ax.set_title('computed solution')
ax.set_zlim(0.00, 0.07)
ax = fig.add_subplot(1, 2, 2, projection='3d')
ax.plot_surface(nodes[:,:,0],nodes[:,:,1],u_fabric_nodes,cmap='viridis', edgecolor='none')
ax.set_title('exact solution')
ax.set_zlim(0.00, 0.07)
plt.show()