WENO-схема — захоплення ударних хвиль високого порядку для гіперболічних рівнянь

1. Гіперболічні закони збереження та проблема Гіббса

Скалярний гіперболічний закон збереження має вигляд ∂u/∂t + ∂f(u)/∂x = 0. Навіть коли потік f(u) і початкові дані гладкі, характеристики можуть перетнутися за скінченний час, породжуючи справжній розрив — удар. Системи на кшталт стисливих рівнянь Ейлера поводяться так само: гладкий початковий імпульс крутішає до фронту удару зі стрибком густини, швидкості й тиску.

Стандартні реконструкції високого порядку скінченних різниць чи скінченних об'ємів (центровані, або поліноміальна інтерполяція з фіксованим шаблоном) розроблені у припущенні локальної гладкості. Поблизу розриву це припущення порушується, і реконструйований поліном перевищує та занижує справжній розв'язок — числовий аналог явища Гіббса в рядах Фур'є. Схеми низького порядку, але стійкі (проти течії першого порядку) уникають коливань, розмазуючи удар на багато комірок, руйнуючи роздільну здатність. WENO-схема претендує на золоту середину: високий порядок там, де розв'язок гладкий, стійкість схожу на першопорядкову — там, де ні — без ручного налаштування детектора удару.

2. Від ENO до WENO — ідея реконструкції

Попередником WENO є ENO (Essentially Non-Oscillatory, Хартен і Ошер, 1987). Для кожної комірки ENO будує кілька шаблонів-кандидатів однакового порядку і обирає єдиний найгладкіший — за тестом гладкості на основі розділених різниць — повністю відкидаючи решту. Це уникає коливань, але реконструкція не є гладкою за коефіцієнтами: крихітна зміна даних може змінити, який шаблон "обрано", що шкодить точності та збіжності до стаціонарного стану.

WENO (Лю, Ошер і Чан, 1994; Цзян і Шу, 1996) замінює бінарний вибір опуклою комбінацією всіх шаблонів-кандидатів, кожен із яких зважений відповідно до своєї гладкості. У гладких областях усі шаблони роблять внесок, близький до їхніх "ідеальних" (лінійних) вагів, і комбінація відновлює максимально можливий порядок (2r−1 для r шаблонів-кандидатів по r точок кожен). Поблизу розриву вага будь-якого шаблону, що перетинає стрибок, прямує до нуля, і схема плавно пониженням переходить до гладких шаблонів низького порядку, які оминають стрибок.

Для WENO5 (п'ятого порядку, найпоширеніший варіант) 5-точковий шаблон вище розбивається на три перекривні 3-точкові шаблони-кандидати: S₀ = {i−2, i−1, i}, S₁ = {i−1, i, i+1}, S₂ = {i, i+1, i+2}. Кожен дає поліноміальну реконструкцію третього порядку для значення на межі комірки i+½; WENO змішує всі три.

3. Індикатори гладкості

Індикатор гладкості Цзяна-Шу β_k для шаблону S_k вимірює повну варіацію похідних локального поліному реконструкції p_k над коміркою:

β_k велике там, де шаблон-кандидат охоплює стрибок чи крутий градієнт, і мале (в ідеалі нуль для лінійної функції) там, де дані гладкі. β_k — та єдина величина, яка дозволяє WENO розрізняти "гладке" від "розривне" без явного прапорця виявлення удару.

4. Нелінійні ваги — WENO5

Три реконструкції-кандидати u на межі i+½:

У поєднанні з лінійними ("ідеальними") вагами d₀ = 1/10, d₁ = 6/10, d₂ = 3/10 вони дають оптимальну центровану схему п'ятого порядку в гладких областях. WENO замінює d_k нелінійними вагами ω_k, що залежать від β_k:

Там, де розв'язок гладкий, β₀ ≈ β₁ ≈ β₂ і ω_k → d_k, відновлюючи повну точність п'ятого порядку. Поблизу удару β_k будь-якого шаблону, що перетинає розрив, зростає невимірно, його ω_k → 0, і реконструкцію домінують решта гладких шаблонів — автоматично й неперервно, без порогового налаштування.

5. Реконструкція потоку та розв'язувачі Рімана

Для системи законів збереження (наприклад, 1D рівнянь Ейлера) WENO реконструює не саме поле, а сам потік, попередньо розщепивши його на праву та ліву частини (розщеплення потоку, наприклад Лакса-Фрідріхса: f± = ½(f(u) ± α u), де α = max|характеристична швидкість|). WENO реконструює f⁺ за шаблоном, зміщеним проти течії (вліво), а f⁻ — за дзеркальним шаблоном, і числовий потік на межі — F_{i+1/2} = f⁺_{i+1/2} + f⁻_{i+1/2}.

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

6. Інтегрування за часом — SSP-RK3

Просторова реконструкція високого порядку марна, якщо поєднана з інтегратором за часом низького порядку чи нестійким. WENO майже завжди поєднують зі схемою Рунге-Кутти, що зберігає сильну стійкість (Shu-Osher SSP-RK3) — опуклою комбінацією підкроків прямого Ейлера, тому вона успадковує властивість TVD (спадання повної варіації) просторового оператора:

Крок часу обмежений звичайною умовою CFL Δt ≤ CFL · Δx / max|λ|, з CFL зазвичай 0.3–0.6 для WENO5 у поєднанні з SSP-RK3.

7. Реалізація на JavaScript — 1D адвекція та Ейлер

// ── Реконструкція WENO5 на межі (зсув вліво) ─────────────────────
function weno5(um2, um1, u0, up1, up2) {
  const EPS = 1e-6;

  // Реконструкції шаблонів-кандидатів
  const u0r = (1/3)*um2 - (7/6)*um1 + (11/6)*u0;
  const u1r = -(1/6)*um1 + (5/6)*u0  + (1/3)*up1;
  const u2r = (1/3)*u0  + (5/6)*up1 - (1/6)*up2;

  // Індикатори гладкості (Цзян-Шу)
  const b0 = (13/12)*(um2-2*um1+u0)**2 + 0.25*(um2-4*um1+3*u0)**2;
  const b1 = (13/12)*(um1-2*u0+up1)**2  + 0.25*(um1-up1)**2;
  const b2 = (13/12)*(u0-2*up1+up2)**2  + 0.25*(3*u0-4*up1+up2)**2;

  // Нелінійні ваги (лінійні ваги d = 0.1, 0.6, 0.3)
  const a0 = 0.1 / (EPS+b0)**2;
  const a1 = 0.6 / (EPS+b1)**2;
  const a2 = 0.3 / (EPS+b2)**2;
  const sum = a0 + a1 + a2;

  return (a0*u0r + a1*u1r + a2*u2r) / sum;
}

// ── Розщеплення потоку Лакса-Фрідріхса + WENO5 + SSP-RK3 ─────────
function wenoRHS(u, flux, alpha, dx, N) {
  const fp = new Float64Array(N), fm = new Float64Array(N);
  for (let i = 0; i < N; i++) {
    const f = flux(u[i]);
    fp[i] = 0.5 * (f + alpha * u[i]);  // права частина
    fm[i] = 0.5 * (f - alpha * u[i]);  // ліва частина
  }
  const rhs = new Float64Array(N);
  for (let i = 2; i < N - 2; i++) {
    const Fp = weno5(fp[i-2], fp[i-1], fp[i], fp[i+1], fp[i+2]);
    const Fm = weno5(fm[i+2], fm[i+1], fm[i], fm[i-1], fm[i-2]);
    const FpL = weno5(fp[i-3]??fp[i-2], fp[i-2], fp[i-1], fp[i], fp[i+1]);
    const FmL = weno5(fm[i+1], fm[i], fm[i-1], fm[i-2], fm[i-3]??fm[i-2]);
    const F_ip = Fp + Fm;      // числовий потік на i+1/2
    const F_im = FpL + FmL;   // числовий потік на i-1/2
    rhs[i] = -(F_ip - F_im) / dx;
  }
  return rhs;
}

// ── Крок часу SSP-RK3 ─────────────────────────────────────────────
function sspRK3Step(u, dt, rhsFn) {
  const N = u.length;
  const u1 = new Float64Array(N), u2 = new Float64Array(N), out = new Float64Array(N);
  const L0 = rhsFn(u);
  for (let i = 0; i < N; i++) u1[i] = u[i] + dt * L0[i];
  const L1 = rhsFn(u1);
  for (let i = 0; i < N; i++) u2[i] = 0.75*u[i] + 0.25*u1[i] + 0.25*dt*L1[i];
  const L2 = rhsFn(u2);
  for (let i = 0; i < N; i++) out[i] = (1/3)*u[i] + (2/3)*u2[i] + (2/3)*dt*L2[i];
  return out;
}

8. Застосування та варіанти

• Ударна труба Сода: канонічний тест 1D рівнянь Ейлера — діафрагма, що розділяє два постійні стани, розривається, породжуючи удар, контактний розрив і віяло розрідження. WENO5 розв'язує всі три з чіткими, безколиваючими профілями.
• WENO-Z та відображена WENO (WENO-M): варіанти, що зменшують надлишкову числову дисипацію в гладких екстремумах, зберігаючи стійкість захоплення удару класичного WENO5.
• Компактні / гібридні схеми: поєднують дешеву центральну компактну схему в гладких областях із WENO-реконструкцією лише поблизу виявлених розривів, скорочуючи обчислювальні витрати для переважно гладких потоків (наприклад, турбулентність із рідкісними ударами).
• Багатовимірні системи: WENO застосовується вимір-за-виміром на структурованих сітках (скінченно-різницева WENO) або через реконструкцію на неструктурованих сітках (скінченно-об'ємна WENO) для складних геометрій, як-от надзвукові профілі крила.
• Зв'язок із LBM та розв'язувачами Нав'є-Стокса на цьому сайті: на відміну від Латтіс-Больцмана чи розв'язувача Нав'є-Стокса з проекцією тиску, WENO орієнтований на повністю стисливий, ударотвірний режим — природний наступний крок, коли потік уже не можна безпечно вважати нестисливим.