是什么导致此 2D 流体模拟挂起？

Question

我试图根据 Jos Stam 的书流体动画艺术 的架构构建一个流体模拟器，只是我使用的是 node.js 和 canvas-素描。我想我可以使它工作，只是当我运行它如下所示时，它会挂起。没有错误，但它只记录了一轮 'Rolling'，动画循环中的输出。我认为我的数组声明可能有些有趣，我尝试将它们分隔到不同的行以防万一它正在复制自己而不是唯一的数组，但没有改变。在草图函数中调用新的 Array 函数之前，我遇到了错误。

大小似乎不是问题，因为当 NUM 非常低时它仍然挂起。该论文通过了很多参考资料，我认为我可以通过使用全局定义的变量并使用大量 void 函数来更改它们的值来解决它。这是一种有缺陷的方法吗？有什么建议吗？

const canvasSketch = require('canvas-sketch');

const NUM=1078;
const SIZE=(NUM+2)*(NUM+2);
const DT = 0.1, DIFF = 0.0, VIS = 0.0;

let U, V, uPREV, vPREV;
let DENS, dPREV;

//Call these macros, okay?
function IX(i, j) {return (i + (NUM + 2) * j);}
function SWAP(x0,x) {let tmp = x0; x0 = x; x = tmp;}

const settings = {
  dimensions: [NUM+2, NUM+2],
  animate:true
};

const animate=()=>{
  requestAnimationFrame(animate); //forgot that
}; //forgot this too

const sketch = () => {
  U = new Array(SIZE);
  V = new Array(SIZE);
  uPREV = new Array(SIZE);
  vPREV = new Array(SIZE);
  DENS = new Array(SIZE);
  dPREV = new Array(SIZE);

  console.log('Init');

  return ({ context, width, height }) => {
    console.log('Rolling'); //only logged once
    context.fillStyle = 'white';
    context.fillRect(0, 0, width, height);

    velStep(NUM, U, V, uPREV, vPREV, VIS, DT);
    densStep(NUM, DENS, dPREV, U, V, DIFF, DT);

    drawDens(context);
    drawVel(context);
  };
};

canvasSketch(sketch, settings);

function addSource(N, x, s, dt){
  for (let i=0; i<SIZE; i++){
    x[i] += dt * s[i];
  }
}

function setBound(N, b, x){
  for (let i=1; i<=N; i++){
    x[IX(0, i)]=b===1 ? -x[IX(1, i)] : x[IX(1, i)];
    x[IX(N+1, i)]=b===1 ? -x[IX(N, i)] : x[IX(N, i)];
    x[IX(i, 0)]=b===2 ? -x[IX(i, 1)] : x[IX(i, 1)];
    x[IX(i, N+1)]=b===2 ? -x[IX(i, N)] : x[IX(i, N)];
  }
  x[IX(0,0)] = 0.5 * (x[IX(1,0)] + x[IX(0 ,1)]);
  x[IX(0,N+1)] = 0.5 * (x[IX(1,N+1)] + x[IX(0 ,N)]);
  x[IX(N+1,0)] = 0.5 * (x[IX(N,0)] + x[IX(N+1 ,1)]);
  x[IX(N+1,N+1)] = 0.5 * (x[IX(N,N+1)] + x[IX(N+1 ,N)]);
}

function linSolve(N, b, x, x0, a, c){
  for (let n= 0; n< 20; n++){
    for (let j=1; j<=N; j++){
      for (let i=1;i<=N; i++){
        x[IX(i, j)] = (x0[IX(i, j)] + a * (x[IX(i-1, j)] + x[IX(i + 1, j)] + x[IX(i,j - 1)] + x[IX(i,j + 1)])) / c;
      }
    }
    setBound(N,b,x);
  }
}

function diffuse(N, b, x, x0, diff, dt){
  let a = dt * diff * N * N;
  linSolve(N, b, x, x0, a, 1 + 4 * a);
}

function advect(N, b, d, d0, u, v, dt){
  let i0, j0, i1, j1;
  let x, y, s0, t0, s1, t1, dt0;
  dt0 = dt*N;
  for (let j=1; j<=N; j++) {
    for (let i = 1; i <= N; i++) {
      x = i - dt0 * u[IX(i,j)];
      y = j - dt0 * v[IX(i,j)];

      if (x<0.5){x=0.5;}
      if (x>N+0.5){x=N+0.5;}
      i0=x; i1 = i0 + 1;

      if (y<0.5){y=0.5;}
      if (y>N+0.5){y=N+0.5;}
      j0=x; j1 = j0 + 1;

      s1 = x - i0;
      s0 = 1 - s1;
      t1 = y - j0;
      t0 = 1 - t1;

      d[IX(i, j)] = s0 * (t0 * d0[IX(i0, j0)] + t1 * d0[IX(i0, j1)]) + s1 * (t0 * d0[IX(i1, j0)] + t1 * d0[IX(i1, j1)]);
    }
  }
  setBound(N ,b, d);
}

function project(N, u, v, p, div){
  for (let j=1; j<=N; j++) {
    for (let i = 1; i <= N; i++) {
      div[IX(i, j)] = - 0.5 * (u[IX(i + 1, j)] - u[IX(i - 1, j)] + v[IX(i, j + 1)] - v[IX(i, j - 1)]) / N;
      p[IX(i, j)] = 0;
    }
  }
  setBound(N, 0, div);
  setBound(N, 0, p);

  linSolve(N, 0, p, div, 1, 4);

  for (let j=1; j<=N; j++) {
    for (let i = 1; i <= N; i++) {
      u[IX(i, j)] -= 0.5 * N * (p[IX(i + 1, j)] - p[IX(i - 1, j)]);
      v[IX(i, j)] -= 0.5 * N * (p[IX(i, j + 1)] - p[IX(i, j - 1)]);
    }
  }
  setBound(N, 1, u);
  setBound(N, 2, v);
}

function densStep(N, x, x0, u, v, diff, dt){
  addSource(N, x, x0, dt);
  SWAP (x0, x);
  diffuse (N, 0, x, x0, diff, dt);
  SWAP (x0, x);
  advect (N, 0, x, x0, u, v, dt);
}

function velStep(N, u, v, u0, v0, vis, dt){
  addSource(N, u, u0, dt);
  addSource(N, v, v0, dt);

  SWAP (u0, u);
  diffuse(N, 1, u, u0, vis, dt);
  SWAP (v0, v);
  diffuse(N, 2, v, v0, vis, dt);

  project(N, u, v, u0, v0);
  SWAP (u0, u);
  SWAP (v0, v);

  advect(N, 1, u, u0, u0, v0, dt);
  advect(N, 2, v, v0, u0, v0, dt);

  project(N, u, v, u0, v0);
}

function drawVel(context){
  let x,y,h;
  h = 1.0/NUM;

  context.save();
  for (let i=1; i<=NUM; i++){
    x = (i - 0.5) * h;
    for (let j=1; j<=NUM; i++){ // BOOM, mistake right here caused i to iterate double and j not at all
      y = (j - 0.5) * h;

      context.beginPath();
      context.moveTo(x,y);
      context.lineTo(x + U[IX(i, j)], y + V[IX(i, j)]);
      context.strokeStyle='black';
      context.stroke();
    }
  }

  context.restore();
}

function drawDens(context){
  let x, y, h, d00, d01, d10, d11;

  h =1.0/NUM;

  context.save();
  for (let i=1; i<=NUM; i++) {
    x = (i - 0.5) * h;
    for (let j = 1; j <= NUM; i++) { // AGAIN, right here. Proofread before you post :\
      y = (j - 0.5) * h;

      d00 = DENS[IX(i, j)];
      d01 = DENS[IX(i, j + 1)];
      d10 = DENS[IX(i + 1, j)];
      d11 = DENS[IX(i + 1, j + 1)];

      context.beginPath();
      context.fillStyle = 'blue';
      context.arc(x, y, d00, 0, Math.PI * 2);
      context.arc(x + h, y, d01, 0, Math.PI * 2);
      context.arc(x + h, y + h, d10, 0, Math.PI * 2);
      context.arc(x, y + h, d11, 0, Math.PI * 2);
      context.fill();
    }
  }
  context.restore();
}

直接来自本书的代码（适用于 c）...缩进编辑：

#define IX(i,j) ((i)+(N+2)*(j))
#define SWAP(x0,x) {float * tmp=x0;x0=x;x=tmp;}
#define FOR_EACH_CELL for ( j=1 ; j<=N ; j++ ) {\
for ( i=1 ; i<=N ; i++ ) {
#define END_FOR }}
void add_source(int N, float *x, float *s, float dt)
{
  int i, size=(N+2)*(N+2);
  for ( i=0 ; i<size ; i++ ) x[i] += dt*s[i];
}
void set_bnd(int N, int b, float *x)
{
  int i;
  for ( i=1 ; i<=N ; i++ ) {
    x[IX(0 ,i)] = b==1 ? -x[IX(1,i)] : x[IX(1,i)];
    x[IX(N+1,i)] = b==1 ? -x[IX(N,i)] : x[IX(N,i)];
    x[IX(i,0 )] = b==2 ? -x[IX(i,1)] : x[IX(i,1)];
    x[IX(i,N+1)] = b==2 ? -x[IX(i,N)] : x[IX(i,N)];
  }
  x[IX(0 ,0 )] = 0.5f*(x[IX(1,0 )]+x[IX(0 ,1)]);
  x[IX(0 ,N+1)] = 0.5f*(x[IX(1,N+1)]+x[IX(0 ,N)]);
  x[IX(N+1,0 )] = 0.5f*(x[IX(N,0 )]+x[IX(N+1,1)]);
  x[IX(N+1,N+1)] = 0.5f*(x[IX(N,N+1)]+x[IX(N+1,N)]);
}
void lin_solve(int N, int b, float *x, float *x0, float a, float c)
{
  int i, j, n;
  for ( n=0 ; n<20 ; n++ ) {
    FOR_EACH_CELL
      x[IX(i,j)] = (x0[IX(i,j)]+a*(x[IX(i-1,j)]+
      x[IX(i+1,j)]+x[IX(i,j-1)]+x[IX(i,j+1)]))/c;
    END_FOR
    set_bnd ( N, b, x );
  }
}
void diffuse(int N, int b, float *x, float *x0, float diff, float dt)
{
  float a=dt*diff*N*N;
  lin_solve ( N, b, x, x0, a, 1+4*a );
}
  void advect(int N, int b, float *d, float *d0, float *u, float *v, float dt)
{
  int i, j, i0, j0, i1, j1;
  float x, y, s0, t0, s1, t1, dt0;
  dt0 = dt*N;
  FOR_EACH_CELL
    x = i-dt0*u[IX(i,j)]; y = j-dt0*v[IX(i,j)];
    if (x<0.5f) x=0.5f; if (x>N+0.5f) x=N+0.5f; i0=(int)x; i1=i0+1;
    if (y<0.5f) y=0.5f; if (y>N+0.5f) y=N+0.5f; j0=(int)y; j1=j0+1;
    s1 = x-i0; s0 = 1-s1; t1 = y-j0; t0 = 1-t1;
    d[IX(i,j)] = s0*(t0*d0[IX(i0,j0)]+t1*d0[IX(i0,j1)])+
    s1*(t0*d0[IX(i1,j0)]+t1*d0[IX(i1,j1)]);
  END_FOR
  set_bnd ( N, b, d );
}
void project(int N, float * u, float * v, float * p, float * div)
{
  int i, j;
  FOR_EACH_CELL
    div[IX(i,j)] = -0.5f*(u[IX(i+1,j)]-u[IX(i-1,j)]+
    v[IX(i,j+1)]-v[IX(i,j-1)])/N;
    p[IX(i,j)] = 0;
  END_FOR
  set_bnd ( N, 0, div ); set_bnd ( N, 0, p );
  lin_solve ( N, 0, p, div, 1, 4 );
  FOR_EACH_CELL
    u[IX(i,j)] -= 0.5f*N*(p[IX(i+1,j)]-p[IX(i-1,j)]);
    v[IX(i,j)] -= 0.5f*N*(p[IX(i,j+1)]-p[IX(i,j-1)]);
  END_FOR
  set_bnd ( N, 1, u ); set_bnd ( N, 2, v );
}
void dens_step(int N, float *x, float *x0, float *u, float *v, float diff, float dt)
{
  add_source ( N, x, x0, dt );
  SWAP ( x0, x ); diffuse ( N, 0, x, x0, diff, dt );
  SWAP ( x0, x ); advect ( N, 0, x, x0, u, v, dt );
}
void vel_step(int N, float *u, float *v, float *u0, float *v0, float visc, float dt)
{
  add_source ( N, u, u0, dt ); add_source ( N, v, v0, dt );
  SWAP ( u0, u ); diffuse ( N, 1, u, u0, visc, dt );
  SWAP ( v0, v ); diffuse ( N, 2, v, v0, visc, dt );
  project ( N, u, v, u0, v0 );
  SWAP ( u0, u ); SWAP ( v0, v );
  advect ( N, 1, u, u0, u0, v0, dt ); advect ( N, 2, v, v0, u0, v0, dt );
  project ( N, u, v, u0, v0 );
}

论文中给出的伪代码，游戏的实时流体动力学：

while ( simulating )
{
  //get_from_UI ( dens_prev, u_prev, v_prev ); //mouse input, can be ignored for the moment
  velStep ( N, U, V, uPREV, vPREV, VIS, DT );
  densStep ( N, DENS, dPREV, U, V, DIFF, DT );
  drawDens ( N, DENS );
}

Answer 1

好吧，对 discord 的敏锐观察让我注意到我在两个函数的 j 的嵌套 for 循环中放置了一个 i。编辑该问题不会产生预期的效果，但现在会显示 canvas 并记录其动画帧。但最初的问题是“为什么会挂起？”已经解决了。我不会在 OP 中删除我的 dumdums，但我会用评论标记它们。

是什么导致此 2D 流体模拟挂起？

What is causing this 2D fluid simulation to hang?

javascript

simulation

canvas

node.js

fluid-dynamics