ML: 代理模型-数据驱动LES (下)
本文代码进行小改，就可以修修补补做一些工作，或许有机会刊发在流体力学顶刊JFM（因为JFM在2023年，2024年发表过完全一样的工作）。本文将JFM的训练工作进行了复现。
训练神经网络的思想主要是，给定一个LES算例，摘取其中某些部分的网格点，提取相应的速度场或者相关数据构建输入参数，输入参数即可以为湍流粘度。对大量的这种一对一的对应关系进行拟合，既可以训练出一个神经网络来替代湍流模型。
步骤1，书接上文，首先准备好一个算例，算例请点击下载。本算例在trainingProperties里面有CFD以及training关键词，在准备数据阶段，CFD关键词要为true，training关键词要为false。首先需要用标准求解器对该算例使用pisoFoam进行计算，获取大量的时间步数据结果。在进入下一步之前，确保算例下有足够的时间步数据。
步骤2，讲下面的代码，复制到pisoFoam求解器中，重命名为lesMLFoam。下面的代码主要可以用于神经网络的训练：
Foam::instantList timeDirs = Foam::timeSelector::select0(runTime, args);
// 0 folder is not included
int BATCH = timeDirs.size() - 1;
Info<< "Batch size: " << BATCH << nl;
scalarList timeU;
forAll(timeDirs, timei)
{
if (timei != 0)
{
runTime.setTime(timeDirs[timei], timei);
timeU.append(runTime.time().value());
}
}
Info<< "Training time list: " << timeU << nl;
runTime.setTime(0, 0);
labelList cellID;
const meshCellZones& cellZones = mesh.cellZones();
if (cellZones.size() == 0)
{
FatalErrorInFunction
<< "No cell zones found, running topoSet??" << abort(FatalError);
}
else
{
forAll(cellZones[0], i)
{
label realID = cellZones[0][i];
cellID.append(realID);
}
}
auto x = torch::full({BATCH, cellID.size(), 9}, 0.0);
auto y = torch::full({BATCH, cellID.size(), 1}, 0.0);
std::cout<< "\nx shape is " << x.sizes() << std::endl;
//Info<< "cell ID is " << cellID << nl;
// input
PtrList<tmp<volTensorField>> gradUlist(BATCH);
PtrList<tmp<volScalarField>> ux(BATCH);
PtrList<tmp<volScalarField>> uy(BATCH);
PtrList<tmp<volScalarField>> uz(BATCH);
PtrList<tmp<volScalarField>> vx(BATCH);
PtrList<tmp<volScalarField>> vy(BATCH);
PtrList<tmp<volScalarField>> vz(BATCH);
PtrList<tmp<volScalarField>> wx(BATCH);
PtrList<tmp<volScalarField>> wy(BATCH);
PtrList<tmp<volScalarField>> wz(BATCH);
// output
PtrList<tmp<volScalarField>> nutlist(BATCH);
forAll(gradUlist, i)
{
gradUlist.set
(
i,
new tmp<volTensorField>
(
new volTensorField
(
IOobject
(
"gradU",
name(timeU[i]),
mesh,
IOobject::MUST_READ,
IOobject::NO_WRITE
),
mesh
)
)
);
nutlist.set
(
i,
new tmp<volScalarField>
(
new volScalarField
(
IOobject
(
"nut",
name(timeU[i]),
mesh,
IOobject::MUST_READ,
IOobject::NO_WRITE
),
mesh
)
)
);
ux.set
(
i,
new tmp<volScalarField>
{
new volScalarField
(
IOobject
(
"ux",
name(timeU[i]),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
gradUlist[i]().component(tensor::XX)
)
}
);
uy.set
(
i,
new tmp<volScalarField>
{
new volScalarField
(
IOobject
(
"uy",
name(timeU[i]),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
gradUlist[i]().component(tensor::XY)
)
}
);
uz.set
(
i,
new tmp<volScalarField>
{
new volScalarField
(
IOobject
(
"uz",
name(timeU[i]),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
gradUlist[i]().component(tensor::XZ)
)
}
);
vx.set
(
i,
new tmp<volScalarField>
{
new volScalarField
(
IOobject
(
"vx",
name(timeU[i]),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
gradUlist[i]().component(tensor::YX)
)
}
);
vy.set
(
i,
new tmp<volScalarField>
{
new volScalarField
(
IOobject
(
"vy",
name(timeU[i]),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
gradUlist[i]().component(tensor::YY)
)
}
);
vz.set
(
i,
new tmp<volScalarField>
{
new volScalarField
(
IOobject
(
"vz",
name(timeU[i]),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
gradUlist[i]().component(tensor::YZ)
)
}
);
wx.set
(
i,
new tmp<volScalarField>
{
new volScalarField
(
IOobject
(
"wx",
name(timeU[i]),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
gradUlist[i]().component(tensor::ZX)
)
}
);
wy.set
(
i,
new tmp<volScalarField>
{
new volScalarField
(
IOobject
(
"wy",
name(timeU[i]),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
gradUlist[i]().component(tensor::ZY)
)
}
);
wz.set
(
i,
new tmp<volScalarField>
{
new volScalarField
(
IOobject
(
"wz",
name(timeU[i]),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
gradUlist[i]().component(tensor::ZZ)
)
}
);
Info<< "Assign values for ML at t: " << timeU[i] << nl;
forAll(cellID, j)
{
// i is time, j is line, k is 9
label C = cellZones[0][j];
x[i][j][0] = ux[i]()[C];
x[i][j][1] = uy[i]()[C];
x[i][j][2] = uz[i]()[C];
x[i][j][3] = vx[i]()[C];
x[i][j][4] = vy[i]()[C];
x[i][j][5] = vz[i]()[C];
x[i][j][6] = wx[i]()[C];
x[i][j][7] = wy[i]()[C];
x[i][j][8] = wz[i]()[C];
y[i][j][0] = 1e5*nutlist[i]()[C];
}
gradUlist[i].clear();
ux[i].clear();
uy[i].clear();
uz[i].clear();
vx[i].clear();
vy[i].clear();
vz[i].clear();
wx[i].clear();
wy[i].clear();
wz[i].clear();
}
for (int epoch = 0; epoch < iter; epoch++)
{
auto yPred = net->forward(x);
auto loss = crit(yPred, y);
opti.zero_grad();
loss.backward();
opti.step();
if (epoch % 20 == 0)
{
std::cout << "Iter: " << epoch << ", Loss: "
<< loss.item<float>() << std::endl;
torch::save(net, "net.pth");
}
}
上述的代码需要嵌入在pisoFoam求解器的CFD求解代码之前。同时将CFD求解部分的代码注释掉，这样在步骤1里面的算例文件下，运行lesMLFoam（记得运行前运行topoSet），就可以进行训练，同时输出net.pth文件。
步骤3，上述的net.pth文件可以用在第一步进行包含，在后续的使用中，参考ML: 代理模型-数据驱动LES (上)。
至此，在OpenFOAM调用数据驱动LES的工作介绍完毕。总结如下：
ML: 代理模型-数据驱动LES (上) 创建了一个求解器框架，调用训练好的net.pth文件来作为代理模型，替换LES模型；
ML: 代理模型-数据驱动LES (下) 讨论了如何对现有的算例进行训练，来获得net.pth文件；
下图是采用训练好的神经网络进行的推理。是一个全新的几何（周期边界的方柱扰流）。数据驱动LES模型计算的结果也与结果高度吻合。
如果对源代码理解有困难，可以在OpenFOAM全序列系统里面下载虚拟机获得完整版本。