从GPU流体力学算例中理解平流有困难

Question

因此，我正在尝试实现流体模拟，如GPU宝石第38章中所描述的。

我还大量地借鉴了帕维尔的出色实现。

我被困住的地方是有意见的。

所以据我所知，平流是随着流体运动而移动的过程。例如，如果在一个时间戳上，在位置p有一个速度D3，那么在下一个时间戳中，在p + v * dt位置有速度V，至少根据平流。

因此，正如它在GPU中所说的，在GPU上实现此逻辑存在困难：

我们在片段程序中实现了我们的模拟，它不能改变他们正在编写的片段的位置。这种前向集成方法需要能够“移动”粒子，因此不能在当前的GPU上实现。解决方法是将问题转化为隐式方法(Stam，1999)。我们没有通过计算粒子在当前时间步长上的移动来提升粒子的数量，而是将粒子的轨迹从每个网格单元追溯到以前的位置，然后将该位置的量复制到起始网格单元。

换句话说，我们不把数值投影到网格中，而是根据当前时间戳中的速度值，用当前的速度在网格中追溯，看看下一个值应该是什么，它是向我们当前网格单元的方向来的。

这种方法的问题似乎发生在速度为零时。例如，假设我们有一个一维的速度网格：

[3|1|2|1]

如果我们在下一个时间戳中根据反向查找来推进这个网格，我们将得到以下内容：

[0|3|3|2]

因此，这是正常工作的预期。

无论如何，如果我们在任何位置上都有一个零，那么就无法获得另一个值来替换零：

t+0: [3|1|2|0]

t+1: [0|3|3|0]

t+2: [0|0|0|0]

因此，基本上，如果单元格的值为零，它将永远不会查找除自身之外的任何值，因此，没有其他量能进入它的方法。

事实上，这正是我在模拟中所看到的行为：

这是我使用的着色器：

#define WORKGROUP_SIZE 16
layout (local_size_x = WORKGROUP_SIZE, local_size_y = WORKGROUP_SIZE, local_size_z = 1 ) in;

precision highp float;
precision highp sampler2D;

layout(push_constant) uniform PushConstant {
    vec2 viewport;
    vec2 texelSize;
    vec2 dyeTexelSize;
    float dt;
    float dissipation;
} pushConstants;


layout(binding = 0, rgba8) uniform writeonly image2D outputImage;


layout(binding = 1) uniform sampler2D uVelocity;
layout(binding = 2) uniform sampler2D uSource;


vec4 bilerp(sampler2D sam, vec2 uv, vec2 tsize) {
    vec2 st = uv / tsize - 0.5;
    vec2 iuv = floor(st);
    vec2 fuv = fract(st);
    vec4 a = texture(sam, (iuv + vec2(0.5, 0.5)) * tsize);
    vec4 b = texture(sam, (iuv + vec2(1.5, 0.5)) * tsize);
    vec4 c = texture(sam, (iuv + vec2(0.5, 1.5)) * tsize);
    vec4 d = texture(sam, (iuv + vec2(1.5, 1.5)) * tsize);
    return mix(mix(a, b, fuv.x), mix(c, d, fuv.x), fuv.y);
}

void main() {
    
    vec2 viewport = pushConstants.viewport;
    vec2 texelSize = pushConstants.texelSize;
    vec2 dyeTexelSize = pushConstants.dyeTexelSize;
    float dt = pushConstants.dt;
    float dissipation = pushConstants.dissipation;
    vec2 position = vec2(gl_GlobalInvocationID.x / viewport.x, gl_GlobalInvocationID.y / viewport.y);
    if (position.x > 1.0 || position.y > 1.0) {
        return;
    }
    
    vec2 uv = position + texelSize/2;
    vec2 vUv = params.uv;
    vec2 coord = vUv - dt * bilerp(uVelocity, uv, texelSize).xy * texelSize;
    vec4 result = bilerp(uSource, coord, dyeTexelSize);
    float decay = 1.0 + dissipation * dt;
    imageStore(outputImage, ivec2(gl_GlobalInvocationID.xy), result / decay);
}

所以我的问题是，怎样才能把数量提升到一个空值呢？我是不是错过了这个过程的一部分？

Answer 1

零速度单元是模拟域的边界(壁)。您需要显式地处理边界条件。平流本身对你没有多大帮助。

从上下文来看，您似乎是在模拟不可压缩的流。直观地说，如果流体流向墙壁，那么流体就会被“推”到墙壁的两侧。这是由压力投影处理的。具体来说，它围绕着你链接的页面上的等式16。

方程16是压力项的离散化。它将引导你进入一个线性系统。您可以迭代地解决它(Stam99)。这个系统的解会给你自由发散的速度，它会像预期的那样把你的流体推进到壁面的两侧。

关于你的问题:如果边界是一道墙，那么任何东西都不会进入它。不可压缩性(这里的压力)将强制流体流动到其他地方。