Explore real-time fluid solvers on the GPU using Smoothed Particle Hydrodynamics (SPH). Learn how to implement multi-pass compute pipelines for density, pressure, and viscosity force calculations, and understand spatial grid partition concepts.

This lesson is currently only available in Vietnamese. Please switch the language toggle in the menu to Vietnamese to read the full guide and take the interactive quiz.

Trong bài học trước, chúng ta đã xây dựng một hệ thống hạt độc lập tự do chịu tác động cơ bản của trọng lực. Để mô phỏng chất lỏng (như dòng nước chảy, dầu loang), các hạt không thể chỉ chuyển động riêng rẽ mà bắt buộc phải có sự tương tác vật lý liên tục với nhau để duy trì thể tích không đổi của khối chất lỏng.

Bài học này sẽ hướng dẫn bạn thuật toán SPH (Smoothed Particle Hydrodynamics) – giải pháp mô phỏng chất lỏng dựa trên hạt Lagrange phổ biến nhất trong đồ họa thời gian thực, cách tổ chức hai lượt Compute Shader lồng nhau (Multi-pass Compute) và tối ưu hóa phân vùng lưới tìm lân cận.

1. Giải thuật Smoothed Particle Hydrodynamics (SPH)

Ý tưởng cốt lõi của SPH là coi chất lỏng là một tập hợp các hạt rời rạc. Mỗi hạt mang một khối lượng nhất định, các đại lượng vật lý (như mật độ, áp suất) tại một điểm bất kỳ trong không gian sẽ được nội suy trơn tru từ các hạt xung quanh thông qua một hàm hạt nhân làm mịn (Smoothing Kernel Function - $W(r, h)$).

Khoảng cách ảnh hưởng tối đa của hạt nhân làm mịn được gọi là bán kính làm mịn ($h$). Nếu khoảng cách $r$ giữa hai hạt vượt quá $h$, tương tác vật lý giữa chúng bằng 0.

2. Thiết lập Multi-pass Compute Pipeline

Vì bước 2 (tính toán lực) yêu cầu giá trị áp suất của toàn bộ các hạt đã được tính toán xong ở bước 1, ta không thể viết chung cả hai bước này trong cùng một Compute Shader duy nhất nếu không có cơ chế đồng bộ hóa toàn cục.

Giải pháp là thiết lập lượt chạy Compute đa pass (Multi-pass Compute) trong mã JavaScript điều phối:

const commandEncoder = device.createCommandEncoder();

// Lượt 1: Chạy pipeline tính mật độ và áp suất
const pass1 = commandEncoder.beginComputePass();
pass1.setPipeline(densityPipeline);
pass1.setBindGroup(0, computeBindGroup);
pass1.dispatchWorkgroups(Math.ceil(particleCount / 64));
pass1.end();

// Lượt 2: Chạy pipeline tính lực và cập nhật tọa độ
const pass2 = commandEncoder.beginComputePass();
pass2.setPipeline(forcePipeline);
pass2.setBindGroup(0, computeBindGroup);
pass2.dispatchWorkgroups(Math.ceil(particleCount / 64));
pass2.end();

3. Khái niệm Spatial Hashing tối ưu hóa $O(N^2) \to O(N)$

Nếu thực hiện duyệt trực tiếp vòng lặp lồng nhau (Double loop), mỗi hạt phải kiểm tra khoảng cách với $N-1$ hạt còn lại, tạo ra độ phức tạp tính toán rất lớn là \(O(N^2)\). Kỹ thuật Spatial Hashing (Phân vùng lưới không gian) chia toàn bộ không gian mô phỏng thành một lưới các ô có kích thước bằng bán kính làm mịn $h$.

Mỗi hạt được ánh xạ (hash) vào một ô lưới dựa trên tọa độ vị trí. Khi cần tính toán mật độ hay lực của hạt $i$, ta chỉ cần duyệt qua các hạt nằm trong cùng ô lưới và 8 ô lưới lân cận xung quanh, giúp giảm số lượng phép tính khoảng cách thừa và hạ độ phức tạp trung bình xuống còn \(O(N)\).

Sân chơi tương tác: GPU SPH Fluid Lab

Dưới đây là mô phỏng động lực học chất lỏng SPH chạy thực tế trên WebGPU. Hãy nhấn chuột kéo rê trực tiếp trên màn hình nước để tạo ra các cơn sóng, và thử thay đổi hệ số độ cứng (áp suất) hoặc độ nhớt để quan sát sự biến chuyển trạng thái từ nước lỏng sang dạng keo dẻo:

struct Particle {
  pos : vec2<f32>,
  vel : vec2<f32>,
  force : vec2<f32>,
  density : f32,
  pressure : f32,
};

struct Params {
  gravity : f32,
  viscosity : f32,
  pressureCoeff : f32,
  h : f32,
  targetDensity : f32,
  particleRadius : f32,
  deltaTime : f32,
  mousePos : vec2<f32>,
  mouseActive : u32,
};

@group(0) @binding(0) var<uniform> u : Params;
@group(0) @binding(1) var<storage, read_write> particles : array<Particle>;

// LƯỢT 1 COMPUTE: Tính mật độ và áp suất của từng hạt
@compute @workgroup_size(64)
fn calcDensity(@builtin(global_invocation_id) id : vec3<u32>) {
  let idx = id.x;
  if (idx >= arrayLength(&particles)) { return; }

  var p = particles[idx];
  var density = 0.0;
  let h2 = u.h * u.h;

  for (var i = 0u; i < arrayLength(&particles); i = i + 1u) {
    let p_j = particles[i];
    let diff = p.pos - p_j.pos;
    let r2 = dot(diff, diff);

    if (r2 < h2) {
      // Smoothing kernel polynomial
      let term = h2 - r2;
      density = density + term * term * term;
    }
  }

  // Chuẩn hóa mật độ dựa trên h9 kernel constant
  p.density = density * (315.0 / (64.0 * 3.141592 * pow(u.h, 9.0)));
  // Tait equation tính áp suất nước chống nén
  p.pressure = u.pressureCoeff * (p.density - u.targetDensity);

  particles[idx] = p;
}

// LƯỢT 2 COMPUTE: Tính lực tương tác và tích hợp tích lũy vị trí
@compute @workgroup_size(64)
fn calcForces(@builtin(global_invocation_id) id : vec3<u32>) {
  let idx = id.x;
  if (idx >= arrayLength(&particles)) { return; }

  var p = particles[idx];
  var force = vec2<f32>(0.0, -u.gravity * p.density);
  let h = u.h;

  for (var i = 0u; i < arrayLength(&particles); i = i + 1u) {
    if (i == idx) { continue; }
    let p_j = particles[i];
    let diff = p.pos - p_j.pos;
    let r = length(diff);

    if (r < h && r > 0.001) {
      let dir = diff / r;
      // Lực áp suất đẩy hạt ra xa nhau (Pressure Force)
      let pForceVal = -(p.pressure + p_j.pressure) / (2.0 * p_j.density) * (45.0 / (3.141592 * pow(h, 6.0))) * pow(h - r, 2.0);
      force = force + dir * pForceVal;

      // Lực nhớt giữ hạt dính nhau (Viscosity Force)
      let vForceVal = u.viscosity * (45.0 / (3.141592 * pow(h, 6.0))) * (h - r);
      force = force + (p_j.vel - p.vel) * (vForceVal / p_j.density);
    }
  }

  // Tương tác chuột tạo sóng nước đẩy hạt
  if (u.mouseActive == 1u) {
    let toMouse = u.mousePos - p.pos;
    let mouseDist = length(toMouse);
    if (mouseDist < 0.6 && mouseDist > 0.01) {
      force = force - (toMouse / mouseDist) * (0.6 - mouseDist) * 12.0 * p.density;
    }
  }

  // Tích phân Euler cập nhật gia tốc, vận tốc, vị trí
  let accel = force / p.density;
  p.vel = p.vel + accel * u.deltaTime;
  p.pos = p.pos + p.vel * u.deltaTime;

  // Va chạm biên hộp mềm mại
  let margin = 0.95;
  if (abs(p.pos.x) > margin) {
    p.pos.x = sign(p.pos.x) * margin;
    p.vel.x = -p.vel.x * 0.4;
  }
  if (abs(p.pos.y) > margin) {
    p.pos.y = sign(p.pos.y) * margin;
    p.vel.y = -p.vel.y * 0.4;
  }

  particles[idx] = p;
}

// Chạy Compute Pass 1: Tính Mật độ & Áp suất
const compPass1 = commandEncoder.beginComputePass();
compPass1.setPipeline(densityPipeline);
compPass1.setBindGroup(0, computeBindGroup);
compPass1.dispatchWorkgroups(Math.ceil(particleCount / 64));
compPass1.end();

// Chạy Compute Pass 2: Tính Lực & Cập nhật Euler Integration
const compPass2 = commandEncoder.beginComputePass();
compPass2.setPipeline(forcePipeline);
compPass2.setBindGroup(0, computeBindGroup);
compPass2.dispatchWorkgroups(Math.ceil(particleCount / 64));
compPass2.end();

// Render Pass vẽ các giọt nước
const renderPass = commandEncoder.beginRenderPass(renderPassDesc);
renderPass.setPipeline(renderPipeline);
renderPass.setBindGroup(0, renderBindGroup);
renderPass.draw(6, particleCount, 0, 0); // Vẽ các instanced quads đại diện hạt
renderPass.end();