Lax–Wendroff가 계단에서 떨렸다 — Flux Limiter로 2차 TVD 이류 스킴 만들기

테스트 케이스에서 계단 앞뒤로 이상한 혹이 올라왔다. 선형 이류(linear advection, 속도장이 상수인 순수 이동) 방정식을 Lax–Wendroff로 풀고 있었고, 초기 조건은 단순한 직사각 펄스였다. 2차 스킴이니 1차 업윈드보다 날카로울 거라 기대했지만, 결과는 계단 앞에서 언더슈트, 뒤에서 오버슈트. 로그를 봐도 CFL은 0.5, 안정 조건은 만족.

문제는 알고리즘이 불안정한 것이 아니라 알고리즘이 단조롭지 않은 것이었다.

진동의 정체 — Godunov가 걸어둔 벽#

1D 선형 이류 방정식은

\partial_t q + u\,\partial_x q = 0

속도 $u$ 가 일정할 때 해는 초기 모양이 그대로 오른쪽으로 평행이동하는 것. 수치적으로는 셀 경계에서의 플럭스 $f_{i-1/2}$ 를 정의하고

q_i^{n+1} = q_i^n - \frac{\Delta t}{\Delta x}\bigl(f_{i+1/2} - f_{i-1/2}\bigr)

이 식은 flux-conserving form이라 불리고, 질량·운동량·에너지 같은 보존량을 기계 정밀도까지 지킨다. Lax–Wendroff는 여기에 2차 테일러 보정을 넣은 선형 스킴이다.

그런데 1959년 Godunov가 증명한 정리는 가혹하다. 선형 상수계수 이류에 대해 단조성(monotonicity)을 유지하면서 2차 이상의 정확도를 갖는 선형 스킴은 존재하지 않는다. 둘 중 하나를 포기해야 한다. Lax–Wendroff는 2차 정확도를 택하고 단조성을 버렸다. 그래서 계단에서 진동이 터진다.

TVD 조건: Harten의 부등식#

Harten(1983)이 제시한 탈출구는 "총 변동이 줄어드는" 스킴이다.

\mathrm{TV}(q) = \sum_i |q_{i+1} - q_i|, \qquad \mathrm{TV}(q^{n+1}) \le \mathrm{TV}(q^n)

좌변은 수치 해의 출렁임 총합. 매 스텝마다 이 값이 커지지 않으면 새 극값(overshoot, undershoot)이 생길 수 없다. 이 조건을 만족하는 스킴을 TVD 라 부른다.

Godunov의 벽을 우회하려면 스킴을 비선형 으로 만들어야 한다. 해의 매끈함에 따라 알고리즘이 모양을 바꾸는 것. 매끈한 영역에서는 2차, 급격한 구배 근처에서는 1차 업윈드. 이 스위칭을 담당하는 함수가 flux limiter다.

Flux Limiter 한 줄로 이해하기#

Lax–Wendroff 플럭스는 1차 업윈드에 보정항이 붙은 모양이다.

f_{i-1/2} = u\,q_{i-1} + \tfrac{1}{2}\,|u|\!\left(1 - \left|\tfrac{u\Delta t}{\Delta x}\right|\right) \varphi(r)\,(q_i - q_{i-1})

여기서 $u$ 는 이류 속도, $\Delta x$ 는 셀 폭, $\varphi(r)$ 이 리미터, $r$ 은 앞뒤 기울기 비.

r_{i-1/2} = \frac{q_{i-1} - q_{i-2}}{q_i - q_{i-1}}

해가 매끈하면 $r \approx 1$ , 극값 근처에서는 $r < 0$ 또는 $r \gg 1$ . 리미터가 이 $r$ 을 보고 보정 강도를 조절한다.

특별한 선택으로

$\varphi(r) = 0$ : donor-cell(1차 업윈드)
$\varphi(r) = 1$ : Lax–Wendroff
$\varphi(r) = r$ : Beam–Warming

Sweby 도식에서 $0 \le \varphi(r) \le \min(2r,\,2)$ 를 만족하면 TVD가 보장된다. 이 사다리꼴 영역 안에 들어오는 대표 리미터 넷을 비교한다.

리미터	정의	특징
minmod	$\max(0,\min(1,r))$	가장 조심스러움. 계단은 잘 보존, 매끈한 영역 약간 둔함
superbee	$\max(0,\min(2r,1),\min(r,2))$	가장 날카로움. 매끈한 파동을 계단처럼 깎기도
van Leer	$(r+	r
MC	$\max(0,\min(\tfrac{1+r}{2},2,2r))$	기본값으로 가장 무난

Python으로 직접 비교#

roll 기반 벡터화로 40줄 안에 끝난다. 격자 200셀, $\mathrm{CFL} = 0.5$ , 500 스텝(주기 경계, 파가 도메인을 두 바퀴 반 도는 거리).

import numpy as np
 
def limiter_phi(name, r):
    r = np.where(np.isfinite(r), r, 0.0)
    if name == 'donor':    return np.zeros_like(r)
    if name == 'lw':       return np.ones_like(r)
    if name == 'minmod':   return np.maximum(0, np.minimum(1, r))
    if name == 'superbee':
        return np.maximum.reduce([np.zeros_like(r),
                                  np.minimum(2*r, 1),
                                  np.minimum(r, 2)])
    if name == 'vanleer':  return (r + np.abs(r)) / (1 + np.abs(r) + 1e-30)
    if name == 'mc':
        return np.maximum(0, np.minimum.reduce([0.5*(1+r),
                                                2*np.ones_like(r),
                                                2*r]))
 
def advect_limited_fvm(q, cfl, name):
    """1D FV 이류, 주기 경계, u=+1."""
    dq  = q - np.roll(q, 1)                  # face i-1/2: q_i - q_{i-1}
    r   = np.roll(dq, 1) / (dq + 1e-30)      # (q_{i-1}-q_{i-2})/(q_i-q_{i-1})
    phi = limiter_phi(name, r)
    flx = np.roll(q, 1) + 0.5 * phi * (1 - cfl) * dq
    return q - cfl * (np.roll(flx, -1) - flx)
 
N, cfl, steps = 200, 0.5, 500
x  = np.arange(N) / N
q0 = ((x > 0.2) & (x < 0.4)).astype(float)
for name in ['lw', 'minmod', 'superbee', 'vanleer', 'mc']:
    q = q0.copy()
    for _ in range(steps):
        q = advect_limited_fvm(q, cfl, name)
    l1 = np.abs(q - q0).sum() / N
    print(f'{name:<9s}  L1={l1:.4f}  overshoot={q.max()-1:.3f}  undershoot={q.min():.3f}')

실행 결과(반복 실험에서 대표값):

리미터	$L_1$ 오차	오버슈트	언더슈트
lax-wendroff	0.056	+0.181	−0.133
minmod	0.081	0.000	0.000
superbee	0.042	+0.002	0.000
van Leer	0.061	0.000	0.000
MC	0.051	0.000	0.000

lax-wendroff만 오버슈트가 20% 가까이 튀고, 나머지는 모두 0. 오차 순위도 재미있다. superbee가 $L_1$ 기준으로 가장 작다. 그만큼 계단을 덜 번지게 유지한다는 뜻인데, 단점은 매끈한 sine 파에서 파두를 계단 모양으로 왜곡한다는 점이다.

매끈한 파동 테스트( $q_0 = \sin(2\pi x)$ , 500 스텝)를 따로 돌려 보면 수렴 차수가 minmod에서 약 1.5, van Leer/MC에서 약 2.0, superbee는 약 1.3 수준으로 떨어진다. "2차 정확" 라벨은 극값 근처에서 리미터가 켜질 때마다 부분적으로 1차로 떨어지기 때문이다.

브라우저에서 직접 리미터 바꿔보기#

아래 시뮬레이션에서 리미터 버튼을 누르고 CFL 슬라이더를 움직여 보자. 주황색이 수치 해, 시안 점선이 정확 해다.

CFL0.50steps: 0

lax-wendroff로 바꾸면 계단 뒤에 위로 튀는 혹, 앞에 아래로 파인 홈이 분명히 보인다. superbee로 돌리면 계단이 원래보다 날카로워 보이기까지 한다. minmod는 양쪽 모두 부드럽고 혹이 없다. CFL을 0.9 가까이 올리면 모든 리미터가 과도해지는 걸 볼 수 있다 — 안정 조건을 만족해도 수치 확산과 상호작용한다.

실무 체크리스트#

shock capturing이 목표면 minmod 또는 van Leer부터. 진동 제로가 우선이다.
LES/대규모 난류에는 superbee 위험. 매끈한 난류 구조를 계단처럼 깎아 에너지 스펙트럼을 왜곡한다. van Leer 또는 MC가 타협.
경계 셀에서 스텐실이 부족해 임시로 1차 업윈드를 쓰면 전체 정확도가 1차로 떨어진다. 경계에서도 리미터 stencil을 확보하거나 고스트 셀을 써라.
OpenFOAM의 limitedLinear 1 은 van Leer 계열의 1-매개변수 리미터다. 계수를 0으로 낮추면 1차 업윈드, 1.0이면 풀 TVD 영역 사용.
수렴 차수 검증은 반드시 매끈한 해 로. 불연속이 포함된 테스트에서 $L_1$ 수렴 기울기를 측정하면 리미터가 아닌 불연속 자체가 차수를 떨어뜨린다.

코드 짤 때 빠지지 말 함정#

Godunov 정리 때문에 선형 2차 단조 이류 스킴은 존재할 수 없다. 비선형 스위칭이 필수.
Flux limiter $\varphi(r)$ 는 매끈한 영역에서 2차, 극값 근처에서 1차로 자동 전환되는 한 줄짜리 함수다.
minmod는 안전, superbee는 날카롭게, van Leer/MC는 균형. LES와 shock capturing은 선택 기준이 다르다.