实战:基于纳维-斯托克斯方程的水墨流体模拟)
计算管线Compute Pipeline实战基于纳维-斯托克斯方程的水墨流体模拟在计算机图形学与东方视觉艺术的交叉领域如何让水墨在屏幕上呈现出宛若真实流体般的灵动扩散一直是最高难度的技术挑战之一。早期的网页水墨特效往往采用简化的粒子贴图或者基于 CPU 的二维网格插值。这种做法在粒子数量较少时勉强可行但一旦要模拟具有真实黏性Viscosity、涡流缠绕Vorticity Confinement以及遇到障碍物时自然分流的复杂水墨流动计算量就会呈指数级暴增一个 $512 \times 512$ 的速度场网格每秒需要执行数千万次偏微分方程数值解算。在单线程的 JavaScript 甚至是 WebAssembly 中这都会直接把主线程卡死到不足 5 帧。WebGPU 计算管线Compute Pipeline的诞生为彻底攻克这一物理算力难题提供了完美的武器。与只能处理顶点和片元的传统图形管线不同WebGPU 的计算管线允许我们直接调用 GPU 底层的成千上万个通用算力核心ALU以原生并发的姿态解算著名的不可压缩流体纳维-斯托克斯Navier-Stokes方程。本文将直接深入到 WGSL 计算着色器的数学世界带你在浏览器中手写一套运行在显卡上的原生水墨流体求解器。物理建模不可压缩流体的四大算子拆解不可压缩流体的运动规律由著名的纳维-斯托克斯方程刻画$$\frac{\partial \mathbf{u}}{\partial t} -(\mathbf{u} \cdot \nabla)\mathbf{u} - \frac{1}{\rho}\nabla p \nu \nabla^2 \mathbf{u} \mathbf{F}$$$$\nabla \cdot \mathbf{u} 0$$其中 $\mathbf{u}$ 是速度矢量场$p$ 是压力场$\nu$ 是动黏度系数$\mathbf{F}$ 是外部驱动力如鼠标划过的外力。在图形学实用流体模拟中Jos Stam 提出的稳定流体算法Stable Fluids将这一复杂的偏微分方程在每个时间步内拆解为四个连续的离散求解阶段[初始速度场与墨汁浓度] │ ▼ 1. 平流阶段 (Advection) - 流体沿着自身速度场迁移半拉格朗日法 │ ▼ 2. 外力注入阶段 (External Forces) - 将用户手势/鼠标动力注入速度场 │ ▼ 3. 扩散阶段 (Diffusion) - 黏度与墨汁分子的自然弥散雅可比迭代求解 │ ▼ 4. 投影与无散度修正 (Projection) - 泊松压力方程求解强制流体不可压缩 │ [下一帧自洽的水墨流体场]在 WebGPU 中这四个阶段被分别映射为独立的Compute Pipeline各阶段之间通过双缓冲离屏纹理Ping-Pong Textures在显存内实现纳秒级的数据传递完全不需要经过 CPU 进行任何数据中转核心 WGSL 实现流体平流与发散度解算计算着色器下面是负责平流阶段Advection的核心 WGSL 计算着色器实现// fluid_advection.compute.wgsl struct FluidUniforms { gridSize: vec2f32, deltaTime: f32, decayRate: f32, }; group(0) binding(0) varuniform u: FluidUniforms; group(0) binding(1) var velocitySampler: sampler; // 上一帧的速度场纹理RG 通道分别存储 X 和 Y 速度 group(0) binding(2) var sourceVelocityTexture: texture_2df32; // 上一帧的墨汁浓度纹理R 通道存储墨浓度 group(0) binding(3) var sourceDyeTexture: texture_2df32; // 本帧平流计算后的输出目标纹理可写存储纹理 Storage Texture group(0) binding(4) var outputDyeTexture: texture_storage_2drgba16float, write; compute workgroup_size(16, 16) fn cs_advection(builtin(global_invocation_id) id: vec3u32) { let coords vec2i32(id.xy); if (f32(coords.x) u.gridSize.x || f32(coords.y) u.gridSize.y) { return; } let texCoord (vec2f32(coords) 0.5) / u.gridSize; let texelSize 1.0 / u.gridSize; // 1. 采样当前位置的速度矢量 let velocity textureSampleLevel(sourceVelocityTexture, velocitySampler, texCoord, 0.0).xy; // 2. 半拉格朗日反向追踪Backtracking寻找在 deltaTime 之前流体质点所在的位置 let previousPos texCoord - velocity * u.deltaTime * texelSize; // 3. 双线性插值采样历史位置的墨汁浓度 let sampledDye textureSampleLevel(sourceDyeTexture, velocitySampler, previousPos, 0.0); // 4. 模拟水墨在宣纸上的微弱衰减与渗透沉淀 let nextDye sampledDye * u.decayRate; // 5. 直接写入显存中的目标存储纹理 textureStore(outputDyeTexture, coords, nextDye); }在这段计算着色器中compute workgroup_size(16, 16)声明了一个包含 256 个并发线程的执行工作组Workgroup。当处理一个 $512 \times 512$ 的流体网格时GPU 硬件会自动并发调度 1024 个工作组将每一个网格像素的微积分求解精准分配到显卡的物理流处理器上利用半拉格朗日反向追踪法算法天然具备绝对的数值稳定性Unconditionally Stable即使外界注入极其剧烈的高速碰撞外力流体系统也绝对不会发生数值发散或爆炸采用 16 位浮点存储纹理rgba16float在保证数值微小精度的同时相比 32 位浮点节约了整整一半的显存总线带宽。前端装配WebGPU 计算管线的编排与调度在 TypeScript 中我们通过GPUCommandEncoder将流体求解的各个 Compute Pass 串联成一条闭环指令包export class FluidSimulationController { private device: GPUDevice; private advectionPipeline!: GPUComputePipeline; private bindGroup!: GPUBindGroup; private gridSize 512; constructor(device: GPUDevice) { this.device device; this.initPipelines(); } private initPipelines() { // 编译计算着色器并创建不可变的 ComputePipeline const module this.device.createShaderModule({ code: fluidAdvectionShaderCode, }); this.advectionPipeline this.device.createComputePipeline({ layout: auto, compute: { module, entryPoint: cs_advection, }, }); } public stepSimulation(commandEncoder: GPUCommandEncoder) { // 开启专属的计算通道Compute Pass const computePass commandEncoder.beginComputePass({ label: FluidAdvectionStep, }); computePass.setPipeline(this.advectionPipeline); computePass.setBindGroup(0, this.bindGroup); // 按照 16x16 的工作组尺寸下发并行分发指令 const workgroupsX Math.ceil(this.gridSize / 16); const workgroupsY Math.ceil(this.gridSize / 16); computePass.dispatchWorkgroups(workgroupsX, workgroupsY); computePass.end(); } }注意整个调度过程的纯洁性在整个物理模拟与动画刷新循环中CPU 没有搬运哪怕一个字节的流体像素数据。CPU 仅仅负责向显卡队列提交一条轻量的dispatchWorkgroups调度指令所有密集的数学解算、双缓冲纹理读写以及最终的水墨渲染全部在显卡的物理显存内部自闭环完成视觉升华工笔意境与流体动力学的交融当这套流体求解器以稳定的 120 帧在浏览器中运转起来时用户的鼠标或触控笔划过屏幕不再是简单地画出一道生硬的几何线条。随着鼠标的移动运动轨迹被转化为局部的动量冲击波速度场在 GPU 内部激荡出优美而复杂的微小涡流浓黑如漆的墨滴顺着流体涡旋自然翻滚、拉伸、交缠遇到之前生成的宣纸纤维阻力场时墨汁在边缘自然受阻、散开呈现出水墨画大师泼墨挥毫时特有的那种苍润、雄浑与灵动。用人类最顶尖的流体力学数学方程驱动现代显卡的最强算力只为在方寸数字屏幕之上重现中华千年来那一份独步天下的水墨气韵。这不仅是工程与美学的完美融合更是技术探索者心中永不熄灭的浪漫。