Płyn SPH — od Wikipedii do 60 FPS

Dlaczego SPH?

Po Boids chciałem czegoś, co wydawało się bardziej fizyczne. Płyn. Woda, którą można zamieszać. SPH (Smoothed Particle Hydrodynamics — hydrodynamika cząstek wygładzanych) to standardowa metoda cząsteczkowa symulacji płynów — używana zarówno w filmach, grach, jak i obliczeniach naukowych.

Metoda reprezentuje płyn jako chmurę cząstek. Każda cząstka niesie masę, pozycję, prędkość i ciśnienie. W każdym kroku czasowym sumujesz wkłady wszystkich pobliskich cząstek, by obliczyć siły, a następnie całkujesz te siły, by przesunąć cząstki.

Proste w teorii. Diabelnie kosztowne w praktyce.

Matematyka (wersja skrócona)

Sercem SPH jest interpolacja ważona jądrem. Dla dowolnego pola ciągłego A wartość w pozycji r jest przybliżana jako:

Gdzie W to jądro wygładzające, h to promień wygładzania, mⱼ to masa cząstki, a ρⱼ to jej gęstość.

Siła ciśnienia między cząstkami wygląda tak:

Użyłem jądra Poly6 dla gęstości i jądra Spiky dla gradientów ciśnienia — niezerowy gradient jądra Spiky w centrum zapobiega zbijaniu się cząstek w gromady.

Naiwna implementacja (i dlaczego jest wolna)

Moja pierwsza implementacja była banalnie prosta: dla każdej cząstki przejdź przez wszystkie pozostałe, sprawdź, czy odległość mieści się w promieniu wygładzania h, zsumuj wkłady.

Przy 1000 cząstek to 1 000 000 porównań na klatkę. Przy 60 FPS to 60 milionów porównań na sekundę. W JavaScript. Symulacja działała przy ~3 FPS.

Hasz przestrzenny — kluczowa optymalizacja

Rozwiązaniem jest haszowanie przestrzenne: podziel przestrzeń 3D na siatkę komórek o rozmiarze h × h × h. Każda cząstka jest przypisana do komórki na podstawie swojej pozycji. Przy obliczaniu sąsiadów sprawdzasz tylko 27 komórek otaczających komórkę cząstki.

// Haszowanie pozycji siatki 3D do wiadra całkowitego
function hashCell(ix, iy, iz) {
  return (ix * 73856093 ^ iy * 19349663 ^ iz * 83492791) % TABLE_SIZE;
}

// Dla pozycji w świecie, znajdź komórkę
function cellIndex(x, y, z, h) {
  return hashCell(
    Math.floor(x / h),
    Math.floor(y / h),
    Math.floor(z / h)
  );
}

Dzięki haszowi przestrzennemu przeciętna cząstka sprawdza tylko ~30 sąsiadów zamiast 999. Symulacja skoczyła z 3 FPS do 58 FPS z dnia na dzień. To jest wartość odpowiedniej struktury danych.

Napięcie powierzchniowe i lepkość

Same siły ciśnienia dają galaretowate bryły. Dwa dodatkowe człony sprawiają, że wygląda to jak prawdziwa woda:

• Lepkość — tłumi względne prędkości między cząstkami, sprawiając, że płyn opiera się ścinaniu. Dodana jako człon Laplace'a: F_visc = μ · Σⱼ mⱼ · (vⱼ − vᵢ) / ρⱼ · ∇²W
• Napięcie powierzchniowe — przyciąga cząstki na powierzchniach do wewnątrz, tworząc gładkie kropelki. Obliczane przez normalną powierzchni i krzywiznę.

Dobranie właściwego współczynnika lepkości zajęło większość popołudnia strojenia. Za niski i woda pryska wszędzie jak kulki. Za wysoki i porusza się jak miód.

Wnioski

1. Profiluj najpierw. Wąskim gardłem zawsze było wyszukiwanie sąsiadów, nie obliczanie sił.
2. Float32 to twój przyjaciel. Użycie Float32Array zamiast tablic obiektów zmniejszyło zużycie przepustowości pamięci o połowę.
3. Jądro ma znaczenie. Użycie złego gradientu jądra (na chwilę spróbowałem Poly6 dla ciśnienia) powodowało zbijanie się cząstek i eksplozje.
4. Stabilność kroku czasowego. SPH ma maksymalny stabilny krok czasowy. Przekroczenie go powoduje eksplozję symulacji. Dodałem sprawdzenie warunku CFL.