SPH kernel functions: a deep-dive into Poly6, Spiky, Viscosity and Cubic-Spline

1. The role of a kernel in SPH

У SPH будь-яка фізична величина (густина, тиск, швидкість) обчислюється не в точці, а як зважена сума значень сусідніх частинок. Ядро W(r, h) — це вагова функція, яка визначає, наскільки сильно частинка на відстані r впливає на центральну частинку в межах радіуса згладжування h (smoothing length).

Ключове спостереження: ядро — це не просто «вага за відстанню». Різні фізичні величини потребують різних похідних ядра: густина використовує саме W, тиск — градієнт ∇W, а в'язкість — лапласіан ∇²W. Погано підібране ядро для градієнта призводить до нестабільності навіть якщо густина обчислюється коректно.

2. Mathematical requirements

Щоб ядро було коректним для SPH, воно повинно задовольняти чотирьом умовам:

• Нормалізація: ∫ W(r, h) dr = 1 у всьому просторі — інакше інтерполяція зміщує величини систематично.
• Компактний носій (compact support): W(r, h) = 0 при r > h — це дозволяє обмежити пошук сусідів радіусом h і використовувати просторове хешування замість O(N²) перебору.
• Додатність: W(r, h) ≥ 0 всюди — фізичні величини на кшталт густини не можуть бути від'ємними.
• Зменшення з відстанню (monotonic decay): найближчі сусіди повинні впливати сильніше за дальні, а на межі носія (r = h) ядро і, бажано, його похідна плавно йдуть до нуля (C¹-неперервність), щоб уникнути розривів сили.

Нормувальна константа σ (наприклад, 315/(64πh⁹) для Poly6 у 3D) виводиться саме з цієї умови і залежить від розмірності простору (1D/2D/3D мають різні константи).

3. Poly6 kernel

Poly6 (запропоновано Müller et al., 2003) — гладка поліноміальна функція, ідеальна для обчислення густини, де важлива гладкість, а не точність градієнта біля центру.

Проблема Poly6 — його градієнт прямує до нуля при r → 0:

Коли r → 0, множник r обнуляє весь градієнт. Це означає, що для двох частинок, які злиплись впритул, сила відштовхування, обчислена через градієнт Poly6, зникає — саме тому це ядро не можна використовувати для тиску. Використання Poly6 лише для густини — стандартна практика.

4. Spiky kernel

Spiky-ядро вирішує проблему нульового градієнта біля центру за рахунок «гострої» вершини (звідси назва) — похідна лишається ненульовою навіть при r → 0.

При r → 0 градієнт прямує до −45/(πh⁶) · h² — ненульова константа, яка забезпечує сильне відштовхування частинок, що злиплись. Це ключова властивість, яка запобігає взаємопроникненню частинок під високим тиском.

5. Viscosity kernel

В'язкість вимагає лапласіана поля швидкості (∇²W), а не градієнта. Ні Poly6, ні Spiky не дають монотонного, стабільного лапласіана — Müller et al. запропонували окреме ядро саме для цієї мети.

Лапласіан цього ядра лінійний і завжди додатний у межах носія — це гарантує, що в'язка сила завжди діє в напрямку вирівнювання швидкостей (як фізична в'язкість і повинна), без паразитних осциляцій, які виникали б із лапласіана Poly6 чи Spiky.

6. Cubic-spline (B-spline) kernel

Ядро Мональхана-Ласерсона (кубічний B-сплайн) — стандарт у науковій астрофізичній SPH (GADGET, PHANTOM), на відміну від ігрового набору Poly6/Spiky/Viscosity. Воно кусочно-визначене на двох інтервалах і дає один набір формул для густини, градієнта і лапласіана — з коректною поведінкою в усіх трьох випадках:

На відміну від Poly6/Spiky/Viscosity, кубічний сплайн C²-гладкий (неперервна друга похідна) — важливо для довгих астрофізичних симуляцій, де накопичення похибки на мільярдах кроків критичне. Носій охоплює 2h замість h — тому потребує більше сусідів (≈ 50 у 3D) для тієї ж роздільної здатності.

7. Kernel comparison

8. Common pitfalls

Using one kernel for everything

Найпоширеніша помилка новачків — обчислювати і густину, і тиск, і в'язкість через градієнт одного й того ж Poly6. Результат: частинки, злипаючись, не відштовхуються (градієнт → 0), і симуляція «замерзає» в купки замість плавної течії рідини.

Wrong normalization constant for dimensionality

Константи 315/(64πh⁹) для Poly6 і 15/(πh⁶) для Spiky коректні лише в 3D. Для 2D-симуляцій (площинна рідина в canvas/WebGL side-view) потрібні окремі нормувальні множники — інакше загальна маса системи «пливе».

Smoothing length too small or too large

Занадто мале h → мало сусідів → шумна, нестабільна апроксимація (недостатньо статистики). Занадто велике h → надмірне згладжування, втрата деталей, вища вартість O(N·k). Практичне правило: h таке, щоб у носії було 20–40 сусідів у 2D, 30–60 у 3D.

9. Псевдокод

Три ядра в одному об'єкті (JS)

// Ядра SPH — 2D варіант (canvas / WebGL side simulations)
function makeKernels(h) {
  const h2 = h * h, h6 = h**6, h9 = h**9;

  return {
    // Poly6 — тільки для густини
    poly6(r) {
      if (r > h) return 0;
      const diff = h2 − r * r;
      return (315 / (64 * Math.PI * h9)) * diff**3;
    },

    // Spiky gradient — сила тиску (ненульовий градієнт при r→0)
    spikyGrad(r, dir) {
      if (r > h || r === 0) return [0, 0];
      const coeff = −45 / (Math.PI * h6) * (h − r)**2;
      return [coeff * dir[0], coeff * dir[1]];
    },

    // Viscosity Laplacian — гарантовано монотонний, стабільний
    viscLaplacian(r) {
      if (r > h) return 0;
      return (45 / (Math.PI * h6)) * (h − r);
    }
  };
}