基于CUDA的火焰模拟

满家巨，邹有，陈传淼

(湖南师范大学数学与计算机科学学院高性能计算与随机信息处理省部共建教育部重点实验室, 中国 长沙 410081)

CUDA(General Unified Device Architecture)是NVIDIA公司在2007年正式发布的GPU通用计算开发环境和软件体系，它不但为用户提供了一系列操作GPU数据读写和命令流控制的API，让用户应用GPU进行通用计算(GPGPU)，还提供了OpenGL和Direct 3D进行互操作的接口，使计算之后的结果可以快速地绘制出来．随着廉价图形卡的处理能力日渐强大，许多行业已经开始使用CUDA进行加速计算，如音视频处理、石油勘探、医学成像、科学计算、物理模拟等．

在计算机的计算能力有限时，研究者们只能通过参数建模的方法来进行流体模拟，如将波浪函数的线性波组合来模拟浪花[1]、使用拉格朗日粒子模拟波浪[2]等，但这些方法不能模拟比较复杂的、细节更为丰富的流体效果，于是研究者们开始把注意力转向基于物理的流体模型．流体模拟使用的基于物理的方法主要有：晶格波尔兹曼方法(LBM)[3]、拉格朗日方法[4]以及基于N-S方程的欧拉方法[5]．其中欧拉方法是最常使用的一种计算方法，能够非常容易地重构光滑的流体表面，借助文献[6]提出的半拉格朗日法进行求解，可以选取非常大的时间步长并能保持计算的稳定，但缺点在于计算量比较大，需要非常高效的方程组求解器才能实现实时模拟．

本文以欧拉方法进行火焰模拟为线索，对模型的计算复杂度进行分析，实现了基于CUDA的线性方程组的求解器；针对火焰模拟中的颜色、漩涡、运动形态等内容，引入漩涡约束因子和随机风力来制造火焰的摇曳效果，并提出一种简便的烟雾和焰心效果绘制策略来增加模拟的真实性，最终计算结果通过基于CUDA的光线投射方法实现可视化．

1 CUDA编程模型与图形API

相较于CPU而言，CUDA主要结构为计算逻辑单元而不是控制逻辑单元，因此其浮点运算能力可达到CPU的几十倍甚至上百倍．CUDA将对GPU硬件的操作予以屏蔽，让用户可以使用简单的C语言接口实现通用计算．

1.1 CUDA编程模式

一个CUDA程序由主机端程序和设备端程序组成.主机端程序是在CPU上执行的串行代码，而设备端程序是在GPU上执行的并行代码，即kernel程序．使用CUDA进行编程时，首先准备好主机端数据，然后将其从主机的内存复制到设备内存中，调用kernel程序进行高密度的并行计算.待计算完成再将计算结果从设备端复制回主机内存．如图1，在CUDA架构中，线程被分成许多彼此独立线程块(Block)，同一个线程块内的所有线程使用共享内存协调计算任务，一个GPU上的所有线程块组织成网格(Grid)．

图1 线程块与线程的组织形式[7] 图2 CUDA存储模型[7]

1.2 CUDA存储模式

CUDA为满足数据的不同需求提供了6种存储器(图2)：寄存器(Registers)、本地存储器(Local Memory)、共享存储器(Shared Memory)、全局存储器(Global Memory)、常数存储器(Constant Memory)和纹理存储器(Texture Memory)．其各自特点如下：

(1)寄存器的访问效率最高，每个线程都拥有专属的寄存器，但它的数量比较少；

(2)本地存储器访问速度较慢，也是线程特有的存储单元，主要在寄存器空间不足时使用；

(3)共享存储器访问速度快，是线程块所有线程公有的存储单元，其大小也比较有限.常见的硬件每个线程块可使用的共享存储器大小只有16 KB．

(4)全局存储器的大小与显卡的显存大小相当，它可通过主机端程序和设备端程序访问，充当主机端和设备端的数据交互媒介．

(5)常数存储器是显存中的只读存储单元，它拥有缓存加速的功能，可以节约带宽而加快访问速度，因此适合存储程序中需要频繁访问的只读参数．

(6)纹理存储器也是只读存储器，它由GPU用于纹理渲染的图形专用单元发展而来，具有地址映射、数据滤波、缓存等功能，适合实现图像处理、查找表或大量数据随机访问等操作．

1.3 CUDA与图形API互操作

CUDA为图形开发用户提供了与OpenGL和Direct 3D这2种主流的图形API互操作接口．一些OpenGL和Direct 3D的资源可以被映射到CUDA的地址空间，使得计算数据可以通过OpenGL或Direct 3D来显示．本文中使用CUDA与OpenGL的互操作API来实现可视化，将OpenGL缓冲对象映射到CUDA的地址空间，这样就可以在CUDA中读取OpenGL写入的数据，也可以用CUDA写入数据待OpenGL使用．

2 基本方程及求解

2.1 基本方程

火焰的燃烧模型通常被当成不可压缩流来处理[8]，其速度场所对应的动量方程和连续性方程(Navier-Stoke方程)为(1)和(2)所示，式(3)为火焰的温度场方程．

(1)

(2)

(3)

根据Helmholtz-Hodge分解定理对基本方程进行投影操作，可将方程分解为多个步骤分别求解[9]：

(4)

(5)

(6)

(7)

(8)

2.2 计算复杂度分析

根据半拉格朗日方法，式(4)～(8)的求解过程可分为4步：添加源、对流、扩散和投影．使用边长为N个划分的立方体网格进行计算，总的节点数为N3．

(1)添加源操作中，需要将速度场和温度场进行初始化和用户交互的输入．给定一个速度场(三维)和温度场，需对每个节点的速度和温度进行赋值操作，共需4N3次乘-加运算．

(2)扩散步中，速度场三维数据和温度场一维数据需全部更新，而每节点数据都由其周边6个位置的数据插值所得，假设最大迭代次数为MAX_ITER，则扩散过程共需要MAX_ITER×4×6N3次浮点运算．

(3)投影由散度计算、线性方程求解和压力修正3部分组成．首先对每一个网格单元进行速度场的散度计算，散度的计算只针对速度场进行，计算涉及每单元的周边6个数值的乘-加操作；线性方程组求解的计算方法和复杂度与扩散步相同；压力修正即对速度场的3个分量使用对应的2个压力值进行调整，各需2次乘加操作．整个投影步需要进行的浮点操作数为(57+MAX_ITER×6)N3次．

(4)对流步计算也是速度场和温度场共同参与的计算过程，经插值和越界判断，共需45N3次浮点运算．

综上所述，在使用稳定流方法进行三维的火焰模型计算中，总共需要进行的浮点运算为(106+30MAX_ITER)N3次．

2.3 基于CUDA的线性方程组求解器

由上一小节可知，模型的求解过程中涉及大量线性方程组的求解，因此加速线性方程组的求解是提高模拟速度的关键．分析方程可知，方程组系数矩阵均为大规模的稀疏矩阵，其非零元呈对角分布且主对角占优，因而非常适合迭代法求解，并可采用一维数组来存储以达到高效访存．

根据CUDA编程模型，将三维数据按高度所在方向分层分配成线程块，并使用循环进行遍历，每一层的网格对应的值可以用线程的索引来获取：

确定对数值的访问模式后，就可以开始设计迭代法求解了．本文实现了3种比较常见的迭代法：Jacobi迭代、Gauss-Seidel和共轭梯度法(CG)，并分别进行求解实验，表1～3分别展示了3种迭代法在不同问题规模下相对于纯CPU运算的加速情况．本文实验使用的配置为：CPU：Intel 四核2.34 Hz；内存：4.0 GB；显卡型号：Geforce GTX 285；显存大小：2 GB；操作系统：Windows 7．

表1 Jacobi迭代(迭代20次)

表2 Gauss-Seidel迭代(迭代20次)

表3 共轭梯度法(迭代20次)

通过以上3个表的数据可以看出，在网格规模较小时，GPU的加速比相对较小，并行加速的优势不能完全发挥，这主要是因为此时线程数目较少，GPU的计算和存储资源不能被充分利用．随着网格数目的增大，加速比开始呈上升趋势．在迭代同样多的次数时，Jacobi迭代的计算速度更胜一筹．

3 燃烧效果控制

3.1 漩涡限制因子

基于N-S方程模拟火焰时，通常将其看作是低粘性的流体，其运动会产生一定的旋涡现象．然而，采用半拉格朗日方法进行计算时会引起一定程度的数值耗散[8]，使得漩涡效果减弱．为了使模拟效果更加逼真，本文使用文献[10]提出的漩涡限制因子方法来恢复这种细节，通过捕捉速度场中每一个计算单元上的旋转运动，计算出加强或维持这种旋转所需要的粘滞力，以此来减少由数值耗散引起的细节丢失．图3展示了加入不同大小人工漩涡得到的效果.可以看出，加入人工旋涡后，模拟的结果表现出更多的形态变化效果，从而更富有表现力．

(a)无人工漩涡 (b)少量人工漩涡 (c)较多人工漩涡图3 人工漩涡效果

3.2 颜色控制与焰心绘制

真实情况下的火焰是一个非常复杂的系统，火焰的颜色也受到许多因素的影响，如燃料的种类，燃烧的程度，燃烧时的气压等，每一个条件的变化都可能引起火焰颜色的变化．本文从火焰颜色的物理模型出发，使用黑体辐射表来确定火焰的颜色，黑体辐射可以通过文献[11]中的工具进行计算．

在通常的模拟中，火焰燃烧产生的烟雾场和火焰是作为2个独立的模型来处理的，这种方法可以很方便地控制烟雾外观和数量，然而计算量也大大增加．本文采用一种简便方法，把烟雾的颜色放到火焰的颜色场之中处理，将烟雾看作是火焰的外延．其基本过程描述如下：

(1)给定温度范围[T0,T1]，计算出对应的颜色表；

(2)根据燃料的性质定义一种烟雾的颜色，如木材燃烧一般冒黑烟，可以把烟雾定义成黑色．把定义的烟雾颜色加入颜色表尾部．

(3)对于计算出的温度场中所有温度小于T0的点，都使用烟雾的颜色进行绘制．

图5展示了采用这一方案进行绘制的烟雾效果.

图4 将烟雾看成是主焰的延伸

(a)无烟雾效果 (b)增加烟雾效果图5 烟雾效果对比

对于焰心的绘制，本文摒弃了理论和计算都颇为复杂的燃烧模型，直接利用OpenGL进行图形绘制时的Alpha通道，使处于火焰中心位置的人造焰心与火焰重叠出比较逼真的焰心效果．如图6所示，对于一个给定的三维火焰效果，从其中心挖去一个平截头体作为人工焰心，对挖掉的部分填充焰心颜色．在绘制启用Alpha通道后，视点位置就可以透过外层火焰看到内部的焰心，内部的焰心颜色会和外层不断变化的火焰颜色叠加在一起，形成逼真的燃心效果．更进一步，对挖掉的平截头体的大小和形状增加一定的随机性，就可以模拟出随机跳动的焰心效果了．图7展示了使用本文方法加入了焰心的绘制效果．

图6 焰心绘制示意图

(a)无焰心效果 (b)加入焰心效果图7 焰心绘制效果

4 模拟效果可视化

4.1 体绘制

体绘制是一个从三维数据集生成其在二维平面的投影过程，实现体绘制最常见的方法为光线投射(Ray Marching)．光线投射基础步骤为：

图8 绘制流程

(1)投射：视点出发向所有的像素发射一条射线；

(2)采样：沿着射线方向，使用相同的间隔选取采样点；

(3)插值：对每一个采样点，使用其周围的点(通常有8个点)进行颜色值和Alpha值的插值运算，求出该采样点的颜色值和Alpha值；

(4)累加：把在射线上的所有采样点的颜色值和Alpha值进行由后向前的合成，即可获得每条光线对应的像素点处的颜色值和Alpha值．

显然，要通过光线投射获得最终的画面，需要对窗口上的每一个像素点都进行一次计算和颜色累加，对于规模较大的模拟，其计算量相当大．

4.2 基于CUDA的光线投射算法

光线投射的计算量虽然庞大，但任意2条光线的投射过程没有任何关联，非常适合使用GPU来并行计算．使用CUDA实现光线投射大致包括绑定像素缓存、光线投射计算、解决缓存绑定、绘制结果4个步骤[12]，其中最重要的部分为光线投射计算．其基本实现思路如下：

(1)设定绘制窗口大小，尽量使用16的整数倍；

(2)将绘制窗口中的像素映射到16×16的线程块中，每个线程块处理256个像素点；

(3)计算每个像素点在世界坐标中的位置，计算从视点指向像素点的射线方向；

(4)求解每条射线与计算空间的6个面的交点；

(5)在最远的交点和最近的交点之间以某固定距离采样，计算采样点对应的颜色值；

(6)将结果写入像素缓存中用于最终的显示．

4.3 可视化效果展示

到此为止，本文已经完成了一整套火焰模拟的计算和绘制方案．通过图9，可以明显地看出原始模拟在依次加入人工漩涡、烟雾和焰心绘制后得到的不同效果．图10展示了加入场景绘制的效果图，其可视化效果基本可以满足一般游戏或普通动画要求．

(a)原始效果 (b)加入漩涡 (c)绘制烟雾 (d)绘制焰心图9 本文实现的效果

图10 燃烧的火把和蜡烛

5 结语

本文基于CUDA实现了火焰模拟的计算、形状颜色控制以及可视化全过程，工作包括线性方程组求解器的实现、燃烧效果控制和可视化3部分，针对方程组求解和可视化，实现了基于CUDA的3种线性求解器和光线投射算法，针对燃烧效果的控制提出了简化的烟雾和焰心绘制方法，总体达到了比较可观的模拟结果．由于本文的立足点并不是模型改进或算法改良方面的工作，因此后续的工作还有很多，比如网格改进、模拟对象改变等都是比较好的研究方向．

参考文献：

[1] DARWYN R, PEACHEY. Modeling waves and surf [J]. Computer Graphics, 1986,20(4): 65-74.

[2] ALAIN F, WILLIAM T. A simple model of ocean waves [J]. Computer Graphics,1986,20(4):75-84.

[3] THÜREY N, RÜDE U. Free surface lattice-Boltzmann fluid simulations with and without level sets[C]//Proceedings of Workshop on Vision, Modeling, and Visualization (Stanford, CA). Amsterdam:IOS Press, 2004,199-208.

[4] REEVES W T. Particle systems: a technique for modeling a class of fuzzy objects[C]//Proceeding of the 10th annual conference on computer graphics and interactive techniques. New York, NY: ACM, USA, 1983.

[5] SCHLATTER B. A pedagogical tool using smoothed particle hydrodynamics to model fluid flow past a system of cylinders[R]. Oregon:Oregon State University, 1999.

[6] STAM J. Stable fluids [C]//Proceedings of SIGGRAPH. Los Angeles, CA:ACM Press, 1999,121-128.

[7] NVIDIA C. Programming Guid: version 2.1[C]//Proceeding of Symposium on Interactive Ray Tracing. Los Angeles, USA: IEEE Press, 2008:185-190.

[8] 李建明，吴云龙，迟忠先，等.基于流体模型和GPU 加速的火焰实时仿真[J].系统仿真学报, 2007,19(19):4382-4385.

[9] BRIDSON R. Fluid simulation for computer graphics [M].Wellesley, MA：A K Peters, Ltd, 2008.

[10] FEDKIW R, STAM J. Visual simulation of smoke [C]//Proceedings of SIGGRAPH 2001, New York: ACM, 2011:15-22.

[11] Black body spectrum[EB/OL].http://www.mi.infm.it/manini/dida/BlackBody.html.

[12] MARSALEK L, HAUBER A. High-speed volume ray casting with CUDA[C]//Interactive Ray Tracing, 9-10 Aug., 2008, Los Angeles, CA. IEEE Symposium on, IEEE: IEEE Press, 2008, 185.