从零开始的 Position Based Fluid

| 0 | 总字数 4.5k | 期望阅读时间 17 min
  1. 1. Part 1 - Infrastructure
  2. 2. Part 2 - Mathematics
  3. 3. Part 3 - Implementation
  4. 4. Part 4 - Optimization
  5. 5. Reference

这个其实是 GAMES101 的大作业

Rendered by blender, cycles-x

(其实这是这篇 report 写的第三次了,前两次写的多多少少都有问题)

说是 report,我当然没打算把它当成 report 来写,就随意一点,写一些自己写这玩意儿中的心得体会。

作为一个工程和算法融合的项目,我真的是花了挺大精力来做这玩意儿的。

认认真真研究了数学部分(从七月中旬就开始了!),认真写了 CMake,认真组织了项目结构,认认真真优化代码,认认真真做抽象。
于是,目前终于接近尾声了。

Part 1 - Infrastructure

这个 Infrastructure 是一个比较粗略的概念,我在开始这个项目的时候,甚至连是使用 taichi 还是使用 C++ with OpenGL 都没有决定好。

我有一丢丢使用 taichi 的经验,从那一丢丢体验上来说,我还是不喜欢写 Python,于是我选择了写 C++,当然,也遇到了很多困难。

说起写 Python,我之前没有使用过任何的 Python 静态查错的工具,之前写的 btrfs auto backup 就因为一些类型一致性的问题出了 bug.

于是我设计了一套最最最最简单的可配置的 GUI,专门为 PBF 准备的。其能够显示每个 particle 的位置,与颜色(用一个 color ramp 从粒子密度映射到一些比较骚的颜色,方便之后的观察,事实证明很有用)

GUI 的接口非常简单,

class RTGUI_particles : public GUI {
public:
  RTGUI_particles(int WIDTH, int HEIGHT);
  ~RTGUI_particles() override = default;

  void set_particles(const std::vector<SPHParticle> &_p);
  void main_loop(const std::function<void()> &callback) override;
  void del();
  // ...
};

至于为什么单独引一个 del 函数而不是放在析构函数里,其实也是因为 OpenGL 的一些神秘特性,要求 glDelete.. 的时间需要提前一些。

其中,set_particles 提供了向外的最简单最简单的接口,这个简单的接口设计之后帮了我很大的忙。

main_loop 则以显示器的特性以特定的 interval 调用 callback 函数,callback 会阻塞 main_loop,反过来也会被刷新率阻塞。

这本是 glfw 相比于 glut 的一个“缺陷”,我使用它也只是因为习惯。但是这个设计莫名的有意思,而且很完美地解决了许多的问题。

还有一些可提升的空间,比如 GUI 端给 callback 传输 GUI 端的部分信息,以便实现 adaptive calculating 的功能。

当然,实践下来,这套无论是系统还是代码都可以复用了。
OpenGL 的细节我就不说了,反正这玩意儿贼难用,就写着写着就写完了。
实际上来说,之后可能要在这上面添加更多的 feature,比如不是以 glPoints 的形式,而是以小球体的形式来展示,但是我暂时不想写了。

Part 2 - Mathematics

这整个项目可能有两个难点,第一个就在这里了。对于这部分,我可能会写得更认真一些,毕竟不同于以前的其他工程项目,光实现细节都能写上万字。

从最开始讲起吧,一整个 PBF 要解决的问题很简单也很显然,但是解决的出发点实在让人觉得很神奇,更神奇的是仅仅使用这么简单的思路居然,就模拟出了液体?

思路是,液体的密度具有一致性。就是这么“简单”的思路。如果不是写过代码,我是不会相信这么一个思路可以运作的。

思路有了,将其数学化。使用 Lagrange view,我们需要使用离散的粒子去拟合出连续的场。基于 SPH 的方法,我们使用这个公式:

$$
\rho_i = \sum_{j}{m_j \mathrm{W}(\left |\boldsymbol{p}_i - \boldsymbol{p}_j \right |, h)}.
$$

其中 $\mathrm{W}$ 是核函数 poly6(很帅的名字),$h$ 是计算的区间。现在对于 $\rho_i$ 的定义清晰了,我们只需要抹去 $\rho_i$ 中的下标变为 $\rho$,同时抹去 $\boldsymbol{p}_i$ 的下标,就有了关于 $\boldsymbol{p}$ 的连续场 $\rho$.

虽然我们还是需要保留离散的形式。我们定义一个 constraint function

$$
C_i(\boldsymbol{p}_1, \cdots, \boldsymbol{p}_n) = \dfrac{\rho_i}{\rho_0} - 1,
$$

这里有一个需要注意的地方,$C_i$ 是关于 $\boldsymbol{p}_k$ 的函数,其变量是变化的所有点的位置,我第一次就被这里搞混了。

所以 $C_i$ 的梯度需要考虑一个具体的关于某个点,$\nabla C_i$ 本身并没有显然的意义(有意义的其实是 $\nabla_{\boldsymbol{p_k}}C_i$)
在考虑 $\nabla_{\boldsymbol{p}_k}C_i$ 时,应该以固定其他所有点,若只移动点 $k$,整体下降最快的方向是——来思考。

还需要深入理解的一点是 $C_i$ 的下标 $i$,它的含义也不是那么显然。对于静态的 $C_i$ 其实不太好考虑,我们从动态的视角去考虑。

首先,$C_i$ 其实是关于所有点位置的函数,所以考虑移动其中某一个点,移动了之后,$C_i$ 的值可能会发生变化。而这个值的意义,其实可以简单考虑为 $\rho_i$ 的变化,即这个粒子所代表的密度变化。

回到更高层的抽象,我们目前有每个 particle 的位置,即 $\boldsymbol{p}$。假设 $C(\boldsymbol{p})$ 的输出是一个列向量,每一行使用 $i$ 进行索引。我们需要一个 $C(\boldsymbol{p} + \Delta \boldsymbol{p})$,使得这个值最小化。其中这个 $\Delta \boldsymbol{p}$ 是所有点的位置向量组成的张量,是一个全局解。

但是,很显然,通过 $n$ 个方程解 $3n$ 个值是不可能的,我们需要考虑如何使用 $n$ 个值来替代这 $3n$ 个值。

于是我们给 $\Delta \boldsymbol{p}$ 限制了方向,$\Delta \boldsymbol{p} = \boldsymbol{\lambda} \nabla C(\boldsymbol{p})$,即这个向量场只能向最快下降方向生长。这样我们只需要 $n$ 个 $\lambda$。这部分具体的实现其实是比较模糊的,之后我还会阐释。
网络上的部分资料解释了这里为什么将 $\Delta \boldsymbol{p}$ 限制为 $\nabla C$ 方向,有动量守恒的考量。

好,那么我们将 $C$ 进行展开为一阶泰勒级数,

$$
\begin{aligned}
C(\boldsymbol{p} + \Delta \boldsymbol{p}) & \approx C(\boldsymbol{p}) + \nabla C^T \Delta \boldsymbol{p} \\
& \approx C(\boldsymbol{p}) + \boldsymbol{\lambda} \nabla C^T \nabla C \\
& \approx 0,
\end{aligned}
$$

这样我们只需迭代求解 $\lambda$ 就可以了.. 对吗?

解出来的每一个 $\lambda_i$,都是针对这行里所有的向量的。即,每一个 $C_i$,都会对 $n$ 个点的运动均产生限制,而最后每个点的运动的限制,就是这 $n$ 个 $C_i$ 产生的运动限制之向量和。
注意,下面这个公式把力拆分成了自己“需要”的作用与其他粒子对自己的作用之和,即:

$$
\begin{aligned}
\Delta \boldsymbol{p_i} &= \lambda_i \nabla_{\boldsymbol{p_i}} C_i + \sum_{j} \lambda_j \nabla_{\boldsymbol{p_j}} C_i \\
&= \dfrac{1}{\rho_0}{\sum_{j}{(\lambda_i + \lambda_j) \nabla_{\boldsymbol{p}_j}C_i}}
\end{aligned}
$$

最后需要解决的问题是,$\nabla_{\boldsymbol{p}_k}C_i$ 如何计算?即,在粒子 $i$ 的密度函数(伪)下,如何为 $n$ 个粒子分别计算一个运动方向。

这个公式是整个 PBF 的核心:

$$
\nabla_{\boldsymbol{p_k}}C_i = \dfrac{1}{\rho_0}
\left\{
\begin{aligned}
&\sum_{j}{\nabla_{\boldsymbol{p_k}}\mathrm{W}(\left |\boldsymbol{p_i} - \boldsymbol{p_j} \right |, h)} &\text{ if } k = i \\
&-\nabla_{\boldsymbol{p_k}} \mathrm{W}(\left| \boldsymbol{p_i} - \boldsymbol{p_j} \right|, h) &\text{ if } k = j
\end{aligned}
\right.
$$

这一个公式是一整篇论文的核心,它将前式分成了两个大部分来进行计算,一个是自己该向最稀疏的方向运动,其他所有粒子该向远离自己的方向运动。至于运动多少,就靠解出 $\lambda$

$$
\lambda_i = \dfrac{C_i(\bf{p_1}, \cdots, \bf{p_n})}{\sum_{k}{|\nabla_{\bf{p}_k}} C_i|^2}
$$

对于这个公式,因为存在分式,就需要考虑分母为 0 会导致系统不稳定,于是我们引入论文中提到的 CFM(constraint force mixing),向公式中添加松弛参数。

$$
\lambda_i = \dfrac{C_i(\bf{p_1}, \cdots, \bf{p_n})}{\sum_{k}{|\nabla_{\bf{p}_k}} C_i|^2 + \varepsilon}
$$

这个参数的取值也是个坑,之后再说。

接下来需要解决的是粒子聚集的问题,粒子的聚集可能有几个原因。(尽管凭借我现在的可视化,并无法观察到粒子聚集的情况)

  1. $C_i < 0$ 的情况,在对于粒子距离极其远的情况,比如两个粒子彻底分离的时候,这种这种时候若还添加约束,有可能使得粒子之间存在其他的排斥力。

  2. 在周围粒子不足以支撑当前粒子质量时,可能会导致大量粒子聚集,从而产生负压力。(我都忘了为啥了,但是解决问题的方案我会在下一个部分说)

总而言之,解决方法又是平滑平滑

$$
s_{corr} = -k\left(\dfrac{\mathrm{W}(\boldsymbol{p_i} - \boldsymbol{p_j}, h)}{\mathrm{W}(\Delta \boldsymbol{q}, h)}\right)^n
$$

其中 $\Delta \boldsymbol{q}$ 是一个固定距离,这个公式想要表达的是,在 $\Delta \boldsymbol{q}$ 一定时,两个粒子的距离越远,修正作用越弱。而这个修正被放在了求解 $\Delta \boldsymbol{p}$ 的公式中

$$
\Delta \boldsymbol{p_i} = \dfrac{1}{\rho_0}{\sum_{j}{(\lambda_i + \lambda_j + s_{corr})}} \nabla \mathrm{W}(\boldsymbol{p_i} - \boldsymbol{p_j}, h)
$$

Part 3 - Implementation

目前最终实现的结果在这里

效果如图(我有点懒得搞 gif 了,看静态吧!)

拖 OMP, GCC 的福,我很轻松就将其优化到了 16 threads 15k particles 60fps,单帧计算 13ms+,尽管还有挺大的优化空间,在自己的系统和并行计算知识完备前我就不做 CUDA 和其他多线程优化了,还有一个原因是我懒得去租 intel 的机器,没有一个足够好的 profiler(AMD!你看看人家的 vtune)

这里有一个比较有意思的小插曲,我在使用 spatial hashing[2] 时,不知道是什么原因(可能是巨大的常数,字面意思的常数),我只要开高等级的优化,编译时间就会很离谱(还没跑出来过)。在随手换成 index_sort 后,只换了 hash 函数,性能就极小数值几十 ms 的级别变为基本所有的操作都是 5ms-,应该是访存还有各类计算的问题,论文中给出了一句意味深长的话,大意是 spatial hashing 的思路本身与访存优化的连续访问相悖。

我首先实现了 $O(n^2)$ 的哈希方法,并根据 CompactNSearch 设计了一套简单的接口。

虽然名字取的是 CompactHash,实际上我没写 CompactHash,因为 index sort 的性能也还挺不错 x

class CompactHash {
 public:
  explicit CompactHash(float _radius);
  // ...
  void build();
  void add_particle(const SPHParticle &p);
  std::vector<SPHParticle> &get_data();
  int n_points() const;
  int n_neighbor(uint index) const;
  int neighbor(uint index, uint neighbor_index) const;
  std::vector<uint> &neighbor_vec(uint index);
  // ...
};

很可惜的是,在大常数下,$n^2$ 完全无法满足性能需求。我那时偷懒,一直懒于实现高性能的数据结构,直到后来完全无法预览,我在写出了第一版的加速数据结构后整个代码的结果终于能预览了。

但是!这时候的预览结果可以说是一团糟。理论上来说 PBF 应该能够在 $\Delta t \approx 0.01$ 时稳定的,但是我的这套系统哪怕已经到了 $\Delta t = 0.0016$ 都无法稳定运行,这事让人颇为迷惑。后来才反应过来,是核函数与整体 scale 的问题。如果我将 radius 设置为 0.01~0.1,poly6 在 $x = 0$ 时的值就会极大,导致(我们在求梯度时使用的是与之接近而 spiky kernel)两粒子在靠近时没有缓冲的空间。我将整个坐标系扩大了十倍,这时 poly6 的梯度明显减小,系统终于稳定了下来(困扰了我好几天!),未经优化而情况下,能在 3k particles 时跑到 20fps。

当然… 那些简单的优化我认为不必再提,不过是加几个 #pragma omp 的事情而已。白嫖了将近 50 倍的性能,已经接近别人在 CUDA 上的实现,但是与原论文的实现依然有较大差距。既然我没有能力 profile,就先暂时到此为止吧?

Part 4 - Optimization

回到学校后,我忽然又产生了继续完成这个项目的动力,于是我在其上再做了一些修改,并实现了 CUDA 版本。

目前 (2021/08/29),这份代码依然存在很大的真实性问题。(特别是在 mesh 之后,丢到 blender 里去预览那叫一个密恐福利)

我主要做了架构优化,CUDA 优化,和 mesh 支持三个改动,还有不少杂七杂八的改动我就不说了。

架构优化比较复杂,也还没做完。比如对 radius 和 border 等等参数的设置依然完全不 OOP. 但是把这几个都 OOP 的话,后面某些部分就会有麻烦。

比如 CUDA 的部分,我对 CUDA 那叫做一个入门,我只简单读了 thrust 的文档,看了两三章 CUDA by Example,不过幸亏这个 project 对 CUDA 的要求比较低,只需要按照最套路的 <<<num_threads, threads_per_blocK>>> 来跑即可。然后在有 overlap 的哈希表上处理一下,你就会发现,CUDA 是真他妈的快。

只可惜 __global__ 函数必须位于全局作用域,于是所有的 OOP 使用的参数必须一个一个传给他,一个典型的参数列表大概是这样的。

__global__ void fill_lambda(SPHParticle *dev_data,
                            size_t data_size,
                            int *dev_neighbor_map,
                            const int *dev_n_neighbor_map,
                            size_t pitch_neighbor,
                            float *c_i,
                            float *lambda,
                            float rho_0,
                            float mass)
// ...

然后每一次调用都需要传入全部的.. 你说噩梦不噩梦。我本来想用宏解决,后来还是决定直接写好了。

将原本的搬到 CUDA 也就两个难点,注意内存的原子访问的问题,当然,如果下标互不相同不用担心,在有可能会同时访问的下标加个锁即可。

这里有个坑,他妈的。我后来才注意到。我的迭代过程是这样的:

最开始,CPU 上的代码能够正常运行,分为数据结构实现 nsearch.hpp/cpp 和算法实现 pbd.hpp/cpp 两个文件。先修改其中任意一者都可以,但是不要同时修改两者。

我先重构了 pbd.cpp,其输入稍微复杂一些,原本是一个 vectorvector 的动态二维数组,我将其改成了一个静态的二维数组与一个静态的数组,设置一个最大长度,尝试模拟一个 vectorvector。(我把所有的动态数组都改成了这种形式了)

在极其激进的流体分布下可能会崩溃,大部分时候还好。而且这样保证了显卡内部的内存对齐与连续访问,所以重构后,pbd.cu 的输入是一个 neighbor map,表示每个粒子的临近粒子的下标与个数,还有一个 data 数组表示所有的粒子。

输出依然是 data,当然,为了抽象,我们不用考虑其内容是在什么平台上,先统一拷贝到内存中。经过简短的调试,目前可以运行。但是 cudaMemcpyDtoH 耗时太长,因为 neighbor map 实在不小,完全无法满足 real-time 的需求。于是下一步是把建立 neighbor map 的数据结构搬到 CUDA 上.

这步搬迁的思路和之前一样,只是在 build hash map 的时候,我因为经验不足和个人傻逼狠狠地被坑了一回,代码如下

// ...
  if (i < data_size) {
    int hash_map_index = PBDSolver::hash(dev_data[i].pos, n_grids);
    int *dev_hash_map_item = (int *)((char *)dev_hash_map +
                                     pitch_hash * hash_map_index);
    if (dev_n_hash_map[hash_map_index] < PBDSolver::MAX_NEIGHBOR_SIZE) {
    while (atomicCAS(&hash_map_mutex[hash_map_index], 0, 1) != 0)
      ;
      atomicExch((int *)&(dev_hash_map_item[dev_n_hash_map[hash_map_index]]),
                 (int)i);
      atomicAdd((int *)&(dev_n_hash_map[hash_map_index]), 1);
      atomicExch(&hash_map_mutex[hash_map_index], 0);
    }
  }
}

我一开始看起来没啥问题,盯了半天也不觉得有什么问题,后来我把 lock 调整到了 if 外才发现这个问题。原来我没锁住 dev_n_hash_map 的访问,这整个 if 应该都被锁住。…

就莫名奇妙地被卡住,然后莫名奇妙地解决。

不过说实话我不太习惯 thrust 的玩法,怎么有点别扭呢… 据猫猫所说,thrust 的函数会同时调用 CPU 和 GPU 导致超功耗。不过说实话,你确实用不满 CPU 到 GPU 的带宽(经过 profile,大概只使用了 1/5),用 CPU 的算力优化一下也无可厚非。

不过不管 CPU 怎么 vectorize 加上 parallize 也跑不过 gpu 啊,在这些计算密集型。虽然 nvidia 的 profiler 乱成一锅粥,我都弄不清楚有几个(我用过的已经有三个?至少比 uprof 好用,AMD 骂的就是你)。经过两天的编码与调试,总算是弄完了。

不过目前,我写的都是 basic CUDA,UM 啊,任何一点稍微高级一点的语法和 feature 都没有使用。如果不出意外的话,之后会 bfs 下去。

关于 mesh 的部分我就不多说了,因为我偷懒了用了 OpenVDB,并没有自己去实现 marching cubes. 但是添加了 OpenVDB 之后,编译时间变得梦幻了起来.. 暂时不太清楚如何解决,如果知道了我会添加。(我的头文件不乱啊!)

这个项目的下一步是实现 WCSPH,并添加 caustics 的支持(既然是水面,当然少不了焦散对吧!)

之后再见!

Reference

[1]: Macklin, Miles, and Matthias Müller. “Position Based Fluids.” ACM Transactions on Graphics (TOG) 32.4 (2013): 104.

[2]: Ihmsen, M., Akinci, N., Becker, M., & Teschner, M. (2010). A Parallel SPH Implementation on Multi-Core CPUs. Computer Graphics Forum, 30(1), 99–112.