비정렬 격자에서 2차 재구성이 만든 가짜 극값 — Barth–Jespersen·Venkatakrishnan 리미터

잘 돌던 해석기다. 2차만 켜자 밀도가 음수로 떨어졌다. 충격파 앞의 면 값이 셀 평균보다 위로 솟아 있었다. 1차 업윈드는 같은 격자에서 멀쩡했다. 정확도를 올린 쪽이 계산을 죽였다. 비정렬 셀 안의 직선이, 이웃에 없는 극값을 만들기 때문이다. Barth–Jespersen은 그 극값을 잘라 낸다. Venkatakrishnan은 자른 자리를 둥글린다.

재구성이 없던 극값을 만든다#

유한체적법(FVM, 셀 평균을 미지수로 두는 이산화)에서 2차는 셀 안의 직선이다.

$u_i$는 셀 $i$의 평균이다. $\nabla u_i$는 재구성된 기울기, $\mathbf{x}_f$는 면 중심, $\Phi_i\in[0,1]$은 리미터다.

$\Phi_i=1$이면 리미터는 꺼져 있다. 기울기는 최소제곱이나 Green–Gauss로 이웃을 평균해 만든다. 불연속 옆에서 그 평균은 과하다. 직선을 면까지 뻗으면 이웃에 없던 최대·최소가 생긴다.

기울기 gain을 올리면 unlimited 재구성(Φ=1)이 회색 밴드, 곧 이웃의 최대·최소를 뚫고 빨강으로 바뀐다.

밴드 안에 선이 머물면 그 셀의 Φ는 1이다. 점프 옆에서만 선분이 먼저 짧아진다. 같은 gain에서 리미터를 Barth–Jespersen으로 바꾸면 선은 밴드 밖으로 나가지 않는다.

가짜 극값은 다음 스텝의 플럭스를 오염시킨다. 진동이 자라는 끝은 음의 밀도와 음의 에너지다. 재구성된 면 값은 이웃이 이미 가진 범위를 넘지 않는다. 국소 최대·최소 원리(LMP)다. $\Phi_i$는 그 범위를 지키려고 기울기를 깎는 스칼라다.

Barth–Jespersen: 이웃의 최대·최소로 가둔다#

Barth와 Jespersen(1989)은 허용 범위를 먼저 고정한다. 셀 $i$와 이웃이 만드는 최대·최소다.

$N(i)$는 셀 $i$의 면 이웃 집합이다. 각 면에서 리미터 없는 증분 $u_f-u_i$를 보고, 면마다 계수 $\phi_f$를 구한다.

셀 하나의 리미터는 모든 면 가운데 가장 보수적인 값이다.

범위 안에 있으면 $\phi_f=1$이라 2차가 남는다. 넘치려는 면은 경계에 딱 닿는 기울기만 남긴다. 극값 옆에서는 $\Phi_i\to 0$이 되어 1차 업윈드로 돌아간다.

양의 증분만 읽으면 $\phi_f=\min(1,(u_i^{\max}-u_i)/(u_f-u_i))$다. 분자가 분모보다 작아지는 순간, 면 값은 $u_i^{\max}$에 붙어 멈춘다. 위 그림의 빨간 선이 회색 경계에서 꺾이는 이유가 이 비율이다.

미분 불가능이라는 두 번째 함정#

Barth–Jespersen은 단조성을 지킨다. 정상상태 해석기에 넣으면 잔차가 어느 선에서 더 내려가지 않고 떨린다. 원인은 $\min$과 나눗셈이다. $\phi_f$는 해에 대해 미분되지 않는다. 어느 면이 최소를 주는지가 반복마다 바뀐다. 리미터가 켜졌다 꺼지면 잔차가 튕긴다. 암시적(implicit) 솔버의 야코비안은 이 불연속을 싫어한다.

매끄러운 곳의 Φ는 1이면 된다. 진짜 불연속에서만 켜져야 하고, 켜지는 경계는 부드러워야 한다. 날카로운 꺾임을 둥글리는 식이 다음이다.

Venkatakrishnan: 부드럽게 깎는다#

Venkatakrishnan(1993)은 $\min(1,y)$의 꺾임을 유리함수로 바꿨다. 소스가 언급한 Michalak–Ollivier-Gooch 리미터도 같은 계열이다. 면마다의 계수는 이렇게 쓴다.

$\Delta_-=u_f-u_i$는 리미터 없는 증분이다. $\Delta_+$는 증분의 부호에 따라 $u_i^{\max}-u_i$ 또는 $u_i^{\min}-u_i$다. 손으로 고르는 양은 $\epsilon^2$다.

$K$는 튜닝 파라미터이고 $\Delta x$는 격자 크기다. 변동이 $\epsilon$보다 작으면 $\phi_f\to 1$이라 리미터가 꺼진다. $K$를 키우면 문턱이 올라가고, 리미터가 풀리는 영역이 넓어진다. 너무 크면 충격파 옆 진동을 놓친다. $K$를 0으로 보내면 Barth–Jespersen으로 돌아간다.

코드는 셀 크기를 1로 두고 $\epsilon^2=K^3$으로 계산한다. 기본값 $K=0.3$이면 $\epsilon^2=0.027$이다. $\epsilon=0$이고 $\Delta_+=\Delta_-=\Delta$이면, 재구성이 허용 한계에 정확히 닿는 자리에서 유리함수는 $(3\Delta^2)/(4\Delta^2)=3/4$가 된다. 같은 자리에서 Barth–Jespersen은 $\phi_f=1$이다. Venkatakrishnan은 한계에 닿기 전에 기울기의 4분의 1을 버린다. 모서리가 한 단계 일찍 뭉툭해지는 이유다.

관행으로 적힌 구간은 $K\approx 0.1$에서 $K\approx 5$다. 코드와 같은 $\epsilon^2=K^3$에 이 끝을 넣으면 0.001과 125가 나온다. 문턱의 비는 12만 5천이다. $K$를 조금 올려도 세제곱이라, 매끄럽다고 면제되는 셀이 급히 늘어난다. 반대로 $K$를 0.1 아래에 두면 $\epsilon^2$는 사실상 0이고, Barth–Jespersen의 꺾임이 그대로 돌아온다.

위 그림에서 Venkatakrishnan을 켜고 $K$를 움직이면 Φ가 계단 없이 미끄러진다.

Python — 스칼라 이류로 세 리미터를 겨룬다#

주기 경계다. 사각파와 가우스 봉우리를 함께 이류(advection)시킨다. 무제한(Fromm), Barth–Jespersen, Venkatakrishnan을 나란히 두고 마지막 최솟값을 본다. 0보다 작으면 그만큼 가짜 골짜기가 팬 것이다.

셀 기울기는 $s=(u_{i+1}-u_{i-1})/2$다. 면으로 나가는 증분은 $0.5s$라서, 좌우 이웃의 차이가 1이면 $s=1/2$이고 면 증분은 0.25다. 그 0.25가 허용 폭보다 크면 Φ는 1 아래로 내려간다. 1차원 중앙 차분은 좌우 이웃만 본다. 비정렬 셀에서 같은 기울기의 가중치는 이동 최소제곱이 셀마다 다시 맞추는 방식과 다르다. 리미터는 기울기가 만들어진 뒤에 붙는다.

import numpy as np
 
NX, A, CFL = 200, 1.0, 0.4
dx = 1.0 / NX
dt = CFL * dx / A
 
def init_profile():
    x = (np.arange(NX) + 0.5) / NX
    u = np.where((x > 0.1) & (x < 0.3), 1.0, 0.0)   # 사각파
    u += np.exp(-((x - 0.65) / 0.06) ** 2)          # 가우스 봉우리
    return u
 
def cell_slope(u):
    return (np.roll(u, -1) - np.roll(u, 1)) / 2.0    # 중심 기울기(Fromm)
 
def barth_jespersen_phi(u, s):
    up, um = np.roll(u, -1), np.roll(u, 1)
    umax = np.maximum(u, np.maximum(up, um))
    umin = np.minimum(u, np.minimum(up, um))
    phi = np.ones_like(u)
    for du in (0.5 * s, -0.5 * s):                   # 두 면
        f = np.ones_like(u)
        pos, neg = du > 1e-12, du < -1e-12
        f[pos] = np.minimum(1.0, (umax[pos] - u[pos]) / du[pos])
        f[neg] = np.minimum(1.0, (umin[neg] - u[neg]) / du[neg])
        phi = np.minimum(phi, f)
    return np.clip(phi, 0.0, 1.0)
 
def venkatakrishnan_phi(u, s, K=0.3):
    up, um = np.roll(u, -1), np.roll(u, 1)
    umax = np.maximum(u, np.maximum(up, um))
    umin = np.minimum(u, np.minimum(up, um))
    eps2 = K ** 3                                    # (K*h)^3, 셀 크기 h=1 단위
    phi = np.ones_like(u)
    for du in (0.5 * s, -0.5 * s):
        d = np.where(du > 0, umax - u, umin - u)
        num = (d * d + eps2) * du + 2 * du * du * d
        den = d * d + 2 * du * du + d * du + eps2
        f = np.where(np.abs(du) < 1e-12, 1.0, num / (du * den))
        phi = np.minimum(phi, f)
    return np.clip(phi, 0.0, 1.0)
 
def muscl_rhs(u, limiter):
    s = cell_slope(u)
    phi = limiter(u, s) if limiter else np.ones_like(u)
    uL = u + 0.5 * phi * s              # 면 i+1/2 좌측 상태
    flux = A * uL                       # 업윈드 플럭스 (a > 0)
    return -(flux - np.roll(flux, 1)) / dx
 
def advance_muscl(u, limiter):         # SSP-RK2
    k1 = muscl_rhs(u, limiter)
    u1 = u + dt * k1
    k2 = muscl_rhs(u1, limiter)
    return 0.5 * (u + u1 + dt * k2)
 
def run_limiter_race(steps=160):
    fields = {"none": init_profile(), "bj": init_profile(), "venk": init_profile()}
    lims = {"none": None, "bj": barth_jespersen_phi, "venk": venkatakrishnan_phi}
    for _ in range(steps):
        for k in fields:
            fields[k] = advance_muscl(fields[k], lims[k])
    for k, u in fields.items():
        print(f"{k:5s}  min={u.min():+.4f}  max={u.max():.4f}")
 
run_limiter_race()

출력은 이렇다.

none   min=-0.0417  max=1.0231
bj     min=+0.0000  max=1.0000
venk   min=-0.0004  max=1.0009

무제한의 최솟값 -0.0417은 높이 1인 사각파가 0 아래로 파인 깊이다. 최댓값 1.0231은 같은 높이 위로 넘친 양이다. Barth–Jespersen은 최소·최대를 $[0,1]$에 고정했다. Venkatakrishnan의 최솟값 -0.0004는 그 골짜기의 약 104분의 1이다. 1을 넘는 폭은 0.0231에서 0.0009로 준다.

160스텝의 이동 거리는 $A \cdot dt \cdot 160$이다. $dx=1/200$이고 CFL이 0.4이면 $dt=0.002$라서 거리는 0.32다. 200셀 중 64칸이다. 사각파 구간 $(0.1, 0.3)$은 $(0.42, 0.62)$로 가고, 가우스 중심 0.65는 0.97에 도착한다. 둘은 같은 속도로 움직이므로 상대 위치는 유지된다. 출력의 최소·최대는 두 파가 부딪힌 결과가 아니다. 각 모서리에서 따로 생긴 값이다. 사각파 길이 0.2는 40셀이고, 가우스 식의 스케일 0.06은 12셀이다. 최소·최대를 가른 것은 그 모서리 근처의 몇 칸이다.

사각파 뒤 모서리만 보면 된다. 무제한(빨강)은 0 아래로 잔물결이 자란다. 두 리미터는 단조롭게 남는다. Venkatakrishnan(노랑)은 Barth–Jespersen(시안)보다 모서리가 조금 더 뭉툭하다. $\epsilon=0$에서 기울기가 3/4만 남던 계산과 같은 방향이다. 코드의 $K=0.3$은 그 꺾임을 없애지 않고, 초과만 0.0009 수준으로 줄인다.

현장에서 리미터를 켤 때#

이웃을 어떻게 묶느냐가 허용 범위를 가른다. Barth–Jespersen의 $\max/\min$을 면 이웃으로 잡느냐, 정점(vertex) 이웃으로 잡느냐에 따라 결과가 달라진다. 비정렬에서 정점 이웃, 곧 꼭짓점을 공유하는 모든 셀을 쓰면 방향 편향이 줄어 더 균형 잡힌다.

$K$는 그다음이다. 너무 작으면 Barth–Jespersen처럼 수렴이 정체된다. 너무 크면 충격파에서 진동을 놓친다. $\epsilon^2$에 $\Delta x^3$이 들어가므로, 격자를 잘게 쪼개면 같은 $K$에서도 리미터가 더 자주 켜진다.

벡터계(Euler 방정식)에서는 보존 변수마다 리미터를 따로 걸면 성분 사이가 어긋나 새 진동이 생긴다. 특성 변수(characteristic variable)로 사영한 뒤 거는 편이 안전하다. 비용은 더 든다.

잔차가 멈추면 기울기부터 의심하지 않는다#

2차 해석기가 음의 밀도로 죽으면, 면 값이 이웃 범위 밖인지를 먼저 본다. 밖이면 기울기 루틴이 아니라 Φ를 만든 이웃 집합이 문제다. 면 이웃만 쓰면 대각 방향의 극값을 놓친다. 재구성은 그 틈으로 다시 나간다.

정상상태에서 잔차만 제자리면 $K$를 연다. $K=0$은 Barth–Jespersen의 꺾임이라, 야코비안이 반복마다 점프한다. 점프에서 Φ가 0으로 내려가면 그 셀의 차수는 1차에 붙는다. 매끄러운 영역은 2차 근처에 남는가. 남는다면 $\epsilon^2$가 그 영역을 놓아 준 것이다. 이 갈림을 한 계산에서 동시에 재면 어떤 숫자가 나오는지, 충격파와 매끈한 와류를 붙여 차수를 잰 기록이 아직 열려 있다.