Рідина Латтіс-Больцмана у WebGL
Запустіть справжню симуляцію рідини повністю на GPU. Використовуючи пару текстур із рухомою комою та фрагментний шейдер, кожен піксель обчислює функції розподілу одного вузла ґратки — а потім стрімить, зіштовхує та візуалізує поле швидкості рідини в реальному часі.
1. Теорія LBM — D2Q9
Метод Латтіс-Больцмана не розвʼязує рівняння Нав’є-Стокса напряму. Натомість він симулює газ частинок на сітці, де кожен вузол зберігає 9 функцій розподілу f_i (модель D2Q9), що представляють популяції частинок, які рухаються в 9 дискретних напрямках швидкості.
9 векторів швидкості (плюс спокій) та їхні ваги:
Напрямок: e_i вага w_i спокій: (0,0) 4/9 E,N,W,S: (±1,0),(0,±1) 1/9 NE,NW,SW,SE: (±1,±1) 1/36
Кожен крок має дві фази: стрімінг (популяції переносяться вздовж свого напрямку швидкості) та зіткнення (популяції релаксують до локальної рівноваги, що керується вʼязкістю):
// Зіткнення BGK (єдиний час релаксації) f_eq_i = w_i * ρ * (1 + 3(e_i·u) + 4.5(e_i·u)² − 1.5(u·u)) f_i* = f_i + (f_eq_i − f_i) / τ // результат зіткнення // макроскопічні величини: ρ = Σ f_i // густина ρu = Σ f_i * e_i // імпульс ν = (τ − 0.5) / 3 // кінематична вʼязкість
2. Розкладка текстур
D2Q9 потребує 9 чисел з рухомою комою на вузол. Ми використовуємо дві текстури RGBA_FLOAT: одна містить f0..f3, упаковані як RGBA (спокій + E + N + W), а друга містить f4..f7 (S + NE + NW + SW) плюс f8 в альфа-каналі, а маску перешкод упаковано окремо.
// Текстура A: RGBA float (ширина × висота) // R = f0 (спокій) // G = f1 (E) // B = f2 (N) // A = f3 (W) // Текстура B: RGBA float // R = f4 (S) // G = f5 (NE) // B = f6 (NW) // A = f7 (SW) // Текстура C: f8 (SE) + маска перешкод // R = f8 (SE) // G = перешкода (0 або 1) // B, A не використовуються
OES_texture_float та WEBGL_color_buffer_float. WebGL2 підтримує float-текстури нативно. Для нових проєктів використовуйте WebGL2.3. Налаштування WebGL
- Створіть контекст і розширення
const canvas = document.querySelector('canvas'); const gl = canvas.getContext('webgl2'); if (!gl) throw new Error('WebGL2 required'); // Підтримка float-текстур (у WebGL2 вона є за замовчуванням) // Але MRT (кілька цілей рендерингу) потребує: // gl.getExtension('EXT_color_buffer_float'); - Створіть ping-pong фреймбуфери
function createFBOTexture(gl, w, h) { const tex = gl.createTexture(); gl.bindTexture(gl.TEXTURE_2D, tex); gl.texImage2D(gl.TEXTURE_2D, 0, gl.RGBA32F, w, h, 0, gl.RGBA, gl.FLOAT, null); gl.texParameteri(gl.TEXTURE_2D, gl.TEXTURE_MIN_FILTER, gl.NEAREST); gl.texParameteri(gl.TEXTURE_2D, gl.TEXTURE_MAG_FILTER, gl.NEAREST); gl.texParameteri(gl.TEXTURE_2D, gl.TEXTURE_WRAP_S, gl.CLAMP_TO_EDGE); gl.texParameteri(gl.TEXTURE_2D, gl.TEXTURE_WRAP_T, gl.CLAMP_TO_EDGE); const fbo = gl.createFramebuffer(); gl.bindFramebuffer(gl.FRAMEBUFFER, fbo); gl.framebufferTexture2D(gl.FRAMEBUFFER, gl.COLOR_ATTACHMENT0, gl.TEXTURE_2D, tex, 0); return { tex, fbo }; } // Створюємо два набори текстур A/B/C (поточний і наступний) const stateA = [createFBOTexture(gl, W, H), createFBOTexture(gl, W, H)]; let readIdx = 0, writeIdx = 1; - Повноекранний квад для обчислювальних проходів
// Вершинний шейдер: просто повноекранний трикутник const vsSource = `#version 300 es in vec2 aPos; out vec2 vUV; void main() { vUV = aPos * 0.5 + 0.5; gl_Position = vec4(aPos, 0.0, 1.0); }`; // Два трикутники, що покривають NDC: const positions = new Float32Array([ -1,-1, 1,-1, -1,1, -1, 1, 1,-1, 1,1 ]);
4. Шейдер стрімінгу
Прохід стрімінгу зчитує дані із сусідів. Для напрямку i зі швидкістю (cx, cy) значення f_i у позиції (x,y) береться з позиції (x−cx, y−cy) на попередньому кроці:
#version 300 es
precision highp float;
uniform sampler2D uTexA; // f0..f3
uniform sampler2D uTexB; // f4..f7
uniform sampler2D uTexC; // f8 + перешкода
uniform vec2 uInvRes;
in vec2 vUV;
out vec4 fragColor; // оновлена texA
vec2 e[9] = vec2[](
vec2(0,0), vec2(1,0), vec2(0,1), vec2(-1,0),
vec2(0,-1), vec2(1,1), vec2(-1,1), vec2(-1,-1), vec2(1,-1)
);
float sampleF(int i, vec2 uv) {
vec2 offset = -e[i] * uInvRes; // навпаки: тягнемо проти течії
vec4 a = texture(uTexA, uv + offset);
vec4 b = texture(uTexB, uv + offset);
float c = texture(uTexC, uv + offset).r;
// витягуємо канал для напрямку i
if (i == 0) return a.r;
if (i == 1) return a.g;
if (i == 2) return a.b;
if (i == 3) return a.a;
if (i == 4) return b.r;
if (i == 5) return b.g;
if (i == 6) return b.b;
if (i == 7) return b.a;
return c; // i == 8
}
void main() {
// Стрімінг: f0..f3 у texA
float f0 = sampleF(0, vUV);
float f1 = sampleF(1, vUV);
float f2 = sampleF(2, vUV);
float f3 = sampleF(3, vUV);
fragColor = vec4(f0, f1, f2, f3);
}5. Шейдер зіткнень
Після стрімінгу прохід зіткнень обчислює макроскопічні ρ та u, а потім релаксує до рівноваги:
#version 300 es
precision highp float;
uniform sampler2D uTexA; // f0..f3 після стрімінгу
uniform sampler2D uTexB; // f4..f7 після стрімінгу
uniform sampler2D uTexC; // f8 після стрімінгу
uniform float uTau; // час релаксації (1.5–2.0 для стабільності)
in vec2 vUV;
out vec4 fragColor;
const float W[9] = float[](4./9., 1./9., 1./9., 1./9., 1./9.,
1./36., 1./36., 1./36., 1./36.);
const vec2 E[9] = vec2[](
vec2(0,0),vec2(1,0),vec2(0,1),vec2(-1,0),vec2(0,-1),
vec2(1,1),vec2(-1,1),vec2(-1,-1),vec2(1,-1));
float feq(int i, float rho, vec2 u) {
float eu = dot(E[i], u);
float uu = dot(u, u);
return W[i] * rho * (1.0 + 3.0*eu + 4.5*eu*eu - 1.5*uu);
}
void main() {
vec4 a = texture(uTexA, vUV);
vec4 b = texture(uTexB, vUV);
float f8 = texture(uTexC, vUV).r;
float f[9];
f[0]=a.r; f[1]=a.g; f[2]=a.b; f[3]=a.a;
f[4]=b.r; f[5]=b.g; f[6]=b.b; f[7]=b.a;
f[8]=f8;
// Макроскопічні величини
float rho = 0.0;
vec2 u = vec2(0.0);
for (int i = 0; i < 9; i++) {
rho += f[i];
u += f[i] * E[i];
}
u /= rho;
// Зіткнення BGK
float inv_tau = 1.0 / uTau;
for (int i = 0; i < 9; i++) {
f[i] += (feq(i, rho, u) - f[i]) * inv_tau;
}
fragColor = vec4(f[0], f[1], f[2], f[3]);
// texB та texC записуються аналогічно через розкладку MRT
}6. Bounce-back для перешкод
Перешкоди зберігаються як двійкова маска (1 = тверда) у зеленому каналі texC. Коли вузол твердий, замість звичайного стрімінгу+зіткнення застосовуйте bounce-back з прилипанням: відбивайте кожну функцію розподілу в протилежний напрямок:
// У шейдері стрімінгу, після зчитування стрімленого f:
float isObstacle = texture(uTexC, vUV).g;
if (isObstacle > 0.5) {
// Bounce-back: міняємо місцями протилежні напрямки
// f1(E) ↔ f3(W), f2(N) ↔ f4(S), f5(NE) ↔ f7(SW), f6(NW) ↔ f8(SE)
float tmp;
tmp = f[1]; f[1] = f[3]; f[3] = tmp;
tmp = f[2]; f[2] = f[4]; f[4] = tmp;
tmp = f[5]; f[5] = f[7]; f[7] = tmp;
tmp = f[6]; f[6] = f[8]; f[8] = tmp;
}Малюйте перешкоди інтерактивно: в обробнику mousemove на JavaScript рендеріть заповнені кола у текстуру перешкод. Bounce-back коректно забезпечує граничні умови прилипання за будь-якої форми перешкоди.
7. Візуалізація швидкості
Фінальний прохід відображення зчитує поточний стан рідини, обчислює ρ та u, а потім відображає швидкість у колір:
// Фрагментний шейдер візуалізації
void main() {
vec4 a = texture(uTexA, vUV);
vec4 b = texture(uTexB, vUV);
float f8 = texture(uTexC, vUV).r;
float obs = texture(uTexC, vUV).g;
float f[9];
f[0]=a.r;f[1]=a.g;f[2]=a.b;f[3]=a.a;
f[4]=b.r;f[5]=b.g;f[6]=b.b;f[7]=b.a;f[8]=f8;
vec2 u = vec2(0.0);
float rho = 0.0;
for (int i=0; i<9; i++) { rho+=f[i]; u+=f[i]*E[i]; }
u /= rho;
float speed = length(u);
// HSV: відтінок = напрямок швидкості, насиченість+яскравість = швидкість
float hue = atan(u.y, u.x) / (2.0 * PI) + 0.5;
vec3 col = hsv2rgb(vec3(hue, 0.85, clamp(speed * 20.0, 0.0, 1.0)));
if (obs > 0.5) col = vec3(0.2); // темно-сіра перешкода
outColor = vec4(col, 1.0);
}