[논문 리뷰] GMRES가 1000번 돌고도 안 멈췄다 — 전속도 melt pool을 위한 물리 기반 블록 전처리

속도–압력 Schur 보수 전처리로 Mach가 0으로 가도 반복수가 늘지 않게 만든다

레이저로 금속을 녹이는 melt pool 시뮬레이션을 implicit으로 돌렸다. 첫 시간 스텝에서 GMRES가 1000번을 돌고도 잔차가 1e-2에서 멈춰 있었다. Mach 수를 1e-3으로 낮추자 증상은 더 심해졌다 — 음속이 유속의 1000배라 자코비안이 병들어 있었다. Weston·Nourgaliev·Delplanque·Barker(2019)는 이 정확한 문제를 풀었다. 비결은 속도–압력 Schur 보수를 전처리로 쓰는 것. 오늘은 그 핵심을 2-필드 toy로 직접 재현해, 전처리 하나로 Mach와 무관한 수렴이 어떻게 살아나는지 본다.

이 논문의 좌표#

저자: B. Weston, R. Nourgaliev, J.-P. Delplanque, A. T. Barker (LLNL / UC Davis)
저널: Journal of Computational Physics, 397, 108847 (2019)
제목: Preconditioning a Newton-Krylov solver for all-speed melt pool flow physics
한 줄: 상변화를 포함한 전속도(all-speed) 압축성 Navier–Stokes를 완전 음함수 JFNK로 풀되, 원시 변수 블록 Schur 보수 전처리로 저Mach 강성을 잡는다.

기존 전처리들 — element-block SOR, monolithic AMG, block Gauss–Seidel — 은 모두 시간 스텝이나 Mach가 극단으로 가면 반복수가 폭발했다. 논문의 기여는 vP-vT Schur 보수라는, 물리 결합 구조를 그대로 베낀 전처리다.

M이 0으로 갈 때 자코비안이 병드는 이유#

전속도란 한 코드로 비압축( $M\to 0$ )부터 초음속까지 다 푼다는 뜻이다. 문제는 저Mach 극한이다. 음향파 속도 $c$ 와 유속 $u$ 의 비가 $c/u = 1/M$ 로 벌어진다. $M=10^{-3}$ 이면 음파가 유체보다 1000배 빠르다.

완전 음함수로 가면 이 두 스케일이 한 자코비안 안에 섞인다. 압력–속도 결합 블록이 $1/M^2$ 로 커지며 행렬의 조건수가 폭발한다. Krylov 솔버는 고유값이 넓게 흩어진 행렬에서 기어간다. 이게 GMRES가 안 멈춘 이유다.

JFNK(Jacobian-Free Newton–Krylov)는 자코비안을 명시적으로 만들지 않고, 행렬–벡터 곱만 Fréchet 차분으로 근사한다.

J\boldsymbol{\kappa} \approx \frac{F(x + h\boldsymbol{\kappa}) - F(x)}{h}

$F$ 는 비선형 잔차, $\boldsymbol{\kappa}$ 는 Krylov 벡터, $h$ 는 작은 유한수다. 행렬을 저장 안 하니 메모리는 아끼지만, 전처리 없이는 저Mach에서 수렴이 죽는다. 전처리는 선택이 아니라 필수다.

보존 변수를 버리고 원시 변수 $(P,v,T)$ 로#

첫 수는 변수 교체다. 보존 변수 $U=(\rho, \rho v, E)$ 대신 원시 변수 $W=(P, v, T)$ 로 푼다. 저Mach에서 원시 변수가 훨씬 잘 조건화되기 때문이다.

\frac{\partial F(x)}{\partial U}\frac{\partial U}{\partial W}\,\delta W = -F(x)

$\partial U/\partial W$ 는 변수 변환 야코비안이다. 잔차 $F$ 는 여전히 보존 법칙을 만족하게 적혀 있어, 변수만 바꿔도 질량·운동량·에너지 보존은 깨지지 않는다. 이제 자유도를 필드별로 묶어 $3\times3$ 블록 행렬로 정렬한다.

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

각 블록은 한 원시 필드(속도·압력·온도)에 대응한다. 이 구조가 다음 수의 발판이다.

블록 LU와 속도–압력 Schur 보수#

압력–온도 결합( $M_{PT}, M_{TP}$ )은 약하다. 둘을 버리면 블록 LU 분해가 깔끔해지고, 핵심은 속도–압력 Schur 보수로 응축된다.

S_{vP} = M_{PP} - M_{Pv}\,M_{vv}^{-1}\,M_{vP}

$S_{vP}$ 는 속도를 소거한 뒤 압력만 남긴 유효 방정식이다 — SIMPLE 계열의 압력 보정 방정식과 같은 뼈대다. 전처리 $P$ 를 이 블록 상삼각으로 잡으면, $P^{-1}A$ 의 고유값이 1 근처로 뭉친다. 이론상 정확한 $S$ 를 쓰면 $A=LU$ 의 $U$ 가 곧 $P$ 가 되어, $P^{-1}A$ 는 멱영(unipotent) 하삼각과 닮아 모든 고유값이 정확히 1이다. GMRES가 두세 번에 끝난다.

Python — Mach를 낮추며 GMRES 반복수를 센다#

온도를 빼고 2-필드 $(v,P)$ toy로 핵심만 재현한다. 블록계는 $A=\begin{bmatrix} I & \beta G \\ D_v & I\end{bmatrix}$ , $\beta=1/M^2$ . 전처리는 정확한 Schur 보수 $S=I-\beta D_v G$ 다.

import numpy as np
from scipy.sparse import diags, eye, bmat
from scipy.sparse.linalg import gmres, LinearOperator, splu
 
n = 64
dx = 1.0 / n
 
def fwd_diff():
    # 주기 경계 전진차분 (압력 경사 연산자)
    D = diags([-1.0, 1.0], [0, 1], shape=(n, n)).tolil()
    D[-1, 0] = 1.0
    return (D / dx).tocsr()
 
G = fwd_diff()
Dv = (-G.T).tocsr()          # 발산 = -G^T → Dv·G ≈ 라플라시안
 
def all_speed_system(M):
    beta = 1.0 / M**2        # 저Mach 강성 계수
    A = bmat([[eye(n), beta * G], [Dv, eye(n)]]).tocsr()
    return A, beta
 
def schur_preconditioner(beta):
    # 속도-압력 Schur 보수 S = M_PP - M_Pv M_vv^{-1} M_vP = I - β Dv G
    S = (eye(n) - beta * (Dv @ G)).tocsc()
    luS = splu(S)
    def apply(r):
        ru, rp = r[:n], r[n:]
        zp = luS.solve(rp)               # S z_p = r_p
        zu = ru - beta * (G @ zp)        # z_u = r_u - β G z_p
        return np.concatenate([zu, zp])
    return LinearOperator((2 * n, 2 * n), matvec=apply)
 
def count_iters(M, use_prec):
    A, beta = all_speed_system(M)
    b = np.random.default_rng(0).standard_normal(2 * n)
    P = schur_preconditioner(beta) if use_prec else None
    it = {"k": 0}
    x, info = gmres(A, b, M=P, rtol=1e-8, maxiter=500, restart=2 * n,
                    callback=lambda xk: it.__setitem__("k", it["k"] + 1))
    return it["k"]
 
print(f"{'Mach':>8}  {'no-prec':>8}  {'Schur':>6}")
for M in [1.0, 1e-1, 1e-2, 1e-3]:
    print(f"{M:>8.0e}  {count_iters(M, False):>8d}  {count_iters(M, True):>6d}")

결과는 명백하다. 전처리 없이는 $M$ 이 작아질수록 반복수가 수십·수백으로 치솟고 $M=10^{-3}$ 에서는 한도 안에 수렴조차 못 한다. Schur 전처리는 $M$ 과 무관하게 두세 번이면 끝난다. 논문이 보여준 "Mach 독립 수렴"이 이 한 줄짜리 전처리에서 나온다.

인터랙티브 — 음향과 대류, 벌어지는 두 속도#

왜 강성이 생기는지부터 눈으로 보자. 아래 시뮬레이션에서 Mach를 직접 낮춰보자. 위 트랙의 음향파는 $c=1/M$ 로, 아래 대류 추적자는 $u=1$ 로 달린다.

Mach number M = 3.2e-2

explicit acoustic CFL forces Δt ∝ M·Δx — the smaller M, the stiffer the coupled system.

$M=10^{-2}$ 로 내리면 음향파가 트랙을 100바퀴 도는 동안 대류 추적자는 한 칸 움직인다. 이 100:1의 속도 간극이 자코비안의 조건수로 그대로 옮겨붙는다. 명시적이라면 음향 CFL이 $\Delta t \propto M\Delta x$ 로 시간 스텝을 짓밟는다.

인터랙티브 — 전처리가 잔차를 무너뜨리는 자리#

이제 실제 GMRES 잔차를 보자. 아래에서 Mach 슬라이더를 움직이면 전처리 유무 두 곡선이 다시 그려진다.

Mach number M = 1.0e-2 (β = 1/M² = 1.0e+4)

unpreconditioned: 0 iters to 1.0e+0 · Schur: 0 iters to 1.0e+0

주황색(전처리 없음)은 $M$ 을 낮출수록 평평해지며 멈춘다. 청록색(Schur 전처리)은 $M$ 이 어디든 두세 번 만에 1e-11로 떨어진다. $M=10^{-3}$ 에서 두 곡선의 격차가 가장 극적이다 — 물리 결합을 베낀 전처리가 강성을 통째로 흡수한 결과다.

직접 짜며 의심스러웠던 세 가지#

첫째, 정확한 $S$ 는 toy에서만 공짜다. 여기선 $S=I-\beta D_v G$ 를 LU로 정확히 풀어 두 번 만에 수렴했다. 실전에서 $S_{vP}=M_{PP}-M_{Pv}M_{vv}^{-1}M_{vP}$ 의 $M_{vv}^{-1}$ 은 밀집 역행렬이라 만들 수 없다. 논문은 $S$ 를 근사로 풀고(AMG 한 V-cycle 등), 그 근사 품질이 곧 외부 반복수를 좌우한다. "Mach 독립"의 진짜 비용은 여기 숨어 있다.

둘째, 결합을 어디까지 버릴지. 논문은 $M_{PT}, M_{TP}$ 를 버려도 강건성이 안 깨진다고 한다. 이건 melt pool 물성에 특화된 가정이다. 강한 압축성 열결합(예: 폭연·연소)에서는 같은 절단이 수렴을 죽일 수 있다. 전처리의 단순화는 물리에 종속적이다.

셋째, JFNK의 Fréchet 차분 $h$ . $h$ 가 너무 작으면 반올림이, 크면 절단 오차가 행렬–벡터 곱을 오염시킨다. 전처리가 좋아도 이 $h$ 튜닝이 어긋나면 Newton이 정체한다. OpenFOAM의 분리형(SIMPLE) 대신 완전 결합 음함수를 쓸 때 가장 먼저 만나는 함정이다.

닫으며 — 전처리는 공짜가 아니다#

GMRES가 안 멈췄던 그 밤의 교훈은 단순했다. 저Mach 강성은 자코비안을 더 세게 푼다고 풀리지 않는다. 결합 구조 자체를 베낀 전처리가 있어야 한다. Schur 보수는 압력을 위한 유효 방정식이고, 그것을 정확히 풀면 고유값이 1로 모인다. 다만 실전의 $S$ 는 근사이고, 그 근사를 얼마나 싸고 정확하게 푸느냐가 다음 싸움이다. 다음엔 $S$ 를 AMG로 근사했을 때 외부 반복수가 어떻게 변하는지, 같은 toy에 한 단계 더 얹어 볼 차례다.