Głębokie zanurzenie #1 — Dynamika płynów SPH: od Naviera-Stokesa do 60 FPS w przeglądarce

Dlaczego SPH?

Równania Naviera-Stokesa opisują ruch płynu, ale ich rozwiązywanie na stałej siatce eulerowskiej wymaga kosztownych kroków projekcji ciśnienia i złożonego śledzenia swobodnej powierzchni. Dla interaktywnej, wizualnej symulacji, w której użytkownik nalewa, rozpryskuje i miesza płyn, podejście lagranżowskie oparte na cząstkach pasuje znacznie lepiej. Wygładzona hydrodynamika cząstek (SPH, po raz pierwszy opisana przez Gingolda i Monaghana w 1977 roku) reprezentuje płyn jako chmurę cząstek niosących masę, prędkość, gęstość i ciśnienie. Równania ruchu są obliczane w całości w pozycjach cząstek — bez siatki, bez modelu siatkowego.

Kompromis: wyszukiwanie sąsiadów jest naiwnie O(N²), co jest katastrofalne przy 2000+ cząstkach potrzebnych do satysfakcjonującej wizualizacji. Niemal cały wysiłek inżynierski w tej symulacji poszedł na przyspieszenie tych wyszukiwań sąsiadów.

Fizyka: SPH w pięciu równaniach

Każda wielkość A przy cząstce i jest szacowana przez sumowanie wkładów sąsiednich cząstek j, ważonych jądrem wygładzającym W(r, h):

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

gdzie:
  mⱼ  = masa cząstki j
  ρⱼ  = gęstość przy cząstce j
  h   = promień wygładzania (nośnik jądra)
  W   = jądro sześciennej funkcji sklejanej

Z tego pojedynczego wzoru interpolacyjnego wynikają kluczowe wielkości płynu:

• Gęstość — suma ważonych mas sąsiadów przy każdej cząstce.
• Ciśnienie — równanie stanu: p = k(ρ − ρ₀), gdzie k to parametr sztywności, a ρ₀ to gęstość spoczynkowa. Słabo ściśliwy SPH (WCSPH) toleruje ~1% odchylenia gęstości, co jest wizualnie niezauważalne.
• Siła ciśnienia — symetryczny gradient: Fᵢ_ciśnienie = −Σⱼ mⱼ(pᵢ+pⱼ)/(2ρⱼ) · ∇W
• Siła lepkości — na bazie laplasjanu: Fᵢ_lepkość = μ Σⱼ mⱼ(vⱼ−vᵢ)/ρⱼ · ∇²W
• Napięcie powierzchniowe — podejście krzywizny pola koloru Müllera — kosztowne, ale niezbędne do realistycznych kropel.

Siatka haszowania przestrzennego

Naiwne wyszukiwanie sąsiadów O(N²) — porównanie każdej cząstki z każdą inną — jest natychmiast zabójcze w skali. Przy 2000 cząstkach to 4 miliony porównań na klatkę przy 60 FPS: 240 milionów porównań na sekundę w JavaScript. Nie ma mowy.

Rozwiązaniem jest siatka haszowania przestrzennego: płaska mapa haszująca ze współrzędnych komórek 2D do list indeksów cząstek. Każda komórka ma długość boku równą promieniowi jądra h, więc dowolne dwie cząstki, które mogłyby być sąsiadami, muszą znajdować się w tej samej lub sąsiedniej komórce. Wyszukiwanie sąsiadów wygląda wtedy tak:

// Faza budowy — O(N)
for each particle i:
    cell = floor(position[i] / h)
    hashGrid[hash(cell)].push(i)

// Faza zapytania — O(1) zamortyzowane na cząstkę
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)

Funkcja haszująca hash(cx, cy) = (cx * 73856093 ^ cy * 19349663) % TABLE_SIZE to standardowy hasz lokalności przestrzennej. TABLE_SIZE jest ustawiane na kolejną liczbę pierwszą powyżej 2× oczekiwanej liczby cząstek, aby utrzymać współczynnik obciążenia poniżej 0,5.

Tablice typizowane: trzymanie GC z dala od gorącej ścieżki

Stan każdej cząstki — pozycja, prędkość, gęstość, ciśnienie — żyje w Float32Array o stałym kroku 8 liczb zmiennoprzecinkowych na cząstkę: [x, y, vx, vy, density, pressure, ax, ay]. Trzymanie danych w tablicach typizowanych eliminuje alokację obiektów i presję garbage collectora w wewnętrznej pętli całkowania. Układ przyjazny dla pamięci podręcznej oznacza, że prefetcher CPU może przewidywać wzorce dostępu.

Całkowanie: Leapfrog, nie Euler

Proste całkowanie Eulera (v += a·dt; x += v·dt) nie zachowuje energii — powoli dodaje energię do systemu, powodując, że cząstki przyspieszają i ostatecznie eksplodują. Używamy metody Leapfrog (Störmera-Verleta), która jest symplektyczna (dokładnie zachowuje dyskretny odpowiednik energii) i drugiego rzędu dokładności przy tym samym koszcie co Euler:

// Leapfrog: pozycje i prędkości są przesunięte o pół kroku czasowego
v_half  = v + 0.5 * a_prev * dt     // kopnięcie
x_new   = x + v_half * dt            // dryf
a_new   = compute_forces(x_new)
v_new   = v_half + 0.5 * a_new * dt // kopnięcie

// Zapisz na następną klatkę
v = v_new;  x = x_new;  a_prev = a_new;

Krok czasowy Δt jest ograniczany warunkiem CFL: dt = 0.4 * h / v_max. Jeśli cząstki poruszają się szybciej niż jeden promień jądra na krok, wyszukiwanie sąsiedztwa się załamuje.

Obsługa granic: bez cząstek widmowych

Wiele implementacji SPH używa „cząstek widmowych” do wymuszania warunków brzegowych — dodatkowych cząstek umieszczonych poza domeną, które zapewniają odpychanie. Cząstki widmowe podwajają liczbę cząstek i komplikują siatkę haszującą. Zamiast tego używamy prostszej siły kary: gdy cząstka przekracza ścianę, sprężynopodobne odpychanie proporcjonalne do głębokości penetracji pcha ją z powrotem. W połączeniu z tłumieniem prędkości przy odbiciu daje to czyste kolizje ze ścianami bez tunelowania przy rozmiarach cząstek, których używamy.

Renderowanie: od punktów do metaball przez przestrzeń ekranu

Renderowanie pojedynczych cząstek jako kółek wygląda jak worek koralików, a nie płyn. Sztuczka to renderowanie płynu w przestrzeni ekranu:

1. Renderuj każdą cząstkę jako sprite punktowy WebGL z gładkim profilem alfa Gaussa, zapisując głębokość do bufora pozaekranowego.
2. Zastosuj rozmycie dwustronne (zachowujące krawędzie), aby wygładzić pole głębokości między sąsiednimi cząstkami, zachowując ostre granice.
3. Zrekonstruuj normalne powierzchni na piksel z rozmytego gradientu głębokości (dFdx, dFdy w GLSL).
4. Zastosuj cieniowanie Phonga z kolorem wody opartym na Fresnelu — niebieski połysk blisko odbicia zwierciadlanego, przezroczystość pod kątami ślizgowymi.

Cały potok renderowania to dwa wywołania rysowania: cząstki → FBO głębokości, następnie pełnoekranowy quad → kompozyt końcowy. Na GPU średniej klasy dodaje to mniej niż 0,3 ms do klatki.

Liczby wydajnościowe

Co zrobiłbym dziś inaczej

Symulacja działa dobrze, ale wracając do niej teraz, wprowadziłbym trzy zmiany:

• Shadery obliczeniowe WebGPU — przenieść całe całkowanie SPH do shadera obliczeniowego. Wyszukiwanie sąsiadów pięknie mapuje się na redukcję równoległą na GPU. Potencjał osiągnięcia 10 000+ cząstek na nowoczesnym sprzęcie.
• Solwer ciśnienia PCISPH — obecny WCSPH wymaga małych kroków czasowych do kontroli błędu gęstości. Predykcyjno- korekcyjny SPH (PCISPH) iteruje do zbieżności ciśnienia, umożliwiając większe kroki czasowe i sztywniejszą nieściśliwość przy tej samej liczbie klatek.
• Adaptacyjny promień jądra — obecnie h jest stałe, co oznacza niedostateczne próbkowanie w rzadkich regionach i marnowanie mocy obliczeniowej w gęstych. Adaptacyjne h oparte na lokalnej gęstości pozwoliłoby symulacji obsługiwać zarówno swobodnie latające kropelki, jak i gęste baseny przy jednej liczbie cząstek.

Możesz zbadać symulację już teraz pod adresem mysimulator.uk/fluid/. Upuszczaj obiekty, zmieniaj lepkość, zwiększaj liczbę cząstek — wszystkie suwaki udostępniają parametry opisane powyżej.