简单案例演示传统CFD和PINN

以一个简单的例子来演示传统 CFD 与 PINN 在求解一维泊松方程的基本过程。

要求解的方程为：

方程右侧边界为 Neumman 边界：

左侧边界为Dirichlet 边界条件

原方程足够简单，很容易求得解析解：

1 传统离散方法

以有限差分法为例展示传统方法求解此方程的基本过程。有限体积法离散方式稍有差别，但最终得到的方程组是一样的。

有限差分法需要构造网格系统，本案例计算区域为  ，构造如下图所示的网格系统，其包含6个网格节点，其中节点1和节点6为网格边界，节点2、3、4、5为内部节点。这里采用均匀网格，网格间距为  。

利用中心差分离散二阶导数项  ，方程可离散为：

因此可得到离散方程组：

对于边界节点1：

对于边界节点6：

然后可以整理成代数方程组：

求解此方程可得到离散解：

利用python代码很容易实现：

import numpy as np
import matplotlib.pyplot as plt

L = 2
N = 5
deltaX = 2/N

# 构造两个空矩阵，用来放系数
A = np.zeros((N+1,N+1))
b = np.zeros((N+1,1))

# 左侧边界
A[0,0] = 1

# 右侧边界
A[N,N] = 1
A[N,N-1] =-1

# 中间节点
for i in range(1,N):
    A[i,i-1] = 1
    A[i,i] = -2
    A[i,i+1] = 1

# 对于b矩阵赋值
b[0] = 0
b[N] = 4*deltaX
for i in range(1,N):
    b[i] = 2 * deltaX * deltaX

# 直接求解矩阵
u_solve = np.linalg.solve(A,b)

# 图形显示结果
x = np.linspace(-1,1,N+1)
# 解析解
y = (x + 1 ) **2

plt.figure(figsize=(8,5))
plt.plot(x,y,label='true',color = 'red',zorder = 1)
plt.scatter(x,u_solve.reshape(1,-1)[0],color = 'blue',label = 'numerical',zorder =2)
plt.xlabel('x')
plt.ylabel('u')
plt.legend()

显示结果为：

看到精度比较差，可以尝试增大网格数量，如将N改为30，可以看到精度好了不少，如下图所示。当然也可以体会到，计算结果依赖于计算网格的数量。但如果想要改善计算精度，则需要更改网格系统重新求解一遍。并不能在粗糙网格结果的基础上直接通过简单的代数运算得到精细网格的结果。

2 PINN求解

PINN是一种处理思路，原则上是不依赖于任何框架或计算机语言的。这里简单起见，使用DeepXDE框架进行处理。DeepXDE属于一种前端框架，其可以支持的后端框架包括Pytorch、Tensorflow、Jax、PaddlePaddle等。

PINN原理啥的我也是新手，讲不明白。直接上代码（这个是DeepXDE文档自带的案例，有小幅改动）：

import deepxde as dde
import numpy as np
import matplotlib.pyplot as plt

# 定义PDE方程损失
def pde(x, y):
    dy_xx = dde.grad.hessian(y, x)
    return dy_xx - 2

# 左边界
def boundary_l(x, on_boundary):
    return on_boundary and dde.utils.isclose(x[0], -1)

# 右边界
def boundary_r(x, on_boundary):
    return on_boundary and dde.utils.isclose(x[0], 1)

# 定义计算区域[-1,1]
geom = dde.geometry.Interval(-1, 1)

# 定义左右两侧的边界条件
bc_l = dde.icbc.DirichletBC(geom, lambda x:0, boundary_l)
bc_r = dde.icbc.NeumannBC(geom, lambda x:4, boundary_r)

# 定义PDE方程
data = dde.data.PDE(geom, pde, [bc_l, bc_r], 16, 2, num_test=100)

# 定义全连接神经网络，包括1个输入层，3个隐藏层，每层50个神经元，以及1个输出层
layer_size = [1] + [50] * 3 + [1]
activation = "tanh"
initializer = "Glorot uniform"
net = dde.nn.FNN(layer_size, activation, initializer)

# 定义Model，采用adam优化器，学习率使用0.001
model = dde.Model(data, net)
model.compile("adam", lr=0.001)

# 训练网络
losshistory, train_state = model.train(iterations=5000)

# 查看结果
x  = np.linspace(-1,1,30).reshape(-1,1)
y_ture = (x + 1) ** 2
u_predict = model.predict(x)

u_predict = u_predict.reshape(1,-1)[0]

plt.figure(figsize=(8,5))
plt.plot(x, y_ture, 'b-', label='True')
plt.scatter(x,u_predict,marker='v',color = 'red',label= 'PINN')
plt.xlabel('x')
plt.ylabel('u')
plt.legend()

得到结果如下：

 

结果还是不错的。当然前面网络训练5000轮实在是有点儿多了，实际上对于这个简单的问题，训练个千八百次精度完全够用。比如训练1000次的结果如下图所示。与上图相比看不出来有啥区别。需要注意的是，这里显示的点数是后处理添加的，并不会对计算精度产生影响。

 

对模型精度产生影响的地方是程序中data定义里面指定的精度点数。如下面代码中的16和2。

data = dde.data.PDE(geom, pde, [bc_l, bc_r], 16, 2, num_test=100)

其他包括网络结构、初始化方法、激活函数、优化方法等也会对计算结果造成影响。当然这里的案例方程足够简单，如果要求解复杂的控制方程（比如NS方程），这里面估计就很烧脑了。这个有兴趣的道友可以自行寻找相关资料查看。

PINN得到的结果有个好处，其获得的网络在推理计算时可以当做是连续的函数，如需要更细密的结果，不需要重新训练模型，如下面获取[-1,1]之间50个点的结果，预测结果依然与解析解吻合良好。而传统方法如果想要得到更加精细的结果，就得重新划分网格，重新进行计算。

 

PINN除了可以嵌入物理信息外，还可以嵌入实验数据（包括CFD仿真数据），嵌入的真实数据有利于提高模型计算精度及收敛速度，也有利于缩短模型训练时间。

随着技术的发展，也许在未来，PINN能够做到在条件工况（如计算区域、边界条件、初始条件等）发生改变后，不需要对网络进行完全训练，而是像今天使用的迁移学习那样，仅需在基础网络模型的基础上进行小幅调整即可适应于新的问题。这问题目前瞅着还挺难的，但并非不可能。而传统方法基本上已经断绝了可能，输入条件变了就只能重新算了。

PINN的问题其实还蛮多的，当然这里的案例过于简单，麻烦的地方都没有体现出来。当然PINN属于新玩意儿，未来的走向谁也说不清楚，也许过些年又会有更好的方法冒出来也说不定。传统CFD方法还能深挖的领域基本上都是硬骨头。

工程应用的话，个人觉得有限体积法在相当长的一段时间里依然会是CFD计算的主流方法，短期内根本看不到PINN能够取代有限体积法的希望。

（完）