ブラウザで流体シミュレーション — 100×60グリッドで実現するリアルタイム演算

const N = 100
const M = 60
const size = (N + 2) * (M + 2) // 境界セル含む

const density = new Float32Array(size)
const velocityX = new Float32Array(size)
const velocityY = new Float32Array(size)

function diffuse(
  b: number,
  x: Float32Array,
  x0: Float32Array,
  diff: number,
  dt: number,
) {
  const a = dt * diff * N * M
  for (let k = 0; k < 20; k++) {
    for (let j = 1; j <= M; j++) {
      for (let i = 1; i <= N; i++) {
        const idx = i + (N + 2) * j
        x[idx] =
          (x0[idx] + a * (x[idx - 1] + x[idx + 1] +
            x[idx - (N + 2)] + x[idx + (N + 2)])) /
          (1 + 4 * a)
      }
    }
    setBoundary(b, x)
  }
}

function advect(
  b: number,
  d: Float32Array,
  d0: Float32Array,
  u: Float32Array,
  v: Float32Array,
  dt: number,
) {
  for (let j = 1; j <= M; j++) {
    for (let i = 1; i <= N; i++) {
      const idx = i + (N + 2) * j

      // 逆追跡
      let x = i - dt * N * u[idx]
      let y = j - dt * M * v[idx]

      // クランプ
      x = Math.max(0.5, Math.min(N + 0.5, x))
      y = Math.max(0.5, Math.min(M + 0.5, y))

      // バイリニア補間
      const i0 = Math.floor(x)
      const j0 = Math.floor(y)
      const s1 = x - i0
      const t1 = y - j0
      const s0 = 1 - s1
      const t0 = 1 - t1

      d[idx] =
        s0 * (t0 * d0[i0 + (N + 2) * j0] + t1 * d0[i0 + (N + 2) * (j0 + 1)]) +
        s1 * (t0 * d0[i0 + 1 + (N + 2) * j0] + t1 * d0[i0 + 1 + (N + 2) * (j0 + 1)])
    }
  }
  setBoundary(b, d)
}

function project(
  u: Float32Array,
  v: Float32Array,
  p: Float32Array,
  div: Float32Array,
) {
  for (let j = 1; j <= M; j++) {
    for (let i = 1; i <= N; i++) {
      const idx = i + (N + 2) * j
      div[idx] = -0.5 * (
        (u[idx + 1] - u[idx - 1]) / N +
        (v[idx + (N + 2)] - v[idx - (N + 2)]) / M
      )
      p[idx] = 0
    }
  }
  setBoundary(0, div)
  setBoundary(0, p)

  // ポアソン方程式をガウス・ザイデルで解く
  for (let k = 0; k < 20; k++) {
    for (let j = 1; j <= M; j++) {
      for (let i = 1; i <= N; i++) {
        const idx = i + (N + 2) * j
        p[idx] = (div[idx] + p[idx - 1] + p[idx + 1] +
          p[idx - (N + 2)] + p[idx + (N + 2)]) / 4
      }
    }
    setBoundary(0, p)
  }

  // 速度場から勾配を引く
  for (let j = 1; j <= M; j++) {
    for (let i = 1; i <= N; i++) {
      const idx = i + (N + 2) * j
      u[idx] -= 0.5 * N * (p[idx + 1] - p[idx - 1])
      v[idx] -= 0.5 * M * (p[idx + (N + 2)] - p[idx - (N + 2)])
    }
  }
  setBoundary(1, u)
  setBoundary(2, v)
}

function renderFluid(
  ctx: CanvasRenderingContext2D,
  density: Float32Array,
  width: number,
  height: number,
) {
  const imageData = ctx.createImageData(width, height)
  const cellW = width / N
  const cellH = height / M

  for (let j = 1; j <= M; j++) {
    for (let i = 1; i <= N; i++) {
      const d = Math.min(255, density[i + (N + 2) * j] * 255)
      const px = Math.floor((i - 1) * cellW)
      const py = Math.floor((j - 1) * cellH)
      // セル内のピクセルを塗る（簡略版）
      const idx = (py * width + px) * 4
      imageData.data[idx] = d * 0.3     // R
      imageData.data[idx + 1] = d * 0.6 // G
      imageData.data[idx + 2] = d        // B
      imageData.data[idx + 3] = 255      // A
    }
  }
  ctx.putImageData(imageData, 0, 0)
}