-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathpressure_utility.py
146 lines (128 loc) · 5.42 KB
/
pressure_utility.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
143
144
145
146
import scipy.io
import numpy as np
import matplotlib.pyplot as plt
def read_data_at_snapshot(filename, num_snap, resolution):
data_mat = scipy.io.loadmat(filename)
x = data_mat['x'] # N*1
y = data_mat['y'] # N*1
t = data_mat['t'][num_snap, 0] # T*1
u = data_mat['u'][:, num_snap] # N*T
v = data_mat['v'][:, num_snap] # N*T
x_unique = np.unique(x).reshape(-1, 1)
y_unique = np.unique(y).reshape(-1, 1)
t = t.item()
u_mat = u.reshape(resolution, resolution)
v_mat = v.reshape(resolution, resolution)
return x_unique, y_unique, t, u_mat, v_mat
def read_data_all(filename, resolution):
data_mat = scipy.io.loadmat(filename)
x = data_mat['x'] # N*1
y = data_mat['y'] # N*1
t = data_mat['t'][:, 0] # T*1
u = data_mat['u'][:, :] # N*T
v = data_mat['v'][:, :] # N*T
x_unique = np.unique(x).reshape(-1, 1)
y_unique = np.unique(y).reshape(-1, 1)
return x_unique, y_unique, t, u, v
# generate grid with one extra virtual point along each direction
def generate_extend_grid(original_mat):
resolution = original_mat.shape[0]
extended_mat = np.ones((resolution+2, resolution+2))
extended_mat[1: 1+resolution, 1:1+resolution] = original_mat
extended_mat[1: 1+resolution, 0] = original_mat[:, -1]
extended_mat[1: 1 + resolution, -1] = original_mat[:, 0]
extended_mat[0, 1: 1 + resolution] = original_mat[-1, :]
extended_mat[-1, 1: 1 + resolution] = original_mat[0, :]
return extended_mat
# get the numerical derivative with 2 order central difference
def get_derivative(extended_mat, delta):
resolution = extended_mat.shape[0]-2
dx = (extended_mat[2:, 1:1+resolution] - extended_mat[:-2, 1:1+resolution]) / (2 * delta)
dy = (extended_mat[1:1+resolution, 2:] - extended_mat[1:1+resolution, :-2]) / (2 * delta)
return dx, dy
def calculate_source(u_mat, v_mat, delta):
# generate grid for u and v
u_extended = generate_extend_grid(u_mat)
v_extended = generate_extend_grid(v_mat)
ux, uy = get_derivative(u_extended, delta)
vx, vy = get_derivative(v_extended, delta)
f = -(ux*ux+2*uy*vx+vy*vy)
return f
def assign_periodic_condition(extend_mat):
extend_mat[0, :] = extend_mat[-2, :]
extend_mat[-1, :] = extend_mat[1, :]
extend_mat[:, 0] = extend_mat[:, -2]
extend_mat[:, -1] = extend_mat[:, 1]
return extend_mat
def plot_solution(x, y, p):
mesh_x, mesh_y = np.meshgrid(x, y)
plt.figure(figsize=(8, 6))
plt.contourf(mesh_x, mesh_y, p, levels=200, cmap='jet')
plt.show()
return
def solve_poisson_sgs(p0, f_extend, delta, tol=1e-5, max_iter=10000):
"""
solve poisson equation with source term, with periodic boundary condition
:param p0: initial guess
:param f_extend: source term
:param delta: delta=dx=dy
:param tol: the minimum tolerance
:param max_iter: max acceptable number of iterations
:return: return the solution (without extension of virtual grid), shape[resolution, resolution]
"""
resolution = p0.shape[0]-2
for it in range(max_iter):
p_old = np.copy(p0)
# upward scan
for j in range(1, 1 + resolution):
for i in range(1, 1 + resolution):
p0[i, j] = 0.25 * (p0[i + 1, j] + p0[i - 1, j] + p0[i, j + 1] + p0[i, j - 1] - (delta ** 2) * f_extend[i, j])
# downward scan
for j in range(1, 1 + resolution):
jr = resolution - j + 1
for i in range(1, 1 + resolution):
ir = resolution - i + 1
p0[ir, jr] = 0.25 * (p0[ir+1, jr] + p0[ir-1, jr] + p0[ir, jr+1] + p0[ir, jr-1] - (delta ** 2) * f_extend[ir, jr])
# assign periodic boundary condition
p0 = assign_periodic_condition(p0)
# Exit condition
res = np.max(np.abs(p0 - p_old))
if res < tol:
break
else:
print("Current residual: ", res, " Current iter num:", it)
p_inside = p0[1:1+resolution, 1:1+resolution]
return p_inside
def solve_poisson_lu_sgs(p0, f_extend, delta, tol=1e-6, max_iter=1000):
"""
solve poisson equation with source term, with periodic boundary condition
:param p0: initial guess
:param f_extend: source term
:param delta: delta=dx=dy
:param tol: the minimum tolerance
:param max_iter: max acceptable number of iterations
:return: return the solution (without extension of virtual grid), shape[resolution, resolution]
"""
resolution = p0.shape[0]-2
for it in range(max_iter):
p_old = np.copy(p0)
# upward scan
for j in range(1, 1 + resolution):
for i in range(1, 1 + resolution):
p0[i, j] = 0.25 * (p0[i - 1, j] + p0[i, j - 1] - (delta ** 2) * f_extend[i, j])
# downward scan
for j in range(1, 1 + resolution):
jr = resolution - j + 1
for i in range(1, 1 + resolution):
ir = resolution - i + 1
p0[ir, jr] = 0.25 * (p0[ir+1, jr] + p0[ir, jr+1]) + p0[ir, jr]
# assign periodic boundary condition
p0 = assign_periodic_condition(p0)
# Exit condition
res = np.max(np.abs(p0 - p_old))
if res < tol:
break
else:
print("Current residual: ", res, " Current iter num:", it)
p_inside = p0[1:1+resolution, 1:1+resolution]
return p_inside, p0