WENO 도식

1. 개요[편집]

WENO 도식 Weighted Essentially Non-Oscillatory Scheme
원조	Liu · Osher · Chan (1994), J. Comput. Phys. 115
표준형	Jiang · Shu (1996), J. Comput. Phys. 126 — WENO-JS
핵심	후보 스텐실을 고르지 않고 가중 결합
차수	매끄러운 영역 $2r-1$ (표준 $r=3$ → 5차)
스텐실	폭 $2r-1$ 점, 후보 $r$ 개
주요 변종	WENO-M · WENO-Z · TENO · MPWENO

ENO는 스텐실 하나를 고르고 나머지를 버렸다. WENO는 “버릴 거면 왜 계산했냐”고 물었다.

WENO 도식(Weighted Essentially Non-Oscillatory scheme)은 ENO 도식의 후보 스텐실 중 하나를 골라 쓰는 대신, 각 후보의 매끄러움을 숫자로 재서 그 값에 따라 전부 가중 평균하는 고차 재구성 기법이다. 불연속을 물고 있는 후보에는 거의 0에 가까운 가중치를, 매끈한 후보에는 큰 가중치를 준다.

이 한 번의 발상 전환이 두 가지를 동시에 해결했다. 첫째, 매끄러운 영역에서는 모든 후보가 살아 있으므로 $r$ 개의 3점 스텐실을 합친 $2r-1$ 점짜리 스텐실의 정확도가 통째로 나온다 — ENO가 $r$ 차에 머물던 것을 $2r-1$ 차로 끌어올린 것이다. 둘째, 가중치가 자료의 매끄러운 함수라서 스텐실이 딸깍딸깍 바뀌며 생기던 ENO의 고질병(정상상태 수렴 정체, 반올림 민감성)이 사라진다.

류·오셔·챈(1994)이 원형을 내놓고, 지앙·슈(1996)가 지금 쓰는 매끄러움 지표를 정의하며 5차를 실현했다. 이후 30년간 고차 충격파 포착의 사실상 표준으로 군림 중이며, 압축성 난류 직접수치모사부터 천체물리 자기유체역학 코드까지 “일단 WENO5”가 국룰이다.

2. 5차 재구성 — 숫자를 전부 적는다[편집]

$r=3$, 즉 WENO5의 경우를 끝까지 적어 둔다. 셀 $i$ 를 중심으로 왼쪽 편향 재구성값 $v_{i+1/2}^{-}$ 를 구한다고 하자. 후보 스텐실은 셋이다.

$$S_0 = \{i-2,\,i-1,\,i\},\qquad S_1 = \{i-1,\,i,\,i+1\},\qquad S_2 = \{i,\,i+1,\,i+2\}$$

각 스텐실이 주는 3차 정확도 재구성값은

$$v^{(0)} = \tfrac{1}{6}\bigl(2v_{i-2} - 7v_{i-1} + 11v_{i}\bigr),\quad v^{(1)} = \tfrac{1}{6}\bigl(-v_{i-1} + 5v_{i} + 2v_{i+1}\bigr),\quad v^{(2)} = \tfrac{1}{6}\bigl(2v_{i} + 5v_{i+1} - v_{i+2}\bigr)$$

이고, 이 셋을 선형 가중치

$$d_0 = \tfrac{1}{10},\qquad d_1 = \tfrac{6}{10},\qquad d_2 = \tfrac{3}{10}$$

로 섞으면 5점 스텐실 전체가 주는 5차 재구성이 정확히 복원된다. 이 $d_k$ 는 “최적 가중치”라고도 불리는데, 매끄러운 영역에서 도달하고 싶은 목표가 바로 이 조합이기 때문이다.

매끄러움 지표(smoothness indicator) $\beta_k$ 는 각 후보 다항식 $p_k(x)$ 의 도함수 제곱을 셀 위에서 적분한 것이다. 일반형은

$$\beta_k = \sum_{l=1}^{r-1}\int_{x_{i-1/2}}^{x_{i+1/2}} \Delta x^{\,2l-1}\left(\frac{d^{\,l}p_k}{dx^{\,l}}\right)^{2} dx$$

이고, $\Delta x^{2l-1}$ 이 붙어 있어 $\beta_k$ 는 격자 폭에 대해 균질하다. $r=3$ 에서 이 적분을 손으로 하면 다음 세 줄이 나온다.¹

$$\beta_0 = \tfrac{13}{12}\bigl(v_{i-2}-2v_{i-1}+v_{i}\bigr)^2 + \tfrac{1}{4}\bigl(v_{i-2}-4v_{i-1}+3v_{i}\bigr)^2$$ $$\beta_1 = \tfrac{13}{12}\bigl(v_{i-1}-2v_{i}+v_{i+1}\bigr)^2 + \tfrac{1}{4}\bigl(v_{i-1}-v_{i+1}\bigr)^2$$ $$\beta_2 = \tfrac{13}{12}\bigl(v_{i}-2v_{i+1}+v_{i+2}\bigr)^2 + \tfrac{1}{4}\bigl(3v_{i}-4v_{i+1}+v_{i+2}\bigr)^2$$

마지막으로 비정규화 가중치와 최종 가중치는

$$\alpha_k = \frac{d_k}{(\varepsilon + \beta_k)^2},\qquad \omega_k = \frac{\alpha_k}{\alpha_0+\alpha_1+\alpha_2},\qquad v^{-}_{i+1/2} = \sum_{k=0}^{2}\omega_k\, v^{(k)}$$

이다. 오른쪽 편향 $v^{+}_{i-1/2}$ 는 격자에 대해 좌우 대칭으로 뒤집어 그대로 쓴다. 코드에서 이 대칭을 손으로 다시 유도하다 부호를 틀리는 것이 입문자의 첫 관문이다.

3. 지표가 실제로 하는 일[편집]

$\beta_k$ 를 매끄러운 영역에서 전개하면 $\beta_k = f'^2\Delta x^2 + O(\Delta x^4)$ 로, 세 지표의 선행항이 전부 같다. 차이는 $O(\Delta x^4)$ 부터 나오므로 상대 편차가 $O(\Delta x^2)$ 이고, 따라서 $\omega_k = d_k + O(\Delta x^2)$ 가 되어 가중치가 선형 가중치로 수렴한다. 매끄러운 곳에서 5차가 나오는 이유가 이것이다.

불연속을 물고 있는 후보에서는 이야기가 달라진다. 그 스텐실의 $\beta_k$ 는 점프 크기의 제곱, 즉 $O(1)$ 이 되고 매끈한 후보의 $\beta$ 는 $O(\Delta x^2)$ 이므로, 제곱이 붙은 $\alpha_k$ 의 비는 $O(\Delta x^{4})$ 수준으로 벌어진다. 나쁜 후보의 가중치가 사실상 소멸하는 것이다. “고르지 않는다”고 했지만 불연속 앞에서는 결과적으로 고른 것과 거의 같아진다 — 다만 그 전환이 계단이 아니라 매끄러운 곡선이라는 점이 다르다.²

여기서 지수 2를 쓰는 것($\alpha_k \propto \beta_k^{-2}$)에 깊은 이유가 있는 것은 아니다. 지앙과 슈가 실험적으로 고른 값이고, 지수를 키우면 불연속 판별이 날카로워지는 대신 매끄러운 영역에서 가중치가 예민해진다. 이 자유도를 조절하는 것이 이후 변종들의 공통된 손잡이다.

4. 임계점 — WENO-JS의 아픈 곳[편집]

$f'(x_i) = 0$ 인 임계점(critical point)에서는 위 전개의 선행항이 사라진다. $\beta_k = \tfrac{13}{12}f''^2\Delta x^4 + O(\Delta x^5)$ 가 되고, 이번에는 $O(\Delta x^5)$ 항이 세 후보에서 서로 달라 상대 편차가 $O(\Delta x)$ 로 커진다. 결과적으로 $\omega_k - d_k = O(\Delta x)$ 이고, 플럭스 차분을 취하면 도식의 차수가 5차에서 3차로 떨어진다. $f'=f''=0$ 이면 2차까지 내려간다.

헨릭·아슬람·파워스(2005)가 이 현상을 정리했고, 5차를 보장하는 충분조건이 $\omega_k - d_k = O(\Delta x^{3})$ 임을 보인 뒤 두 가지 처방을 제시했다.

• WENO-M (mapped WENO). 계산된 $\omega_k$ 를 사상함수

  $$g_k(\omega) = \frac{\omega\bigl(d_k + d_k^2 - 3d_k\omega + \omega^2\bigr)}{d_k^2 + \omega(1-2d_k)}$$

  에 한 번 통과시킨 뒤 다시 정규화한다. 이 함수는 $g_k(d_k)=d_k$, $g_k'(d_k)=g_k''(d_k)=0$ 이라 $d_k$ 근방에서 납작하다. 그래서 $\omega_k - d_k = O(\Delta x)$ 였던 것이 $g_k(\omega_k) - d_k = O(\Delta x^3)$ 으로 눌리고 5차가 회복된다. 대가는 후보마다 유리함수 평가가 하나씩 추가되는 비용.

• WENO-Z (Borges 외, 2008). 사상 대신 대역 매끄러움 지표 $\tau_5 = \lvert \beta_0 - \beta_2 \rvert$ 를 하나 더 만들어

  $$\alpha_k = d_k\left[\,1 + \left(\frac{\tau_5}{\beta_k + \varepsilon}\right)^{p}\,\right]$$

  로 쓴다. $\tau_5$ 는 매끄러운 영역에서 $\beta_k$ 보다 몇 차수 작아 괄호가 1에 붙고, 불연속에서는 $\tau_5/\beta_k$ 가 폭발해 나쁜 후보를 죽인다. $p=1$ 이 원 논문의 기본값이고 $p=2$ 를 쓰면 1차 임계점에서도 형식 차수가 살아난다. 곱셈 몇 번으로 WENO-M과 비슷한 효과를 내면서 소산이 더 낮다는 것이 정설이라, 요즘 새로 짜는 코드는 대개 WENO-Z로 간다.

5. $\varepsilon$ 은 0으로 나누지 말라는 뜻이 아니다[편집]

원 논문의 $\varepsilon = 10^{-6}$ 은 ”$\beta_k$ 가 0일 때 터지지 말라”는 안전장치로 소개되지만, 실제로는 도식의 성격을 정하는 파라미터다. $\varepsilon \gg \beta_k$ 인 영역에서는 $\alpha_k \approx d_k/\varepsilon^2$ 가 되어 가중치가 선형 가중치로 굳고, 도식은 그냥 5차 선형 상류 도식이 된다 — 정확하지만 진동에는 무방비다. 반대로 $\varepsilon \to 0$ 이면 판별은 날카롭지만 매끄러운 영역에서 가중치가 반올림 잡음에 흔들린다.

더 고약한 것은 $\varepsilon$ 이 차원을 가진다는 사실이다. $\beta_k$ 의 단위는 $v^2$ 이므로 $\varepsilon = 10^{-6}$ 은 $v = O(1)$ 로 정규화된 문제에서만 의도대로 작동한다. 밀도가 $10^{-10}$ 인 천체물리 문제나 압력이 $10^{5}$ 인 연소 문제에 그대로 들고 가면 도식이 선형 쪽이나 ENO 쪽 극단으로 붙어 버린다. 실무 처방은 문제의 대표 규모 $v_{\rm ref}$ 로 $\varepsilon = 10^{-6}v_{\rm ref}^2$ 처럼 무차원화하거나, $\varepsilon \propto \Delta x^2$ 로 격자에 묶는 것이다. 후자는 임계점에서 차수까지 회복시켜 주지만 불연속 근처의 판별력을 조금 내준다.³

6. 유한체적 재구성 vs 유한차분 플럭스 — 헷갈리면 안 되는 구분[편집]

WENO를 처음 구현할 때 가장 많이 나는 사고가 이 둘을 섞는 것이다. 위에 적은 $v^{(k)}$·$d_k$ 숫자가 양쪽에서 똑같이 나오기 때문에 더 헷갈린다.

유한체적(FV) WENO. 자료는 셀 평균 $\bar u_i$ 다. WENO 재구성으로 면 양쪽의 점값 $u^{\pm}_{i+1/2}$ 를 복원한 뒤, 리만 솔버로 수치 플럭스를 만든다. 유한체적법의 정통 경로이고 임의 격자·비정렬 격자에 쓸 수 있지만, 다차원에서 비용이 폭발한다. 면마다 고차 면적분을 위해 가우스점을 여러 개 두어야 하고, 비정렬 격자에서는 선형 가중치가 음수로 나오는 경우가 생겨 슈의 양·음 분해 같은 별도 처리가 필요하다.

유한차분(FD) WENO. 자료는 점값 $u_i$ 다. 핵심은 플럭스 함수 $f$ 에 대해

$$f(x) = \frac{1}{\Delta x}\int_{x-\Delta x/2}^{x+\Delta x/2} h(\xi)\,d\xi \qquad\Longrightarrow\qquad \left.\frac{\partial f}{\partial x}\right|_{x_i} = \frac{h_{i+1/2} - h_{i-1/2}}{\Delta x}\ \ (\text{정확})$$

를 만족하는 보조함수 $h$ 를 도입하는 것이다. 즉 점값 $f(u_j)$ 를 $h$ 의 셀 평균으로 간주해 위의 재구성 공식을 그대로 먹이면 $h_{i+1/2}$ 가 나오고, 그 차분이 미분의 고차 근사가 된다. 재구성 연산자가 FV와 동일한 이유가 이것이다.

상류 편향을 주려면 플럭스를 나눠야 한다. 표준은 락스-프리드리히스 분해다.

$$f^{\pm}(u) = \tfrac{1}{2}\bigl(f(u) \pm \alpha u\bigr),\qquad \alpha = \max_u \lvert f'(u)\rvert$$

이러면 $df^{+}/du \ge 0$, $df^{-}/du \le 0$ 이 보장되어 $f^{+}$ 는 왼쪽 편향, $f^{-}$ 는 오른쪽 편향 재구성을 쓰면 된다. $\alpha$ 를 전역 최대로 잡으면 전역 LF, 면 근방의 국소 최대로 잡으면 국소 LF이고 후자가 덜 소산적이다.

정리하면 이렇다.

	유한차분 WENO	유한체적 WENO
자료	점값 $u_i$	셀 평균 $\bar u_i$
재구성 대상	분해된 플럭스 $f^{\pm}$	보존 변수 $u$
면에서 하는 일	플럭스 분해 (LF 등)	리만 솔버 호출
격자 제약	균등 또는 매끄러운 사상 필수	임의 격자 가능
다차원 비용	방향별 1차원 반복, 싸다	면당 가우스점 다수, 비싸다
실무 점유율	압도적	비정렬·AMR 특수 상황

“WENO5를 썼다”고 하면 열에 아홉은 FD 판본이다. FD가 균등 격자를 요구하는 이유도 위 유도에서 바로 보인다 — $h$ 의 정의가 셀 폭이 일정하다는 가정 위에 서 있기 때문이다. 곡선 격자에서는 계산 좌표계로 사상한 뒤 자유흐름 보존(free-stream preservation)을 만족하도록 계량항을 같은 차수로 이산화해야 하고, 여기서 어긋나면 균일류가 저절로 일그러진다.

7. 계로 갈 때 — 성분별이냐 특성별이냐[편집]

오일러 방정식 같은 계에 WENO를 먹이는 방법은 둘이다. 보존 변수마다 독립적으로 재구성하는 성분별(component-wise) 방식은 싸지만 접촉 불연속과 충격파 근처에서 압력·속도에 진동이 새어 나온다. 정석은 면마다 로 평균(Roe average)으로 국소 고유벡터 행렬 $R$, $R^{-1}$ 을 만들어 특성 변수로 사영한 뒤 성분별로 WENO를 돌리고 다시 되돌리는 특성별(characteristic-wise) 방식이다.

비용 차이가 만만치 않다. $m$ 개 방정식이면 면마다 $O(m^2)$ 짜리 행렬-벡터 곱이 네 번씩 추가되고, 실제 프로파일링에서 3차원 오일러 WENO5 계산 시간의 절반 이상이 여기로 간다. 그래도 강한 충격파 문제에서는 이 비용을 내는 게 맞다 — 성분별로 갔다가 접촉면에서 압력 진동이 나 결국 다시 짜는 경우가 흔하다.

8. 비용, 그리고 왜 아직도 국룰인가[편집]

먼저 흔한 오해 하나. WENO5의 스텐실 폭은 5점으로, 같은 5차 선형 상류 도식과 똑같다. 세 배가 되는 것은 스텐실이 아니라 면당 산술 연산량이다 — 후보 다항식 3개, 매끄러움 지표 3개, 나눗셈과 정규화가 추가된다. 대략적인 청구서는 이렇다.

• 면당 부동소수점 연산이 2차 MUSCL 대비 5~10배.
• 시간 적분은 보통 3단 SSP 룽게-쿠타라 스텝당 룽게-쿠타법 우변 평가가 3회.
• 헤일로 폭이 방향마다 3셀. MUSCL이 1~2셀이면 되는 것과 비교하면 이게 병렬·적응 격자 세분화 환경에서 진짜 비싼 부분이다. 통신량과 경계 처리 코드가 여기서 불어난다.
• 나눗셈과 분기가 많아 벡터화 효율이 선형 도식보다 낮다.

그럼에도 30년째 표준인 이유는 단순하다.

• 문제마다 맞춰야 하는 상수가 없다. 전변분 감소 계열의 TVB 상수 $M$ 처럼 사용자가 감으로 정해야 하는 값이 없다는 것이 실무에서 생각보다 크다.
• 차수를 $r$ 로 자유롭게 올릴 수 있다. 7차·9차·11차 WENO가 같은 틀에서 그대로 나온다.
• FD 판본이 방향별로 분리돼 구현이 단순하다. 1차원 커널 하나 짜서 세 방향에 돌리면 끝이라 GPU 이식도 쉽다.
• 양수성 보존 제한자(Zhang–Shu)와 궁합이 좋아, 극단적 마하수·저밀도 문제에서 밀도·압력이 음수로 가는 것을 사후에 막을 수 있다.

경쟁자가 없는 것은 아니다. 복잡 형상에서는 불연속 갤러킨법이 강하고, 소산을 더 줄이려면 TENO(targeted ENO)처럼 후보를 문턱으로 완전히 켜고 끄는 계열이나 단조 보존(MP) 계열이 있다. 다만 이들 모두 “WENO 대비”로 논문을 쓰는 처지라는 것이 WENO의 위치를 말해 준다.

9. 실무에서 밟는 지뢰[편집]

• 차수 검증은 매끄러운 문제로, 진동 검증은 충격파 문제로 따로 한다. 그리고 차수 시험 함수에 임계점이 있으면(예: $\sin$ 파) WENO-JS는 5차가 안 나온다. 그건 버그가 아니라 이 문서 §임계점의 그 현상이다.
• $\varepsilon$ 을 안 만졌다면 그건 튜닝을 안 한 게 아니라 기본값을 선택한 것이다. 문제 규모가 1에서 멀면 반드시 확인한다.
• 매끄러운 영역에서도 WENO는 선형 도식보다 소산적이다. 저해상도 대와류 모사에서 WENO를 그대로 쓰면 수치 소산이 아격자 모형보다 커지는 일이 흔하다. 이 문제는 수치 소산과 분산 쪽 이야기.
• 양수성이 깨지면 원인은 대개 재구성이 아니라 시간 적분이다. 스텝을 줄이거나 SSP 계수를 확인한다.
• 가중치를 출력해 보는 것이 최고의 디버깅이다. 균일류에서 $\omega_k$ 가 $(0.1, 0.6, 0.3)$ 에 붙어 있지 않으면 지표 공식이나 편향 방향이 틀린 것이다.

10. 관련 문서[편집]

• ENO 도식 · 전변분 감소 · MUSCL · TENO
• 유한체적법 · 유한차분법 · 리만 솔버 · 불연속 갤러킨법
• 충격파 · 압축성 유동 · 깁스 현상 · 수치 소산과 분산
• 룽게-쿠타법 · CFL 조건 · 차분 도식 · 수정 방정식
• 적응 격자 세분화 · 직접수치모사 · 자기유체역학
• 오일러 방정식 · 특성곡선법 · 웰밸런스드 도식

11. Footnotes[편집]

1. $13/12$ 라는 계수가 대체 어디서 왔냐는 질문을 다들 한 번씩 한다. 2계 도함수 제곱 항의 셀 적분에서 나오는 값으로, $\int_{-1/2}^{1/2}(\text{2계항})^2$ 을 정직하게 계산하면 $1 + 1/12$ 이 나온다. 마법도 아니고 오타도 아니다. ↩

2. 후보를 “고르지 않는다”는 WENO의 슬로건은 사실 절반만 맞다. 불연속 근처에서 가중치가 $10^{-8}$ 까지 떨어지면 그건 부동소수점 관점에서 고른 것과 구별되지 않는다. 다만 그 전환이 미분가능하다는 점 하나로 정상상태 잔차가 기계 정밀도까지 떨어지느냐 마느냐가 갈리므로, 절반이 아니라 결정적인 절반이다. ↩

3. 그래서 논문에 ”$\varepsilon = 10^{-6}$ 을 썼다”고만 적힌 것을 보면 재현이 안 되는 경우가 생긴다. 무차원화 방식을 안 적어 놨기 때문이다. 검증 및 확인이 잔소리하는 “모든 수치 파라미터를 명시하라”가 정확히 이런 곳을 겨냥한 것이다. ↩