我正在使用带有TypeScript的Electron来使用Lattice Boltzmann算法对一些流体模拟代码进行原型设计,该模型最终将进入游戏。到目前为止,我一直在使用静态边界条件(仅在网格内部进行模拟计算,并且边界单元的值保持固定),在这种情况下,一切似乎都可以正常工作。特别是,我可以通过手动设置每个单元格之间的单元格值来施加内部边界条件(例如,强制在每帧上始终以一定密度的流体离开某个帧上的某个晶格位置,以模拟软管/火箭喷嘴/任何物体)。模拟步骤。
但是,如果我转为使用周期性边界条件(即环绕的,环形拓扑Asteroids世界),则整个模拟将变为静态。我一直到处都获得恒定的流体密度,这就像擦除所有边界条件一样,无论在模拟周期中的任何位置(在流之前还是在碰撞之前),我都选择声明它们。我不确定周期性边界条件是否最终会与游戏相关,但是这种失败使我认为模拟中的某处一定存在一些细微的错误。
完整的代码可在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++)
和for (let x=1; x<xdim-1; x++)
循环的边界更改为for (let y=0; y<ydim; y++)
,它应该使周期性边界随时可用。和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值移动到晶格的内部(实际上必须每帧重新设置边界条件)确实符合您的期望-会使流体在特定位置出现或消失。但是,切换到周期性边界条件会导致该代码不再具有任何可见的效果。
所以...有人知道我可能做错了吗?
我不完全确定为什么这个特定的错误会产生这种奇怪的效果,但是事实证明,问题出在我使用%
运算符时-已签名。因此,当输入负晶格索引时,%
的天真用法不会执行您希望从适当的模运算符得到的环绕操作;相反,它只给您返回相同的负值,并导致超出范围的数组访问。
[在获取余数之前添加数组维数可确保所有值均为正,并且我们获得了必要的环绕行为。
顺便说一句,能够在不打扰边缘的情况下在整个晶格范围内进行调整,特别允许将嵌套循环折叠为整个晶格的单个线性扫描,这消除了对主索引计算功能的需求,并极大地简化了碰撞-流偏移索引函数cIndex
,现在看起来像const cIndex = (i: number, s: -1|1, j: number) => 9*((i+s*(cxs[j]+cys[j]*xdim)+max)%max)+j;
,只需要一个模数,而不是每个维一个。简化字符串的结果是极大地提高了代码速度,并提高了帧速率。