javascript node.js canvas simulation fluid-dynamics

What is causing this 2D fluid simulation to hang?

I attempted to build a fluid simulator based off the architecture of Jos Stam's book, The Art of Fluid Animation, only I'm using node.js and canvas-sketch. I think I can finagle it to work, it's just that when I run it as shown below, it hangs. No errors, but it only logs one bout of 'Rolling', the output in the animation loop. I think there may be something funny with my array declarations, I've tried spacing them out to different lines just in case it was making copies of itself instead of unique arrays, but no change. I was getting errors before I called the new Array functions in the sketch function.

Size doesn't seem to be the problem, since it still hangs when NUM is super low. The paper passes a lot of references, I thought I could get around it by using globally defined variables and using a lot of void functions to change their values. Is this a flawed approach? Any suggestions?

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

Code directly from the book (meant for c)... Edited with indents:

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

Pseudo code given in the paper, Real-Time Fluid Dynamics for Games:

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

Solution

Okay, a keen eye on discord brought to my attention that I had put an i in a nested for loop for j in two functions. Editing that problem does not create the desired effects, but now the canvas is shown and it logs its animation frames. but the original question of "Why is it Hanging?" has been solved. I won't edit out my dumdums in the OP, but I will mark them with comments.