음속 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인지가 이 글의 내용이다. 전처리기 자체의 구성은 같은 논문을 다룬 이전 글에 적어 두었고, 여기서는 버려지는 블록 하나만 본다. CFL_a를 1 위로 올리면 블록 Gauss-Seidel의 잔차 기울기가 감소에서 증가로 뒤집힌다.

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.

붉게 끊긴 화살표가 하삼각이 버린 $A$ 다. 세 전처리기를 번갈아 누르면, 같은 슬라이더에서 잔차가 줄어드는 쪽과 불어나는 쪽이 갈린다.

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

논문은 자코비안을 보존변수가 아니라 원시변수 $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의 제곱이다. 이 숫자 하나가 아래 모든 이야기를 지배한다. 위에서 잔차 기울기가 뒤집힌 지점이 $\mathrm{CFL}_a = 1$ 이다.

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

블록 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했던 순서가 여기서 나온다. 음속 CFL이 시간간격을 잡아먹지 않게 음향과 수송을 가르는 법은 Peluchon의 음향-수송 IMEX에 있다.

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

Mach 슬라이더를 $10^{-2}$ 에서 $10^{-6}$ 으로 내리면 파란 입자는 0.5칸인데 주황 음파 전선은 도메인을 여러 바퀴 돈다.

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.

주황 전선이 한 스텝에 도는 바퀴 수가 아래 막대의 $\mathrm{CFL}_a$ 다. Mach가 내려갈수록 Gauss-Seidel 막대부터 빨간 선을 넘는다.

측정값이 손으로 얻은 식과 자릿수까지 맞는다. 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 근처까지 밀린다. 밀리는 위치가 다를 뿐, 밀어 올리는 것이 무엇인지는 같다.

CFL_a가 0.1이면 $\mathrm{CFL}_a^2 = 0.01$ 이라 GS 이득은 0.01, sweep은 3번이다. 0.5에서는 이득 0.25와 sweep 10, 1.0에서는 이득 1과 diverge다. 코드는 r이 0.999 이상이면 diverge로 친다. 이득 1은 오차가 줄지 않는다는 뜻이다. CFL 2에서 GS는 4, SOR은 1이라 둘 다 diverge다. CFL 10에서 GS는 100, SOR은 4.045. CFL 100에서 GS는 1e+04, SOR은 40이다. $\omega = 0.4$ 이면 $(1-\omega)^2 = 0.36$ , $\omega^2 = 0.16$ 이라 SOR 이득은 $\sqrt{0.36 + 0.16\,\mathrm{CFL}_a^2}$ 다. CFL 0.1에서 루트 안은 0.36에 0.0016을 더한 0.3616이고, 표의 0.601과 맞는다. CFL 10이면 루트 안은 0.36에 16을 더한 16.36이고, 표의 4.045와 맞는다. 물질 CFL 0.5를 고정하면 $\mathrm{CFL}_a = 0.5/M$ 이다. $M = 10^{-2}$ 에서 CFL_a는 50, GS 이득 2500, CG 9회다. $M = 10^{-3}$ 에서 CFL_a 500, GS 2.5e+05, SOR 200.3, CG는 아직 9다. $M = 10^{-4}$ 에서 CG가 12로 늘고, $10^{-5}$ 와 $10^{-6}$ 에서 CG는 17에서 멈춘다. GS 이득은 2.5e+09, 그다음 2.5e+11까지 간다. Schur의 CG가 4에서 9, 그리고 17까지인 것이 논문이 말한 CPU 시간의 약한 증가와 같은 자리다.

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

논문의 두 번째 실험은 시간간격을 고정하고 음속만 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을 받는다. 저마하 해석의 어려움은 물리가 아니라 이 숫자에 있다. 음파를 명시적으로 풀지 않으려고 음향과 대류를 나누는 접근이 나온 것도 같은 이유다. 위에서 주황 전선이 한 스텝에 도는 바퀴가 바로 이 $\mathrm{CFL}_a = \mathrm{CFL}_u/M$ 이다.

위 출력의 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가 여기서는 제 영역을 만난다. 바깥 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 보수 전처리기는 중간~높은 마하수에서도 동작하지만, 비용 대비 효과가 있는 것은 저마하 영역이라고 논문은 명시한다.

이득이 1을 넘으면 그 전처리기는 느린 게 아니다#

"physics-based preconditioner"라는 말은 보통 모호하게 쓰인다. 이 논문에서는 뜻이 좁다. 어떤 블록 결합이 시간간격이나 마하수와 함께 커지는지를 알고, 그 결합만은 근사하지 말라는 것이다. 음향 결합의 이득은 $\mathrm{CFL}_a^2$ 다. 표에서 그 제곱은 CFL 10에 100, CFL 100에 1e+04, $M = 10^{-6}$ 에 2.5e+11까지 간다. 1을 넘는 순간 오차는 sweep마다 불어난다. 반복 횟수를 늘려도 부호가 반대인 이득은 되돌아오지 않는다.

Schur 쪽 CG는 같은 표에서 4, 9, 그리고 17에서 멈춘다. 바깥 반복이 길어지는 자리는 여기다. $\omega = 0.4$ 의 SOR이 CFL 2에서 이득 1, CFL 10에서 4.045인 것도 같은 문턱이다. 문턱 아래에서만 그 전처리기를 쓴다. 음속이 커져도 시간간격을 살리는 다른 경로는 Orlando–Bonaventura의 AP-IMEX 쪽에 있다.

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