1001배 큰 시간스텝의 청구서 — 반음해 압력 스킴이 L∞ 안정이 아닌 이유

음속을 CFL 조건에서 지우면 시간스텝은 풀리지만, 불연속면의 진동은 감쇠되지 않고 그대로 남는다.

1971년, 압력을 음속에서 떼어낸 계산법#

Francis Harlow와 Anthony Amsden이 Los Alamos에서 ICE 방법을 발표한 해다. 논문 제목에 "all flow speeds"가 들어갔다. 압축성 코드가 저마하 영역에서 기어가는 문제를, 공간이 아니라 시간스텝 쪽에서 풀겠다는 선언이었다. 그 계보의 최신판이 Boscheri와 Pareschi의 2021년 반음해(semi-implicit) 압력 기반 스킴이다. 3차원 Navier–Stokes를 전 마하수 영역에서 3차 정확도로 푼다. 버는 것과 잃는 것을 표의 숫자로 맞춰 본다. 시간스텝은 최대 1001배를 벌었고, 불연속면의 진동은 네 바퀴를 돌아도 그대로 남았다.

먼저 무엇을 버는지부터 눈으로 보자. 아래 그림에서 마하수 슬라이더를 직접 내려보자.

Mach 0.100CFL 0.50

dt_explicit = 1.136e-2   dt_semi = 1.250e-1   gain = 11.0x   steps so far: 0 vs 0

Both panels cover the same physical time and the same 40 cells. Pull the Mach slider down and the red acoustic characteristics tilt until they are nearly horizontal — the explicit panel has to slice that window into ever thinner layers to keep them inside one cell, while the green panel keeps taking the same eight steps. The step counter is the whole argument for treating pressure implicitly.

위아래 두 패널은 같은 물리 시간, 같은 40개 셀을 덮는다. 마하수를 낮출수록 빨간 음향 특성선이 눕고, 위 패널은 그 선을 셀 하나 안에 가두려고 시간층을 잘게 썰기 시작한다. 아래 패널의 스텝 수는 변하지 않는다.

시간스텝이 마하수에 묶이는 자리#

무차원화된 Navier–Stokes에서 강성(stiffness)을 재는 변수는 $\varepsilon$ , 곧 마하수의 제곱이다. 음속은 $c/\sqrt{\varepsilon}$ 로 나타난다. 완전 양해(explicit) 스킴의 CFL 조건은 이 음속을 그대로 물고 있다.

\Delta t \le \mathrm{CFL}\,\frac{\Delta x}{\max\left(|u| + c/\sqrt{\varepsilon}\right)}

$\Delta x$ 는 격자 간격, $u$ 는 유체 속도, $c$ 는 기준 음속이다. $\varepsilon \to 0$ 이면 분모가 발산하고 $\Delta t$ 는 $\sqrt{\varepsilon}$ 처럼 0으로 간다. 물리적으로 아무 일도 일어나지 않는 저마하 유동에서, 아무도 관심 없는 음파 때문에 계산이 멈춰 서는 것이다.

반음해 스킴은 압력 항만 음해(implicit)로 넘긴다. 그러면 안정성 조건에서 음속이 사라진다.

\Delta t \le \mathrm{CFL}\,\frac{\Delta x}{\max |u|}

이 두 식의 비가 $1 + 1/M$ 이다. $M = 0.001$ 이면 $1 + 1000 = 1001$ 배다. 표는 그 비를 gain 열에 1001.0으로 찍고, 양해 시간스텝 $1.249\times 10^{-3}$ 과 반음해 $1.250$ 을 나란히 둔다. $M=0.1$ 이면 이득은 11.0, 양해 $1.136\times 10^{-3}$ 대 반음해 $1.250\times 10^{-2}$ 다.

코드의 시간간격은 $\mathrm{CFL}\,\Delta x/|U_0|$ 다. $\Delta x = 1/400 = 0.0025$ , CFL $= 0.5$ , $U_0 = 0.1$ 이면 $dt = 0.0125$ 다. 한 바퀴 시간 $1/(U_0 + c) = 1/1.1$ 을 0.0125로 나누면 72.7이고, 출력은 73스텝이다. 음향 CFL은 $0.5 \times 1 / 0.1 = 5$ 다. 음속을 CFL에서 빼는 분할을 시간적분 쪽에서 다시 자른 글이 음향과 수송을 나누는 Peluchon IMEX에 있다.

압력이 타원형 방정식으로 바뀌는 과정#

논문의 1차 반음해 반이산화는 세 줄이다. 연속 방정식은 완전 양해, 운동량 방정식은 압력 구배만 음해, 에너지 방정식은 엔탈피 플럭스의 운동량을 음해로 둔다.

\frac{\rho^{n+1} - \rho^{n}}{\Delta t} + \nabla\cdot(\rho\mathbf{u})^{n} = 0

\frac{(\rho\mathbf{u})^{n+1} - (\rho\mathbf{u})^{n}}{\Delta t} + \nabla\cdot\left((\rho\mathbf{u})^{n}\otimes\mathbf{u}^{n}\right) + \frac{1}{\varepsilon}\nabla p^{n+1} = 0

$n$ 은 현재 시간층, $n+1$ 은 다음 시간층이다. 여기서 운동량 식을 에너지 식에 대입하면 미지수가 $p^{n+1}$ 하나뿐인 방정식이 남는다.

\varepsilon(\rho e)^{n+1} + \varepsilon\frac{\Delta t}{2}\frac{(\rho\mathbf{u})^{n}}{\rho^{n}}\nabla p^{n+1} - \Delta t^{2}\,\nabla\cdot\left(h^{n}\nabla p^{n+1}\right) = b^{n}

$h = (\rho e + p)/\rho$ 는 비엔탈피, $b^n$ 은 양해 항만으로 만들어진 우변이다. $\Delta t^2 \nabla\cdot(h\nabla p)$ 항이 이 식을 타원형으로 만든다. 이상기체라면 $(\rho e)^{n+1} = p^{n+1}/(\gamma-1)$ 이므로 선형 시스템이고, Redlich–Kwong 같은 비선형 상태방정식이면 Newton 반복이 붙는다. 논문은 허용오차 $10^{-10}$ 에서 최대 4회로 수렴했다고 보고한다.

압력이 풀리면 운동량을 되돌려 계산하고, 총에너지는 운동에너지의 음해 부분을 붙여 마무리한다. 요약하면 질량 → 압력(타원형) → 운동량 → 에너지 순서다. 비압축성 투영법과 구조가 같지만, 여기서는 상태방정식이 살아 있어서 압축성이 그대로 남는다.

저마하 극한에서 발산이 정말 0으로 가는가#

점근보존(asymptotic preserving)은 $\varepsilon \to 0$ 에서 이산 스킴이 저절로 비압축 극한을 재현한다는 성질이다. 모든 변수를 $\varepsilon$ 의 멱급수로 전개해 선행 항만 남기면 압력의 선행 항이 공간에 대해 상수가 된다.

\nabla p_{(0)}^{n+1} = 0

이것을 에너지 방정식의 극한에 넣으면 이상기체에서 다음이 나온다.

h^{n}\,\nabla\cdot\mathbf{u}^{n+1} = O(\Delta t)

$O(\Delta t)$ 의 상수는 $\varepsilon$ 에 의존하지 않는다. 마하수를 아무리 낮춰도 속도장의 발산이 $\Delta t$ 수준으로 유지된다는 뜻이다. 격자를 바꾸지 않고 저마하로 내려가도 압력이 체커보드로 갈라지지 않는 이유가 여기 있다.

Python으로 같은 펄스를 네 바퀴 돌렸다#

논문 3.6절은 불편한 문장을 하나 적어둔다. 고차 음해 시간이산화는 강한 불연속 앞에서 오버슈트를 제거하지 못하며, 그 진동은 커지지도 사라지지도 않는다는 것이다. 음해니까 안정하다는 말은 $L_2$ 안정이라는 뜻이지 $L_\infty$ 안정이라는 뜻이 아니다.

이 주장만 따로 떼어 확인하려면 비선형 솔버 전체가 필요하지 않다. 선형화된 1차원 Euler 방정식에 같은 분해를 걸면 된다. 물질 이송은 양해, 음향 결합은 음해로 두고, 음해 가중치 $\theta$ 를 0.5(2차, 무감쇠)에서 1.0(1차, 감쇠)까지 움직인다. 엇갈림 격자에 압력을 셀 중심, 속도를 면에 두면 압력 방정식이 정확히 3점 Helmholtz 시스템이 된다. 주기 경계라서 순환 삼중대각 solver를 쓴다.

초기 조건은 오른쪽으로 달리는 사각 음향 펄스다. 정확해는 그 펄스를 $u_0 + c$ 만큼 평행이동한 것이므로 오차를 해석적으로 잴 수 있다. 물질 CFL을 0.5로 잡으면 음향 CFL은 5가 된다.

import numpy as np
 
GAMMA, RHO0, C0 = 1.4, 1.0, 1.0
P0 = C0 ** 2 * RHO0 / GAMMA          # so that c = sqrt(gamma p0 / rho0) = 1
N, CFL, U0 = 400, 0.5, 0.1
DX = 1.0 / N
 
 
def cyclic_thomas(a, b, c, d):
    """주기 삼중대각 풀이: Thomas + Sherman-Morrison 보정."""
    n = len(b)
    al, be, g = a[0], c[-1], -b[0]
    bb = b.copy()
    bb[0] -= g
    bb[-1] -= al * be / g
 
    def tri(rhs):
        cp, dp = np.zeros(n), np.zeros(n)
        cp[0], dp[0] = c[0] / bb[0], rhs[0] / bb[0]
        for i in range(1, n):
            m = bb[i] - a[i] * cp[i - 1]
            cp[i], dp[i] = c[i] / m, (rhs[i] - a[i] * dp[i - 1]) / m
        x = np.zeros(n)
        x[-1] = dp[-1]
        for i in range(n - 2, -1, -1):
            x[i] = dp[i] - cp[i] * x[i + 1]
        return x
 
    v = np.zeros(n)
    v[0], v[-1] = g, be
    y, z = tri(d), tri(v)
    fac = (y[0] + al * y[-1] / g) / (1.0 + z[0] + al * z[-1] / g)
    return y - fac * z
 
 
def grad_face(p):
    return (np.roll(p, -1) - p) / DX
 
 
def div_cell(u):
    """div_cell(grad_face(.))가 정확히 3점 라플라시안이 된다."""
    return (u - np.roll(u, 1)) / DX
 
 
def upwind(f):
    return (f - np.roll(f, 1)) / DX
 
 
def si_update(u, p, dt, theta):
    """물질 이송은 양해, 음향 결합은 음해 가중치 theta로."""
    ps = p - dt * U0 * upwind(p)
    ue = u - dt * U0 * upwind(u) - dt * (1 - theta) / RHO0 * grad_face(p)
    k = C0 ** 2 * (theta * dt) ** 2 / DX ** 2
    rhs = ps - dt * GAMMA * P0 * (theta * div_cell(ue) + (1 - theta) * div_cell(u))
    p_n = cyclic_thomas(np.full(N, -k), np.full(N, 1 + 2 * k), np.full(N, -k), rhs)
    return ue - dt * theta / RHO0 * grad_face(p_n), p_n
 
 
def chi_flattener(u, k1):
    """식 (104): 압축률이 k1 * c_min을 넘는 셀에 표시."""
    divu = div_cell(u)
    return np.minimum(1.0, np.maximum(0.0, -(divu + k1 * C0) / (k1 * C0)))
 
 
def exact_wave(x, t):
    s = np.mod(x - (U0 + C0) * t, 1.0)
    return np.where((s >= 0.3) & (s < 0.6), 1.0, 0.0)
 
 
def solve_case(n_transit, theta, k1=None):
    xc, xf = (np.arange(N) + 0.5) * DX, (np.arange(N) + 1.0) * DX
    p, u = exact_wave(xc, 0.0), exact_wave(xf, 0.0) / (RHO0 * C0)
    t_end, dt = n_transit / (U0 + C0), CFL * DX / abs(U0)
    t, nstep, nflag = 0.0, 0, 0
    while t < t_end - 1e-14:
        h = min(dt, t_end - t)
        u2, p2 = si_update(u, p, h, theta)
        if k1 is not None:
            chi = chi_flattener(u, k1)
            nflag += int(np.sum(chi > 0))
            u1, p1 = si_update(u, p, h, 1.0)
            cf = np.maximum(chi, np.roll(chi, -1))
            u2, p2 = cf * u1 + (1 - cf) * u2, chi * p1 + (1 - chi) * p2
        u, p, t, nstep = u2, p2, t + h, nstep + 1
    ex = exact_wave(xc, t_end)
    return dict(l1=float(np.mean(np.abs(p - ex))), over=float(np.max(p) - 1.0),
                under=float(-np.min(p)), l2=float(np.sqrt(np.mean(p ** 2))),
                steps=nstep, flagged=nflag, acfl=CFL * C0 / abs(U0))
 
 
print("=== (1) time step from Eq.(20) vs Eq.(106), CFL=0.5, 400 cells ===")
print(f"{'Mach':>8}{'dt explicit':>14}{'dt semi-impl':>15}{'gain':>8}")
for M in [0.5, 0.2, 0.1, 0.02, 0.005, 0.001]:
    de, di = CFL * DX / (M * C0 + C0), CFL * DX / (M * C0)
    print(f"{M:8.3f}{de:14.3e}{di:15.3e}{di / de:8.1f}")
 
b = solve_case(1, 1.0)
print(f"\n=== (2) square acoustic pulse, M={U0}, material CFL={CFL}")
print(f"    -> acoustic CFL {b['acfl']:.1f}, {b['steps']} steps per transit ===")
print(f"{'theta':>7}{'over@1':>10}{'over@2':>10}{'over@4':>10}{'L1@4':>9}{'|p|2@4':>9}")
for th in (0.50, 0.52, 0.55, 0.60, 0.70, 1.00):
    r = [solve_case(n, th) for n in (1, 2, 4)]
    ov = [f"{x['over']:10.4f}" if abs(x['over']) < 1e3 else f"{'blow-up':>10}" for x in r]
    l1 = f"{r[2]['l1']:9.4f}" if r[2]['l1'] < 1e3 else f"{'-':>9}"
    l2 = f"{r[2]['l2']:9.4f}" if r[2]['l2'] < 1e3 else f"{'-':>9}"
    print(f"{th:7.2f}{ov[0]}{ov[1]}{ov[2]}{l1}{l2}")
 
print("\n=== (3) flattener blend, theta=0.55 against theta=1, 4 transits ===")
print(f"{'k1':>9}{'flagged/step':>14}{'over':>9}{'under':>9}{'L1':>9}")
base = solve_case(4, 0.55)
print(f"{'off':>9}{0:14.1f}{base['over']:9.4f}{base['under']:9.4f}{base['l1']:9.4f}")
for k1 in [1e-3, 1e-1, 1e0, 1e1, 1e2]:
    r = solve_case(4, 0.55, k1)
    print(f"{k1:9.0e}{r['flagged'] / r['steps']:14.1f}{r['over']:9.4f}"
          f"{r['under']:9.4f}{r['l1']:9.4f}")
mono = solve_case(4, 1.0)
print(f"{'theta=1':>9}{'(all cells)':>14}{mono['over']:9.4f}"
      f"{mono['under']:9.4f}{mono['l1']:9.4f}")

=== (1) time step from Eq.(20) vs Eq.(106), CFL=0.5, 400 cells ===
    Mach   dt explicit   dt semi-impl    gain
   0.500     8.333e-04      2.500e-03     3.0
   0.200     1.042e-03      6.250e-03     6.0
   0.100     1.136e-03      1.250e-02    11.0
   0.020     1.225e-03      6.250e-02    51.0
   0.005     1.244e-03      2.500e-01   201.0
   0.001     1.249e-03      1.250e+00  1001.0
 
=== (2) square acoustic pulse, M=0.1, material CFL=0.5
    -> acoustic CFL 5.0, 73 steps per transit ===
  theta    over@1    over@2    over@4     L1@4   |p|2@4
   0.50  488.1288   blow-up   blow-up        -        -
   0.52    4.5456  212.3962   blow-up        -        -
   0.55    0.2285    0.2023    0.1825   0.1348   0.5169
   0.60    0.0538    0.0286   -0.0346   0.1270   0.4630
   0.70   -0.0163   -0.0900   -0.2320   0.1992   0.4064
   1.00   -0.1610   -0.3210   -0.4823   0.3098   0.3347
 
=== (3) flattener blend, theta=0.55 against theta=1, 4 transits ===
       k1  flagged/step     over    under       L1
      off           0.0   0.1825   0.1663   0.1348
    1e-03         302.9  -0.3611   0.4734   0.2870
    1e-01         280.0  -0.3544   0.4998   0.2888
    1e+00         196.0  -0.1857   0.5381   0.2448
    1e+01           8.2   0.1275   0.3135   0.1469
    1e+02           0.0   0.1801   0.1726   0.1355
  theta=1   (all cells)  -0.4823  -0.0985   0.3098

$\theta$ 슬라이더를 0.55에 두면 네 바퀴 오버슈트가 0.1825에서 멈추고, 0.50으로 내리면 한 바퀴에 488에서 터진다.

theta 0.55k1 = 1e1

transits 0.00   overshoot 0.0000   undershoot 0.0000   L1 0.0000   flagged 0 / 128

The acoustic CFL here is 5, so the implicit solve is doing all the work. Leave theta at 0.55 and watch the overshoot readout while the pulse laps the domain: it settles near 0.18 and stays there — bounded, but never shrinking. Push theta to 1.00 and the wiggles vanish while the pulse melts. Turning on the flattener with k1 = 1e-3 lights up the yellow ticks along almost the whole axis, which is the scheme quietly dropping to first order everywhere.

1.00까지 밀면 오버슈트는 음수가 되고 $L_{1}$ 은 0.3098까지 커진다. 0.55의 $L_{1}$ 0.1348과 그 끝이 한 슬라이더 위에 있다.

0.55와 0.50 사이에서 갈라지는 것#

두 번째 표가 갈림이다. $\theta = 0.5$ 는 Crank–Nicolson, 곧 형식상 2차이고 음향 부분만 보면 에너지를 정확히 보존한다. 그런데 양해 이송과 합치자 한 바퀴 만에 오버슈트가 488까지 올라가고 두 바퀴째에 터진다. 음해 반쪽이 감쇠를 전혀 제공하지 않으면, 양해 반쪽의 풍상차분 감쇠만으로는 분해 오차를 못 덮는다. Casulli 계열 스킴이 $\theta$ 를 0.5보다 조금 크게 잡으라고 권하는 이유가 이것이다.

$\theta = 0.55$ 에서 상황이 바뀐다. 오버슈트가 0.2285 → 0.2023 → 0.1825로 간다. 0.2285에서 0.2023까지는 0.0262, 그다음 0.1825까지는 0.0198이다. 네 바퀴의 0.1825는 높이 1인 펄스의 18%다. 커지지 않으므로 $L_{2}$ 로는 안정이고, 줄지도 않으므로 $L_{\infty}$ 로는 안정이 아니다. 논문이 적어둔 문장이 그대로 재현된다. 반대편 끝의 $\theta = 1.0$ 은 오버슈트가 음수다. 펄스가 제 높이에 닿지도 못했다는 뜻이고, $L_{1}$ 은 0.1348 대 0.3098이다. $0.3098/0.1348 \approx 2.30$ 이라 2.3배와 맞다. $|p|_{2}$ 는 오히려 0.5169에서 0.3347로 준다. $L_{2}$ 로 보면 1.0이 더 조용하고, 봉우리만 보면 0.55가 덜 녹는다.

flattener의 $k_1$ 은 임계값이 아니었다#

논문의 해법은 감지 후 혼합이다. 속도장의 발산을 이웃 셀 최소 음속과 비교해 충격파 셀을 찾는다.

\chi^{n} = \min\left[1,\ \max\left(0,\ -\frac{\nabla\cdot\mathbf{u}^{n} + k_1 c^{n}_{\min}}{k_1 c^{n}_{\min}}\right)\right]

$\chi$ 가 0이면 고차해를 그대로 쓰고, 1이면 1차해로 바꾼다. 그 사이는 볼록결합이다. 논문은 모든 계산에 $k_1 = 10^{-3}$ 을 썼고 "표시되는 셀은 매우 적다"고 적었다.

세 번째 표가 그 상수를 옮겨 쓸 때 일어나는 일이다. $k_1 = 10^{-3}$ 이면 400개 중 302.9개가 걸린다. 오버슈트는 0.1825에서 $-0.3611$ 로 사라지지만 $L_1$ 오차가 0.1348에서 0.2870으로 2.1배가 된다. 이건 감지가 아니라 스킴 전체를 1차로 내린 것이다. 실제로 $\theta=1$ 결과(0.3098)와 거의 같은 자리에 도착한다.

이유는 식 자체에 있다. $\chi$ 는 $\nabla\cdot\mathbf{u}$ 가 $-k_1 c$ 에서 $-2k_1 c$ 로 가는 사이에 0에서 1로 올라간다. 그 폭은 $k_1 c$ 하나다. $k_1 c = 10^{-3}$ 이면 0.001이다. $k_1 = 10^{-3}$ 에서 302.9/400 = 0.757, 셀의 약 76%다. $L_{1}$ 은 $0.2870/0.1348 \approx 2.13$ 으로 2.1배와 맞다. 이 문제에서 발산은 $O(1)$ 크기로 요동하므로, 조금이라도 압축 중인 셀은 전부 $\chi = 1$ 이 된다. 임계값이 아니라 부호 판정기로 동작하는 것이다.

셀을 실제로 골라내려면 $k_1$ 을 네 자릿수 올려야 했다. $k_1 = 10$ 에서 스텝당 8.2개만 걸리고, 오버슈트는 0.1825에서 0.1275로 간다. $0.1275/0.1825 \approx 0.70$ 이라 약 30%가 줄고, $0.1469/0.1348 \approx 1.09$ 라 $L_{1}$ 은 9%만 커진다. $k_1 = 10^2$ 이면 아무것도 걸리지 않는다. 논문의 값이 틀렸다는 뜻이 아니다. 그 값은 매끄러운 영역에서 발산이 사실상 0인 비선형 충격관 계산에 맞춰져 있고, 진동이 이미 깔려 있는 문제로 옮기면 기준이 의미를 잃는다.

음속을 지운 자리에 무엇이 남는가#

음속을 CFL에서 지운 만큼, 시간적분은 그 음파를 못 따라간다. 저마하에서는 그게 맞다. 마하수가 올라가면 그 음파가 충격파다. $\theta$ 는 정확도와 단조성을 맞바꾸는 손잡이다. 0.5는 이 분해에서 한 바퀴에 488까지 올라 안정하지도 않았다.

감지기 상수는 문제의 발산 스케일에 붙어 있다. 논문의 $10^{-3}$ 을 그대로 복사하면 400개 중 302.9개가 걸려 스위치가 된다. 압력 부분계로 음속을 옮기는 다른 분할은 Toro–Vázquez가 강성을 옮기는 자리에 있다. 거기도 음속은 사라지지 않고 다른 방정식으로 이사한다.