Fluid Engine Development PIC/FLIP 代码分析

本文主要是介绍Fluid Engine Development PIC/FLIP 代码分析，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！

把 Fluid Engine Development 看完了，但是仍然感觉不懂

https://github.com/doyubkim/fluid-engine-dev

感觉还是应该了解整体代码怎么写的，所以做个总结

看着看着，感觉还是从底层开始看起

从底层重新开始看的时候，感觉就来了

而且作者也有很多注释，感觉能够体会到别人的思路

他这里也有很多内容，我选择从 PIC/FLIP 开始看起

然后我选择它的 hybrid_liquid_sim 工程

它的 main 文件就很容易找到了，在 src\examples\hybrid_liquid_sim\main.cpp

Animation

Animation 类里面提供了 update 函数，update 函数里面调用了虚函数 onUpdate

虚函数 onUpdate 用于给子类提供具体的实现的

这就完成了第一层的抽象，就是把所有的动画都视为一个可以随时间更新的东西

那么我们用户统一地使用 update 函数去更新它

具体怎么更新，那就是程序员怎么在子类中实现虚函数 onUpdate 的事情了

也就是下层（怎么实现虚函数 onUpdate）对上层（用户调用 update）透明

例如，在 src\examples\hybrid_liquid_sim\main.cpp 中，有不同的求解器和情景可供选择，PIC，FLIP，APIC 等

最后这些求解器都进入函数 runSimulation，在其中，求解器的更新都是只需要一句 solver->update(frame);

就很优雅

PhysicsAnimation

继承 Animation

这个类主要实现了把某一帧拆分成若干个子步骤的功能

具体实现方法就是它 final 重写了虚函数 onUpdate，限定死了怎么处理某一个帧的

因为他一定要保证他自己这个拆分子步骤的功能嘛，所以是 final

具体内容就是，只有帧数实参大于类中存储的当前帧号时，才说明这个类要沿时间推进

if (frame.index > _currentFrame.index) {...

初始化

推进的具体办法是，帧号 < 0 时调用函数 initialize 进行初始化，

if (_currentFrame.index < 0) {initialize();
}...

这里涉及到了初始化函数，是因为初始化函数的调用跟帧号有关，所以必须在这里给一个钩子，才能在这个正确的时间调用初始化

每个算法的初始化肯定是不一样的，所以如果我们对初始化的时间没有要求的话，我们肯定不用这么麻烦，但是就是因为他是跟帧号相关，所以才要提早到这里

那么自然地，这个初始化函数应该是一个虚函数

但是它这里的 initialize 还不是虚函数。反而，initialize 里面调用了虚函数 onInitialize，它的设计是让子类去重写这个虚函数

可能这是它自己的设计哲学吧，没懂为什么要额外一层

如果不用初始化，那么就统计要前进的帧数，对每一帧都调用 advanceTimeStep。也就是说，对于每一帧，都尝试拆分子步骤

其实这个函数名和 update 也很像，其实不用纠结，他就是为了完成拆分子步骤

...int32_t numberOfFrames = frame.index - _currentFrame.index;for (int32_t i = 0; i < numberOfFrames; ++i) {advanceTimeStep(frame.timeIntervalInSeconds);
}_currentFrame = frame;

固定子步骤数量

具体他是怎么实现拆分子步骤的呢？也很简单

如果你设置了固定子步骤数量 n，那么就把这帧的时间平均分为 frame time/n 份，对每份都执行虚函数 onAdvanceTimeStep 也就是相当于把实际每一步的流体求解又推迟到这个 onAdvanceTimeStep 函数中了

const double actualTimeInterval =timeIntervalInSeconds /static_cast<double>(_numberOfFixedSubTimeSteps);for (unsigned int i = 0; i < _numberOfFixedSubTimeSteps; ++i) {onAdvanceTimeStep(actualTimeInterval);_currentTime += actualTimeInterval;...

为什么说是推迟呢？想象一下如果我们始终不知道要拆分子步骤这件事，那么流体求解部分应该在子类中重写 Animation 类的虚函数 onUpdate

但是现在的话，我们的流体求解部分要放到在 PhysicsAnimation 的子类对虚函数 onAdvanceTimeStep 的重写中了，所以是，我觉得这就像推迟了一个函数

自适应子步骤数量

如果你没有设置固定子步骤数量，那就是你要是用自适应子步骤数量。但是自适应算法应该是有不同的实现的，所以这里给出一个虚函数 numberOfSubTimeSteps 允许子类实现这个自适应的算法

虚函数 numberOfSubTimeSteps 接受剩余的时间作为参数，这应该是自适应算法的常规输入吧

获取了自适应的子步骤数量 n 之后，当前子步骤的时间步是 remaining time/n

剩余的时间 remaining time 最初是 frame time

执行了一步之后是 remaining time -= remaining time/n

执行完一步之后，下一次将 remaining time 输入到虚函数 numberOfSubTimeSteps 获得下次的子步骤数量 n2，下一次的子步骤的时间步就是 remaining time/n2

所以每一个子步骤的时间步可能是不一样的

double remainingTime = timeIntervalInSeconds;
while (remainingTime > kEpsilonD) {unsigned int numSteps = numberOfSubTimeSteps(remainingTime);double actualTimeInterval =remainingTime / static_cast<double>(numSteps);onAdvanceTimeStep(actualTimeInterval);remainingTime -= actualTimeInterval;_currentTime += actualTimeInterval;...

GridFluidSolver3

这里开始加入了很多东西

我觉得应该从他对虚函数的重写开始看起

其他的他自己定义的函数，看上去像是为了实现它对虚函数的重写中的逻辑，而创建的

初始化

初始化函数中，要初始化边界条件，边界条件就包括了碰撞体和粒子发射器

void GridFluidSolver3::onInitialize() {updateCollider(0.0);updateEmitter(0.0);

这两个更新函数，就很简单，就是调用各自的 update 而已

_collider->update(...);
_emitter->update(...);

我觉得具体的可以之后再看，比如碰撞体那部分的类是怎么样的

还是回到 GridFluidSolver3

对初始化函数的重写就是这些了，那么下一个是步进函数，对虚函数 onAdvanceTimeStep 的重写

步进 求解框架

void GridFluidSolver3::onAdvanceTimeStep(double timeIntervalInSeconds) {beginAdvanceTimeStep(timeIntervalInSeconds);computeExternalForces(timeIntervalInSeconds);computeViscosity(timeIntervalInSeconds);computePressure(timeIntervalInSeconds);computeAdvection(timeIntervalInSeconds);endAdvanceTimeStep(timeIntervalInSeconds);
}	

这个结构就很清晰，这就是求解 NS 方程要做的事情

前处理 后处理

前后的 beginAdvanceTimeStep endAdvanceTimeStep 分别包括了对虚函数 onBeginAdvanceTimeStep onEndAdvanceTimeStep 的调用，为更具体的流体算法子类提供虚函数接口

额外的是，beginAdvanceTimeStep 里面还有更新碰撞体和发射器，然后还要调用 GridBoundaryConditionSolver 的更新 updateCollider 和速度约束函数 constrainVelocity

当你想要看一下这个 GridBoundaryConditionSolver 是啥的时候，跳到定义，你就会看到他是又是一个类的指针，顺便也看到求解平流、扩散（粘性力）、压力的方法也是调用别的类的指针

AdvectionSolver3Ptr _advectionSolver;
GridDiffusionSolver3Ptr _diffusionSolver;
GridPressureSolver3Ptr _pressureSolver;
GridBoundaryConditionSolver3Ptr _boundaryConditionSolver;

这和碰撞体和发射器也是一样的

GridSystemData3Ptr _grids;
Collider3Ptr _collider;
GridEmitter3Ptr _emitter;

所以现在对于 GridFluidSolver3 就有一个大概的印象了，它定义了求解 NS 方程的大概框架，分解了求解步骤（外部力、粘性力、压力）但是在每个力的具体计算，实际上还是依赖更进一步的定义。这通过提供工具类的示例来实现。比如确定怎么求解压力的话，需要给 GridPressureSolver3Ptr _pressureSolver 赋值。或者说这是一种依赖注入的思路。

有点突然理解了依赖注入！主要是，概念很容易理解，就是调用方需要的东西，由别人来提供。主要是我觉得我似乎理解了为什么要这么做，为什么调用方要的东西要别人来给？因为调用方只知道算法流程，但是具体怎么实现还是要看子类的，或者别的类之类……？总之就是有种老板提方案员工负责执行的感觉

外部力

然后是 computeExternalForces。里面只有一个处理重力的函数 computeGravity

这个函数是第一个涉及到具体怎么计算物理量的函数

因为所有的变量都是 auto，所以理论上应该不用管类型也很容易理解代码

实际上也是这样的，获取速度场，从速度场获取访问器

速度场是有一个遍历所有元素的函数 forEachUIndex，这个函数可以输入一个 lambda 函数，在这个匿名函数里，访问器就相当于物理量，用它来实现物理公式就好了

auto vel = _grids->velocity();
auto u = vel->uAccessor();
auto v = vel->vAccessor();
auto w = vel->wAccessor();vel->forEachUIndex([&](size_t i, size_t j, size_t k) {u(i, j, k) += timeIntervalInSeconds * _gravity.x;
});vel->forEachVIndex([&](size_t i, size_t j, size_t k) {v(i, j, k) += timeIntervalInSeconds * _gravity.y;
});vel->forEachWIndex([&](size_t i, size_t j, size_t k) {w(i, j, k) += timeIntervalInSeconds * _gravity.z;
});...

这里确实是需要知道一下“场”这个类的概念，才有这种直观的认识。但是既然都过了一遍书本了，基本的印象还是有的。

然后就是要调用一下 applyBoundaryCondition 函数

观察发现，其他计算力的函数后面都有这个

因为 applyBoundaryCondition 里面是使用约束的那个工具类来约束速度，所以我们可以知道它的设计目的就是，每一次更新速度之后，都要约束一下速度。

粘性力 压力

在 computeViscosity 和 computePressure 里面，单纯是对 _diffusionSolver 和 _pressureSolver 的 solve 调用

那么它的设计目的就是，粘性力和压力的计算也是根据具体算法来的

平流

在 computeAdvection 也是单纯对各个变量的 solve 的调用

他看上去复杂，只是因为他要考虑所有标量，所有定义在网格中心的向量 CollocatedVectorGrid3 和定义在网格面上的向量 FaceCenteredGrid3

它把所有向量场都保存在一个列表中，以向量场的基类保存，所以你在遍历这个列表的时候不知道每个元素是 CollocatedVectorGrid3 还是 FaceCenteredGrid3，所以他这里做了两次 cast，把指向基类的指针 cast 到指向子类的指针，指针不为 null 表示这个示例确实是这个子类。

那么现在对这两个子类有不同的处理情况，平流求解器 _advectionSolver 是需要知道场的类型的（表现在函数重载上），但是我们可以把这两种处理情况串行写在一起。因为如果是子类 A 就不会是子类 B，所以实际上只会执行其中一个

auto collocated =std::dynamic_pointer_cast<CollocatedVectorGrid3>(grid);
auto collocated0 =std::dynamic_pointer_cast<CollocatedVectorGrid3>(grid0);
if (collocated != nullptr) {_advectionSolver->advect(*collocated0, *vel,...
}auto faceCentered =std::dynamic_pointer_cast<FaceCenteredGrid3>(grid);
auto faceCentered0 =std::dynamic_pointer_cast<FaceCenteredGrid3>(grid0);
if (faceCentered != nullptr && faceCentered0 != nullptr) {_advectionSolver->advect(*faceCentered0, *vel,...
}

最后还有速度要平流自己

这里还有另外一个设计就是，如果更新了速度，比如这里平流了速度，那么之后就要跟上对速度的约束 applyBoundaryCondition；如果是更新了标量场或者是向量场，那么之后就要跟上对这个场的外插 extrapolateIntoCollider，这都是为了保证每一步之后，每个场，包括标量场、向量场、速度场都是正确的（不是指没有耗散，我指的是，逻辑上是正确的，更新后的值，不会出现那种新值与旧值并存的情况）

自适应子步骤数量

这里的逻辑也很简单，就是算当前的 CFL 条件数

$$(\max V)\Delta t/\Delta x$$

这个东西最直观的意义就是速度最快的粒子单位时间走过的格子的数量，希望他不要太大

所以如果他太大的话，就用一个约定的最大值去除，比如你速度最大会走 10 个格子，但是我约定只能走 5 个，所以我就拆成两个子步骤

SDF

传给 _diffusionSolver _pressureSolver _advectionSolver 的碰撞体的 SDF 和碰撞体的速度都是从边界条件工具类 _boundaryConditionSolver 中得来

也就是这个工具类来负责

而流体的 SDF，在这个 GridFluidSolver3 里面是调用了一个虚函数 fluidSdf 这表示子类来负责流体的 SDF 怎么算

边界条件

边界条件这里的怎么约束速度的方法在工具类里面，暂时看不到

但是将值外推到边界外的方法 extrapolateIntoCollider 在这里定义

它定义了三个重载，分别处理不同类型的值，ScalarGrid3 CollocatedVectorGrid3 FaceCenteredGrid3 的

三种情况的逻辑都是一样的，新建一个标记数组 marker，它的大小和要外插值的网格的大小相当

然后对每一个格子，都基于碰撞体的 SDF 判断这个点在不在碰撞体内，将结果记入 marker

然后得到的 marker 和要插值的值都输入到函数 extrapolateToRegion 中，这是一个更一般的工具函数。因为 marker 一旦构建了，他就可以表示任意边界，不一定是碰撞体的边界，还可以是流体的边界等等，所以需要这么一个工具函数

void GridFluidSolver3::extrapolateIntoCollider(CollocatedVectorGrid3* grid) {Array3<char> marker(grid->dataSize());auto pos = grid->dataPosition();marker.parallelForEachIndex([&](size_t i, size_t j, size_t k) {if (isInsideSdf(colliderSdf()->sample(pos(i, j, k)))) {marker(i, j, k) = 0;} else {marker(i, j, k) = 1;}});unsigned int depth = static_cast<unsigned int>(std::ceil(_maxCfl));extrapolateToRegion(grid->constDataAccessor(), marker, depth,grid->dataAccessor());
}

在这个工具函数中，指定迭代次数

每一次迭代，都遍历 marker 中的所有元素，对于那些 flag 为 0 的位置，尝试从他上下左右前后 6 个邻居中找 flag 为 1 的元素，如果找到了就把邻居的值统计平均赋给当前位置，并且置当前位置的 flag 为 1

这样的话，刚好在边界外一格远的位置，有可能获得插值，如果是在边界外很远的地方，它的邻居都不在边界内，那么这次迭代之后他依然也不会获得插值

也就是，每一次迭代，都会把边界整体往外拓展一格（如果空间允许的话）

现在我们看完了 extrapolateToRegion。我们知道了 extrapolateToRegion 只需要知道怎么访问数据就行了，也就是得到 dataAccessor 就行了，dataAccessor 的类型与 ScalarGrid3 CollocatedVectorGrid3 FaceCenteredGrid3 无关，那么这里为什么要分三种情况呢

因为实际上不同类型的 Grid 的尺寸是不一样的，

Size3 CellCenteredScalarGrid3::dataSize() const {// The size of the data should be the same as the grid resolution.return resolution();
}Size3 VertexCenteredScalarGrid3::dataSize() const {if (resolution() != Size3(0, 0, 0)) {return resolution() + Size3(1, 1, 1);} else {return Size3(0, 0, 0);}
}

这篇关于Fluid Engine Development PIC/FLIP 代码分析的文章就介绍到这儿，希望我们推荐的文章对编程师们有所帮助！

---