질량은 정확히 보존됐는데 밀도가 음수다 — 양수 보존 flux 리미터

보존형 스킴은 질량을 한 톨도 잃지 않는다. 면을 통해 나간 양이 곧 이웃이 받은 양이니, 도메인 전체 합은 기계 정밀도까지 일정하다. 그런데 그 스킴이 밀도 −0.003을 만든다. 총합은 맞는데 개별 셀이 음수다. 오늘은 왜 보존성이 양수성을 보장하지 않는지, 그리고 고차 정확도를 거의 그대로 유지한 채 부호만 지켜내는 flux 리미터를 어떻게 짜는지 다룬다. 실제로 돌려서 몇 %의 면이 실제로 손을 타는지까지 세어 본다.

발산이 아니라 정의역 밖으로 나간 것이다#

로그가 NaN으로 끝났을 때 가장 먼저 의심하는 것은 시간 스텝이다. CFL을 절반으로 줄인다. 그래도 죽는다. 격자를 세분화한다. 더 빨리 죽는다.

이때는 안정성 문제가 아닐 가능성이 높다. 음속을 구하는 줄을 보자.

$\gamma$는 비열비, $p$는 압력, $\rho$는 밀도다. $p$나 $\rho$ 중 하나가 음수가 되는 순간 $c$는 NaN이 된다. 다음 스텝의 시간 간격도 NaN, 그다음 flux도 전부 NaN. 실제 사고는 NaN이 찍힌 줄보다 한 스텝 앞에서 이미 끝나 있다.

핵심은 이것이다. Euler 방정식의 해가 살아 있으려면 보존변수 $\mathbf{U} = (\rho,\ \rho u,\ E)$가 허용 집합(admissible set) 안에 있어야 한다.

보존성은 "총합이 유지된다"는 성질이지 "각 셀이 $G$ 안에 남는다"는 성질이 아니다. 둘은 완전히 다른 요구다.

고차 재구성은 부호를 지켜 주지 않는다#

$G$를 벗어나는 대표적인 상황은 진공 근처다. 강한 팽창파, 고고도 재진입 유동, 캐비테이션 기포 내부, 그리고 폭발파 뒤쪽. 이런 곳에서 $\rho$와 $p$는 $10^{-4}$ 수준까지 떨어진다.

여기에 MUSCL이나 WENO 재구성을 얹으면 셀 평균 $\bar\rho_i$가 양수여도 면 값 $\rho_{i+1/2}^{L}$은 음수가 될 수 있다. 기울기가 셀 폭의 절반을 곱해 더해지기 때문이다. 재구성은 셀 평균만 보존할 뿐, 부호에는 아무 약속을 하지 않는다.

압력은 더 나쁘다. $p$는 보존변수의 비선형 함수다. $\rho$와 $E$가 각각 양수여도 운동에너지 $(\rho u)^2/2\rho$가 $E$를 넘으면 $p$는 음수가 된다. 실제로 뒤에서 돌려 볼 이중 희박파에서는 밀도가 0.37로 멀쩡한 상태에서 압력이 먼저 −0.163으로 떨어진다.

Lax–Friedrichs가 CFL 0.5에서 버티는 이유#

그럼 무엇을 믿을 수 있는가. 1차 Lax–Friedrichs(LF) flux다.

$\alpha = \max_j(|u_j| + c_j)$는 최대 특성속도다. 이 flux로 업데이트를 전개하면 다음처럼 정리된다.

$\lambda = \Delta t/\Delta x$다. 세 계수의 합은 정확히 1이고, $\lambda\alpha \le 1$이면 모두 음수가 아니다. 즉 볼록 결합이다.

$G$는 볼록 집합이므로, 결합에 들어가는 세 항이 모두 $G$ 안이면 결과도 $G$ 안이다. 문제는 $\mathbf{U} \pm \mathbf{F}/\alpha$가 다시 $G$에 들어오느냐인데, Perthame–Shu가 이 조건이 $\lambda\alpha \le 1/2$에서 성립함을 보였다. 흔히 쓰는 "LF는 CFL 0.5에서 양수 보존"이라는 문장의 출처가 여기다.

정리하면 우리 손에는 두 개의 flux가 있다. 정확하지만 부호를 안 지키는 고차 flux $\hat{\mathbf{F}}^{H}$, 그리고 부정확하지만 부호를 지키는 $\hat{\mathbf{F}}^{LF}$.

두 flux를 잇는 선분 위에서 θ 찾기#

Hu, Adams, Shu(2013)의 아이디어는 단순하다. 둘을 섞되, 섞는 비율을 면마다 따로 정한다.

$\theta = 1$이면 순수 고차 스킴, $\theta = 0$이면 LF로 후퇴한다. 어떤 $\theta$를 써도 flux 형식이므로 보존성은 자동으로 유지된다. 리미터가 질량을 만들거나 없애지 않는다는 뜻이다. 이 점이 인공 점성을 국소적으로 들이붓는 방식과 결정적으로 다르다.

$\Delta\mathbf{F}_{i+1/2} \equiv \hat{\mathbf{F}}^{H} - \hat{\mathbf{F}}^{LF}$로 쓰면 업데이트는 이렇게 갈라진다.

첫 항은 CFL 조건만 지키면 $G$ 안에 있음이 보장된다. 나머지는 그 안전한 점에서 고차 해 쪽으로 뻗는 선분이다. 우리가 할 일은 선분이 $G$를 벗어나기 직전에서 멈추는 것뿐이다.

바닥값은 0이 아니라 아주 작은 양수 $\varepsilon$으로 잡는다.

$\rho_j^0$, $p_j^0$는 초기 상태의 밀도·압력이다. 정확히 0을 목표로 하면 반올림 한 번에 다시 음수가 된다.

예산을 반씩 나눠 쓴다 — 면 하나가 셀 둘을 건드린다#

밀도는 보존변수 자체라서 $\rho_i^{n+1}$이 $\theta$에 대해 1차 함수다. 그래서 조건을 대수적으로 풀 수 있다.

셀 $i$가 가진 여유분을 예산으로 부르자. $b_i = \rho_i^{LF} - \varepsilon \ge 0$이다. 면 $i+1/2$에서 $\Delta F^\rho > 0$이면 이 면은 왼쪽 셀 $i$의 밀도를 깎는다. 반대면 오른쪽 셀을 깎는다.

여기서 함정이 하나 있다. 내부 셀은 좌우 두 면에서 동시에 깎인다. 각 면이 자기가 깎는 셀의 예산 전부를 쓸 수 있다고 계산하면, 면 둘이 각각은 합법이면서 합쳐서 예산의 두 배를 써 버린다. 그래서 각 면에 절반씩만 허용한다.

아래 도식에서 직접 확인해 보자.

|dF| scale을 올리면 셀 2 양옆의 두 $\theta$가 노란색으로 내려간다. 여기서 "full budget per face"를 누르면 두 $\theta$가 다시 1 근처로 튀어 오르는데, 아래쪽 셀 2 막대는 초록 바닥선을 뚫고 내려간다. 면 단위로만 보면 아무 잘못이 없는데 셀 단위로는 위반이다.

압력은 밀도 다음에, 그리고 이분법으로#

밀도가 정리되면 압력을 본다. 순서가 중요하다. $p$를 계산하려면 $\rho$로 나눠야 하므로, $\rho > 0$이 먼저 확보되어야 한다.

$p$는 $\mathbf{U}$의 비선형 함수라 1차식으로 풀리지 않는다. 대신 좋은 성질이 하나 있다. $p(\mathbf{U})$는 $G$ 위에서 오목 함수(concave)다. $\mathbf{U}(\theta)$가 $\theta$에 대해 선형이므로 $p(\theta)$도 오목하다. 오목 함수의 상위 준위 집합 $\{\theta : p(\theta) \ge \varepsilon\}$은 구간이고, $\theta = 0$(LF 상태)에서 이미 조건을 만족하므로 그 구간은 0을 포함한다.

즉 근이 하나뿐이다. 이분법이 안전하게 작동한다. 20~40회면 배정밀도 한계까지 좁혀진다.

한 가지 주의할 점. 셀 $i$ 때문에 $\theta_{i\pm1/2}$를 줄이면 이웃 셀 $i\pm1$의 계산이 달라진다. 한 번 훑고 끝내면 드문 경우에 위반이 남는다. 위반 셀이 사라질 때까지 스윕을 반복해야 한다. $\theta \to 0$이면 LF 상태로 수렴하므로 반복은 반드시 끝난다.

Python으로 되살린 이중 희박파#

토이 문제는 이중 희박파다. 관 가운데를 기준으로 좌우가 $\mp u_0$로 서로 멀어진다. 가운데에 진공에 가까운 구멍이 열리고, 그 구멍이 고차 스킴을 죽인다.

import numpy as np
 
GAMMA = 1.4
 
 
def to_primitive(U):
    rho = U[0]
    u = U[1] / rho
    p = (GAMMA - 1.0) * (U[2] - 0.5 * rho * u * u)
    return rho, u, p
 
 
def euler_flux(U):
    rho, u, p = to_primitive(U)
    return np.array([rho * u, rho * u * u + p, (U[2] + p) * u])
 
 
def lf_face_flux(UL, UR, alpha):
    return 0.5 * (euler_flux(UL) + euler_flux(UR) - alpha * (UR - UL))
 
 
def density_theta(rho_lf, dF_rho, lam, eps_rho):
    """면마다 theta 결정. 각 면은 자기가 깎는 셀의 예산을 절반까지만 쓴다."""
    n = rho_lf.size
    budget = np.maximum(rho_lf - eps_rho, 0.0)
    theta = np.ones(dF_rho.size)
    for f in range(1, n):                      # 내부 면만
        d = dF_rho[f]
        if d > 0.0:                            # 왼쪽 셀을 깎는다
            cap = 0.5 * budget[f - 1] / (lam * d)
        elif d < 0.0:                          # 오른쪽 셀을 깎는다
            cap = 0.5 * budget[f] / (lam * (-d))
        else:
            cap = 1.0
        theta[f] = min(1.0, cap)
    return theta
 
 
def pressure_at(U_lf, dFl, dFr, tl, tr, lam):
    U = U_lf - lam * (tr * dFr - tl * dFl)
    return to_primitive(U)[2]
 
 
def pressure_theta(U_lf, dF, theta, lam, eps_p):
    """p(theta)는 오목하므로 안전 구간이 하나. 이분법으로 경계를 찾는다.
    한 셀의 theta를 줄이면 이웃이 영향을 받으므로 위반이 없어질 때까지 반복."""
    n = U_lf.shape[1]
    for _ in range(20):
        dirty = False
        for i in range(n):
            tl, tr = theta[i], theta[i + 1]
            if pressure_at(U_lf[:, i], dF[:, i], dF[:, i + 1], tl, tr, lam) >= eps_p:
                continue
            dirty = True
            lo, hi = 0.0, 1.0
            for _ in range(40):
                mid = 0.5 * (lo + hi)
                ok = pressure_at(U_lf[:, i], dF[:, i], dF[:, i + 1],
                                 tl * mid, tr * mid, lam) >= eps_p
                lo, hi = (mid, hi) if ok else (lo, mid)
            theta[i] *= lo
            theta[i + 1] *= lo
        if not dirty:
            return theta
    return theta

시간 전진 루프는 매 스텝 두 개의 flux를 만들고, $\theta$를 구하고, 섞어서 업데이트한다.

def minmod(a, b):
    return np.where(a * b <= 0.0, 0.0, np.where(np.abs(a) < np.abs(b), a, b))
 
 
def face_states(U):
    """MUSCL-minmod 재구성 -> 각 면의 좌/우 상태"""
    d = minmod(U[:, 1:-1] - U[:, :-2], U[:, 2:] - U[:, 1:-1])
    s = np.zeros_like(U)
    s[:, 1:-1] = d
    return U[:, :-1] + 0.5 * s[:, :-1], U[:, 1:] - 0.5 * s[:, 1:]
 
 
def march_double_rarefaction(n=200, cfl=0.45, u0=4.0, t_end=0.15, limiter=True):
    dx = 1.0 / n
    x = (np.arange(n) + 0.5) * dx
    rho = np.ones(n)
    u = np.where(x < 0.5, -u0, u0)
    p = np.full(n, 0.4)
    U = np.vstack([rho, rho * u, p / (GAMMA - 1.0) + 0.5 * rho * u * u])
 
    eps = min(1e-13, rho.min(), p.min())
    t, step, clipped, total = 0.0, 0, 0, 0
 
    while t < t_end:
        r, v, pr = to_primitive(U)
        if r.min() <= 0.0 or pr.min() <= 0.0:                 # 허용 집합 이탈
            return dict(crashed=True, step=step,
                        rho_min=r.min(), p_min=pr.min())
        a = np.sqrt(GAMMA * pr / r)
        alpha = float(np.max(np.abs(v) + a))
        dt = min(cfl * dx / alpha, t_end - t)
        lam = dt / dx
 
        Ug = np.hstack([U[:, :1], U, U[:, -1:]])              # zero-gradient ghost
        UL, UR = face_states(Ug)
        Flow = np.zeros((3, n + 1))
        Fhigh = np.zeros((3, n + 1))
        for f in range(n + 1):
            Flow[:, f] = lf_face_flux(Ug[:, f], Ug[:, f + 1], alpha)
            Fhigh[:, f] = lf_face_flux(UL[:, f], UR[:, f], alpha)
 
        dF = Fhigh - Flow
        U_lf = U - lam * (Flow[:, 1:] - Flow[:, :-1])          # 안전한 기준점
 
        if limiter:
            th = density_theta(U_lf[0], dF[0], lam, eps)
            th = pressure_theta(U_lf, dF, th, lam, eps)
            clipped += int(np.sum(th[1:n] < 1.0 - 1e-12))
            total += n - 1
        else:
            th = np.ones(n + 1)
 
        F = Flow + th * dF                                     # 섞인 flux
        U = U - lam * (F[:, 1:] - F[:, :-1])
        t += dt
        step += 1
 
    r, _, pr = to_primitive(U)
    return dict(crashed=False, step=step, rho_min=float(r.min()),
                p_min=float(pr.min()), clipped=100.0 * clipped / max(total, 1))
 
 
for lim in (False, True):
    o = march_double_rarefaction(limiter=lim)
    tag = "limiter ON " if lim else "limiter OFF"
    if o["crashed"]:
        print(f"{tag}: crashed at step {o['step']}  "
              f"rho_min={o['rho_min']:.4f}  p_min={o['p_min']:+.4f}")
    else:
        print(f"{tag}: reached t=0.15 in {o['step']} steps  "
              f"rho_min={o['rho_min']:.3e}  p_min={o['p_min']:.3e}  "
              f"theta<1 on {o['clipped']:.3f}% of faces")

$u_0 = 4$, 격자 200개, CFL 0.45에서 출력은 이렇다.

limiter OFF: crashed at step 2  rho_min=0.3698  p_min=-0.1630
limiter ON : reached t=0.15 in 300 steps  rho_min=9.560e-04  p_min=5.140e-04  theta<1 on 0.027% of faces

리미터 없이는 두 번째 스텝에서 끝난다. 주목할 것은 그 순간의 밀도가 0.3698이라는 점이다. 밀도는 전혀 위험해 보이지 않는데 압력이 먼저 음수로 떨어졌다. 밀도만 감시하는 코드는 이 사고를 놓친다.

아래 시뮬레이션에서 직접 조작해 보자.

limiter OFF를 누르고 pull-apart u0를 3.5 이상으로 올리면, 위쪽 밀도 곡선이 아직 멀쩡한 상태에서 가운데 압력 곡선이 빨간 선을 뚫는다. limiter ON으로 되돌리면 맨 아래 $\theta$ 막대 중 극소수만 초록에서 내려오고 계산은 끝까지 간다.

θ가 1보다 작아진 면은 전체의 몇 %인가#

측정값은 0.027%다. 200셀 × 300스텝 = 약 6만 개의 면 중 열여섯 개 정도만 손을 탔다. $u_0$를 8까지 올려도 0.093%에 그친다.

이 숫자가 이 기법의 핵심이다. 리미터는 사실상 잠들어 있다. 진공이 열리는 몇 개의 셀, 몇 개의 스텝에서만 깨어나서 그 면의 flux를 LF 쪽으로 살짝 당긴다. 나머지 99.97%의 면에서는 원래 고차 스킴이 그대로 돌아간다.

인공 점성을 전역으로 올려서 문제를 덮는 방식과 비교하면 차이가 분명하다. 그쪽은 매끄러운 영역의 정확도까지 함께 깎는다. 이쪽은 $\theta = 1$이 유지되는 한 원래 스킴과 비트 단위로 동일하다. 수렴 차수를 재는 격자 수렴 시험에서도 리미터가 차수를 떨어뜨리지 않는 이유가 이것이다.

코드에 넣기 전 확인할 세 가지#

첫째, CFL 상한을 실제로 지키고 있는가. 이 기법 전체가 "$\mathbf{U}^{LF}$는 안전하다"는 전제 위에 서 있다. LF의 양수 보존은 $\lambda\alpha \le 1/2$에서만 성립한다. 평소 CFL 0.8로 돌리던 코드라면 기준점부터 이미 무너져 있어서 $\theta$를 0으로 내려도 살아나지 않는다. 리미터가 안 듣는다면 이것부터 의심하자.

둘째, Runge–Kutta 각 스테이지마다 적용하고 있는가. SSP-RK는 각 스테이지가 전진 Euler의 볼록 결합이다. 스테이지 하나하나가 $G$ 안에 있어야 최종 결과도 $G$ 안에 들어온다. 마지막에 한 번만 검사하면 중간 스테이지에서 이미 $\sqrt{}$ 안이 음수가 된 뒤다.

셋째, $\varepsilon$을 0으로 두지 않았는가. 정확히 0을 목표로 이분법을 돌리면 마지막 반올림에서 $-10^{-17}$이 나온다. 초기 최소값 기준의 작은 양수를 바닥으로 깔아야 한다.

다시 진공을 만났을 때 꺼낼 것#

보존성과 양수성은 다른 성질이다. 전자는 flux 형식이 공짜로 주고, 후자는 따로 강제해야 한다.

flux를 섞는 방식은 이 강제를 보존성을 깨지 않고 해낸다. $\theta$가 무엇이든 여전히 flux 차분이기 때문이다.

그리고 실제로 개입하는 면은 0.1% 미만이다. 안전장치의 비용이 이 정도면 안 달 이유가 없다.

참고 문헌 X.Y. Hu, N.A. Adams, C.-W. Shu, "Positivity-preserving method for high-order conservative schemes solving compressible Euler equations", Journal of Computational Physics 242 (2013) 169–180.