格子Boltzmann数值计算的性能优化

Performance Optimization of Lattice Boltzmann Numerical Calculation

摘要: 格子Boltzmann方法是科学与工程计算的重要模型与数值方法，该方法能够从底层描述和计算细尺度特征，同时计算量也较大。为了提高计算效率，本文主要从指令级别对格子Boltzmann数值计算进行性能优化。数值实验结果表明，优化后计算性能显著提高。

Abstract: Lattice Boltzmann method is an important model and numerical method in scientific and engi-neering calculation. It can describe and calculate the fine scale characteristics from the bottom layer, and the calculation cost is also large. In order to improve the computational efficiency, this paper mainly optimizes the performance of lattice Boltzmann numerical calculation from the in-struction level. The results of numerical experiments show that the computational performance is improved significantly after optimization.

1. 引言

随着数值计算方法成为研究科学与工程问题的重要手段，对计算效率的要求也越来越高。为了确保在尽可能短的时间内完成一项计算任务，在单机中我们主要有两个选择，第一个选择是设计一个更好的算法，可以用更少的步骤达到相同的结果。第二种方法是在指令级别执行优化，该方法是在需要较多时间完成的指令的位置上使用一个或多个更快的指令。

格子Boltzmann方法是模拟微观模型的数值计算方法。格子Boltzmann方法除了被应用于一般的流体力学问题之外，还在湍流 [1]、多孔介质流 [2] [3]、粒子悬浮流 [4]、磁流体力学 [5]、多相流 [6] 等相关领域也取得了比较成功的应用。但格子Boltzmann方法的计算量较大，为提高计算效率，本文将探讨泰勒–格林涡流(Taylor-Green Vortex)格子Boltzmann数值计算的性能优化。

2. 格子Boltzmann模型

格子Boltzmann方法的一个主要优点是基于Boltzmann方程而不是连续方程和动量方程，它的实现比传统方法简单。格子Boltzmann方法的基本量是离散速度分布函数 $f_i\left(x,t\right)$，通常称为粒子分布。可以通过此分布函数来求解质量密度和动量密度：

$\rho \left(x,t\right)={\displaystyle \underset{i}{\sum }{f}_i\left(x,t\right)}$ (1)

$\rho u\left(x,t\right)={\displaystyle \underset{i}{\sum }{c}_i{f}_i\left(x,t\right)}$ (2)

$f_i$ 函数表示所有的参数变量都是离散的，下标i代表一个离散速度集，定义 $f_i$ 在空间中为正方形格子各个方向的分布。速度集通常用DdQq表示，d是速度集覆盖的空间维数，q是速度方向的个数。最常用的速度集是D1Q3、D2Q9、D3Q15、D3Q19和D3Q27。本文采用D2Q9模型，i取值范围是0，1，2，3，4，5，6，7，8九个方向。

通过在物理空间、速度空间和时间上离散Boltzmann方程，得到了格子Boltzmann方程：

$f_i\left(x+c_i\Delta t,t+\Delta t\right)=f_i\left(x,t\right)+{\Omega }_i\left(x,t\right)$ (3)

同时，粒子受到碰撞算符 ${\Omega }_i\left(x,t\right)$ 的影响，该算符通过在每个位置的分布中重新分配粒子来模拟粒子碰撞。虽然有许多不同的碰撞算符可用，但最常用的一个算符是Bhatnagar-Gross-Krook (BGK)运算符 [7]：

${\Omega }_i\left(f\right)=-\frac{f_i-f_i^{eq}}{\tau }\cdot \Delta t$ (4)

它使分布 $f_i$ 以由松弛时间 $\tau$ 决定的速率趋于平衡分布 $f_i^{eq}$，平衡分布为

$f_i^{eq}\left(x,t\right)=w_i\rho \left(1+\frac{u\cdot c_i}{c_s^2}+\frac{{\left(u\cdot c_i\right)}^2}{2c_s^4}-\frac{u\cdot u}{2c_s^2}\right)$ (5)

格子BGK (LBGK)方程(即用BGK碰撞算符完全离散的Boltzmann方程)可以表示为：

$f_i\left(x+c_i\Delta t,t+\Delta t\right)=f_i\left(x,t\right)-\frac{\Delta t}{\tau }\left(f_i\left(x,t\right)-f_i^{eq}\left(x,t\right)\right)$ (6)

总的来说，格子Boltzmann方程由碰撞和迁移两部分组成。第一部分是碰撞(或松弛)：

$f_i^{\ast }\left(x,t\right)=f_i\left(x,t\right)-\frac{\Delta t}{\tau }\left(f_i\left(x,t\right)-f_i^{eq}\left(x,t\right)\right)$ (7)

第二部分是迁移(或传播)：

$f_i\left(x+c_i\Delta t,t+\Delta t\right)=f_i^{\ast }\left(x,t\right)$ (8)

碰撞只是一个代数局部运算，首先要计算密度，找到平衡分布后的宏观速度，碰撞后，将得到相邻节点的分布 $f_i^{\ast }$，当碰撞和迁移两个操作完成时，经过一个时间步，再重复这些操作。

格子Boltzmann数值计算的实现过程可分为碰撞和迁移两个子过程，通常格子Boltzmann方法的程序结构有两种形式：碰撞–迁移结构和迁移–碰撞结构。本文程序设计采用迁移–碰撞结构，这种程序结构可以看作是求解离散速度方程

${\partial }_t{f}_i+c_i\cdot \nabla f_i={\Omega }_i\left(f\right)$ (9)

的时间分裂方法 [8]，具体步骤如下：

(I) 初始化分布函数 $f_i\left(x,0\right)$

(II) 执行迁移(stream)

$f_i\left(x+c_i\Delta t,t+\Delta t\right)=f_i^{\ast }\left(x,t\right)$ (10)

(III) 计算宏观量(computeRhoU)

$\rho \left(x,t\right)={\displaystyle \underset{i}{\sum }{f}_i\left(x,t\right)}$， $\rho u\left(x,t\right)={\displaystyle \underset{i}{\sum }{c}_i{f}_i\left(x,t\right)}$ (11)

(IV) 执行碰撞(collide)

$f_i^{\ast }\left(x,t\right)=f_i\left(x,t\right)-\frac{\Delta t}{\tau }\left(f_i\left(x,t\right)-f_{ieq}\left(x,t\right)\right)$ (12)

(V) 保存数据(save) (每saveN步保存一次)

(VI) 重复(II)~(V)直到满足终止条件。

3. 计算与性能优化

本文计算泰勒–格林涡流问题，速度和压力的初始状态解析地设定，泰勒–格林涡流在 ${\mathcal{l}}_x\times {\mathcal{l}}_y$ 区域内是非定常的全周期流动，其速度场和压力场在二维空间上表示为

$u\left(x,t\right)=u_0\left[\begin{array}{c}-\sqrt{k_y/k_x}\cos \left(k_{x}x\right)\sin \left(k_{y}y\right)\\ \sqrt{k_x/k_y}\sin \left(k_{x}x\right)\cos \left(k_{y}y\right)\end{array}\right]{\text{e}}^{-t/t_d}$ (13)

$t_d=\frac{1}{\upsilon \left(k_x^2+k_y^2\right)}$ (14)

$p\left(x,t\right)=p_0-\rho \frac{u_0^2}{4}\left(\frac{k_y}{k_x}\cos \left(2k_{x}x\right)+\frac{k_x}{k_y}\cos \left(2k_{y}y\right)\right)$ (15)

其中 $u_0$ 为初始速度大小， $k_{x,y}=2\pi /{\mathcal{l}}_{x,y}$ 为波矢量分量， $t_d$ 为涡流衰减时间，平均压力 $p_0$ 可以是任意的，计算区域为 ${\mathcal{l}}_x\times {\mathcal{l}}_y=640\times 640$ (格子单位)。计算采用标准平衡分布，平衡分布初始化采用Mei [9] 的初始化方案。泰勒–格林涡流的速度 $u_x$ 和 $u_y$ 如图1所示。

Figure 1. Taylor-green vortex velocity: (a) x directional velocity $u_x$, (b) y directional velocity $u_y$

图1. 泰勒–格林涡流速度：(a) x方向速度 $u_x$，(b) y方向速度 $u_y$

上面的计算是根据格子Boltzmann方法典型计算步骤(I)~(VI)进行设计程序的，下面对该程序进行一些优化处理。我们首先注意到一个相对较小的优化：迁移步不需要复制0方向分布f₀。为了避免这种不必要的内存访问，我们对0方向分布使用单独的变量f₀，对其余八个方向的分布f_1-8保留两个变量f₁和f₂。

第二个优化是注意到在迁移期间为每个节点存储的分布的值是在碰撞步期间读取的那些值。因此，四个函数stream、computeRhoU、collide和save可以组合成一个函数stream_computeRhoU_collide_save，该函数访问内存的频率显著降低。

第三个优化是循环展开。当循环的内部块计算很快完成时，递增计数器并检查其是否超过界限的开销构成了循环执行时间的一个显著部分。循环展开可以避免这些开销，这里将九个分布上迭代的所有循环展开，其展开公式如下：

$\rho =f_0+f_1+f_2+f_3+f_4+f_5+f_6+f_7+f_8$ (16)

$u_x=\left[\left(f_1+f_5+f_8\right)-\left(f_3+f_6+f_7\right)\right]/\rho$ (17)

$u_y=\left[\left(f_2+f_5+f_6\right)-\left(f_4+f_7+f_8\right)\right]/\rho$ (18)

$f_0^{eq}=\frac{2\rho }{9}\left(2-3u^2\right)$ (19)

$f_1^{eq}=\frac{\rho }{18}\left(2+6u_x+9u_x^2-3u^2\right)$ (20)

$f_2^{eq}=\frac{\rho }{18}\left(2+6u_y+9u_y^2-3u^2\right)$ (21)

$f_3^{eq}=\frac{\rho }{18}\left(2-6u_x+9u_x^2-3u^2\right)$ (22)

$f_4^{eq}=\frac{\rho }{18}\left(2-6u_y+9u_y^2-3u^2\right)$ (23)

$f_5^{eq}=\frac{\rho }{36}\left[1+3\left(u_x+u_y\right)+9u_x{u}_y+3u^2\right]$ (24)

$f_6^{eq}=\frac{\rho }{36}\left[1-3\left(u_x-u_y\right)-9u_x{u}_y+3u^2\right]$ (25)

$f_7^{eq}=\frac{\rho }{36}\left[1-3\left(u_x+u_y\right)+9u_x{u}_y+3u^2\right]$ (26)

$f_8^{eq}=\frac{\rho }{36}\left[1+3\left(u_x-u_y\right)-9u_x{u}_y+3u^2\right]$ (27)

这里 $u={\left(u_x+u_y\right)}^{{\rm T}}$， $u^2=u_x^2+u_y^2$。

还有其它的一些优化处理：将频繁调用的节点索引函数定义为内联函数，省去了函数调用的开销，从而提高函数的执行效率。另外，程序中包含一个布尔参数saveN，用于指示是否将该时刻的数据写入内存，通过选择saveN最优的值，定期保存中间密度场和速度场的数据，避免过多的内存写入。

表1是程序在未进行优化和优化后的执行的时间和速率。从表中可以看出，程序经过优化处理后，计算规模为64 × 64、160 × 160、160 × 160和1920 × 1920的计算时间分别减少到未优化的22.98%、22.21%、20.81%和20.74%，计算速率分别是未优化的4.35倍、4.50倍、4.80倍和4.82倍。

Table 1. Comparison of implementation results between unoptimized and optimized

表1. 优化前后执行结果对比

4. 总结

本文以泰勒–格林涡流作为算例给出了格子Boltzmann数值模拟的算法和计算实现，以及讨论了如何优化代码提高计算效率。计算结果表明，程序通过优化处理后，计算时间明显减少，计算速率明显提高，而且计算规模越大，这种优势越显著。

基金项目

广东省教育厅青年创新人才项目(2017GkQNCX110)；国家自然科学基金(11804355, 31800083)。

NOTES

^*通讯作者。