Глибоке занурення #1 — Гідродинаміка SPH: від Нав'є-Стокса…

Чому SPH?

Рівняння Нав'є-Стокса описують рух рідини, але їхнє розв'язання на фіксованій ейлеровій сітці вимагає дорогих кроків проєкції тиску та складного відстеження вільної поверхні. Для інтерактивної, візуальної симуляції, де користувач наливає, розбризкує та перемішує рідину, лагранжів підхід на основі частинок набагато краще підходить. Згладжена гідродинаміка частинок (SPH, вперше описана Гінголдом і Монаганом у 1977 році) представляє рідину як хмару частинок, що несуть масу, швидкість, густину і тиск. Рівняння руху обчислюються повністю в позиціях частинок — без сітки, без сітчастої моделі.

Компроміс: пошук сусідів наївно має складність O(N²), що катастрофічно при 2000+ частинок, необхідних для задовільного візуалу. Майже всі інженерні зусилля в симуляції були спрямовані на прискорення цих пошуків сусідів.

Фізика: SPH у п'яти рівняннях

Кожна величина A у частинці i оцінюється підсумовуванням внесків сусідніх частинок j, зважених функцією згладжування W(r, h):

A(xᵢ) ≈ Σⱼ  mⱼ · (Aⱼ / ρⱼ) · W(|xᵢ - xⱼ|, h)

де:
  mⱼ  = маса частинки j
  ρⱼ  = густина в частинці j
  h   = радіус згладжування (носій ядра)
  W   = кубічне сплайнове ядро

З цієї єдиної формули інтерполяції випливають ключові величини рідини:

Густина — сума зважених мас сусідів у кожній частинці.
Тиск — рівняння стану: p = k(ρ − ρ₀), де k — параметр жорсткості, а ρ₀ — густина спокою. Слабко стисливий SPH (WCSPH) допускає відхилення густини ~1%, що візуально непомітно.
Сила тиску — симетричний градієнт: Fᵢ_тиск = −Σⱼ mⱼ(pᵢ+pⱼ)/(2ρⱼ) · ∇W
Сила в'язкості — на основі лапласіана: Fᵢ_в'язкість = μ Σⱼ mⱼ(vⱼ−vᵢ)/ρⱼ · ∇²W
Поверхневий натяг — підхід кривизни кольорового поля Мюллера — дорогий, але необхідний для реалістичних крапель.

Сітка просторового хешування

Наївний пошук сусідів O(N²) — порівняння кожної частинки з кожною іншою — миттєво летальний у масштабі. При 2000 частинках це 4 мільйони порівнянь за кадр на 60 FPS: 240 мільйонів порівнянь за секунду на JavaScript. Це нереально.

Рішення — сітка просторового хешування: плоска хеш-мапа з 2D-координат комірок до списків індексів частинок. Кожна комірка має довжину сторони, що дорівнює радіусу ядра h, тому будь-які дві частинки, які можуть бути сусідами, мають бути в тій самій або сусідній комірці. Пошук сусідів тоді виглядає так:

// Фаза побудови — O(N)
for each particle i:
    cell = floor(position[i] / h)
    hashGrid[hash(cell)].push(i)

// Фаза запиту — O(1) амортизовано на частинку
for each particle i:
    for dx in [-1, 0, 1]:
        for dy in [-1, 0, 1]:
            cell = floor(position[i] / h) + (dx, dy)
            for j in hashGrid[hash(cell)]:
                if |xᵢ - xⱼ| < h:
                    process_neighbour(i, j)

Хеш-функція hash(cx, cy) = (cx * 73856093 ^ cy * 19349663) % TABLE_SIZE — стандартний хеш просторової локальності. TABLE_SIZE встановлюється рівним наступному простому числу вище за подвоєну очікувану кількість частинок, щоб утримувати коефіцієнт завантаження нижче 0,5.

Типізовані масиви: тримаємо GC подалі від гарячого шляху

Стан кожної частинки — позиція, швидкість, густина, тиск — зберігається в Float32Array з фіксованим кроком у 8 чисел з плаваючою комою на частинку: [x, y, vx, vy, density, pressure, ax, ay]. Зберігання даних у типізованих масивах усуває виділення об'єктів і тиск збирача сміття у внутрішньому циклі інтегрування. Кеш-дружнє розташування означає, що процесорний попереджувач може прогнозувати шаблони доступу.

Інтегрування: Leapfrog, а не Ейлер

Просте інтегрування Ейлера (v += a·dt; x += v·dt) не зберігає енергію — воно повільно додає енергію в систему, змушуючи частинки прискорюватися і врешті-решт вибухати. Ми використовуємо Leapfrog (Штермер-Верле), який симплектичний (точно зберігає дискретний аналог енергії) і другого порядку точності за ту саму ціну, що й Ейлер:

// Leapfrog: позиції та швидкості зміщені на половину кроку часу
v_half  = v + 0.5 * a_prev * dt     // поштовх
x_new   = x + v_half * dt            // зсув
a_new   = compute_forces(x_new)
v_new   = v_half + 0.5 * a_new * dt // поштовх

// Зберегти для наступного кадру
v = v_new;  x = x_new;  a_prev = a_new;

Крок часу Δt обмежується за умовою CFL: dt = 0.4 * h / v_max. Якщо частинки рухаються швидше, ніж один радіус ядра за крок, пошук околиці ламається.

Обробка меж: без «привидних» частинок

Багато реалізацій SPH використовують «привидні частинки» для застосування граничних умов — додаткові частинки, розміщені поза межами домену, що забезпечують відштовхування. Привидні частинки подвоюють кількість частинок і ускладнюють хеш-сітку. Натомість ми використовуємо простішу силу штрафу: коли частинка перетинає стіну, пружиноподібне відштовхування, пропорційне глибині проникнення, штовхає її назад. У поєднанні з демпфуванням швидкості при відскоку це дає чисті зіткнення зі стінами без тунелювання за розмірів частинок, які ми використовуємо.

Рендеринг: від точок до метакуль через екранний простір

Рендеринг окремих частинок у вигляді кіл виглядає як мішок намиста, а не рідина. Хитрість — рендеринг рідини в екранному просторі:

Рендерити кожну частинку як точковий спрайт WebGL з плавним гаусовим профілем альфи, записуючи глибину в позаекранний буфер.
Застосувати білатеральне розмиття (зі збереженням країв), щоб згладити поле глибини між сусідніми частинками, зберігаючи різкі межі.
Реконструювати нормалі поверхні на піксель з розмитого градієнта глибини (dFdx, dFdy у GLSL).
Застосувати тіньування за Фонгом з кольором води на основі Френеля — синій блиск біля дзеркального відблиску, прозорість під пологими кутами.

Весь конвеєр рендерингу — це два виклики рендерингу: частинки → FBO глибини, потім повноекранний квад → фінальна композиція. На GPU середнього класу це додає менше 0,3 мс до кадру.

Показники продуктивності

Оптимізація	До	Після	Виграш
O(N²) → просторовий хеш	4 FPS (500 ч.)	60 FPS (500 ч.)	15×
Об'єкти → TypedArrays	60 FPS (500 ч.)	60 FPS (2000 ч.)	4× частинок
Ейлер → Leapfrog	Вибухає при dt > 0,008	Стабільно при dt = 0,016	2× dt
Прохід розмиття в екранному просторі	Вигляд намиста	Суцільна поверхня рідини	Візуал

Ключовий висновок: сітка просторового хешування дає ~90% усього виграшу продуктивності. Якщо ви будуєте будь-яку симуляцію частинок O(N²), сітка хешування — перша й найвпливовіша оптимізація, до якої варто звернутися.

Що я зробив би інакше сьогодні

Симуляція працює добре, але, переглядаючи її зараз, я б зробив три зміни:

Обчислювальні шейдери WebGPU — перенести все інтегрування SPH у обчислювальний шейдер. Пошук сусідів чудово відображається на паралельну редукцію на GPU. Потенціал досягти 10 000+ частинок на сучасному обладнанні.
Розв'язувач тиску PCISPH — поточний WCSPH вимагає малих кроків часу для контролю похибки густини. Прогнозно-коригувальний SPH (PCISPH) ітеративно збігається до тиску, дозволяючи більші кроки часу і жорсткішу нестисливість за тієї самої частоти кадрів.
Адаптивний радіус ядра — наразі h фіксований, що означає недостатню вибірку в розріджених регіонах і марну обчислювальну потужність у щільних. Адаптивний h на основі локальної густини дозволив би симуляції обробляти як окремі краплі, що вільно летять, так і щільні басейни з однією кількістю частинок.

Ви можете дослідити симуляцію прямо зараз за адресою mysimulator.uk/fluid/. Кидайте об'єкти, змінюйте в'язкість, збільшуйте кількість частинок — усі повзунки відкривають параметри, описані вище.