采用基于 mipmap（多级纹理映射）的方法实现多尺度流体动力学模拟- 1 尺度计算无重复或依赖读取，仅在平流过程中存在相关读取操作，减少数据冗余处理。
- 2 无需随机或随机采样，避免因随机因素导致的模拟结果波动，提升模拟稳定性。
- 3 在保证与其他方法相当保真度的前提下，大幅减少总读取次数，降低计算资源消耗。
- 4 可轻松实现单次遍历采样整个缓冲区（平均情况下），简化数据采样流程，提升5 处理效率。
- 6 计算复杂度转移至 mipmap 生成过程（虽系数较大），优化核心算法计算压力。
- 固定尺度数量下，算法时间复杂度为 O (n)，保证大规模流体模拟时的高效运行。
算法核心围绕纳维 - 斯托克斯方程（Navier-Stokes Equations）展开，同时结合了涡量、旋度、泊松方程等流体力学相关的数学工具👇Shadertoy 中流体模拟的核心是「缓冲区循环迭代」：每个缓冲区（Buffer A-D）对应方程的一个分步，通过 texture(iChannelX, uv)读取上一步结果，计算后输出到当前缓冲区。整体流程对应 GLSL 主函数结构：
void mainImage( out vec4 fragColor, in vec2 fragCoord ) { vec2 uv = fragCoord / iResolution.xy; // 归一化纹理坐标（0-1） vec2 texel = 1.0 / iResolution.xy; // 纹理像素步长（用于梯度/拉普拉斯计算） // 1. 读取上一帧/其他缓冲区数据（速度场、涡量场等） vec4 dataA = texture(iChannel0, uv); // 假设 iChannel0 绑定 Buffer A vec2 vel = dataA.xy; // 速度场（u, v）存在纹理 xy 通道 // 2. 分步计算（平流、扩散、压强投影、涡量限制等） vel = advection(vel, uv, texel); // 平流项 vel = diffusion(vel, uv, texel); // 扩散项 vel = pressureProjection(vel, uv, texel); // 压强投影（满足不可压缩） vel = vorticityConfinement(vel, uv, texel); // 涡量限制（外力项） // 3. 输出结果到当前缓冲区 fragColor = vec4(vel, 0.0, 1.0);}
2 平流项（半拉格朗日方法）
对应公式：u(x,t+Δt)=u(x−u(x,t)⋅Δt,t)
功能：Buffer A 负责子采样平流，模拟流体 “携带” 速度场运动实现思路：回溯粒子过去的位置，采样历史速度场（避免数值不稳定）
float ADVECTION_STEPS = 4.0; // 项目可配置参数（原项目通用选项卡）float ADVECTION_SCALE = 0.9; // 平流尺度系数vec2 advection(vec2 vel, vec2 uv, vec2 texel) { vec2 dt = iTimeDelta * ADVECTION_SCALE; // 时间步长（乘以尺度系数） vec2 pos = uv - vel * dt; // 回溯粒子位置（核心公式） // 多步采样优化（减少数值扩散，对应 ADVECTION_STEPS） vec2 advectedVel = vec2(0.0); for (float i = 0.0; i < 1.0; i += 1.0/ADVECTION_STEPS) { vec2 samplePos = mix(uv, pos, i); // 插值采样位置 advectedVel += texture(iChannel1, samplePos).xy; // 读取 Buffer B 的速度场 } advectedVel /= ADVECTION_STEPS; return advectedVel;}
3 扩散项（拉普拉斯算子 + mipmap 多尺度）对应公式：ν∇2u（二维拉普拉斯：∇2u=∂x2∂2u+∂y2∂2u）
float VISCOSITY = 0.01; // 运动粘度系数（可配置）int MIP_LEVELS = 3; // mipmap 层级数（多尺度）vec2 diffusion(vec2 vel, vec2 uv, vec2 texel) { // 1. 基础拉普拉斯计算（3x3 邻域） vec2 laplacian = vec2(0.0); laplacian += texture(iChannel0, uv + vec2(texel.x, 0.0)).xy; // 右 laplacian += texture(iChannel0, uv - vec2(texel.x, 0.0)).xy; // 左 laplacian += texture(iChannel0, uv + vec2(0.0, texel.y)).xy; // 上 laplacian += texture(iChannel0, uv - vec2(0.0, texel.y)).xy; // 下 laplacian -= 4.0 * vel; // 拉普拉斯核心公式：4u - 邻域和 // 2. 多尺度 mipmap 采样（融合不同尺度扩散） vec2 multiScaleDiff = vec2(0.0); float weightSum = 0.0; for (int k = 0; k < MIP_LEVELS; k++) { float lod = float(k); float weight = 1.0 / (1.0 + lod); // 尺度权重（小尺度权重更高） vec2 mipVel = textureLod(iChannel0, uv, lod).xy; // 采样第 k 层 mipmap multiScaleDiff += mipVel * weight; weightSum += weight; } multiScaleDiff /= weightSum; // 3. 扩散更新（粘性项：nu * laplacian * dt） vec2 diffusedVel = multiScaleDiff + VISCOSITY * laplacian * iTimeDelta; return diffusedVel;}
4 压强投影（多尺度泊松求解器）
对应公式：∇2p=∇⋅u∗（泊松方程）+ u=u∗−∇p（速度修正）
功能：Buffer D 核心功能，满足不可压缩条件 ∇⋅u=0
int POISSON_ITERATIONS = 8; // 泊松迭代次数（平衡精度与效率）// 辅助函数：计算散度（∇·vel = du/dx + dv/dy）float divergence(vec2 uv, vec2 texel) { float du = texture(iChannel3, uv + vec2(texel.x, 0.0)).x - texture(iChannel3, uv - vec2(texel.x, 0.0)).x; float dv = texture(iChannel3, uv + vec2(0.0, texel.y)).y - texture(iChannel3, uv - vec2(0.0, texel.y)).y; return (du + dv) / (2.0 * texel.x); // 除以 2*texel 归一化}// 辅助函数：计算压强梯度（∇p = (dp/dx, dp/dy)）vec2 pressureGradient(vec2 uv, vec2 texel, sampler2D pressureTex) { float dpdx = texture(pressureTex, uv + vec2(texel.x, 0.0)).r - texture(pressureTex, uv - vec2(texel.x, 0.0)).r; float dpdy = texture(pressureTex, uv + vec2(0.0, texel.y)).r - texture(pressureTex, uv - vec2(0.0, texel.y)).r; return vec2(dpdx, dpdy) / (2.0 * texel.x);}vec2 pressureProjection(vec2 velStar, vec2 uv, vec2 texel) { // 1. 计算散度（泊松方程右边项） float div = divergence(uv, texel); // 2. 多尺度泊松求解压强 p（迭代优化） vec4 pressure = vec4(0.0); // 初始压强为 0 for (int iter = 0; iter < POISSON_ITERATIONS; iter++) { for (int k = 0; k < MIP_LEVELS; k++) { // 多尺度迭代 float lod = float(k); vec2 mipUV = uv * (1.0 / pow(2.0, lod)); // mipmap 层级对应 UV 缩放 vec2 mipTexel = texel * pow(2.0, lod); // mipmap 层级对应 texel 缩放 // 泊松迭代核心公式：p_new = (p_up + p_down + p_left + p_right - div * texel²) / 4 float pNeighbors = 0.0; pNeighbors += textureLod(iChannel3, mipUV + vec2(mipTexel.x, 0.0), lod).r; pNeighbors += textureLod(iChannel3, mipUV - vec2(mipTexel.x, 0.0), lod).r; pNeighbors += textureLod(iChannel3, mipUV + vec2(0.0, mipTexel.y), lod).r; pNeighbors += textureLod(iChannel3, mipUV - vec2(0.0, mipTexel.y), lod).r; float pNew = (pNeighbors - div * dot(mipTexel, mipTexel)) / 4.0; pressure.r = mix(pressure.r, pNew, 0.5); // 平滑更新 } } // 3. 用压强梯度修正速度场（满足不可压缩） vec2 gradP = pressureGradient(uv, texel, iChannel3); vec2 vel = velStar - gradP; return vel;}
5 涡量与涡量限制
对应公式：ω=∂x∂v−∂y∂u（涡量）+ Fc=ϵ⋅∇(∣∇ω∣+δω)（涡量限制力）
功能：Buffer B 计算涡量，Buffer C 施加限制力，增强涡流细节
float VORTICITY_STRENGTH = 0.1; // 涡量限制强度（可配置）float EPS = 1e-6; // 防止分母为 0// 辅助函数：计算涡量（ω = dv/dx - du/dy）float vorticity(vec2 uv, vec2 texel) { float dvdx = texture(iChannel1, hannel1, uv - vec2(texel.x, 0.0)).y; float dudy = texture(iChannel1, uv + vec2(0.0, texel.y)).x - texture(iChannel1, uv - vec2(0.0, texel.y)).x; return (dvdx - dudy) / (2.0 * texel.x); // 中心差分归一化}// 辅助函数：计算标量场的梯度（∇s = (ds/dx, ds/dy)）vec2 gradient(float s, vec2 uv, vec2 texel, sampler2D tex) { float dsdx = texture(tex, uv + vec2(texel.x, 0.0)).r - texture(tex, uv - vec2(texel.x, 0.0)).r; float dsdy = texture(tex, uv + vec2(0.0, texel.y)).r - texture(tex, uv - vec2(0.0, texel.y)).r; return vec2(dsdx, dsdy) / (2.0 * texel.x);}vec2 vorticityConfinement(vec2 vel, vec2 uv, vec2 texel) { // 1. 计算当前位置涡量 ω float omega = vorticity(uv, texel); // 2. 计算涡量的梯度 ∇ω vec2 gradOmega = gradient(omega, uv, texel, iChannel2); // iChannel2 绑定 Buffer B（涡量场） // 3. 计算涡量限制力 Fc（公式核心） vec2 Fc = VORTICITY_STRENGTH * normalize(gradOmega + EPS) * sign(omega); // 4. 叠加限制力到速度场（外力项） return vel + Fc * iTimeDelta;}
我们的工作是探索编程和数学在艺术和设计中一切可能性，形式不限，艺术创作，跨界合作，展览展示，品牌活动等等。位于大洋彼岸的洛杉矶，我们的使命是建设和维护全球华人创意编程和算法艺术社群，形式有开展线上教学，举办社群活动，策划举办艺术展览，挖掘发现本地社群里值得被传播的人和事通过自媒体进行传播。
课程和业务介绍请扫码查看👇
算法艺术实验室
探索数学与编程在设计与艺术中一切之可能
用运算和美学让你变更酷
主营业务
科学艺术咨询 | 公共艺术 | 未来舞台美术 | 大数据可视化 | 设计人才猎头 | 品牌Event | 教育
