版本对应：
本网页讨论的机器学习内容，为特定的机器学习内容（并未与CFD结合）。在《无痛苦NS方程笔记》中， 特定的将CFD与机器学习互相结合起来，无普适性机器学习内容。
ML: PINN代码实例-1D槽道流
CFD中的一维管道流可以通过1D传输方程来进行。在OpenFOAM中，可以通过boundaryFoam求解器来进行，更详细的理论背景请参考boundaryFoam解析。本网页通过OpenFOAM libtorch环境，调用PINN来对其进行求解。
首先回顾以下基本的控制方程：
在一个充分发展的管道里面，速度型线的分布可以用方程(1)来表示。boundaryFoam解析里曾讨论过，若已知\(\frac{\partial p}{ \partial x}\)的值，则可以获得\(u_1\)的解。本文将其假定为\(1.20179\times 10^{-5}\)。因此方程(1)可以理解为一个泊松方程，上边界为零法向梯度，下边界的\(u_1=0\)（固定值边界条件）。同时假定粘度为定值：\(\nu_t=1\times 10^{-8}\)。因此，求解方程可以写为：
至此，需要给定到PINN的参数即为方程(2)，以及无滑移壁面边界条件。对方程(2)如果感觉不太理解，请参考boundaryFoam解析。然而也可以不理解方程(2)的流体力学背景，仅仅将其当做一个具有一个未知数\(u_1\)的泊松方程。
PINN与有限体积法FVM
PINN与FVM是完全两条路。二者都是求解PDE的方法。FVM是从数学上严格推导出离散的数值解。PINN更像是通过神经网络去拟合出一个数值解。
PINN网格划分
在PINN里面同样需要进行网格划分。然而PINN领域的人一般不将其称之为网格，而叫做配点（Collocation Points）。PINN在求解过程中，需要保证配点上满足：
方程约束，也即方程(2)；
边界约束，也即本文给定的速度为0的边界条件；
在OpenFOAM libtorch环境下，可以通过下面几行来生成网格：
// Computational domain
double dy = 0.00125;
auto init = torch::full({}, 1.0);
auto mesh =
torch::arange(-0.05, 0, dy, torch::requires_grad()).unsqueeze(1);
int cellN = static_cast<int>(0.05/dy);
上述几行命令生成的网格mesh为从-0.05至0，网格单元大小为0.00125dy，网格数量为0.05/0.00125cellN。有关OpenFOAM libtorch环境下的tensor使用，请参考ML: libtorch张量（一）。
Warning
机器学习里面的tensor与计算物理里面的tensor不太一样。一个很明显的区别就是，计算物理里面的tensor大部分情况下是正方形的，比如剪切应力是\(3\times 3\)的。但是机器学习里面的tensor什么样的都有，很少有\(3\times 3\)的。详情请参考计算物理的tensor，机器学习的tensor，两个一样么？
神经网络搭建
在OpenFOAM libtorch环境下，可以通过下述代码搭建神经网络：
class NN
:
public torch::nn::Module
{
torch::nn::Sequential net_;
public:
NN()
{
net_ = register_module
(
"net",
torch::nn::Sequential
(
torch::nn::Linear(1,8),
torch::nn::Tanh(),
torch::nn::Linear(8,1)
)
);
}
auto forward(torch::Tensor x)
{
return net_->forward(x);
}
};
在这里主要的是，输入层为1，隐藏层具有8个神经元，输出层为1。其可以理解为，输入一个网格坐标，就会输出一个速度值。在这里我们使用的是Tanh激活函数。
超参数以及自动微分
PINN在训练网络的时候，需要指定梯度下降的方法。这部分理论请参考ML: 手算反向传播与自动微分。具体的，本算例训练10000次，采用AdamW梯度更新方法。学习率为0.001。损失计算的方法采用均方差最小化。在OpenFOAM libtorch环境下的代码可以写为
int iter = 1;
int IterationNum = 100000; //训练10000次
double learning_rate = 0.001; //学习率为0.001
auto crit = torch::nn::MSELoss(); //损失函数均方差最小化
auto model = std::make_shared<NN>();
auto opti =
std::make_shared<torch::optim::AdamW> //AdamW梯度更新方法
(
model->parameters(),
torch::optim::AdamWOptions(learning_rate)
);
在训练的过程中，每次都需要将网格信息输入到神经网络，然后其会预测出一个速度值，这可以通过下述代码来实现：
auto upred = model->forward(mesh);
不同于FVM，在PINN中，导数的计算是通过自动微分来实现。这部分理论请参考ML: 手算反向传播与自动微分。在OpenFOAM libtorch环境下的代码可以写为
auto dudy =
torch::autograd::grad
(
{upred},
{mesh},
{torch::ones_like(upred)},
true,
true
)[0];
auto dudyy =
torch::autograd::grad
(
{dudy},
{mesh},
{torch::ones_like(upred)},
true,
true
)[0];
上面的代码声明的dudy以及dudyy分别是\(\frac{\p u_1}{\p y}\)以及\(\frac{\p}{\p y}\left(\frac{\p u_1}{\p y}\right)\)。
损失
在计算了导数之后，要定义损失。本算例定义两类损失：
方程损失；
边界损失；
首先看方程损失。下列代码定义方程(2)中的源项：
auto dpdxByNu = torch::full({cellN}, 1.2e-05/1e-8);
结合刚才计算的dudyy\(\frac{\p}{\p y}\left(\frac{\p u_1}{\p y}\right)\)，可以对比二者的差值（理想情况下需要是\(0\)）。
auto loss_pde = crit(-dudyy.reshape({cellN}), dpdxByNu);
下面定义边界损失，首先定义下边界处的速度值\(u_1=0.0375\)：
auto ubottom = torch::full_like(upred[0], 0.0375);
Warning
之所以定义\(u_1=0.0375\)是为了与OpenFOAM的输出数据做一个完全的对应。在后续我们将对比PINN与OpenFOAM网格点速度的差异。OpenFOAM网格点的速度值并不是网格边界的速度值。\(0.0375\)为OpenFOAM第一层网格的速度值。
我们希望预测的upred在第0个网格点上的速度，与给定的0.0375无差异，那么损失就可以定义为：
auto loss_bottom = crit(upred[0], ubottom);
类似的，我们可以定义上边界零法向梯度的损失：
auto dudyTop = dudy[cellN - 1][0];
auto loss_top = crit(dudyTop, 0.0*dudyTop);
在有了方程损失以及边界损失之后，我们可以计算总损失：
auto loss = loss_pde + loss_bd;
计算总损失后，可以进行反向传播。在迭代的过程中，总体思想就是通过预测的upred，来计算各项损失，最终要求损失最小。在这个情况下，即同时满足方程损失以及边界损失。也即满足偏微分方程的解。
运行算例
请首先按照OpenFOAM libtorch tutorial step by step来配置OpenFOAM libtorch环境。也可以在这里下载配置好的虚拟机。本算例的源代码请参考本文文末。
将源代码直接创立一个.C文件后，即可通过wmake来编译。编译后可以进行训练。训练过程的log如下：
dyfluid@dyfluid-virtual-machine:~/OKS202409/boundary$ boundary
250 loss:1.43925e+06 pde: 1.43903e+06 bd: 209.787 ini: 8.42688
500 loss:1.43611e+06 pde: 1.4351e+06 bd: 878.38 ini: 138.162
750 loss:1.42121e+06 pde: 1.41882e+06 bd: 1588.68 ini: 795.64
1000 loss:1.38646e+06 pde: 1.38089e+06 bd: 3019.29 ini: 2548.81
1250 loss:1.32955e+06 pde: 1.31894e+06 bd: 4920.52 ini: 5679.98
1500 loss:1.25027e+06 pde: 1.23302e+06 bd: 6992.6 ini: 10250.6
1750 loss:1.15032e+06 pde: 1.12526e+06 bd: 8942.17 ini: 16121.1
2000 loss:1.03329e+06 pde: 999892 bd: 10473.4 ini: 22929.6
2250 loss:904337 pde: 862790 bd: 11413.6 ini: 30133.7
2500 loss:769703 pde: 721042 bd: 11529.2 ini: 37131.9
2750 loss:636071 pde: 581951 bd: 10917.2 ini: 43202.3
3000 loss:509803 pde: 452502 bd: 9669.03 ini: 47632.3
最终求解器会输出相关的信息，如dudyy\(\frac{\p}{\p y}\left(\frac{\p u_1}{\p y}\right)\)等。在这里我们将PINN预测的速度与OpenFOAM的boundaryFoam预测的速度进行对比，如下图（通过gnuplot绘制）：
源代码：
#include <torch/torch.h>
using namespace std;
class NN
:
public torch::nn::Module
{
torch::nn::Sequential net_;
public:
NN()
{
net_ = register_module
(
"net",
torch::nn::Sequential
(
torch::nn::Linear(1,8),
torch::nn::Tanh(),
torch::nn::Linear(8,1)
)
);
}
auto forward(torch::Tensor x)
{
return net_->forward(x);
}
};
int main()
{
int iter = 1;
int IterationNum = 1000000;
double learning_rate = 0.001;
auto crit = torch::nn::MSELoss();
auto model = std::make_shared<NN>();
auto opti =
std::make_shared<torch::optim::AdamW>
(
model->parameters(),
torch::optim::AdamWOptions(learning_rate)
);
// Computational domain
double dy = 0.00125;
auto init = torch::full({}, 1.0);
auto mesh =
torch::arange(-0.05, 0, dy, torch::requires_grad()).unsqueeze(1);
int cellN = static_cast<int>(0.05/dy);
auto dpdxByNu = torch::ones({cellN, 1}, torch::kFloat);
auto nut = torch::zeros({cellN});
for (int i = 0; i < IterationNum; i++)
{
opti->zero_grad();
auto upred = model->forward(mesh);
auto dudy =
torch::autograd::grad
(
{upred},
{mesh},
{torch::ones_like(upred)},
true,
true
)[0];
auto dudyy =
torch::autograd::grad
(
{dudy},
{mesh},
{torch::ones_like(upred)},
true,
true
)[0];
auto dudyBottom = dudy[0][0];
auto dudyTop = dudy[cellN - 1][0];
auto meanU = torch::mean(upred);
auto ubottom = torch::full_like(upred[0], 0.0375);
auto dpdxByNu = torch::full({cellN}, 1.2e-05/1e-8);
//#include "nut40.H"
//auto dpdxByNu = 0.018/nut;
//auto ubottom = torch::full_like(upred[0], 0.77);
auto loss_bottom = crit(upred[0], ubottom);
auto loss_top = crit(dudyTop, 0.0*dudyTop);
auto loss_bd = loss_bottom + 1000*loss_top;
auto loss_ini = 1000*crit(meanU, init);
auto loss_pde = crit(-dudyy.reshape({cellN}), dpdxByNu);
auto loss = loss_pde + loss_bd + loss_ini;
//auto loss = loss_pde + loss_bd;
loss.backward();
opti->step();
#include "output.H"
iter++;
}
torch::save(model, "model.pth");
std::cout<< "Done!" << std::endl;
return 0;
}
boundaryFoam预测的速度原始数据：
0.0375
0.1107
0.1821
0.2516
0.3192
0.3849
0.4487
0.5106
0.5706
0.6288
0.6851
0.7394
0.7919
0.8425
0.8912
0.9380
0.9830
1.0260
1.0672
1.1064
1.1438
1.1793
1.2129
1.2446
1.2745
1.3025
1.3286
1.3528
1.3752
1.3957
1.4143
1.4310
1.4459
1.4589
1.4701
1.4794
1.4868
1.4924
1.4961
1.4980