Palabos 源码解读 | 浸入边界法Immersed Boundary Method

源码中单独的浸入边界法类

首先我们需要了解ib-lbm计算中，先（首次需initial）计算Lagrangian点的力g，再将g spreading on Eulerian格点，接下来执行碰撞和流动，然后通过interpolation得到Lagrangian的速度，最后更新Lagrangian格点位置。
中英穿插的描述是因为方便读者理解论文中用到的英语。

迭代过程中边界力的计算

> template<typename T> void
> ComputeImmersedBoundaryForce3D<T>::processGenericBlocks (
>         Box3D domain, std::vector<AtomicBlock3D*> blocks ) {
   
>     PLB_PRECONDITION( blocks.size()==2 );
>     TensorField3D<T,3>* force = dynamic_cast<TensorField3D<T,3>*>(blocks[0]);
>     AtomicContainerBlock3D* container = dynamic_cast<AtomicContainerBlock3D*>(blocks[1]);
>     PLB_ASSERT(force);
>     PLB_ASSERT(container);
>     Dot3D location = force->getLocation();
>     ImmersedWallData3D<T>* wallData = 
>         dynamic_cast<ImmersedWallData3D<T>*>( container->getData() );
>     PLB_ASSERT(wallData);
>     std::vector< Array<T,3> > const& vertices = wallData->vertices;
>     std::vector<Array<T,3> >& g = wallData->g;
>     PLB_ASSERT( vertices.size()==g.size() );
> 
>     for (plint iX = 0; iX < force->getNx(); iX++) {
   
>         for (plint iY = 0; iY < force->getNy(); iY++) {
   
>             for (plint iZ = 0; iZ < force->getNz(); iZ++) {
   
>                 force->get(iX, iY, iZ).resetToZero();
>             }
>         }
>     }
> 
>     for (pluint i=0; i<vertices.size(); ++i) {
   
>         Array<T,3> const& vertex = vertices[i];//遍历所有固体格点
>         Array<plint,3> intPos((plint) vertex[0] - location.x, (plint) vertex[1] - location.y, (plint) vertex[2] - location.z);
>         const Array<plint,2> xLim((vertex[0] < (T) 0 ? Array<plint,2>(-2, 1) : Array<plint,2>(-1, 2)));//判断距离，并设定样条
>         const Array<plint,2> yLim((vertex[1] < (T) 0 ? Array<plint,2>(-2, 1) : Array<plint,2>(-1, 2)));
>         const Array<plint,2> zLim((vertex[2] < (T) 0 ? Array<plint,2>(-2, 1) : Array<plint,2>(-1, 2)));
>         const Array<T,3> fracPos(util::frac(vertex[0]), util::frac(vertex[1]), util::frac(vertex[2]));
>         for (plint dx = xLim[0]; dx <= xLim[1]; dx++) {
   
>             for (plint dy = yLim[0]; dy <= yLim[1]; dy++) {
   
>                 for (plint dz = zLim[0]; dz <= zLim[1]; dz++) {
   
>                     Array<plint,3> pos(intPos+Array<plint,3>(dx,dy,dz));
>                     Array<T,3> r((T)dx-fracPos[0],(T)dy-fracPos[1],(T)dz-fracPos[2]);
>                     T W = inamuroDeltaFunction<T>().W(r);
>                     force->get(pos[0], pos[1], pos[2]) += W*g[i];
>                 }
>             }
>         }
>     } }

主要思路是以张量场的方式收集和计算区块的力。在IBM实现中，已定义rhobar，j，和container，定义的g即为力。

定义force为3D的张量场，定义location为force的位置，定义wallData为container里数据，定义vertices为wallData里的顶点数据（即流域内固体顶点位置），定义g为wallData里的力。

force是基于Eulerian mesh，而g基于Lagrangian mesh。

目测这个force就是作用在eulerian网格上的力，location也是eulerian网格坐标，vertices的顶点坐标vertex和力g都是lagragian网格上的参数。

> template<typename T> struct ImmersedWallData3D : public
> ContainerBlockData {
   
>     std::vector< Array<T,3> > vertices;
>     std::vector<T> areas;
>     std::vector< Array<T,3> > normals;
>     std::vector< Array<T,3> > g;
>     std::vector<int> flags; // Flag for each vertex used to distinguish between vertices for conditional reduction operations.
>     std::vector<pluint> globalVertexIds;
>     virtual ImmersedWallData3D<T>* clone() const {
   
>         return new ImmersedWallData3D<T>(*this);
>     } };

循环的第一步将force重置为0，此处我猜测是否了避免整个计算过程中力的累加。

第二个循环则依次在每个顶点上计算。并定义vertex为单个顶点，该顶点由0,1，2表示三个轴坐标。

<表达式1>?<表达式2>:<表达式3> ;

若表达式1为真，则执行表达式2；若表达式1的值为假，则执行表达式3。所以xLim，yLim，zLim的结果只有两种:

1. xLim[0]为-2，xLim[1]为1。
2. xLim[0]为-1，xLim[1]为2。

循环中的dx也是从-2，-1,0，1,2或者-1,0,1,2,3循环。

Array<plint,3> intPos((plint) vertex[0] - location.x, (plint) vertex[1] - location.y, (plint) vertex[2] - location.z);这段代码的作用是得到vertex的Lagrangian坐标，再四舍五入，减去location坐标，intPos得到值是Eulerian坐标，根据IBM作用邻域的概念判断，后续是为了确定计算邻域。

像dx-fracPos[0]即计算Dirac函数中的x方向距离。

Array<plint,3> pos(intPos+Array<plint,3>(dx,dy,dz));dx，dy，dz用于得到周围邻域被施加力的格点Eulerian坐标，由负到正，实际效果则是在这些Eulerian格点pos上作用。

force->get(pos[0], pos[1], pos[2]) += W*g[i];主要是实现下图的g的迭代累加（这里的g是Eulerian点的力）。

梳理一遍

1）得到Eulerian的force，得到边界Lagrangian坐标
2）利用整型四舍五入拉格朗日坐标来获得欧拉坐标，intPos=(plint) vertex[0] - location.x
3）确定邻域范围，利用pos(intPos+Array<plint,3>(dx,dy,dz))得到周围的欧拉坐标点
4）util::frac(vertex[0]得到小数值，(T)dx-fracPos[0]即计算距离r
5）计算得到体积力force->get(pos[0], pos[1], pos[2]) += W*g[i]

补档

文上几处我猜测的地方还未确定，等确定后会回来修改。

浸入边界的生成

> template<typename T>
> InstantiateImmersedWallData3D<T>::InstantiateImmersedWallData3D (//可见这里需要三个输入参数，其中最后一个是normals
>             std::vector< Array<T,3> > const& vertices_,
>             std::vector<T> const& areas_,
>             std::vector< Array<T,3> > const& normals_)//这个是算例中出现的container
>     : vertices(vertices_),
>       areas(areas_),
>       normals(normals_) {
   
>     PLB_ASSERT(vertices.size() == areas.size());
>     PLB_ASSERT(normals.size()==0 || normals.size() == areas.size()); }
// 这个rhobar来源于生成标量场
//		MultiScalarField3D<T> *rhoBar = generateMultiScalarField<T>((MultiBlock3D&) *lattice, 
//		largeEnvelopeWidth).release();
//    	rhoBar->toggleInternalStatistics(false);
//container以rhoBar生成
//		MultiContainerBlock3D container(*rhoBar);
> template<typename T> void
> InstantiateImmersedWallData3D<T>