1 ##### Point-cloud deep learning for prediction of fluid flow fields on irregular geometries (supervised learning) #####
9 import os # 提供与操作系统交互的功能，例如文件和目录操作。
10 import linecache # 提供从文件中读取特定行的方法。
11 import math # 提供数学函数和常量。
12 from operator import itemgetter # 提供获取对象元素的函数，通常用于排序和筛选操作。
13 import numpy as np # 导入NumPy库，并将其命名为np，用于进行数组操作和数学计算。
14 from numpy import zeros # 从NumPy库中导入zeros函数，用于创建全零数组。
15 import matplotlib # 导入Matplotlib库，用于数据可视化和绘图。
16
17
18 matplotlib.use('Agg') # 设置Matplotlib的后端为'Agg'，适用于无图形界面环境，如服务器上的绘图。
19 import matplotlib.pyplot as plt # 导入Matplotlib中的pyplot模块，用于创建各种图表和图形。
20
21
22 plt.rcParams['font.family'] = 'serif'
23 plt.rcParams['font.serif'] = ['Times New Roman'] + plt.rcParams['font.serif']
24 plt.rcParams['font.size'] = '12'
25 from mpl_toolkits.mplot3d import Axes3D
26 import tensorflow as tf
27 from tensorflow.python.keras.layers import Input
28 from tensorflow.python.keras import optimizers
29 from tensorflow.python.keras.models import Model
30 from tensorflow.python.keras.layers import Reshape
31 from tensorflow.python.keras.layers import Convolution1D, MaxPooling1D, BatchNormalization
32 from tensorflow.python.keras.layers import Lambda, concatenate
33
34 ##### Importing data #####
35 # For your convinient, we have already prepared data as a numpy array. You can download it from https://github.com/Ali-Stanford/PointNetCFD/blob/main/Data.npy
36 # The data for this specific test case is the spatial coordinates of the finite volume (or finite element) grids and the values of velocity and pressure fields on those grid points.
37 # The spatial coordinates are the input of PointNet and the velocity (in the *x* and *y* directions) and pressure fields are the output of PointNet.
38 # Here, our focus is on 2D cases.
39
40 Data = np.load('Data.npy') # 加载数据
41 data_number = Data.shape[0] # 获取数据的行数，即数据的数量
42
43 print('Number of data is:')
44 print(data_number)
45
46 point_numbers = 1024 # 每个数据点的数量 采样点
47 space_variable = 2 # 2 in 2D (x,y) and 3 in 3D (x,y,z) # 空间变量的维度
48 cfd_variable = 3 # (u, v, p)变量的维度，通常为（u，v，p） 速度场的x分量、y分量和压力场
49
50 input_data = zeros([data_number, point_numbers, space_variable], dtype='f') # 创建输入数据数组，维度为 [数据数量, 数据点数量, 空间变量维度]
51 output_data = zeros([data_number, point_numbers, cfd_variable], dtype='f') # 创建输出数据数组，维度为 [数据数量, 数据点数量, CFD变量维度]
52
53 for i in range(data_number):
54 input_data[i, :, 0] = Data[i, :, 0] # x coordinate (m)
55 input_data[i, :, 1] = Data[i, :, 1] # y coordinate (m)
56 output_data[i, :, 0] = Data[i, :, 3] # u (m/s)
57 output_data[i, :, 1] = Data[i, :, 4] # v (m/s)
58 output_data[i, :, 2] = Data[i, :, 2] # p (Pa)
59
60 ##### Normalizing data ##### 数据归一化
61 # We normalize the output data (velocity and pressure) in the range of [0, 1] using *Eq. 10* of the [journal paper](https://aip.scitation.org/doi/full/10.1063/5.0033376).
62 # Please note that due to this choice, we set the `sigmoid` activation function in the last layer of PointNet.
63 # Here, we do not normalize the input data (spatial coordinates, i.e., {*x*, *y*, *z*}). However, one may to normalize the input in the range of [-1, 1].
64
65 u_min = np.min(output_data[:, :, 0])
66 u_max = np.max(output_data[:, :, 0])
67 v_min = np.min(output_data[:, :, 1])
68 v_max = np.max(output_data[:, :, 1])
69 p_min = np.min(output_data[:, :, 2])
70 p_max = np.max(output_data[:, :, 2])
71
72 output_data[:, :, 0] = (output_data[:, :, 0] - u_min) / (u_max - u_min)
73 output_data[:, :, 1] = (output_data[:, :, 1] - v_min) / (v_max - v_min)
74 output_data[:, :, 2] = (output_data[:, :, 2] - p_min) / (p_max - p_min)
75
76
77 ##### Data visualization ##### 数据可视化
78 # We plot a few of input/output data.
79 # If you are using Google Colab, you can find the saved figures in the file sections (left part of the webpage).
80
81 # 绘制2D点云图
82 def plot2DPointCloud(x_coord, y_coord, file_name): # file_name为保存文件的名称
83 plt.scatter(x_coord, y_coord, s=2.5) # plt.scatter函数绘制点云图 s=2.5指定散点的大小
84 plt.xlabel('x (m)')
85 plt.ylabel('y (m)')
86 x_upper = np.max(x_coord) + 1 # 坐标轴范围
87 x_lower = np.min(x_coord) - 1
88 y_upper = np.max(y_coord) + 1
89 y_lower = np.min(y_coord) - 1
90 plt.xlim([x_lower, x_upper])
91 plt.ylim([y_lower, y_upper])
92 plt.gca().set_aspect('equal', adjustable='box') # 等比例缩放
93 plt.savefig(file_name + '.png', dpi=300)
94 # plt.savefig(file_name+'.eps') #You can use this line for saving figures in EPS format
95 plt.clf() # 清楚图形缓存，以便下次使用
96 # plt.show()
97
98
99 # 绘制二维结果云图
100 def plotSolution(x_coord, y_coord, solution, file_name, title):
101 plt.scatter(x_coord, y_coord, s=2.5, c=solution, cmap='jet') # c=solution颜色映射的数值来源，颜色映射为 'jet'彩虹色
102 plt.title(title)
103 plt.xlabel('x (m)')
104 plt.ylabel('y (m)')
105 x_upper = np.max(x_coord) + 1
106 x_lower = np.min(x_coord) - 1
107 y_upper = np.max(y_coord) + 1
108 y_lower = np.min(y_coord) - 1
109 plt.xlim([x_lower, x_upper])
110 plt.ylim([y_lower, y_upper])
111 plt.gca().set_aspect('equal', adjustable='box')
112 cbar = plt.colorbar()
113 plt.savefig(file_name + '.png', dpi=300)
114 # plt.savefig(file_name+'.eps') #You can use this line for saving figures in EPS format
115 plt.clf()
116 # plt.show()
117
118 # 数据进行可视化 通过修改number数值
119 number = 0 # It should be less than 'data_number'
120 plot2DPointCloud(input_data[number, :, 0], input_data[number, :, 1], 'PointCloud')
121 plotSolution(input_data[number, :, 0], input_data[number, :, 1], output_data[number, :, 0], 'u_velocity',
122 'u (x-velocity component)')
123 plotSolution(input_data[number, :, 0], input_data[number, :, 1], output_data[number, :, 1], 'v_velocity',
124 'v (y-velocity component)')
125 plotSolution(input_data[number, :, 0], input_data[number, :, 1], output_data[number, :, 2], 'pressure', 'pressure')
126
127 ##### Spliting data #####
128 # We split the data *randomly* into three categories of training, validation, and test sets. 训练集、验证集和测试集
129 # A reasonable partitioning could be 80% for the training, 10% for validation, and 10% for test sets. 8 1 1
130
131 indices = np.random.permutation(input_data.shape[0]) # 创建一个随机排列的数据索引数组
132 # 索引数组 indices 按照比例分成三个部分
133 training_idx, validation_idx, test_idx = indices[:int(0.8 * data_number)], indices[int(0.8 * data_number):int(
134 0.9 * data_number)], indices[int(0.9 * data_number):]
135 #划分输入输出值
136 input_training, input_validation, input_test = input_data[training_idx, :], input_data[validation_idx, :], input_data[
137 test_idx, :]
138 output_training, output_validation, output_test = output_data[training_idx, :], output_data[validation_idx,
139 :], output_data[test_idx, :]
140
141 ##### PointNet architecture ##### pointnet结构
142 # We use the segmentation component of PointNet. One of the most important features of PointNet is its scalability. 使用了 PointNet 的分割组件。PointNet 最重要的特性之一是其可扩展性。
143 # The variable `scaling` in the code allows users to make the network bigger or smaller to control its capacity.
144 # Please note that the `sigmoid` activation function is implemented in the last layer, since we normalize the output data in range of [0, 1].
145 # Additionally, we have removed T-Nets (Input Transforms and Feature Transforms) from the network implemented here.
146 # However, please note that we had embedded T-Nets in the network in our [journal paper](https://aip.scitation.org/doi/figure/10.1063/5.0033376) (please see *Figure 5* of the journal paper).
147
148 scaling = 1.0 # reasonable choices for scaling: 4.0, 2.0, 1.0, 0.25, 0.125
149
150
151 def exp_dim(global_feature, num_points):
152 return tf.tile(global_feature, [1, num_points, 1])
153
154
155 PointNet_input = Input(shape=(point_numbers, space_variable)) # 输入数据的形状 1024*2 二维是2
156 # Shared MLP (64,64)
157 branch1 = Convolution1D(int(64 * scaling), 1, input_shape=(point_numbers, space_variable), activation='relu')(
158 PointNet_input)
159 branch1 = BatchNormalization()(branch1)
160 branch1 = Convolution1D(int(64 * scaling), 1, input_shape=(point_numbers, space_variable), activation='relu')(branch1)
161 branch1 = BatchNormalization()(branch1)
162 Local_Feature = branch1
163 # Shared MLP (64,128,1024)
164 branch1 = Convolution1D(int(64 * scaling), 1, activation='relu')(branch1)
165 branch1 = BatchNormalization()(branch1)
166 branch1 = Convolution1D(int(128 * scaling), 1, activation='relu')(branch1)
167 branch1 = BatchNormalization()(branch1)
168 branch1 = Convolution1D(int(1024 * scaling), 1, activation='relu')(branch1)
169 branch1 = BatchNormalization()(branch1)
170 # Max function
171 Global_Feature = MaxPooling1D(pool_size=int(point_numbers * scaling))(branch1)
172 Global_Feature = Lambda(exp_dim, arguments={'num_points': point_numbers})(Global_Feature)
173 branch2 = concatenate([Local_Feature, Global_Feature])
174 # Shared MLP (512,256,128)
175 branch2 = Convolution1D(int(512 * scaling), 1, activation='relu')(branch2)
176 branch2 = BatchNormalization()(branch2)
177 branch2 = Convolution1D(int(256 * scaling), 1, activation='relu')(branch2)
178 branch2 = BatchNormalization()(branch2)
179 branch2 = Convolution1D(int(128 * scaling), 1, activation='relu')(branch2)
180 branch2 = BatchNormalization()(branch2)
181 # Shared MLP (128, cfd_variable)
182 branch2 = Convolution1D(int(128 * scaling), 1, activation='relu')(branch2)
183 branch2 = BatchNormalization()(branch2)
184 PointNet_output = Convolution1D(cfd_variable, 1, activation='sigmoid')(
185 branch2) # Please note that we use the sigmoid activation function in the last layer.
186
187 ##### Defining and compiling the model ##### 负责定义和编译模型
188 # We use the `Adam` optimizer. Please note to the choice of `learning_rate` and `decaying_rate`. Adam 优化器 学习率和衰减率
189 # The network is also sensitive to the choice of `epsilon` and it has to be set to a non-zero value. epsilon 参数 设置为非零值
190 # We use the `mean_squared_error` as the loss function (please see *Eq. 15* of the [journal paper](https://aip.scitation.org/doi/full/10.1063/5.0033376)). 均方误差（mean_squared_error）作为损失函数
191 # Please note that for your specific application, you need to tune the `learning_rate` and `decaying_rate`. 学习率和衰减率
192
193 learning_rate = 0.0005
194 decaying_rate = 0.0 # 衰减率为零，表示不对学习率进行衰减
195 model = Model(inputs=PointNet_input, outputs=PointNet_output)
196 # compile配置模型
197 model.compile(optimizers.Adam(lr=learning_rate, beta_1=0.9, beta_2=0.999, epsilon=0.000001, decay=decaying_rate)
198 , loss='mean_squared_error', metrics=['mean_squared_error'])
199
200 ##### Training PointNet ##### 训练pointnet
201 # Please be careful about the choice of batch size (`batch`) and number of epochs (`epoch_number`). 批处理大小 (batch) 和时期数
202 # At the beginning, you might observe an increase in the validation loss, but please do not worry, it will eventually decrease.
203 # Please note that this section might take approximately 20 hours to be completed (if you are running the code on Google Colab). So, please be patient.
204 # Alternatively, you can run this section on your cluster computing.
205
206 batch = 32
207 epoch_number = 4000
208 # fit训练模型
209 results = model.fit(input_training, output_training, batch_size=batch, epochs=epoch_number, shuffle=True, verbose=1,
210 validation_split=0.0, validation_data=(input_validation, output_validation))
211
212 ##### Plotting training history ##### 制训练过程中损失值的变化曲线
213 # Trace of loss values over the training and validation set can be seen by plotting the history of training.
214
215 plt.plot(results.history['loss'])
216 plt.plot(results.history['val_loss'])
217 plt.yscale('log')
218 plt.ylabel('Loss')
219 plt.xlabel('Epoch')
220 plt.legend(['Training', 'Validation'], loc='upper left')
221 plt.savefig('Loss_History.png', dpi=300)
222 # plt.savefig('Loss_History.eps') #You can use this line for saving figures in EPS format
223 plt.clf()
224 # plt.show()
225
226 ##### Error analysis ##### 误差分析
227 # We can perform various error analyses. For example, here we are interested in the normalized root mean square error (RMSE).
228 # We compute the normalized root mean square error (pointwise) of the predicted velocity and pressure fields over each geometry of the test set, and then we take an average over all these domains.
229 # We normalize the velocity error by the free stream velocity, which is 1 m/s. Moreover, we normalized the pressure error by the dynamic pressure (without the factor of 0.5) at the free stream.
230
231 # 初始化三个变量
232 average_RMSE_u = 0
233 average_RMSE_v = 0
234 average_RMSE_p = 0
235
236 test_set_number = test_idx.size # 计算测试集的样本数量
237 sample_point_cloud = zeros([1, point_numbers, space_variable], dtype='f')
238 truth_point_cloud = zeros([point_numbers, cfd_variable], dtype='f')
239
240 for s in range(test_set_number):
241 for j in range(point_numbers):
242 for i in range(space_variable):
243 sample_point_cloud[0][j][i] = input_test[s][j][i]
244
245 prediction = model.predict(sample_point_cloud, batch_size=None, verbose=0)
246
247 # Unnormalized
248 prediction[0, :, 0] = prediction[0, :, 0] * (u_max - u_min) + u_min
249 prediction[0, :, 1] = prediction[0, :, 1] * (v_max - v_min) + v_min
250 prediction[0, :, 2] = prediction[0, :, 2] * (p_max - p_min) + p_min
251
252 output_test[s, :, 0] = output_test[s, :, 0] * (u_max - u_min) + u_min
253 output_test[s, :, 1] = output_test[s, :, 1] * (v_max - v_min) + v_min
254 output_test[s, :, 2] = output_test[s, :, 2] * (p_max - p_min) + p_min
255
256 average_RMSE_u += np.sqrt(np.sum(np.square(prediction[0, :, 0] - output_test[s, :, 0]))) / point_numbers
257 average_RMSE_v += np.sqrt(np.sum(np.square(prediction[0, :, 1] - output_test[s, :, 1]))) / point_numbers
258 average_RMSE_p += np.sqrt(np.sum(np.square(prediction[0, :, 2] - output_test[s, :, 2]))) / point_numbers
259
260 average_RMSE_u = average_RMSE_u / test_set_number
261 average_RMSE_v = average_RMSE_v / test_set_number
262 average_RMSE_p = average_RMSE_p / test_set_number
263
264 print('Average normalized RMSE of the x-velocity component (u) over goemtries of the test set:')
265 print(average_RMSE_u)
266 print('Average normalized RMSE of the y-velocity component (v) over goemtries of the test set:')
267 print(average_RMSE_v)
268 print('Average normalized RMSE of the pressure (p) over goemtries of the test set:')
269 print(average_RMSE_p)
270
271 # 创建一个包含数据的字典
272 data = {
273 'Average RMSE of u': [average_RMSE_u],
274 'Average RMSE of v': [average_RMSE_v],
275 'Average RMSE of p': [average_RMSE_p]
276 }
277
278 # 将字典转换为DataFrame
279 df = pd.DataFrame(data)
280
281 # 检查文件是否存在，如果不存在则创建
282 file_path = 'average_RMSE_values.xlsx'
283 if not os.path.isfile(file_path):
284 df.to_excel(file_path, index=False)
285 else:
286 # 如果文件存在，追加数据
287 with pd.ExcelWriter(file_path, engine='openpyxl', mode='a') as writer:
288 df.to_excel(writer, index=False, header=False)
289
290
291 ##### Saving the trained model or weights #####
292 # Save the entire model
293 model.save('pointnet_model.h5')
294
295 # Alternatively, save only the weights
296 model.save_weights('pointnet_weights.h5')
297
298 ##### Visualizing the results #####
299 # For example, let us plot the velocity and pressure fields predicted by the network as well as absolute pointwise error distribution over these fields for one of the geometries of the test set.
300
301 s = 0 # s can change from "0" to "test_set_number - 1"
302
303 # 遍历测试集某一样本中的所有点和空间变量
304 for j in range(point_numbers):
305 for i in range(space_variable):
306 sample_point_cloud[0][j][i] = input_test[s][j][i]
307
308 prediction = model.predict(sample_point_cloud, batch_size=None, verbose=0)
309
310 # Unnormalized 反归一化
311 prediction[0, :, 0] = prediction[0, :, 0] * (u_max - u_min) + u_min
312 prediction[0, :, 1] = prediction[0, :, 1] * (v_max - v_min) + v_min
313 prediction[0, :, 2] = prediction[0, :, 2] * (p_max - p_min) + p_min
314
315 # Please note that we have already unnormalized the 'output_test' in Sect. Error analysis.
316
317 plotSolution(sample_point_cloud[0, :, 0], sample_point_cloud[0, :, 1], prediction[0, :, 0], 'velocity prediction',
318 'velocity prediction (u)')
319 plotSolution(sample_point_cloud[0, :, 0], sample_point_cloud[0, :, 1], output_test[s, :, 0], 'velocity ground truth',
320 'velocity ground truth (u)')
321 plotSolution(sample_point_cloud[0, :, 0], sample_point_cloud[0, :, 1],
322 np.absolute(output_test[s, :, 0] - prediction[0, :, 0]), 'u_point_wise_error',
323 'absolute point-wise error of velocity (u) m/s')
324
325 plotSolution(sample_point_cloud[0, :, 0], sample_point_cloud[0, :, 1], prediction[0, :, 1], 'v_prediction',
326 'velocity prediction (v) m/s')
327 plotSolution(sample_point_cloud[0, :, 0], sample_point_cloud[0, :, 1], output_test[s, :, 1], 'v_ground_truth',
328 'velocity ground truth (v) m/s')
329 plotSolution(sample_point_cloud[0, :, 0], sample_point_cloud[0, :, 1],
330 np.absolute(output_test[s, :, 1] - prediction[0, :, 1]), 'v_point_wise_error',
331 'absolute point-wise error of velocity (v) m/s')
332
333 plotSolution(sample_point_cloud[0, :, 0], sample_point_cloud[0, :, 1], prediction[0, :, 2], 'p_prediction',
334 'pressure prediction (p) Pa')
335 plotSolution(sample_point_cloud[0, :, 0], sample_point_cloud[0, :, 1], output_test[s, :, 2], 'p_ground_truth',
336 'pressure ground truth (p) Pa')
337 plotSolution(sample_point_cloud[0, :, 0], sample_point_cloud[0, :, 1],
338 np.absolute(output_test[s, :, 2] - prediction[0, :, 2]), 'p_point_wise_error',
339 'absolute point-wise error of pressure (p) Pa')
340
341 # Questions?
342 # If you have any questions or need assistance,
343 # please do not hesitate to contact Ali Kashefi ([email protected]) or Davis Rempe ([email protected]) via email.
标签：plt,训练,point,cfd,pointnet,min,test,output,data From： https://www.cnblogs.com/arwen-xu/p/17878503.html