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

φ(r)이 Sweby 영역 안에 있으면 계단의 새 극값은 사라진다. 매끈한 파의 차수는 리미터마다 갈린다.

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

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

리미터 버튼은 lax-wendroff에서 시작한다. CFL 슬라이더를 0.5에 두면 주황 곡선이 계단 앞뒤로 혹을 만들고, superbee로 바꾸면 그 혹이 걷힌다.

CFL0.50steps: 0

주황이 수치 해, 시안 점선이 정확 해다. minmod는 양쪽이 무디고 혹은 없다. superbee는 계단이 원래보다 날카로워 보이기까지 한다. CFL을 0.9 근처로 올리면 혹이 없는 리미터도 더 퍼진다. 안정 조건 안에서도 수치 확산과 맞물린다.

진동의 정체 — 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가 보장된다. 이 사다리꼴 영역 안에 들어오는 대표 리미터 넷을 비교한다. $\varphi(r) = 0$ 으로 끈 자리는 1차 업윈드다. 그 업윈드가 진동 대신 무엇을 무디게 하는지는 FTCS가 떨릴 때의 맞교환에 있다.

리미터	정의	특징
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차로 떨어지기 때문이다.

r=1에서는 네 리미터가 같은 φ를 준다#

표의 정의에 $r=1$ 을 넣으면 minmod, superbee, van Leer, MC가 모두 1이 된다. $\varphi=1$ 은 Lax–Wendroff다. 매끈한 한가운데에서는 네 리미터가 같은 보정을 쓴다. 차이는 $r$ 이 1에서 벗어날 때 난다.

$r=2$ 를 넣어 보자. minmod는 1에서 멈춘다. superbee는 2로 Sweby 상단에 붙는다. van Leer는 $4/3$ 이다. MC는 1.5다. superbee가 계단의 $L_1$ 을 0.042로 가장 작게 유지하는 이유가 이 상단이다. 같은 선택이 sine 테스트의 차수를 약 1.3까지 내린다. van Leer와 MC는 약 2.0에 남고, minmod는 약 1.5다. $r$ 이 0 이하인 극값에서는 넷 다 0으로 떨어져 donor-cell이 된다.

코드의 보정 계수도 이미 정해져 있다. 속도가 1이고 $\mathrm{CFL}=0.5$ 이면 점프 앞에 붙는 계수는 0.25다. $\varphi=1$ 이어도 이웃 점프의 4분의 1만 더한다. 500스텝이 전진하는 거리는 $0.5\times 500=250$ 셀이다. 격자 200칸이면 주기 경계를 1.25바퀴 돈다. 오버슈트 +0.181은 파고 1의 18퍼센트다. 언더슈트 −0.133은 바닥이 음수로 내려앉은 깊이다. minmod의 $L_1$ 0.081은 superbee 0.042의 약 1.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$ 수렴 기울기를 측정하면 리미터가 아닌 불연속 자체가 차수를 떨어뜨린다.

계단에 superbee를 쓰면 sine은 1.3차로 내려앉는다#

표에서 superbee의 $L_1$ 은 0.042로 가장 작다. 오버슈트는 +0.002라 계단에서는 거의 안 보인다. 같은 리미터를 sine에 얹으면 차수가 약 1.3까지 떨어진다. van Leer와 MC는 약 2.0에 남는다. 불연속을 지키는 선택이 매끈한 파의 차수를 깎는다.

스칼라가 밀도라면 lax-wendroff의 언더슈트 −0.133은 음수 밀도의 입구다. 그 구멍을 막는 한정자는 양수가 새는 자리에 따로 있다. 수동 스칼라의 계단이라면 minmod나 van Leer로 충분하고, 밀도라면 그 다음 문을 연다.