1

我正在使用 Electron 和 TypeScript 对一些流体模拟代码进行原型设计，使用 Lattice Boltzmann 算法，最终将进入游戏。到目前为止，我一直在使用静态边界条件（模拟计算仅发生在网格内部，并且边界单元的值保持固定），在该状态下一切似乎都运行良好。特别是，我可以通过手动设置每个之间的单元格值来施加内部边界条件（例如，强制一定密度的流体总是在每一帧上离开特定晶格位置，以模拟软管/火箭喷嘴/其他）模拟步骤。

但是，如果我改用周期性边界条件（即环绕式、环形拓扑的小行星世界），整个模拟就会变成静态的。我只是随时随地获得恒定的流体密度，就像我所有的边界条件都被删除了，无论在模拟循环中的哪个位置（流式传输之前或碰撞之前）我选择断言它们。我不确定周期性边界条件最终是否会与游戏相关，但这次失败让我认为模拟中的某个地方一定存在一些微妙的错误。

完整的代码可在https://github.com/gliese1337/balloon-prototype/tree/deopt获得，但我期望相关部分如下：

    class LatticeBoltzmann {        
      private streamed: Float32Array; // microscopic densities along each lattice direction
      private collided: Float32Array;
      public rho: Float32Array;    // macroscopic density; cached for rendering

      ...

      public stream(barriers: boolean[]) {
        const { xdim, ydim, collided, streamed } = this;
        const index = (x: number, y: number) => (x%xdim)+(y%ydim)*xdim;
        const cIndex = (x: number, y: number, s: -1|1, j: number) =>
                          9*(((x+s*cxs[j])%xdim)+((y+s*cys[j])%ydim)*xdim)+j;

        // Move particles along their directions of motion:
        for (let y=1; y<ydim-1; y++) {
          for (let x=1; x<xdim-1; x++) {
            const i = index(x, y);
            const i9 = i*9;
            for (let j=0;j<9;j++) {
              streamed[i9 + j] = collided[cIndex(x, y, -1, j)];
            }
          }
        }

        // Handle bounce-back from barriers
        for (let y=0; y<ydim; y++) {
          for (let x=0; x<xdim; x++) {
            const i = index(x, y);
            const i9 = i*9;
            if (barriers[i]) {
              for (let j=1;j<9;j++) {
                streamed[cIndex(x, y, 1, j)] = collided[i9 + opp[j]];
              }
            }
          }
        }
      }

      // Set all densities in a cell to their equilibrium values for a given velocity and density:
      public setEquilibrium(x: number, y: number, ux: number, uy: number, rho: number) {
        const { xdim, streamed } = this;
        const i = x + y*xdim;
        this.rho[i] = rho;

        const i9 = i*9;
        const u2 =  1 - 1.5 * (ux * ux + uy * uy);
        for (let j = 0; j < 9; j++) {
          const dir = cxs[j]*ux + cys[j]*uy;
          streamed[i9+j] =  weights[j] * rho * (u2 + 3 * dir + 4.5 * dir * dir);
        }
      }
    }

晶格数据存储在两个平面数组中，collided它们保存碰撞步骤后的最终状态并用作流步骤的输入streamed，而 保存流步骤之后的最终状态并用作下一个碰撞步骤的输入。D2Q9 晶格的 9 个向量分量存储在连续的块中，然后分组为行。请注意，我已经在使用 mod 操作从晶格坐标计算数组索引；只要模拟计算仅在晶格内部范围内，这完全无关紧要，但是一旦for (let y=1; y<ydim-1; y++)andfor (let x=1; x<xdim-1; x++)循环将其边界更改为for (let y=0; y<ydim; y++)and ，它应该使周期性边界准备就绪for (let x=0; x<xdim; x++)， 分别。事实上，我遇到的问题是特定的 6 个字符的变化。

该setEquilibrium方法用于施加边界条件。在驱动程序代码中，它当前被这样调用，每帧一次：

    // Make fluid flow in from the left edge and out through the right edge
    function setBoundaries(LB: LatticeBoltzmann, ux: number) {
        for (let y=0; y<ydim; y++) {
            LB.setEquilibrium(0, y, ux, 0, 1);
            LB.setEquilibrium(xdim-1, y, ux, 0, 1);
        }
    }

对于静态边界条件，每帧调用一次恰好是多余的，因为它只会改变边界晶格位置。然而，将硬编码的 x 值移动到晶格的内部（实际上每帧重新设置一次边界条件是必要的）完全符合您的预期——它使流体在特定位置出现或消失。然而，切换到周期性边界条件会导致该代码不再具有任何可见的效果。

所以...有人知道我可能做错了什么吗？

4

1 回答 1

0

我不完全确定为什么这个特殊错误会产生这种特别奇怪的效果，但事实证明问题出在我对%运算符的使用上——它已签名。因此，当输入一个负格索引时，天真的使用%不会执行一个适当的模运算符想要的环绕；相反，它只是给您返回相同的负值，并导致超出范围的数组访问。

在取余数之前添加数组维度可确保所有值都是正数，并且我们获得了必要的环绕行为。

顺便说一句，能够在整个晶格上进行范围而不需要特别处理边缘允许将嵌套循环折叠为整个晶格上的单个线性扫描，这消除了对主索引计算函数的需要，并极大地简化了碰撞流偏移索引函数 ，cIndex现在看起来像const cIndex = (i: number, s: -1|1, j: number) => 9*((i+s*(cxs[j]+cys[j]*xdim)+max)%max)+j;，只需要一个模数而不是每个维度一个。这一系列简化的结果是代码的巨大加速，以及相关的帧速率提高。

于 2019-10-29T04:55:17.557 回答