Palabos基本的数据类型，来自于文档的翻译。

包含： BlockXD，格子描述符，Dynamic，数据处理器

上周末做汇报，其实材料已经准备好了，没有发上来。这个月要把圆柱扰流的代码搞定，还是任务艰巨。。。

1 BlockXD

BlockXD构造表示常规的2D或3D数组（或矩阵）数据，是保存变量的基本结构。关于BlockXD，在之前的文章中介绍了BlockLattice相关的类。但其实BlockXD的子类非常丰富。除了网格变量之外，还有标量张量等。下图是Block结构的继承图。可以看得很清楚。

Block结构的继承图
Block结构的继承图

由图可知，BlockXD派生出了两个大类，一个是AtomicBlockXD，一个是MultiBlockXD。而事实上，由BlockXD还派生出来了MultiGridXD，这一点可以从开发文档里看到。同时，还派生出了一些背景类（或者说是基类，如LatticeBaseXD等）。可以看到，所有类型为“block-lattice”的块都有一个方法collideAndStream（），无论它们是作为multi-block还是atomic-block实现，这个collideAndStream()方法执行格子的碰撞和迁移。 同样，所有multi-block都有一个方法getComponent（），无论它们是block-lattice, scalar-field, 还是tensor-field。

此图中显示的情况称为“多重继承”，因为最终用户类以两种方式从BlockXD继承：其一是通过两大类atomic-block和multi-block路径，其二是背景类，即通过LatticeBaseXD，ScalarFieldBaseXD，TensorFieldBaseXD而来。比如BlockLattice2D这个类，就集成自BlockLatticeBase2D和AtomicBlock2D，而AtomicBlock2D又集成自Block2D。如下图：

BlockLattice2Ｄ

可见，BlockXD提供了一切关于求解域，求解域的变量存储等相关的类，并且值得注意的是，其中，BlockLattice这些类具有collide和stream函数，来执行参数的推演。是BlockXD里面最为复杂和重要的类。

2 格子描述子

所有BlockXD构造都根据底层数据类型进行模板化。比如下面的代码：

1
2
3
4
5
6
7
// Construct a 100x100 scalar-field with double-precision floating point values.

MultiScalarField2D<double> a(100,100);

// Construct a 100x100 scalar-field with single-precision floating point values.

MultiScalarField2D<float> b(100,100);

但是对于block-lattice来说，格子还具有其他的模板参数，其中最为重要的是格子描述符，其指定格子的一些拓扑性质（粒子群的数量，离散速度，方向的权重和其他格子常数）。因此，很容易在应用程序中尝试不同的格子：

1
2
3
4
5
6
7
// Construct a 100x100 block-lattice using the D3Q19 structure.

MultiBlockLattice2D<double, D3Q19Descriptor> lattice1(100,100);

// Construct a 100x100 block-lattice using the D3Q27 structure.

MultiBlockLattice2D<double, D3Q27Descriptor> lattice2(100,100);

文档中写道：编写新的网格描述符也很容易（目前没有记录，但你可以查看目录src / latticeBoltzmann / nearestNeighborLattices2D.h中的文件，看看它是如何工作的）。 这非常有用，因为这意味着当您切换到新类型的晶格时，不需要重新编写用于实现BlockLatticeXD的冗长代码部分。 这个论点是Palabos非侵入式程序开发部分中描述的一般概念的一部分。最近邻2D格的描述符。 原则上，由于代码的静态通用性，格子的定义和动态的实现是独立的。 当然，有特例。 例如，使用诸如热通量的高阶矩的动力学需要具有扩展邻域的格子。 另一个例子是D3Q13格子，它只能用它自己的动力学。

不难看出，邻近2D网格描述符（nearestNeighborLattices2D）是通常使用的网格描述符。我顺便查了查nearestNeighborLattices2D这个文件，发现其下面包含了常见的网格描述符：D2Q9Constants、D2Q9DescriptorBase、D2Q9Descriptor、ForcedD2Q9Descriptor等。基本上常见的2D格子描述符都在这里面，其中D2Q9Constants作为父结构体而存在。其中包含了成员变量，里面大部分都是静态的成员变量，没有看到方法。

D2Q9Descriptor
D2Q9Descriptor

如图是D2Q9描述符的成员变量，可以看到，都是D2Q9的一些固有属性的变量。因此，如果需要开发新的描述符的话，完全可以在此基础上继承。

不过值得注意的是，网格描述符属于结构体模板，并非类模板。

3 动态类（Dynamics）

动态类代表了采用的LBM模型。比如是LBGK模型还是MRT模型，或者是包含外力的Guo格式模型。

不同的模型采用不同的平衡态分布函数，乃至离散方式。而动态类Dynamic就是

  • 以下大部分来自于文档的翻译：

在格子Boltzmann模拟的时间迭代期间，block-lattice的所有单元执行局部碰撞步骤，然后执行扩散步骤。扩散步骤是固定的，并且只能通过定义网格描述符中（比如D2Q9Descriptor）的离散速度来影响。然而，碰撞步骤可以完全定制，并且可以实现在不同的单元上执行不同的碰撞方式。通过这种方式，可以在局部调整在网格上模拟的物理特性，并且诸如边界的特定区域可以得到单独的处理。

由于上述目的，比如说要针对不同的边界条件或动力学类型（Dynamics）来执行不同的碰撞模型。那么需要在网格的每个单元保持一个指向Dynamics类型的对象的指针。这个对象提供了碰撞步骤的实现。此外，它还提供了计算宏观变量以及更多有关模型的信息的方法。除了模拟的变量之外，block-lattice的单元格还包含指向Dynamics类型的对象的指针。 该对象最提供了冲突步骤的实现。 此外，它还提供了一种计算宏观变量的方法，用于在不同大小的网格之间重新调整变量的信息，以及更多与模型相关的信息。因此，在Dynamics类中声明的方法列表很长，从这个类继承它似乎很乏味。出于这个原因，Palabos提供了一系列继承自Dynamics的类，其中许多方法是为特定目的预定义的。 因此，定义新的动态类本质上包括在继承层次结构中选择正确的位置，然后重写两个或三个方法。 为了便于说明，此继承图的一部分如下图：

Dynamic继承图
Dynamic继承图

Palabos提供了一系列继承自Dynamics的类，其中许多方法是为特定目的预定义的。要学习如何定义新的动力学类，最简单的方法是查看Palabos中定义的一个类。例如，BGK动态定义在文件src / basicDynamics / isoThermalDynamics.hh中。复合动力学（一个修改另一个现有动力学类的行为的类）的一个很好的例子是Smagorinsky动力学，在src / complexDynamics / smagorinskyDynamics.hh中定义。

简单看了看isoThermalDynamics的代码，发现在里面实现了很多动态类的成员方法，比如碰撞，计算宏观变量等。

对于碰撞的实现，动态对象获得对单个单元的引用，因此，碰撞步骤必然是局部的。模拟的非局部成分用数据处理器实现，如下一节所示。

与block-lattice一样，动态对象依赖于两个模板参数，一个浮点参数，一个是网格描述符。通过使用格子描述符中提供的信息，碰撞步骤应以通用的，与网格无关的方式编写。在以通用方式编写算法时显然效率会打折扣，因此通过对给定晶格使用模板特化，可以避免此问题。例如，再次采用BGKdynamics类的实现。通用动态类的碰撞实现通过dynamicsTemplates实现，它在文件src / latticeBoltzmann / dynamicsTemplates.h中定义。然而，针对2D和3D的特殊情况，在文件dynamicsTemplates2D.h和dynamicsTemplates3D.h对碰撞的实现做了重写，这部分代码实现了高效化运行。

3.1 Dynamics的collide方法与BlockLatticeXD的collide方法

值得注意的是，在BlockLatticeXD中有了Collide和Stream这两个成员方法，而在Dynamic里面也有两个这样的成员方法。那么到底执行碰撞的是在BlockLatticeXD里面还是Dynamic里面？

先看看BlockLatticeXD中的collide方法。

BlockLatticeXD有两种形式，看其中一个，后面那个是前面的直接调用，所以看第一个：

1
2
3
4
5
6
7
8
9
10
11
12
13
template<typename T, template<typename U> class Descriptor>
 void ‪BlockLattice2D::collide(‪Box2D domain) {
     ‪PLB_PRECONDITION( (‪plint)Descriptor::q==(‪plint)Descriptor::numPop );
     // Make sure domain is contained within current lattice
     ‪PLB_PRECONDITION( ‪contained(domain, this->getBoundingBox()) );
 
     for (‪plint iX=domain.‪x0; iX<=domain.‪x1; ++ix) {< span>
         for (‪plint iY=domain.‪y0; iY<=domain.‪y1; ++iy) {< span>
             grid[iX][iY].collide(this->getInternalStatistics());
             grid[iX][iY].revert();
         }
     }
 }

其中PLB_PRECONDITION是一个判断语句，括号里面的为真则继续往下执行，为假则报错。

问题的关键在于循环里面的：“grid[iX][iY].collide(this->getInternalStatistics());”这一句。这一句调用了一个collide函数，而grid[iX][iY]是个什么的东西呢？

它的定义是：

1
Cell< T, Descriptor > ** grid

可见，grid是一个二维数组，其里面存储的是Cell的对象。所以上一个语句调用的是Cell类里面的成员函数collide。

那么Cell里面的collide是这样的：

1
2
3
4
5

void ‪collide(‪BlockStatistics& statistics) {
         ‪PLB_PRECONDITION( ‪dynamics );
         ‪dynamics->collide(*this, statistics);
     }

可以看到里面其实调用的是dynamic里面的collide方法。这样就比较清楚了，BlockLatticeXD里面调用的是对整块格子域进行碰撞操作，而dynamic提供的是具体的方法，也就是每个格子进行碰撞的时候都要调用dynamic里面的这个collide方法。

再来看看dynamic里面的collide方法：

1
2
3
4
5
6
7
8
9
10
11
12
13
template<typename T, template<typename U> class Descriptor>
void ‪BGKdynamics::collide (
         ‪Cell& cell,
         ‪BlockStatistics& statistics )
 {
     T rhoBar;
     ‪Array::d> j;
     ‪momentTemplates::get_rhoBar_j(cell, rhoBar, j);
     T uSqr = ‪dynamicsTemplates::bgk_ma2_collision(cell, rhoBar, j, this->getOmega());
     if (cell.‪takesStatistics()) {
         ‪gatherStatistics(statistics, rhoBar, uSqr);
     }
 }

形参里面包括了一个Cell的引用，这说明调用这个方法会改变Cell里面数值。Cell这个类里面最为关键的成员变量就是分布函数f和动态指针dynamics。

第一句对密度进行了声明，第二局声明了一个Array，数组容器。由于“T”代表的是“double”，Descriptor::d是描述子的维度，如果是D2Q9，那么这个Descriptor::d就等于2。所以j就是一个包含两个元素的数组。

第三句调用了一个函数，从函数名可以看到是得到rhoBar_j，看了代码以后，我知道这个东西是代表的是什么，但是我依旧不清楚这个有什么用。

假设采用的是D2Q9描述符，并且其离散速度如下图所示（Palabos里面的D2Q9的速度分布就是这样的）：

D2Q9

其中:

1
2
3
4
5
6
7
8
9
10
lineX_P1  = f[5] + f[6] + f[7];
lineX_0   = f[0] + f[4] + f[8];
lineX_M1  = f[1] + f[2] + f[3];
 
lineY_P1  = f[7] + f[8] + f[1];
lineY_M1  = f[3] + f[4] + f[5];

rhoBar = lineX_P1 + lineX_0 + lineX_M1;
j[0] = (lineX_P1 - lineX_M1);
j[1] = (lineY_P1 - lineY_M1);

j就是这么个东西了，但是我还不知道它用来干啥的。

第四句，也就是最核心的部分，就是执行碰撞了。可以看到，他这里调用dynamicsTemplates类下的bgk_ma2_collision方法。

bgk_ma2_collision碰撞代码如下:

1
2
3
4
5
static T ‪bgk_ma2_collision(‪Cell& cell, T rhoBar, ‪Array::d> const& j, T omega)
 {
     return ‪dynamicsTemplatesImpltypename Descriptor::BaseDescriptor>
         ‪::bgk_ma2_collision(cell.‪getRawPopulations(), rhoBar, j, omega);
 }

其调用了dynamicsTemplatesImpl类下面的bgk_ma2_collision。而这个函数具体为：

1
2
3
4
5
6
7
8
9
10
static T ‪bgk_ma2_collision(‪Array& f, T rhoBar, ‪Array const& j, T omega) {
     T invRho = Descriptor::invRho(rhoBar);
     const T jSqr = ‪VectorTemplateImpl::normSqr(j);
     for (‪plint iPop=0; iPop < Descriptor::q; ++iPop) {
         f[iPop] *= (T)1-omega;
         f[iPop] += omega * ‪dynamicsTemplatesImpl::bgk_ma2_equilibrium (
                                 iPop, rhoBar, invRho, j, jSqr );
     }
     return jSqr*invRho*invRho;
 }

第一句定义了一个invRho，查了一下文档，发现其是一个这样的东西：

1
2
3
static T ‪invRho(T ‪rhoBar) {
         return (T)1 / (‪rhoBar + (T)1);
     }

后面的可以看出来，也不宜过于深究了。其碰撞公式就是：

$f_{k}(x+\Delta x, t+\Delta t)=f_{k}(x, t)[1-\omega]+\omega f_{k}^{\mathrm{eq}}(x, t)$

绕的有点多了。

4 数据处理器

我认为这个部分是比较复杂的，原因是我在看后面的边界条件的时候，根据不同的边界条件大量的应用了数据处理器。单看用户文档还是无法完全了解。目前先按照用户文档的内容记录。

以下主要来自文档翻译。

数据处理器作用在整个求解域，或者是某个块（block）。在block-lattice中，他们用来实现所有无法通过动态类（Dynamic)实现的数据处理操作。这些操作主要是那些非局部操作，比如用于实现边界条件的有限差分模板的评估。在标量场和张量场中，数据处理器提供了在空间扩展域上执行操作的唯一（足够有效且可并行化）的方式。

此外，数据处理器能够在几个相同类型或不同类型的块之间交换信息（例如：块格和标量场之间的耦合）。这通常用于实现物理耦合（多组分流体，具有Boussinesq近似的热流体），初始条件的设置（从矢量场中的值初始化速度），或执行经典基于数组的操作（以元素方式添加两个标量字段）。

最后，数据处理器用于在整个块或子域上执行缩减操作。示例的范围从计算平均动能到计算作用在障碍物上的阻力。

数据处理器的定义和应用采用了两种不同的观点。在应用程序级别，用户指定给定块的区域（可以是矩形或不规则的），在该区域上执行数据处理器。如果块内部具有多块结构，则数据处理器被细分为几个更具体的数据处理器，各自负责atomic-block的数据处理。因此，在执行级别，数据处理器总是作用于atomic-block，即作用于先前通过将原始区域与原子块的域相交来确定的区域。

最后，定义新的数据处理器非常容易。您需要做的就是编写一个函数，该函数接收atomic-block和sub-domain的坐标作为参数，并在该子域上执行算法。所有复杂的操作，例如存在多块的操作的细分，或代码的并行化，都是自动的。

目前看起来比较头大，有很多东西现阶段不是很理解。不过好在用户文档里面也提到了：Palabos通过所谓的数据处理功能提供简化的界面，隐藏技术细节并让您专注于基本部分。也就是说知道怎么用就行了。看看，可能传递的就是一个atomi-block和一个sub-domain的坐标。然后指定一些palabos已经替我们处理好的算法。

文档提到，可以在文件\src\dataProcessors\ latticeInitializerXD.h和.hh中找到更多关于块格的数据处理器的定义，并且在文件\src\dataProcessors\dataFieldInitializerXD中定义作用于标量或张量字段的数据处理器。减少操作的评估示例在文件\src\dataProcessors\dataAnalysisXD.h和.hh中提供。

实际文档中给出的地址路径不对，且latticeInitializerXD查找不到。究其原因，就是因为palabos官方的文档是1.1的，而我下下来的代码是2.0的。

因此我才用doxygen生成了2.0的开发文档。不得不说用起来非常的方便，查找功能很好。再次安利Palabos的2.0版开发文档，如果能对研究Palabos的伙伴产生作用，那是再好不过了。

5 吐个槽

暂时写到这里，有点晚了，本来以为一两个小时可以搞定，没想到从下午写到了晚上。主要就是在研究collide函数花了比较多的时间。

如果只是想要初步的了解，不建议挖的过深。

我就是有时候挖到底层了，让后让自己迷失在代码里，抓不住重点。这是阅读代码所不应该的地方。

以后还是要注意一下，避免陷入代码海洋，而只见树木不见森林。