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

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

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

Grid 的尺寸

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

因为实际上不同类型的 Grid 的尺寸是不一样的，而获取正确的 dataAccessor 需要 Array 具有正确的尺寸

这里展示不同的类定义了不同的尺寸的 getter

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);
    }
}

怎么从这个 getter 到实际存储数据的类的尺寸？

查找 dataSize 的引用，发现它在 onResize 中被用到，并且这个函数是在 CollocatedVectorGrid3 中被定义的

进一步发现这个 onResize 是 CollocatedVectorGrid3 FaceCenteredGrid3 这一个层级的

对比就发现，这两个类的区别就在于 CollocatedVectorGrid3 是直接定义一个 Vector 作为 data 了，但是 FaceCenteredGrid3 因为把速度的不同分量定义在不同的面，所以现在速度的不同分量的位置不同了，就不能组成一个 vector 了（组成 vector 就说明三个分量在同一个位置）所以速度要存在三个不同的数组里面

所以 CollocatedVectorGrid3 就只用 resize 一个数组，而 FaceCenteredGrid3 需要 resize 三个数组

void CollocatedVectorGrid3::onResize(
    const Size3& resolution,
    const Vector3D& gridSpacing,
    const Vector3D& origin,
    const Vector3D& initialValue) {
   
    _data.resize(dataSize(), initialValue);
    ...
}

void FaceCenteredGrid3::onResize(const Size3& resolution,
                                 const Vector3D& gridSpacing,
                                 const Vector3D& origin,
                                 const Vector3D& initialValue) {
   
    if (resolution != Size3(0, 0, 0)) {
   
        _dataU.resize(resolution + Size3(1, 0, 0), initialValue.x);
        _dataV.resize(resolution + Size3(0, 1, 0), initialValue.y);
        _dataW.resize(resolution + Size3(0, 0, 1), initialValue.z);
    } else {
   
        _dataU.resize(Size3(0, 0, 0));
        _dataV.resize(Size3(0, 0, 0));
        _dataW.resize(Size3(0, 0, 0));
    }
    ...
}

这里可能一开始有点陌生，因为我们是顺着 solver 看过来的，没有从数据结构看过来，所以对数据结构可能有点陌生，但是看过了书的话应该有点印象

CollocatedVectorGrid3 和 FaceCenteredGrid3 的 onResize 都是对基类 VectorGrid3 的重写

VectorGrid3 在自己的 resize 中调用了 onResize，说明 onResize 只是一个钩子

之后，VertexCenteredVectorGrid3 CellCenteredVectorGrid3 FaceCenteredGrid3 在自己的构造函数中都会调用 resize，这就实现了最终对 Array 的 resize 的调用

Builder

感觉现在已经知道大概的逻辑了，但是回到 main 开始看的话，他这里创建求解器是使用一个 builder

看了一下，感觉他这思路就是，创建一个 builder 类，这个类提供一些函数用于配置默认参数，最后返回一个实例的共享指针

auto solver =
    PicSolver3::builder()
        .withResolution({
   resolutionX, 2 * resolutionX, resolutionX})
        .withDomainSizeX(1.0)
        .makeShared();

因为 B 和 C 都是继承自 A，所以只要 B 和 C 的 builder 都是继承自 A 的 builder，那么 B 和 C 配置默认参数的方法就会非常类似

这样子设计，应该是为了初始化配置上的相似性？

因为实际上 B 和 C 也没有什么自己的独特的要初始化的参数，要初始化的全都在基类 A 中定义好了

话又说回来，正是因为这样，才需要继承关系，B 和 C 只需要基类 A 定义的初始化函数，应该是理想的完美继承关系吧

构造函数

然后我打算再看 GridFluidSolver 的构造函数

GridFluidSolver3::GridFluidSolver3(const Size3& resolution,
                                   const Vector3D& gridSpacing,
                                   const Vector3D& gridOrigin) {
   
    _grids = std::make_shared<GridSystemData3>();
    _grids->resize(resolution, gridSpacing, gridOrigin);

    setAdvectionSolver(std::make_shared<CubicSemiLagrangian3>());
    setDiffusionSolver(std::make_shared<GridBackwardEulerDiffusionSolver3>());
    setPressureSolver(
        std::make_shared<GridFractionalSinglePhasePressureSolver3>());
    setIsUsingFixedSubTimeSteps(false);
}

因为这里前面看的都是 GridFluidSolver 求解 NS 方程的框架，实际上具体的求解逻辑还是要看工具类，而构造函数这里就是给出了怎么配置这些工具类

之后的类如果是继承 GridFluidSolver，比如 PicSolver，那么因为是继承嘛，所以父类构造函数也会被调用，也就相当于这些工具类的配置在 GridFluidSolver 这里就是默认值了

AdvectionSolver

AdvectionSolver 的虚函数 advect 有三个重载，分别处理 ScalarGrid3 CollocatedVectorGrid3 FaceCenteredGrid3

这个虚函数确定了输入的场，输出的场，速度场，还有边界 sdf

感觉看到这里的话，就知道这个代码里面大部分都用 sdf 来表示边界了

就很统一

然后把边界作为参数，就是说明一定是要考虑边界

给人感觉就是规范了怎么求解

SemiLagrangian

同样地，CollocatedVectorGrid3 和 FaceCenteredGrid3 的区别就是前者只是一个场做向后追踪，而后者是三个分量都要分别做向后追踪

向后追踪的时候，也要拆分若干个子步骤，也是跟之前一样，根据 CFL 的物理意义来拆分，也就是所谓的自适应子步骤。算出了子步骤的数量之后，当前子步骤的 dt 也是用 remaining time / n 来算

// Adaptive time-stepping
Vector3D vel0 = flow.sample(pt0);
double numSubSteps
    = std::max(std::ceil(vel0.length() * remainingT / h), 1.0);
dt = remainingT / numSubSteps;

然后找点用的是中点法

// Mid-point rule
Vector3D midPt = pt0 - 0.5 * dt * vel0;
Vector3D midVel = flow.sample(midPt);
pt1 = pt0 - dt * midVel;

然后要考虑边界问题，通过两个点的 SDF 值相乘是否 < 0 来判断这两个点是否是一个在边界外一个在边界内

如果是的话，那就说明向后追踪的时候追到外面去了，那么就要把点钳制在边界上

这里通过线性插值把点定位在边界上

如果遇到了这种跨过边界的情况，那么钳制之后就直接结束循环了，不用继续执行剩下的子步骤了

double phi0 = boundarySdf.sample(pt0);
double phi1 = boundarySdf.sample(pt1);

if (phi0 * phi1 < 0.0) {
   
    double w = std::fabs(phi1) / (std::fabs(phi0) + std::fabs(phi1));
    pt1 = w * pt0 + (1.0 - w) * pt1;
    break;
}

这就是向后追踪的过程了

最后把向后追踪得到的位置是一个三维世界坐标，不是网格序号，所以要用 sampler，这个 inputSamplerFunc 返回的就是网格的 sampler

auto outputDataPos = output->dataPosition();
auto outputDataAcc = output->dataAccessor();
auto inputDataPos = input.dataPosition();

output->parallelForEachDataPointIndex([&](size_t i, size_t j, size_t k) {
   
    if (boundarySdf.sample(inputDataPos(i, j, k)) > 0.0) {
   
        Vector3D pt = backTrace(
            flow, dt, h, outputDataPos(i, j, k), boundarySdf);
        outputDataAcc(i, j, k) = inputSamplerFunc(pt);
    }
});

这里就可以看出来 sampler 和 accessor 的区别

sampler 是根据世界坐标来获得值的，而 accessor 是根据网格的编号直接访问数组元素

forEachIndex parallelForEachIndex

我突然很好奇的就是，这样子直接 foreach 的话，如果传入的 input 和 output 是同一个的话，又没有做缓存，那么应该会出现数据相关的问题吧

我主要是想确认一下他这里面确实没有更深一层东西来解决数据相关的问题，于是看了一下 forEachIndex

Grid 的 forEachIndex 实际上是调用的 Array 的 ArrayAccessor 的 forEachIndex

最终看 Accessor 的 forEachIndex 发现他这个就是一个单纯的串行的多层循环，用模板把不同层循环都封装成一个了

看到这里的时候也要注意他这个是列优先存储

template <typename T>
template <typename Callback>
void ArrayAccessor<T, 1>::forEachIndex(Callback func) const {
   
    for (size_t i = 0; i < size(); ++i) {
   
        func(i);
    }
}

template <typename T>
template <typename Callback>
void ArrayAccessor<T, 2>::forEachIndex(Callback func) const {
   
    for (size_t j = 0; j < _size.y; ++j) {
   
        for (size_t i = 0; i < _size.x; ++i) {
   
            func(i, j);
        }
    }
}

template <typename T>
template <typename Callback>
void ArrayAccessor<T, 3>::forEachIndex(Callback func) const {
   
    for (size_t k = 0; k < _size.z; ++k) {
   
        for (size_t j = 0; j < _size.y; ++j) {
   
            for (size_t i = 0; i < _size.x; ++i) {
   
                func(i, j, k);
            }
        }
    }
}

那么更进一步，查看 ArrayAccessor 的 parallelForEachIndex 的话，可以看到它们最终都是调用 parallelFor

同样的，parallelForEachIndex 不同的变体只是为了处理不同的维度

最终可以找到 parallelFor 里面有各种方法实现的多线程并行计算

void parallelFor(IndexType start, IndexType end, const Function& func,
                 ExecutionPolicy policy)

至于这个 parallelRangeFor 输入三个维度的范围，最终居然只拆分 z 维度的范围，很神奇

sampler

因为平流那里涉及到了 sampler，所以想看一些具体是怎么做的

CubicSemiLagrangian 与 SemiLagrangian3 的差别也就只在于这个采样器的配置

虽然我知道这里面有插值

记得看书的时候也看到过，一开始书上讲的是线性插值然后讲三次样条插值

ScalarGrid3 CollocatedVectorGrid3 FaceCenteredGrid3 都是默认的 LinearArraySampler

插值方法都要把世界坐标用网格中心和网格单元尺寸进行归一化

这让我突然想到了，其实这里所有的坐标都是在物体坐标系的，如果你这个流体模拟外部还有别的坐标系，跟这个也无关的

所以才没有别的什么额外的归一化操作，比如处理整个 box 的旋转啊什么的，没有，因为 box 自己是一个物体坐标系

然后他这个求重心坐标的方法我还真看不懂

对于线性插值来说，在 i, j, k 和 i+1, j+1, k+1 之间这八个块内三线性插值

对于三次样条插值，在 i-1 i i+1 i+2 之间这 64 个块内做单调 Catmull-Rom 插值

这个单调 Catmull-Rom 在书中也有讲，搜了一下，这个也是可以做文章的，就是说具体怎么设计这个插值函数

但是具体单调 Catmull-Rom 怎么做的还没看懂……看书的时候就没看懂，估计之后要找文章来对应地看

GridDiffusionSolver

GridDiffusionSolver 也就是提供了虚函数，确定了形参

GridForwardEulerDiffusionSolver

GridForwardEulerDiffusionSolver3 更简单，因为是前向的，所以不用算矩阵求解，直接逐项更新就好了

void GridForwardEulerDiffusionSolver3::solve(
    const ScalarGrid3& source,
    double diffusionCoefficient,
    double timeIntervalInSeconds,
    ScalarGrid3* dest,
    const ScalarField3& boundarySdf,
    const ScalarField3& fluidSdf) {
   
    auto src = source.constDataAccessor();
    Vector3D h = source.gridSpacing();
    auto pos = source.dataPosition();

    buildMarkers(source.resolution(), pos, boundarySdf, fluidSdf);

    source.parallelForEachDataPointIndex(
        [&](size_t i, size_t j, size_t k) {
   
            if (_markers(i, j, k) == kFluid) {
   
                (*dest)(i, j, k)
                    = source(i, j, k)
                    + diffusionCoefficient
                    * timeIntervalInSeconds
                    * laplacian(src, _markers, h, i, j, k);
            } else {
   
                (*dest)(i, j, k) = source(i, j, k);
            }
        });
}

GridBackwardEulerDiffusionSolver

GridBackwardEulerDiffusionSolver 用了后向方法，所以就要求解矩阵了

因为是求解矩阵，所以可以看到多了构建系数矩阵和列向量的步骤，也就是 Ax=b 中的 A 和 b

void GridBackwardEulerDiffusionSolver3::solve(
    const ScalarGrid3& source,
    double diffusionCoefficient,
    double timeIntervalInSeconds,
    ScalarGrid3* dest,
    const ScalarField3& boundarySdf,
    const ScalarField3& fluidSdf) {
   
    auto pos = source.dataPosition();
    Vector3D h = source.gridSpacing();
    Vector3D c = timeIntervalInSeconds * diffusionCoefficient / (h * h);

    buildMarkers(source.dataSize(), pos, boundarySdf, fluidSdf);
    buildMatrix(source.dataSize(), c);
    buildVectors(source.constDataAccessor(), c);

    if (_systemSolver != nullptr) {
   
        // Solve the system
        _systemSolver->solve(&_system);

        // Assign the solution
        source.parallelForEachDataPointIndex(
            [&](size_t i, size_t j, size_t k) {
   
                (*dest)(i, j, k) = _system.x(i, j, k);
            });
    }
}

构建矩阵和列向量的思路在书上确实讲了

书上先讲了一个一维的例子，然后再说，推广到多维也很简单，就是把多维线性化。对于 i j k 的循环，最内层会算出线性化的序号，对于每个点，点的中心是 kc+1，点的邻居是 -c

其实一开始我看这个的话，只是对着书上那个二维例子对照着看，但是没有真正理解三维情况下是什么状况

书上的例子：

$$\left[\begin{array}{cccccc}c+1& -c& 0& ...& 0& 0\\ -c& 2c+1& -c& ...& 0& 0\\ ⋮& ⋮& \ddots & ⋮& ⋮& \\ 0& 0& ...& -c& c+1& \end{array}\right]\cdot f^{n+1}=f^n$$

这里可能有很多坑……或者只是我个人的问题

系数矩阵的意义

第一个是没有注意到线性化的话，或者是没有注意看书的话，那么就会想着，这个矩阵是在表示二维情况，会以为这个矩阵和网格坐标有这样的对应关系

$$\left[\begin{array}{cccccc}(0,0)& (0,1)& (0,2)& ...& (0,n-1)& \\ (1,0)& (1,1)& (1,2)& ...& (1,n-1)& \\ ⋮& ⋮& \ddots & ⋮& ⋮& \\ (n-1,0)& (n-1,1)& (n-1,2)& ...& (n-1,n-1)& \end{array}\right]\cdot f^{n+1}=f^n$$

然后自己递推到三维的话，就会以为要用到三维矩阵，比如把一些二维矩阵堆叠起来的这种形式……

比如系数矩阵的第一层是

$$\left[\begin{array}{cccccc}(0,0,0)& (0,0,1)& (0,0,2)& ...& (0,0,n-1)& \\ (0,1,0)& (0,1,1)& (0,1,2)& ...& (0,1,n-1)& \\ ⋮& ⋮& \ddots & ⋮& ⋮& \\ (0,n-1,0)& (0,n-1,1)& (0,n-1,2)& ...& (0,n-1,n-1)& \end{array}\right]$$

第二层是

[ ( 1 , 0 , 0 ) ( 1 , 0 , 1 ) ( 1 , 0 , 2 ) . . . ( 1 , 0 , n − 1 ) ( 1 , 1 , 0 ) ( 1 , 1 , 1 ) ( 1 , 1 , 2 ) . . . ( 1 , 1 , n −