SPH-рідина — від Вікіпедії до 60 FPS

Чому саме SPH?

Після симуляції зграйної поведінки (Boids) мені захотілося чогось більш фізичного. Рідини. Води, яку можна перемішати. SPH (Smoothed Particle Hydrodynamics, гідродинаміка згладжених частинок) — це стандартний частинковий метод моделювання рідин; він використовується і в кіно, і в іграх, і в наукових обчисленнях.

Метод представляє рідину як хмару частинок. Кожна частинка несе масу, положення, швидкість і тиск. На кожному кроці часу треба підсумувати внески всіх сусідніх частинок, щоб обчислити сили, а потім проінтегрувати ці сили, щоб перемістити частинки.

Просто в теорії. Диявольськи дорого на практиці.

Математика (коротка версія)

Ядром SPH є інтерполяція, зважена ядерною функцією. Для будь-якого неперервного поля A значення в точці r наближається так:

де W — це згладжувальне ядро, h — радіус згладжування, mⱼ — маса частинки, а ρⱼ — густина.

Сила тиску між частинками виглядає так:

Я використав ядро Poly6 для обчислення густини та ядро Spiky для градієнтів тиску — ненульовий градієнт ядра Spiky в центрі не дає частинкам злипатися в купу.

Наївна реалізація (і чому вона повільна)

Моя перша реалізація була до непристойності проста: для кожної частинки перебираємо всі інші частинки, перевіряємо, чи відстань менша за радіус згладжування h, і накопичуємо внески.

З 1000 частинками це 1 000 000 порівнянь на кадр. При 60 FPS це 60 мільйонів порівнянь за секунду. На JavaScript. Симуляція працювала зі швидкістю приблизно 3 кадри за секунду.

Просторовий хеш — ключова оптимізація

Рішення — просторове хешування: поділити 3D-простір на сітку комірок розміром h × h × h. Кожна частинка прив'язується до комірки за своїм положенням. При обчисленні сусідів перевіряються лише 27 комірок навколо комірки самої частинки.

// Хешування позиції на 3D-сітці в цілочисельний бакет
function hashCell(ix, iy, iz) {
  return (ix * 73856093 ^ iy * 19349663 ^ iz * 83492791) % TABLE_SIZE;
}

// За світовим положенням знаходимо комірку
function cellIndex(x, y, z, h) {
  return hashCell(
    Math.floor(x / h),
    Math.floor(y / h),
    Math.floor(z / h)
  );
}

З просторовим хешем середня частинка перевіряє лише близько 30 сусідів замість 999. Симуляція за одну ніч підскочила з 3 FPS до 58 FPS. Ось у чому цінність правильної структури даних.

Поверхневий натяг і в'язкість

Самі лише сили тиску дають желеподібні грудки. Ще два доданки роблять рідину схожою на справжню воду:

• В'язкість — гасить відносні швидкості між частинками, змушуючи рідину чинити опір зсуву. Додається як лапласіан: F_visc = μ · Σⱼ mⱼ · (vⱼ − vᵢ) / ρⱼ · ∇²W
• Поверхневий натяг — притягує частинки на поверхні всередину, утворюючи гладкі краплі. Обчислюється через нормаль поверхні та кривину.

Підбір правильного коефіцієнта в'язкості зайняв більшу частину дня. Занадто низький — і вода розбризкується всюди, наче кульки. Занадто високий — і вона рухається, як мед.

Отримані уроки

1. Спочатку профілюй. Вузьким місцем завжди був пошук сусідів, а не обчислення сил.
2. Float32 — твій друг. Використання Float32Array замість масивів об'єктів удвічі скоротило навантаження на пам'ять.
3. Ядро має значення. Використання неправильного градієнта ядра (я ненадовго спробував Poly6 для тиску) спричинило накопичення частинок і вибухи симуляції.
4. Стабільність кроку часу. У SPH є максимальний стабільний крок часу. Перевищення цього значення призводить до «вибуху» симуляції. Я додав перевірку умови CFL.