음속 CFL 10에서 전처리기가 오차를 sweep마다 100배로 키웠다

저마하에서 전처리기가 무너지는 것은 반복이 느려서가 아니라, 버려진 블록 하나의 되먹임 이득이 1을 넘었기 때문이다.

시간간격을 10배씩 올리자 전처리기가 하나씩 죽었다#

Weston 등이 2019년 JCP에 낸 전속도(all-speed) 용융풀 해석기 논문에는 표가 한 장 있다. 뚜껑 구동 공동(lid-driven cavity) 문제에서 시간간격만 네 번 10배씩 올린 기록이다. 그러면 음속 CFL은 10.3에서 10,300까지 간다.

같은 격자, 같은 비선형 반복, 같은 Krylov 솔버(FGMRES)다. 바뀐 것은 전처리기뿐이다. 그런데 결과는 "몇 배 느려졌다"가 아니라 "수렴하느냐 마느냐"로 갈렸다. 전체 결합 시스템에 건 대수적 멀티그리드(AMG)는 음속 CFL이 10을 넘자 아예 수렴하지 못했다. 원시변수 블록 Gauss-Seidel은 100까지 버티다 그 위에서 무너졌다. 요소블록 SOR은 전부 수렴했지만 한 시간스텝에 FGMRES를 수백 번 돌렸다. 반복 횟수가 시간간격과 무관하게 일정했던 것은 Schur 보수 전처리기와 LU 분해뿐이었다.

전처리기는 보통 상수배의 문제다. 여기서는 임계값의 문제였다. 그 임계값이 어디서 오는지, 그리고 왜 하필 음속 CFL인지가 이 글의 내용이다. 전처리기 자체의 구성은 같은 논문을 다룬 이전 글에 적어 두었고, 여기서는 버려지는 블록 하나만 본다.

압력과 속도 사이에는 화살표가 두 개다#

논문은 자코비안을 보존변수가 아니라 원시변수 $W = (P, \mathbf{v}, T)$ 로 조립한다. 같은 물리를 쓰더라도 미지수를 무엇으로 잡느냐에 따라 행렬의 조건수가 달라지기 때문이다. 그러면 자코비안은 미지수 종류별로 묶인 3×3 블록 행렬이 된다.

M=\begin{bmatrix} M_{PP} & M_{Pv} & M_{PT}\\ M_{vP} & M_{vv} & M_{vT}\\ M_{TP} & M_{Tv} & M_{TT}\end{bmatrix}

여기서 $M_{Pv}$ 는 속도가 압력 방정식에 주는 기여, $M_{vP}$ 는 압력이 운동량 방정식에 주는 기여다. 논문은 압력-온도 결합이 약하다는 이유로 $M_{PT}$ 와 $M_{TP}$ 를 버린다. 그러고 나면 남는 뼈대는 압력-속도 2×2다.

후진 오일러로 시간을 이산화한 저마하 압축성 방정식에서 이 두 블록은 대각을 1로 맞춘 뒤 이렇게 생겼다.

A \equiv M_{Pv} \sim \rho c^{2} \Delta t\, \partial_x, \qquad B \equiv M_{vP} \sim \frac{\Delta t}{\rho}\, \partial_x

$\rho$ 는 밀도, $c$ 는 음속, $\Delta t$ 는 시간간격이다. $A$ 는 속도 발산이 압력을 밀어 올리는 항이고, $B$ 는 압력 기울기가 속도를 미는 항이다. 둘은 화살표 두 개짜리 되먹임 고리를 이룬다. 그 고리를 한 바퀴 돌면

AB \sim c^{2}\Delta t^{2}\, \partial_{xx}

즉 음향 연산자다. 격자 간격 $\Delta x$ 로 이산화하면 크기는 $(c\Delta t/\Delta x)^2$ , 바로 음속 CFL의 제곱이다. 이 숫자 하나가 아래 모든 이야기를 지배한다.

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

CFL_a1.00

Drag CFL_a with block Gauss-Seidel selected: the dropped arrow is the whole story — past CFL_a = 1 the loop gain 1.00 exceeds one and the residual curve turns upward. SOR keeps the arrow but only pointwise, so it survives to CFL_a = 2. Schur complement closes the loop and the residual drops to roundoff in a single sweep, at every CFL_a on the slider.

CFL_a 슬라이더를 끌면서 세 전처리기를 번갈아 눌러 보면 된다. 블록 Gauss-Seidel에서 붉게 끊긴 화살표가 이 글의 주인공이고, 그 화살표를 끊은 대가가 오른쪽 잔차 곡선의 기울기 부호로 나타난다.

하삼각만 쓰면 되돌아오는 화살표가 사라진다#

블록 Gauss-Seidel은 위 행렬의 하삼각 부분만 쓴다. 순서대로 풀면 이렇다.

x_P = M_{PP}^{-1} b_P, \qquad x_v = M_{vv}^{-1}(b_v - M_{vP} x_P)

압력을 먼저 푼다. 이때 속도는 전혀 보지 않는다. $A = M_{Pv}$ 가 통째로 빠졌다는 뜻이다. 되먹임 고리의 화살표 하나가 잘린 것이다.

잘린 대가는 정확히 계산된다. $M = L + U$ 로 쪼개면 $U$ 에는 $A$ 만 남는다. 오차 $e$ 에 대한 반복 연산자는

G = -L^{-1}U = \begin{bmatrix} 0 & -A \\ 0 & BA \end{bmatrix}

이 행렬의 고유값은 0과 $BA$ 의 고유값뿐이다. 주기 경계 중심차분 $D$ 의 고유값이 $i\sin\theta/\Delta x$ 이므로 $BA$ 의 고유값은 $-\mathrm{CFL}_a^2 \sin^2\theta$ 가 되고,

\rho(G) = \max_{\theta}\ \mathrm{CFL}_a^{2}\sin^{2}\theta = \mathrm{CFL}_a^{2}

$\rho$ 는 스펙트럼 반경, 즉 sweep 한 번당 오차가 곱해지는 최대 배율이다. 임계값이 $\mathrm{CFL}_a = 1$ 에 정확히 놓인다. 그 위에서는 전처리기가 오차를 줄이는 게 아니라 키운다.

요소블록 SOR은 사정이 조금 낫다. 중심차분은 대각 성분이 0이라 대각 블록이 단위행렬이 되고, 완화계수 $\omega$ 를 쓰면 이득은 $\sqrt{(1-\omega)^2 + \omega^2\,\mathrm{CFL}_a^2}$ 가 된다. $\omega = 0.4$ 면 임계값이 $\mathrm{CFL}_a = 2$ 로 밀린다. CFL에 제곱이 아니라 1차로 붙기 때문에 더 오래 버틴다. 논문에서 SOR이 Gauss-Seidel보다 robust했던 순서가 여기서 나온다.

Python으로 세 전처리기의 이득을 재봤다#

주기 경계 1D 선형 음향계를 후진 오일러로 세우고, 세 전처리기의 오차 연산자를 그대로 반복 적용했다. 외부 라이브러리는 쓰지 않았다.

import math
 
N, L, RHO = 32, 1.0, 1.0
dx = L / N
 
 
def deriv(v):
    """중심차분 1계 미분, 주기 경계"""
    return [(v[(i + 1) % N] - v[(i - 1) % N]) / (2 * dx) for i in range(N)]
 
 
def build_ops(c, dt):
    """M = [[I, A], [B, I]] — A는 속도→압력, B는 압력→속도 블록"""
    A = lambda u: [RHO * c * c * dt * w for w in deriv(u)]
    B = lambda p: [dt / RHO * w for w in deriv(p)]
    return A, B
 
 
def gain(step, warm=200, n=400):
    """오차 연산자를 반복 적용하고 한 sweep당 증폭률의 기하평균을 낸다"""
    e = [math.sin(1.7 * i * i + 0.9 * i + 1.0) for i in range(2 * N)]
    acc = 0.0
    for k in range(warm + n):
        f = step(e)
        r = math.sqrt(sum(x * x for x in f)) / math.sqrt(sum(x * x for x in e))
        if r == 0.0:
            return 0.0
        if k >= warm:
            acc += math.log(r)
        e = [x / r for x in f]
    return math.exp(acc / n)
 
 
def gs_step(A, B):
    """블록 Gauss-Seidel: 하삼각만 쓰므로 A 블록이 통째로 빠진다"""
    def step(e):
        Aeu = A(e[N:])
        return [-x for x in Aeu] + B(Aeu)
    return step
 
 
def sor_step(A, B, w):
    """점블록 SOR: 중심차분은 대각이 0이라 대각 블록이 I가 된다"""
    def step(e):
        ep, eu = e[:N], e[N:]
        jp, ju = A(eu), B(ep)
        return ([(1 - w) * ep[i] - w * jp[i] for i in range(N)]
                + [(1 - w) * eu[i] - w * ju[i] for i in range(N)])
    return step
 
 
def schur_cg(A, B, rhs, tol=1e-10, cap=200):
    """S = I - A B 는 대칭 양정치 — 켤레기울기법 반복 횟수를 돌려준다"""
    S = lambda p: [p[i] - v for i, v in enumerate(A(B(p)))]
    x, r = [0.0] * N, rhs[:]
    d, rr = rhs[:], sum(v * v for v in rhs)
    r0 = math.sqrt(rr)
    for k in range(1, cap + 1):
        Sd = S(d)
        al = rr / sum(d[i] * Sd[i] for i in range(N))
        x = [x[i] + al * d[i] for i in range(N)]
        r = [r[i] - al * Sd[i] for i in range(N)]
        rn = sum(v * v for v in r)
        if math.sqrt(rn) < tol * r0:
            return k
        d = [r[i] + (rn / rr) * d[i] for i in range(N)]
        rr = rn
    return cap
 
 
def sweeps_to(r, drop=1e-6):
    """오차를 백만분의 1로 줄이는 데 필요한 sweep 수"""
    return "diverge" if r >= 0.999 else str(int(math.ceil(math.log(drop) / math.log(r))))
 
 
rhs = [1.0 if N // 3 <= i < 2 * N // 3 else 0.0 for i in range(N)]  # 여러 모드가 섞인 우변
 
print("CFL_a   rho(GS)    rho(SOR)   sweep(GS)  sweep(SOR) CG on S")
for cfl in [0.1, 0.5, 1.0, 2.0, 10.0, 100.0]:
    A, B = build_ops(1.0, cfl * dx)
    rg, rs = gain(gs_step(A, B)), gain(sor_step(A, B, 0.4))
    print("%-7g %-10.4g %-10.4g %-10s %-10s %d"
          % (cfl, rg, rs, sweeps_to(rg), sweeps_to(rs), schur_cg(A, B, rhs)))
 
print()
print("Mach sweep  (material CFL fixed at 0.5)")
print("Mach     CFL_a    rho(GS)    rho(SOR)   CG on S")
for mach in [1e-2, 1e-3, 1e-4, 1e-5, 1e-6]:
    dt = 0.5 * dx / 1.0          # 물질 속도 |u| = 1 이 시간간격을 정한다
    c = 1.0 / mach               # 음속은 마하수가 정한다
    A, B = build_ops(c, dt)
    print("%-8.0e %-8.4g %-10.4g %-10.4g %d"
          % (mach, c * dt / dx, gain(gs_step(A, B)), gain(sor_step(A, B, 0.4)),
             schur_cg(A, B, rhs)))

CFL_a   rho(GS)    rho(SOR)   sweep(GS)  sweep(SOR) CG on S
0.1     0.01       0.601      3          28         4
0.5     0.25       0.6321     10         31         8
1       1          0.721      diverge    43         9
2       4          1          diverge    diverge    9
10      100        4.045      diverge    diverge    9
100     1e+04      40         diverge    diverge    9
 
Mach sweep  (material CFL fixed at 0.5)
Mach     CFL_a    rho(GS)    rho(SOR)   CG on S
1e-02    50       2500       20.06      9
1e-03    500      2.5e+05    200.3      9
1e-04    5000     2.5e+07    2005       12
1e-05    5e+04    2.5e+09    2.006e+04  17
1e-06    5e+05    2.5e+11    2.006e+05  17

측정값이 손으로 얻은 식과 자릿수까지 맞는다. Gauss-Seidel의 이득은 0.01, 0.25, 1, 4, 100, 10,000 — 정확히 $\mathrm{CFL}_a^2$ 다. SOR은 0.601, 0.632, 0.721, 1.0, 4.045, 40 — $\sqrt{0.36 + 0.16\,\mathrm{CFL}_a^2}$ 그대로다. 음속 CFL 10에서 Gauss-Seidel은 sweep 한 번마다 오차를 100배로 키운다.

토이 모델이라 임계값이 1과 2에 딱 떨어진다. 실제 코드는 sweep을 10번씩 돌리고 $\omega = 0.4$ 로 눌러 두기 때문에 실패 지점이 CFL 100 근처까지 밀린다. 밀리는 위치가 다를 뿐, 밀어 올리는 것이 무엇인지는 같다.

마하수를 내리는 것은 시간간격을 올리는 것과 같다#

논문의 두 번째 실험은 시간간격을 고정하고 음속만 10배씩 올린다. 저마하 해석에서는 시간간격을 물질 시간척도로 잡는 것이 상식이다.

\Delta t = \frac{\mathrm{CFL}_u\, \Delta x}{|u|}, \qquad c = \frac{|u|}{M} \ \Longrightarrow\ \mathrm{CFL}_a = \frac{\mathrm{CFL}_u}{M}

$M$ 은 마하수다. 물질 CFL을 0.5로 얌전히 두어도 $M = 10^{-6}$ 이면 선형 솔버는 음속 CFL 500,000을 받는다. 저마하 해석의 어려움은 물리가 아니라 이 숫자에 있다. 음파를 명시적으로 풀지 않으려고 음향·대류를 분리하는 접근이 나온 것도 같은 이유다.

Mach1e-2CFL_u0.5

Pull Mach down one notch at a time and watch the orange front: at 1e-2 it laps the domain about once per step, at 1e-6 it laps it 1.3e+0 times while the blue particle still crawls 0.5 cells. Nothing about the physics changed — only the number the linear solver is handed. The bars below are the same sweep as the article's Mach table: Gauss-Seidel crosses the red line first, SOR next, and the Schur bar never leaves zero.

Mach 슬라이더를 한 칸씩 내리면서 주황색 음파 전선이 한 시간스텝 동안 도메인을 몇 바퀴 도는지 세어 보자. 파란 물질 입자는 계속 0.5칸만 걷는데, 아래 세 막대는 왼쪽부터 차례로 빨간 선을 넘는다.

위 출력의 Mach 스윕이 논문 Fig. 4와 같은 순서를 준다. 블록 Gauss-Seidel은 $M = 10^{-3}$ 아래에서, 요소블록 SOR은 $10^{-4}$ 아래에서 수렴하지 못했다. $M = 10^{-6}$ 까지 간 것은 Schur 보수와 LU뿐이었다.

Schur 보수는 그 화살표를 근사하지 않고 소거한다#

버려진 $A$ 를 되살리는 방법은 그것을 근사하는 것이 아니라 소거하는 것이다. 압력에 대한 Schur 보수(Schur complement, 한 블록을 소거하고 남는 유효 연산자)는 이렇다.

S = M_{PP} - M_{Pv}M_{vv}^{-1}M_{vP} = I - AB

$S$ 로 압력을 풀고 그 값을 운동량에 넣으면 블록 LU 분해가 정확해진다. 반복이 필요 없다. 오차 연산자가 0이고, 위 코드에서도 어떤 CFL에서든 sweep 한 번에 반올림 오차까지 떨어진다.

값을 안 치르는 것은 아니다. 치르는 자리가 바뀔 뿐이다.

S = I - c^{2}\Delta t^{2}\,\partial_{xx}

이것은 헬름홀츠 형태이고, $-\partial_{xx}$ 가 양정치이므로 $S$ 는 대칭 양정치다. 음속 CFL이 커질수록 $I$ 가 묻히고 압력 포아송 방정식에 가까워진다. 위 표의 마지막 열이 그 대가다. CG 반복이 4에서 9로, 저마하 스윕에서는 17까지 늘어난다. 논문이 "Schur 전처리기는 시간간격에 따라 CPU 시간이 약간 증가한다"고 적은 것이 이 증가다.

대신 남은 문제가 대칭 양정치라는 점이 중요하다. 비대칭 결합 시스템에서 힘을 못 쓰던 AMG가 여기서는 제 영역을 만난다. GMRES가 부분공간을 쌓는 방식을 떠올리면, 바깥 FGMRES가 짧게 끝나도록 안쪽에 잘 맞는 solver를 꽂아 주는 구조다.

논문이 실제로 치른 값 — 근사 세 단계와 자코비안 지연#

실전 구현은 위 유도보다 층이 많다. 논문은 전처리를 세 단계로 나눈다. 첫째, 근사 자코비안 위에 어떤 전처리기를 얹을지(AMG, 요소블록 SOR, 블록 Gauss-Seidel, vP-vT Schur 보수, LU). 둘째, Schur 보수 자체를 어떻게 근사할지(세 가지 전략). 셋째, 각 블록을 어떤 스무더로 풀지(다섯 가지 선택지). "AMG (#1)", "AMG-FGMRES (#3)" 같은 표기가 이 조합을 가리킨다.

자코비안 조립도 값이 든다. 유한차분으로 만들고, 섭동 크기는 $\varepsilon_{rel} = 10^{-8}$ , $x_{min} = 10^{-7}$ 이다. PETSc의 그래프 컬러링을 써 봤지만 잔차 평가 횟수가 과하게 늘었고, 고차 스킴과 3D에서 더 나빠졌다. 결국 요소 단위로 국소 섭동을 주어 요소 자코비안을 조립하는 쪽을 택했다. 잔차 평가가 훨씬 적고 근사도 더 정확했다.

그리고 자코비안을 매 Newton 반복마다 다시 만들지 않는다. 얼려 두었다가, 한 Newton 반복 안에서 바깥 FGMRES가 20~50회를 넘기면 그때 재조립을 켠다. 근사 자코비안은 전처리용일 뿐이고 JFNK가 쓰는 실제 자코비안 벡터곱은 항상 최신이라 가능한 절충이다.

한 가지 단서도 붙어 있다. Schur 보수 전처리기는 중간~높은 마하수에서도 동작하지만, 비용 대비 효과가 있는 것은 저마하 영역이라고 논문은 명시한다.

전처리기가 물리를 알아야 한다는 말의 실제 의미#

"physics-based preconditioner"라는 말은 보통 모호하게 쓰인다. 이 논문에서는 뜻이 좁고 분명하다. 어떤 블록 결합이 시간간격이나 마하수와 함께 커지는지를 알고, 그 결합만은 근사하지 말라는 것이다.

그래서 암시적 솔버가 시간간격을 올릴 때 갑자기 수렴하지 않는다면, 반복 횟수나 tolerance를 만지기 전에 물어볼 것이 하나 있다. 지금 전처리기가 버리고 있는 블록이 무엇이고, 그 블록의 이득이 무엇에 비례하는가. 음향 결합이면 $\mathrm{CFL}_a^2$ 이고, 늘어난 경계층 격자라면 종횡비가 그 자리에 들어간다. 답이 1을 넘는 순간, 그 전처리기는 더 이상 느린 것이 아니라 틀린 것이다.

음속 CFL 10에서 전처리기가 오차를 sweep마다 100배로 키웠다 — 블록 전처리기가 버린 화살표 하나