OpenFOAM-v2006重叠网格挖洞问题研究以及overset代码解析(八)walkFront函数

在我之前的博客里最后提到一个解决方案是对walkFront函数进行改动的。一个叫louisgag的前辈在github上贴了他的代码。这次想要学习walkFront函数以及前辈的nPushWalkFront函数,为了弄明白这个函数在整一个重叠网格中担当的角色。这一函数比较重要,会尽量讲解每一行代码的意思(我理解的,如果有误欢迎指正!)

看到walkFront的参数walkFront(layerRelax, allStencil, allCellTypes, allWeight);其函数定义时参数: 

其中,layerRelax是用户定义的标量,起作用的地方在后面的seedCell函数里。allStencil是上一篇博文提到的markDonors输出的下标列表类型变量。allCellType是在markPatchesAsHoles(在我前面五篇博文里有介绍)这个函数里输出的,在这个函数中还会进行进一步的改动。最后是allWeight,是一个标量场,其在调用walkFront前被定义:scalarField allWeight(mesh_.nCells(), Zero);此标量场是储存每一个网格的权重信息的,至于这个权重是用来干什么的,接着看下去。

    // Current front
    bitSet isFront(mesh_.nFaces());

    const fvBoundaryMesh& fvm = mesh_.boundary();

fvm调用了网格的边界mesh.boundary();isFront则是用于标记挖洞之后的HOLE外面一圈的面。bitSet相当于一个只含有真假两种状态的列表。这里IsFront初始化时容量设置成所有面的数量,并初始化为0,对于需要的面会在后面标记为1。Front是一层需要被标记的插值层(buffer interpolation layer(s) around holes)。需要标记的区域有两类,一是OVERSET边界上的面,二是背景网格挖洞之后HOLE外面的区域。 

1. OVERSET PATCH:

    // 'overset' patches

    forAll(fvm, patchI)
    {
        if (isA<oversetFvPatch>(fvm[patchI]))
        {
            const labelList& fc = fvm[patchI].faceCells();
            forAll(fc, i)
            {
                label cellI = fc[i];
                if (allCellTypes[cellI] == INTERPOLATED)
                {
                    // Note that acceptors might have been marked hole if
                    // there are no donors in which case we do not want to
                    // walk this out. This is an extreme situation.
                    isFront.set(fvm[patchI].start()+i);
                }
            }
        }
    }

如果边界上的某个patch为overset类型,那么这个patch上的面会被标记成Front。这个区域是两层网格的交界处,所以即使没有挖洞也必须标记为Front。这里还另外判断面上对应的cell是不是别标记为INTERPOLATED,是为了防止出现极端的情况:如果某一单元没有别哪个单元贡献的话,那么这个单元可能会被标记为HOLE,这种情况可能出现在重叠网格的外边界陷入到wall边界里面时。

2. outside of HOLE region:

    // Outside of 'hole' region
    {
        const labelList& own = mesh_.faceOwner();
        const labelList& nei = mesh_.faceNeighbour();

        for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
        {
            label ownType = allCellTypes[own[faceI]];
            label neiType = allCellTypes[nei[faceI]];
            if
            (
                 (ownType == HOLE && neiType != HOLE)
              || (ownType != HOLE && neiType == HOLE)
            )
            {
                //Pout<< "Front at face:" << faceI
                //    << " at:" << mesh_.faceCentres()[faceI] << endl;
                isFront.set(faceI);
            }
        }

        labelList nbrCellTypes;
        syncTools::swapBoundaryCellList(mesh_, allCellTypes, nbrCellTypes);

        for
        (
            label faceI = mesh_.nInternalFaces();
            faceI < mesh_.nFaces();
            faceI++
        )
        {
            label ownType = allCellTypes[own[faceI]];
            label neiType = nbrCellTypes[faceI-mesh_.nInternalFaces()];

            if
            (
                 (ownType == HOLE && neiType != HOLE)
              || (ownType != HOLE && neiType == HOLE)
            )
            {
                //Pout<< "Front at coupled face:" << faceI
                //    << " at:" << mesh_.faceCentres()[faceI] << endl;
                isFront.set(faceI);
            }
        }
    }

这一块是标记背景网格里面HOLE外面一圈的Front的,其判断主要用:

if((ownType == HOLE && neiType != HOLE)
        || (ownType != HOLE && neiType == HOLE))

来实现。上部分是对背景网格内部面进行遍历的,而下部分则是对背景网格边界面进行遍历,总之每一次都要遍历背景网格的所有的面。这样,isFront就标记完毕,接下来是处理权重。

    // Current interpolation fraction
    scalarField fraction(mesh_.nFaces(), Zero);

    forAll(isFront, faceI)
    {
        if (isFront.test(faceI))
        {
            fraction[faceI] = 1.0;
        }
    }

这一段定义了一个标量场fraction,其大小和网格面元数量相等,且初始化为0,并根据isFront来设置某一个face上的值为1。这里用到的test是bool返回类型函数,返回faceI位置时的位值。在前面的代码中特定的面上isFront是1,所以这里这些面对应的fraction也会设置为1。

while (returnReduce(isFront.any(), orOp<bool>()))
    {
        // Interpolate cells on front
        bitSet newIsFront(mesh_.nFaces());
        scalarField newFraction(fraction);
        forAll(isFront, faceI)
        {
            if (isFront.test(faceI))
            {
                label own = mesh_.faceOwner()[faceI];
                if (allCellTypes[own] != HOLE)
                {
                    if (allWeight[own] < fraction[faceI])
                    {
                        // Cell wants to become interpolated (if sufficient
                        // stencil, otherwise becomes hole)
                        if (allStencil[own].size())
                        {
                            allWeight[own] = fraction[faceI];
                            allCellTypes[own] = INTERPOLATED;
                            // Add faces of cell (with lower weight) as new
                            // front
                            seedCell
                            (
                                own,
                                fraction[faceI]-layerRelax,
                                newIsFront,
                                newFraction
                            );
                        }
                        else
                        {
                            allWeight[own] = 0.0;
                            allCellTypes[own] = HOLE;
                            // Add faces of cell as new front
                            seedCell
                            (
                                own,
                                1.0,
                                newIsFront,
                                newFraction
                            );
                        }
                    }
                }
                if (mesh_.isInternalFace(faceI))
                {
                    label nei = mesh_.faceNeighbour()[faceI];
                    if (allCellTypes[nei] != HOLE)
                    {
                        if (allWeight[nei] < fraction[faceI])
                        {
                            if (allStencil[nei].size())
                            {
                                allWeight[nei] = fraction[faceI];
                                allCellTypes[nei] = INTERPOLATED;
                                seedCell
                                (
                                    nei,
                                    fraction[faceI]-layerRelax,
                                    newIsFront,
                                    newFraction
                                );
                            }
                            else
                            {
                                allWeight[nei] = 0.0;
                                allCellTypes[nei] = HOLE;
                                seedCell
                                (
                                    nei,
                                    1.0,
                                    newIsFront,
                                    newFraction
                                );
                            }
                        }
                    }
                }
            }
        }

        syncTools::syncFaceList(mesh_, newIsFront, orEqOp<unsigned int>());
        syncTools::syncFaceList(mesh_, newFraction, maxEqOp<scalar>());

        isFront.transfer(newIsFront);
        fraction.transfer(newFraction);
    }

这是一个while的循环块:while (returnReduce(isFront.any(), orOp<bool>())),isFront.any()是bool返回函数,当isFront里面只要有一个值是1时就会返回1;orOp<bool>()或运算,而returnReduce则是// Reduce using either linear or tree communication schedule ,实际上还是一个reduce函数。因为本人代码基础不佳,这里不能完全解读开来。但毕竟是一个while循环,我想当isFront所有的位都置0的时候,这一循环就终止了。为什么要有这个循环呢?我们先看里面的内容。它先是定义了两个新的量newIsFront和newFraction。可以推断出这个是用来替代原来的变量的。forAll(isFront, faceI),对所有的面进行遍历,然后进行四层if判断,分别是:

1. 该面上对应的isFront的值是否为1?

2. 面的owner对应的cell的类型是不是非HOLE?或者该面是不是内部面?

3. allWeight是否小于fraction?

4. allStencil是否有大小?

newIsFront和newFraction会在if内的seedCell上进行改动。seedCell代码如下:

void Foam::cellCellStencils::inverseDistance::seedCell
(
    const label cellI,
    const scalar wantedFraction,
    bitSet& isFront,
    scalarField& fraction
) const
{
    const cell& cFaces = mesh_.cells()[cellI];
    forAll(cFaces, i)
    {
        label nbrFacei = cFaces[i];
        if (fraction[nbrFacei] < wantedFraction)
        {
            fraction[nbrFacei] = wantedFraction;
            isFront.set(nbrFacei);
        }
    }
}

结合函数调用时的参数,可以看到,最后对应面的newIsFront会被set,newFraction也会设为wantedFraction,这个值在调用时的第二个参数里。这个值被设置成fraction[faceI]-laterRelax或1。事实上,laterRelax一般默认为一个极小的值。在seedCell函数里,调用了这个体元的所有面,如果其中一个面的fraction低于wantedFraction,那么这个面的fraction就会等于wantedFraction,也就是说取二者较大的值,且此时newIsFront也会被set。

总的来看,这个walkFront函数的主角是面元face,而我们要进行标记的是cell。每一个face都会对应两个cell分别是own和nei,判断其是否要进行插值运算,则是根据allStencil是否有大小。而allStencil是在markDonors进行定义的。在markDonors里,如果两层网格里的两个单元相互有映射关系且其中任何一个都不是HOLE类型,那么就会有一个donor,也就是说allStencil就会有大小,在这里加上isFront的条件,使得这一个单元会被标记为INTERPOLATED类型,作为HOLE外层或者overset类型边界的插值单元,通过各种插值格式进行不同网格层之间的信息传输。

在循环体结尾还有四行代码:

syncTools::syncFaceList(mesh_, newIsFront, orEqOp<unsigned int>());
syncTools::syncFaceList(mesh_, newFraction, maxEqOp<scalar>());

isFront.transfer(newIsFront);
fraction.transfer(newFraction);

syncFaceList是 //- Synchronize values on all mesh faces.,然后transfer是把newIsFront替换掉isFront。至于为什么要这样做,以及为什么要进行循环,我想可能是和程序的优化有关系,这就不深究了。

挖洞的问题在各论坛里反响已久。其中louisgag提出了自己的解决方案:通过修改walkFront函数来使得挖洞可以往外继续进行。本篇博文介绍他的做法。

放上相关的代码:

int iPushFront = nPushFront;

    while (iPushFront > 0 && returnReduce(isFront.any(), orOp<bool>()))
    {
        // Interpolate cells on front
        bitSet newIsFront(mesh_.nFaces());
        //~ scalarField newFraction(fraction);
        forAll(isFront, faceI)
        {
            if (isFront.test(faceI))
            {
                label own = mesh_.faceOwner()[faceI];
                if (allCellTypes[own] != HOLE)
                {
                    //~ if (allWeight[own] < fraction[faceI])
                    //~ {
                        // Cell wants to become interpolated (if sufficient
                        // stencil, otherwise becomes hole)
                        //~ if (allStencil[own].size())
                        //~ {
                            //~ allWeight[own] = fraction[faceI];
                            allCellTypes[own] = HOLE;
                            newIsFront.set(mesh_.cells()[own]);
                            // Add faces of cell (with lower weight) as new
                            // front
                            //~ seedCell
                            //~ (
                                //~ own,
                                //~ fraction[faceI]-layerRelax,
                                //~ newIsFront,
                                //~ newFraction
                            //~ );
                        //~ }
                        //~ else
                        //~ {
                            //~ allWeight[own] = 0.0;
                            //~ allCellTypes[own] = HOLE;
                            //~ // Add faces of cell as new front
                            //~ seedCell
                            //~ (
                                //~ own,
                                //~ 1.0,
                                //~ newIsFront,
                                //~ newFraction
                            //~ );
                        //~ }
                    //~ }
                }
                if (mesh_.isInternalFace(faceI))
                {
                    label nei = mesh_.faceNeighbour()[faceI];
                    if (allCellTypes[nei] != HOLE)
                    {
                        //~ if (allWeight[nei] < fraction[faceI])
                        //~ {
                            //~ if (allStencil[nei].size())
                            //~ {
                                //~ allWeight[nei] = fraction[faceI];
                                allCellTypes[nei] = HOLE;
                                newIsFront.set(mesh_.cells()[nei]);
                                //~ seedCell
                                //~ (
                                    //~ nei,
                                    //~ fraction[faceI]-layerRelax,
                                    //~ newIsFront,
                                    //~ newFraction
                                //~ );
                            //~ }
                            //~ else
                            //~ {
                                //~ allWeight[nei] = 0.0;
                                //~ allCellTypes[nei] = HOLE;
                                //~ seedCell
                                //~ (
                                    //~ nei,
                                    //~ 1.0,
                                    //~ newIsFront,
                                    //~ newFraction
                                //~ );
                            //~ }
                        //~ }
                    }
                }
            }
        }

        syncTools::syncFaceList(mesh_, newIsFront, orEqOp<unsigned int>());
        //~ syncTools::syncFaceList(mesh_, newFraction, maxEqOp<scalar>());

        isFront.transfer(newIsFront);
        //~ fraction.transfer(newFraction);
        iPushFront -= 1;
    }

这个实现起来其实很简单,就是用户设置往外挖的层数iPushFront,在每一次循环时,对HOLE外围的cell继续设置成HOLE,由此实现了向外挖洞,而且不影响其他的东西,十分巧妙。

© 版权声明
THE END
喜欢就支持一下吧
点赞0 分享
评论 抢沙发
头像
欢迎您留下宝贵的见解!
提交
头像

昵称

取消
昵称表情代码图片