学习笔记：Palabos（6）：数据处理器相关

数据处理器用来处理非局部操作。我们知道Dynamic处理局部操作，数据处理器来处理Dynamic完成不了的工作。数据处理器比较难，我之前以及接触到了一些，是一个非常庞大的东西。

鉴于能力有限，本文的讨论仅基于用户文档中关于数据处理的部分，以及参考了Palabos V2.0的开发文档。

1 数据处理器的作用

Dynamic类和数据处理器是Palabos构建程序的基础部分，文档上说是Palabos的螺钉螺母，就是这个意思。Dynamic是处理LBM单元计算，也就是局部操作：比如计算平衡态分布函数，碰撞等操作。而数据处理器是一个更加通用，且具有并行化特征的数据处理方法。只要不涉及边界，细化网格，并行程序或者其他一些Palabos里复杂的东西，Dynamic还是相对简单的。

数据处理器解决除了Dynamic解决不了的一切问题。处理模型中的非局部部分，在标量场和张量场上执行操作，并且处理多块网格之间的耦合，并且为了增加性能，它们必须高效而且满足并行化的特征。

因此数据处理器理解起来会比较困难，种类繁多，层层的迭代调用，并且包含了很多复杂的C++用法。好在不需要全部理解它们。Palabos的封装性做得很好，很大程度上我们需要去知道其调用方法即可以解决一般性的问题。

当在处理标量场和张量场时，数据处理器的重要性尤其明显。与block-lattice不同，这些场没有智能单元，比如lattice-block都有一个引用指向其dynamic。因此，数据处理器是对它们执行集体操作的唯一有效方式。比如考虑三维的速度张量场：TensorField3D velocity(其中每个单元都有Array 类型的速度变量)，可以通过调用函数computeVelocity从LB模拟获得这样的速度场。假使我们需要计算velocity的空间导数，那么只能用数据处理器来进行。因为数据处理器是并行的、可拓展的、高效的。例如computeVorticity()和computeStrainRate()那样。

2 数据处理器的分类

根据用途，数据处理器可分为三大类：

• 将动态对象分配给lattice的子域，初始化分布函数等等，这些方法在密度和速度的初始值和定义边界条件一节中进行了解释，这些函数在附录中中列出。
• 实现特定的物理模型，许多模型是预先设定好的。
• 用于数据评估，也就是后处理。比如computeVelocity这些

3 定义自己的数据处理器

3.1 用于局部操作的数据处理器

用于局部操作的数据处理器通常出现于Palabos内置的数据处理器或者Dynamic无法满足条件的情况。对于局部的数据操作，Palabos提供了一些抽象的接口：OneCellFunctionalXD，OneCellIndexedFunctionalXD。应该在这些接口上去实现你要的功能，也就是继承他们，并实现它的虚方法，而不是去重写这部分代码。

比如通过实现OneCellFunctionalXD可以执行局部化的操作，而不用访问周围的单元。如果操作取决于周围的单元，则可以继承OneCellIndexedFunctionalXD并实现该方法

这样的案例可以在examples / showCases / multiComponent2d中找到相应的示例程序。该示例程序中需要定制的初始化过程来访问lattice的外部标量以进行热仿真。

3.2 用于矩阵中的数据处理器

当处理的对象不再是一个局部的单元而是一个矩阵时，Palabos也提供了相应的接口，类似于3.1节的那样，是一些抽象的接口。为了实现你自己的功能，需要继承他们并且实现里面的方法。

在矩阵上执行操作一般是在矩阵的所有元素执行某个循环来达到目的。如果矩阵的内存被分为多块，就像Palabos里的Multiblock一样，并且这些组件分布在并行机器的节点上，这样你的循环也需要分为几个小的循环。Palabos的数据处理器的目的就是为了自动执行这些细分，然后它提供单个域的坐标，并且要求在这些单个的域上执行操作，而不是整个完整的域。

下面是一个例子，假设要操作一个域中单元格的分布函数，比如将f5和f1的值互换:

template<typename T, template<typename U> class Descriptor>
class Invert_1_5_Functional2D : public BoxProcessingFunctional2D_L {
    public:
    // ... implement other important methods here.
    void process(Box2D domain, BlockLattice2D& lattice)
    {
        for (plint iX=domain.x0; iX<=domain.x1; ++ix) {< span>            for (plint iY=domain.y0; iY<=domain.y1; ++iy) {< span>                Cell& cell = lattice.get(iX,iY);                // Use the function swap from the C++ STL to swap the two values.                std::swap(cell[1], cell[5]);            }        }    }};
            for (plint iY=domain.y0; iY<=domain.y1; ++iy) {< span>                Cell& cell = lattice.get(iX,iY);                // Use the function swap from the C++ STL to swap the two values.                std::swap(cell[1], cell[5]);            }        }    }};
                Cell& cell = lattice.get(iX,iY);
                // Use the function swap from the C++ STL to swap the two values.
                std::swap(cell[1], cell[5]);
            }
        }
    }
};

从上面可以看到，Invert_1_5_Functional2D这个类继承自BoxProcessingFunctional2D_L。函数的第一个参数是该函数作用的域。一般是分块之后的子域。同样第二个参数也是对应的一个格子块，通常是原子块（block-lattice）而不是一个多块（multi-block-lattice）。因为process方法，通常作为一个域处理器。

在这个方法中，值得注意的是内部的循环需要自己写。而cell级别的操作不一样。cell只需要处理一个单元的信息。因此，process的效率更高。

此外，除了局部操作外，你还可以通过BoxProcessingFunctional2D_L进行非局部操作，比如把f1和右边的f5互换：

void process(Box2D domain, BlockLattice2D& lattice)
{
    for (plint iX=domain.x0; iX<=domain.x1; ++ix) {< span>        for (plint iY=domain.y0; iY<=domain.y1; ++iy) {< span>            Cell& cell = lattice.get(iX,iY);            Cell& partner = lattice.get(iX+1,iY);            std::swap(cell[1], partner[5]);        }    }}
        for (plint iY=domain.y0; iY<=domain.y1; ++iy) {< span>            Cell& cell = lattice.get(iX,iY);            Cell& partner = lattice.get(iX+1,iY);            std::swap(cell[1], partner[5]);        }    }}
            Cell& cell = lattice.get(iX,iY);
            Cell& partner = lattice.get(iX+1,iY);
            std::swap(cell[1], partner[5]);
        }
    }
}

在用这种非局部操作的时候，为了保证在操作的时候不要超过lattice的索引范围，要遵守两个规则：

1. 非局部操作不能超越相邻2个格子，即只能与附近的格子进行数据操作。比如你可以写lattice.get(iX+1,iY)这样的，但是不能有lattice.get(iX+2,iY)这样的；
2. 在通信范围内作用的数据处理器中不允许进行非本地操作。

这里的通信范围（communication envelope）指的是一个block外面包含的一层通信用的节点。

总结一下数据处理器的功能：

• 访问一个块方位内的单元，可以读取单元的变量并且改变他们
• 可以根据索引访问单元及邻近的单元
• 对单元内数据的操作可能与工件有关，因此数据处理器也具有空间依赖性，但是需要注意的是索引时需要采用绝对坐标。

4 需要覆盖的其他方法

除了process需要覆盖以外，还有一些其他的方法也是必须要重载的。还是以3.2节中的例子Invert_1_5_Functional2D
为例，以下三个方法需要重载：

Invert_1_5_Functional2D* clone() const
{
    return new Invert_1_5_Functional2D(*this);
}


void getModificationPattern(std::vector<bool>& isWritten) const
{
    isWritten[0] = true;
}

BlockDomain::DomainT appliesTo() const
{
    return BlockDomain::bulk;
}

第一个方法是一个克隆方法的继承。

第二个方法中，isWritten是一个bool类型的向量，这个字段等于所涉及的块的数量。而这个函数的作用就是为这个字段赋值true/false。在例子中的情况是一个block，并且该方法将这个block的isWritten赋予true值。Palabos利用该信息来决定在执行数据处理器之后是否需要块的进程间通信。

第三个方法中，方法applyTo用于决定数据处理器是仅作用于模拟的实际域（BlockDomain :: bulk）还是包括了通信包络和实际域的情况（BlockDomain :: bulkAndEnvelope）。其中BlockDomain是一个枚举类型，其值如下：

值	含义
bulk	Refers to bulk-nodes, without envelope
envelope	Refers to nodes on the envelope, without the bulk
bulkAndEnvelope	Refers to the full domain

这里的bulk应该就是block-lattice内部的节点，envelope是外部的节点，bulkAndEnvelope是两者都包含的全部的节点。通常情况下，选择bulk就可以了因为包含envelope会在一些非局部数据操作上产生索引溢出的问题。不过envelope的存在也在一定程度上能防止产生索引溢出。

不过envelope的设置并不是为了防止非局部操作的索引溢出而存在的，envelope的设置是为了多块之间的通信，以在需要的时候更新envelope这部分的值。只有在两种情况下，会需要到envelope。第一种是在为一些或者所有单元分配新的Dynamic时，会需要用到envelope；第二种是在修改dynamic的内部状态，就像调用函数defineVelocity在Dirichelet边界节点上分配新的速度值那样。在这些情况下，

只有在（1）为一些或所有单元格分配新的动力学对象时才需要包含包络（就像在调用函数时一样）defineDynamics），或者（2）如果修改动态对象的内部状态（就像调用函数 defineVelocity在Dirichlet边界节点上分配新的速度值时那样）。在这些情况下，包括包络是必要的，因为在原子块之间的通信步骤期间不转移动态对象的性质和内容。传递的唯一信息是单元数据（分布函数和外部标量）。

如果决定要把envelope包含在数据处理器的方法中，则必须要遵守两条规则，否则，将出现未定义的行为

1. 数据处理器必须是局部的，没有额外的envelope用于处理非局部数据访问；
2. 数据处理器最多仅可以对一个涉及的块具有写访问权，即从方法getModificationPattern（）返回的向量isWritten最多只可以在一个地方为true。

关于第二条规则，如上面的例子所示，仅仅让isWritten(0)有了真值。不过我看了下文档，发现了在BoxProcessingFunctional2D中的getModificationPattern函数定义如下：

void ‪BoxProcessingFunctional2D::getModificationPattern(std::vector<bool>& isWritten) const {
     std::vector modified(isWritten.size());
     ‪getTypeOfModification(modified);
     ‪PLB_ASSERT(modified.size()==isWritten.size());
     for (‪pluint iBlock=0; iBlock         isWritten[iBlock] = modified[iBlock]==‪modif::nothing ? false:true;     } }
         isWritten[iBlock] = modified[iBlock]==‪modif::nothing ? false:true;
     }
 }

可以看到，其实在这个里面的getModificationPattern()函数可以不止对一个块赋予真值。所以这一点有点不太清楚。

总而言之，数据处理器中的appliesTo函数，通常还是让啊选择bulk即可。

5 绝对和相对坐标

当process作用在一个原子块（atomic-block）的时候，单元在各自块里的索引(iX, iY)只是在做内部循环的时候有用。因为他们只是代表局部的变量。如果要涉及全局坐标的话，需要对这个相对的坐标做一个转换：

// Access the position of the atomic-block inside the multi-block.
Dot2D relativePosition = lattice.getLocation();

// Convert local coordinates to global ones.
plint globalX = iX + relativePosition.x;
plint globaly = iY + relativePosition.x;

根据这个案例可以看到，原子块本身应该是具有坐标的。这个案例可以在examples/showCases/boussinesqThermal2d/中找到。这种转换有些尴尬，所以对单个单元的局部操作，采用OneCellIndexedFunctionalXD这样的接口可能更好。

当一个process作用在多个块的时候，相对坐标没有必要完全一样。方法process的参数domian始终作为第一个原子块的本地坐标提供。要获取其他块的坐标，必须应用相应的转换：

Dot2D offset_0_1 = computeRelativeDisplacement(lattice0, lattice1);
Dot2D offset_0_2 = computeRelativeDisplacement(lattice0, lattice2);

plint iX1 = iX + offset_0_1.x;
plint iY1 = iY + offset_0_1.y;
plint iX2 = iX + offset_0_2.x;
plint iY2 = iY + offset_0_2.y;

这里涉及到了三个原子块，lattice0和lattice1，以及lattice3。而函数computeRelativeDisplacement则是计算两个格子块的相对位移。这样每个块里面的坐标就可以求出来了，值得注意的是，计算位移通常出现在以下三种情况：

1. 应用数据处理器的多块不具有相同的数据分布，因为它们的构造不同。
2. 应用数据处理器的多块不具有相同的数据分布，因为它们的大小不同。例如computeVelocity，它计算子域上的速度。因为这个数据处理器作用在原始域和对应的子域上。
3. 数据处理器包括envelope。在这种情况下，相对位移源于体节点与来自不同原子块的包络节点耦合。这是为什么通常最好不将envelope包括在数据处理器的应用领域中的另一个原因。

6 未完待续

果然比较复杂。

这个用户文档都看的云里雾里的。。。整篇文章不做任何参考。我估计会有很多理解上的错误和偏颇。后续了解了再来填坑。

不能再啃下去了，要结合例子看才会有效果。

这部分工作先暂时这样。包括用户文档的阅读。

下一步要仔细研究几个案例。

先过一遍Tutorials。

接下来就是Guo的案例。