pointnet cfd训练

  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 (kashefi@stanford.edu) or Davis Rempe (drempe@stanford.edu) via email.