Blaze：GPU加速下的并行流体模拟与高效渲染

Blaze 是一款完全在 GPU 上运行的流体解算器和渲染器。该框架通过仅在Device上分配Grid信息来避免主机和设备之间昂贵的内存传输，这在实践中可带来显著的性能提升。Blaze 可以使用极简用户界面读取生成的场景描述，从而实现简单而强大的流体工作流程。该项目是Alexandre Sirois-Vigneu在麦吉尔大学COMP-559的最终项目，前后历时6周完成，项目地址：

一，内容介绍

流体模拟和渲染在VFX领域通常是两个独立的连续过程，其中模拟的内容需要存储在磁盘上，以便稍后由另一个软件进行渲染。这些大内存传输对如何使用图形硬件加速这些应用程序施加了严重的限制。将体积数据存储在磁盘上也会对网络带宽造成明显的压力，并且需要大量存储解决方案，因为其中一些模拟每个版本很容易达到数十TB。同样值得注意的是，这些模拟通常使用图形管道进行可视化，而图形管道并不能准确地表示最终的光线行进图像。这种差异通常是问题的根源，因为流体中的一些伪影可能只有使用高质量渲染引擎才能看到，而这通常发生在整个模拟之后。本文提供的方法跳过了这些中间步骤，直接进入最终的高质量渲染，步骤最少。这样程序的输出就可以始终值得信赖，而不会有任何保真度损失。只有粒子文件和渲染的帧需要存储在磁盘上，与存储体积网格的要求相比，这仅需要少量的存储硬件。

二，与课题有关的工作

1999 年，Stam 发表了一篇关于稳定流体的论文，首次提出了一种无条件稳定模型来模拟复杂的流体流动。这篇论文对计算机图形学的计算流体动力学 (CFD) 做出了重大贡献，并在之后的几年中引发了多篇后续论文的发表 [Stam 1999]。2008 年，Bridson 出版了他的第一版著作《计算机图形学的流体模拟》，该书于 2015 年修订并出版了第二版，涵盖了流体模拟的所有方面，从数学和算法到实现。这将成为许多有兴趣进入 CFD 领域的人的参考 [Bridson 2015]。2009 年，Gomes 等人研究了将非发散投影的稀疏线性系统的求解移植到 GPU 可能带来的加速 [Gomes 2009]他们对 CPU 和 GPU 实现的比较揭示了使用可编程 GPU 进行流体模拟的好处和挑战。

三，方法讨论

解算器和渲染器的实现完全是在GPU上使用CUDA Toolkit 10.2完成的，想法是将整个流体模拟和渲染管道移到显卡上以避免在主机和设备内存之间传输网格信息的昂贵成本。总所周知，这些内存传输是大多数GPU加速应用程序中瓶颈的主要来源。在本文所讨论的框架中，只有源粒子才会被推送到设备，同样，只有包含渲染像素的帧缓冲区才会从 GPU 内存中检索。这种设计选择使我们能够最大限度地提高GPU的占用率，并使其在几乎任何执行点都保持完全饱和。程序受两个开源库的启发，分别是：

[1] Bitterli 2021 的 incremental-fluids:

[2] Allen 2021 的 raytracinginoneweekendincude：

程序执行过程：

四，流体模拟

不再求解N-S方程，而是求解无粘性流体的欧拉方程，方法是从方程中剔除粘度项：

将复杂方程分解为多个较小部分，然后逐个进行数值求解，求解的三个主要步骤依次是：

外力项:

压力项:

平流项:

高效内存使用：为了最大限度地减少GPU内存占用，我们选择将密度和温度模拟为单个字段。温度值的范围在渲染时重新映射，使用全范围作为烟雾密度，使用光谱的高端来模拟火焰发出的光。这确实对我们的求解器可以建模的现象类型施加了一些限制，但也允许它处理更高的网格分辨率，下图详细显示了运行时GPU内存中所需要的主要数据结构：

MAC网格：中心差分，交错网格采样

无散保证：为了保证无散速度场，需要使用以下公式来更新速度：

这导致以下泊松问题：

存在多种方法可以求解此系统，例如“雅可比，高斯-赛德尔或共轭梯度”等。本程序实现了“红黑高斯-赛德尔求解器”，它是原型的一种优化，要由于雅可比法和共轭梯度。与雅可比算法相比，该算法的一个重要优势是它可以就地执行，这可以降低内存消耗，同时利用先前计算的值来加快收敛速度。该技术背后的想法是概念性地在整个网络上应用红黑棋盘格图案，这样相邻地体素就不会共享相同的颜色：

压力体素的高斯-赛德尔更新将涉及自身及其所有用不同颜色标记的直接相邻体素，这意味着在前半个迭代过程中，这些体素是只读的。后半个迭代将更新另一种颜色，同时保持先前更新的颜色不变。这个两步过程消除了高斯-赛德尔算法并行实现所固有的写入风险，从而确保了确定性结果。实际上，颜色变化是在域的每个 2D 切片上按顺序执行的，以尽可能多地利用预先计算的值。

网格平流与插值：速度和温度通过半拉格朗日[Stam 2003]平流和三阶精确龙格-库塔方法[Ralston 1962]在网格中传输。

这样即使应用较大的时间步长，也能达到良好的精度。对中等速度样本使用三线性插值（trilinearinterpolation），对流体量的实际平流使用三次插值（Tricubic interpolation）。三次插值比三线性插值成本高得多，但精度也高两个数量级，由于程序中大量使用三次和三线性内核，因此对其进行了优化以减少对局部变量的使用。这不可避免地会导致代码可读性降低，但事实证明这是一种将这些函数调用速度提高两倍以上的有效方法。三次内核的原始实现有如此多的局部变量声明，以至于它会耗尽 GPU 上所有可用的寄存器，从而限制可以并行启动的线程数量。

粘度数值误差：如前所述，我们选择从N-S方程中删除粘度项。这背后的想法有两个方面。首先，我们希望确保在网格分辨率的限制内捕捉尽可能多的视觉细节。其次，大多数用于模拟流体的数值方法也会不可避免地引入不明确建模粘度来减少每个时间步所需的计算量。

通过替换将粒子转换为网格：在每个时间步开始时，粒子通过GPU内核散布到网格上（P2G过程），其中每个线程散布单个粒子，使用在网格上执行读取-修改-写入原子操作的原子函数解决写入风险。这确保没有线程会覆盖另一个同时运行的线程刚刚设置的值。众所周知，使用这些操作会导致严重的性能下降，但可以通过直接最小化粒子文件中的粒子重叠来减少原子函数的影响，这种优化可以在Houdini中轻松完成，如下图：

执行这一简单步骤可使性能与完全避免使用原子函数相媲美，因为线程很少遇到写入冲突。与常见实现不同，我们的源是通过替换而不是添加来完成的，这在控制方面具有多个优势。首先，更容易模拟突然的变化，例如剧烈的速度变化，而不必注入过大的力来抵消模拟中已经存在的动量。随着模拟的发展，网格中的这种不连续性也往往会产生有趣的细节。其次，网格内的值保持受源的限制，这允许用户根据网格内值的范围轻松定义着色器的ramps以及其他参数。

对于P2G过程分为三个步骤。首先，Clear掉所有网格的back buffer；然后，使用原子函数将粒子散射到Back Buffer; 最后，将back buffer中超过1 x 10^-5阈值的值覆盖到front buffer。此操作比add在计算上稍微昂贵一些，但与求解器的其他步骤相比可以忽略不计。粒子散射过程使用Cubic Pulse Function（立方脉冲函数）[Quilez 2021]，该函数从粒子表面的零平滑过渡到其中心的1。需要注意的一点是，由于每个线程处理单个粒子，因此在散射大粒子时可能会导致该步骤非常慢。因此，为了利用硬件的并行特性，使用数千个球形小粒子而不是一个大粒子来填充相同空间的点云更有效率。

稠密风：与流体模拟的当前趋势相反，本方法是使用密集域而不是稀疏域。在计算机图形学中，忽视空气运动在烟雾和火灾模拟中的贡献是一个常见的错误。虽然有些现象可以用稀疏域正确建模，但大多数现象都缺乏重要的特征，而这些特征往往被过多的程序湍流所取代。本程序允许用户将域的每一侧定义为固体边界或自由表面，自由表面可以与风力耦合，以模拟环境在有限的模拟域上的相互作用。用户可以使用 cudaNois [Lehtinen 2021] 定义程序矢量 Perlin 噪声，并可以在一秒钟内混合在域的两个相对侧的风向。

风从一侧灌入，从另外一侧流出。这可以产生比直接在渲染体素时添加程序湍流更自然的行为，因为它可以让解算器顺利处理程序风与模拟域中已经发生事件的集成。

空气湍流扰动：也可以将空气湍流注入模拟中。可以使用两套掩码来控制添加湍流的位置。第一套掩码使用温度将湍流保持在可渲染体积之外，因为我们通常希望在空气和密度之间的界面添加细节，而不会扰乱模拟的一般运动。第二个掩码来自速度大小，用于向模拟中可能出现强压力锋的快速移动区域添加更多湍流。这也使我们能够以不会影响模拟中非常缓慢移动部分的方式调节湍流。与风力一样，这种湍流也是使用cudaNoise生成的，但我们不是直接使用 Perlin 噪声 ，而是使用 [Bridson 2015] 中提供的公式将其转换为 curlnoise。

curlnoise 在体素中心计算，类似于散度计算，并暂时存储在速度back buffer中，然后散射到front buffer。散射步骤需要依赖原子函数，因为需要使用后向三线性插值对来自体素中心的量进行加权并将其分配到相应的面。由于写入风险，不使用原子函数会导致不确定的行为。为了控制浮点精度的 curl 计算生成的向量的大小，我们使用以下表达式限制 curlnoise 的大小。

涡度约束：很多求解器都引入了额外的涡度约束以防止湍流特征的缺失，该求解器也同样支持这一点。这一项技术背后的想法是恢复非发散投影期间丢失的一些动量。这里在实现上遵循 Bridson 书中的论述，结果在视觉上很吸引人，但MAC网需要进行速度平均才能正确计算旋度作为 上的中心差分，从而使得此操作的计算成本很高。与curlnoise实现类似，我们限制了涡度矢量的幅度，以防止由于数值误差而导致的不稳定性。

五，渲染实现

为了使用高度优化的硬件三线性插值，我们将温度网格复制到CUDA 3D 纹理中。这样做的缺点是增加了内存占用，但在撰写本文时，当前 GPU 硬件不允许对纹理进行原子操作，这解释了为什么必须同时维护一个用于模拟的温度网格和一个用于渲染的温度纹理 [Gao et al. 2018]。

高斯卷积近似多重散射:我们还以温度网格分辨率的1/4做为纹理分辨率创建另外一个CUDA3D纹理用于存储着色器指认的自发光贡献。然后，该纹理将与3x3x3离散高斯核进行迭代卷积，以近似多重散射。该卷积也有5 x 5 x 5和7 x 7 x 7的版本，核规模越大，迭代次数越低，但实验来看3 x 3 x 3 是最好的一个选择。

网格射线行进:通过对体素网格执行RayMarching来进行渲染，渲染器只支持定向光。在执行射线行进时，温度和近似的多重散射都会根据着色器定义进行采样和解释。温度值可以通过黑体辐射过程驱动自发光颜色。

将帧缓冲区写入磁盘:程序支持将多种文件格式，例如PPM，PNG和OpenEXR, 使用RLE压缩将文件写入磁盘以尽可能减少GPU等待时间。

六，结果展示

下面是使用提供的 HDA 创建的六个场景。场景以及用于生成它们的Houdini文件均与GitHub上的完整源代码一起提供。所有场景均使用配备 11 GB VRAM 的 Nvidia GeForce GTX 1080 Ti GPU 运行

具体性能如下表：