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

고유값이 1로 모이는 건 정확한 Schur 보수일 때뿐이다. AMG로 S를 느슨하게 풀면 외부 반복이 다시 벌어진다.

레이저로 금속을 녹이는 melt pool을 implicit으로 돌렸다. 첫 시간 스텝에서 GMRES가 1000번을 돌고도 잔차가 1e-2에 붙어 있었다. Mach를 1e-3으로 내리자 음속이 유속의 1000배가 됐다. 자코비안이 그 비에서 병들었다. Weston·Nourgaliev·Delplanque·Barker(2019)는 속도–압력 Schur 보수를 전처리로 써서 이 정체를 풀었다. 오늘은 그 전처리를 2-필드 toy로 옮겨, 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바퀴 도는 동안 대류 추적자는 한 칸이다. $M=10^{-3}$ 이면 비가 1000:1이 된다. 잔차가 1e-2에 붙어 있던 밤과 같은 간극이다.

이 논문의 좌표#

저자: 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가 안 멈춘 이유다.

$M$ 을 코드의 네 점으로 찍으면 $\beta=1/M^2$ 가 바로 나온다. $M=1$ 이면 $\beta=1$ . $M=10^{-1}$ 이면 $\beta=10^{2}$ . $M=10^{-2}$ 이면 $\beta=10^{4}$ . $M=10^{-3}$ 이면 $\beta=10^{6}$ . 압력–속도 블록의 칸 하나가 백만 배로 커진다. 격자 수는 $n=64$ , 칸 너비는 $dx=1/64$ 로 네 경우 모두 같다. 행렬 크기도 $2n=128$ 로 같다. 달라지는 것은 $\beta$ 를 곱한 스펙트럼뿐이다.

음향 CFL로 시간 스텝을 조이면 $\Delta t \propto M\Delta x$ 다. $M$ 을 1에서 $10^{-3}$ 으로 내리면 허용 스텝이 1000분의 1이 된다. 완전 음함수는 그 제한을 풀고, 대가로 조건수를 $\beta$ 에 넘긴다. 같은 강성을 전처리가 거꾸로 키운 자리는 음속 CFL 10의 블록 전처리에 따로 적혀 있다.

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가 두세 번에 끝난다.

그 두세 번은 정지 조건과도 맞물린다. 아래 코드의 허용 잔차는 rtol=1e-8, 반복 상한은 maxiter=500, restart 길이는 2n=128이다. 고유값이 한 점에 모이면 Krylov 부분공간은 두 걸음에서 닫힌다. 모이지 않으면 500번 상한에 먼저 닿는다.

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}$ 에서는 500번 한도 안에 들어오지 못한다. Schur 쪽은 $S=I-\beta D_v G$ 를 LU로 한 번 풀어 두고, GMRES가 두세 번이면 잔차가 1e-11까지 떨어진다. $M=1$ ( $\beta=1$ )과 $M=10^{-3}$ ( $\beta=10^{6}$ )의 Schur 반복수가 같다는 것이, 고유값 1 뭉침의 수치다.

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

앞에서 트랙은 이미 돌렸다. $M=10^{-2}$ 의 100:1은 $c/u=1/M$ 그대로다. 명시적 CFL이 $\Delta t \propto M\Delta x$ 로 스텝을 짓밟는 자리도 그 화면의 속도 차와 같다. 자코비안이 받는 것은 그 비를 제곱한 $\beta$ 다. $M=10^{-2}$ 에서 $\beta=10^{4}$ , $M=10^{-3}$ 에서 $\beta=10^{6}$ 이다.

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

이제 실제 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}$ 에서 두 곡선의 격차가 가장 크다. $\beta=10^{6}$ 인 칸을, 물리 결합을 베낀 전처리가 통째로 흡수한 결과다.

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

정확한 $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가 1000번, 잔차 1e-2에서 멈춘 밤은 자코비안을 더 세게 푼다고 끝나지 않았다. 결합 구조 자체를 베낀 전처리가 있어야 했다. Schur 보수는 속도를 소거한 압력 방정식이고, 그것을 정확히 풀면 고유값이 1로 모인다. 실전의 $S_{vP}$ 안 $M_{vv}^{-1}$ 은 밀집 역행렬이라 논문은 AMG 한 V-cycle로 근사한다. 그 근사가 거칠면 외부 반복이 두세 번에서 다시 벌어진다.

야코비안을 저장하지 않고 행렬–벡터 곱만 쓰는 갈림은 연산자 오버로딩 자동미분과 맞닿아 있다. JFNK는 그 곱을 Fréchet 차분 $h$ 로 대체한다. $h$ 가 어긋나면 전처리가 좋아도 Newton은 그 자리에서 멈춘다. 같은 toy에서 $S$ 를 AMG로 느슨하게 풀 때, 외부 반복이 두세 번에서 얼마나 벌어지는지가 다음 측정이다.