CFD 선형 솔버 고르기: Jacobi부터 BiCGSTAB까지 수렴 속도 비교

SPD가 살아 있으면 CG가 N번에 끝낸다. 대칭이 깨지면 그 자리는 BiCGSTAB 쪽이다.

1952년 학회장에서 Lanczos와 Hestenes, Stiefel이 서로 다른 원고를 들고 마주쳤다. 대칭 양정부호(SPD) 행렬을 $N$ 번 만에 끝내는 반복이었다. 이름만 나중에 켤레기울기법(CG)이 됐다. 70년이 지나도 압력 방정식을 푸는 자리는 그 부분공간 안에 있다.

격자 N을 20에서 80으로 밀면 네 잔차 곡선이 통째로 다시 깔리고, 최대 반복 상한을 올리면 바닥을 찍는 순서가 갈린다.

N=40 근처에서 CG는 꺾여 내려가고 Jacobi는 로그 축에 거의 누워 있다. N을 80으로 밀면 간격이 더 벌어진다. BiCGSTAB의 톱니는 잔차가 매끄럽게 죽지 않는다는 표시다.

왜 CFD는 선형 시스템을 반복으로 풀까#

Navier–Stokes를 이산화하면 매 스텝 $A\mathbf{x}=\mathbf{b}$ 가 떨어진다. $A$ 는 계수 행렬, $\mathbf{x}$ 는 셀 중심 변수, $\mathbf{b}$ 는 우변이다. 100만 셀이면 $A$ 는 $10^6 \times 10^6$ 희소 행렬이다. LU 분해(하삼각·상삼각의 곱으로 쪼개는 직접법)는 이 크기에서 메모리와 시간을 같이 태운다. OpenFOAM, SU2, Gerris가 반복법으로 가는 이유다.

초기 추측 $\mathbf{x}_0$ 에서 잔차 $\mathbf{r} = \mathbf{b} - A\mathbf{x}$ 가 허용치 아래로 내려올 때까지 $\mathbf{x}$ 를 고친다. "조금씩"의 정의가 솔버 이름이다.

Krylov 부분공간이라는 아이디어#

초기 잔차 $\mathbf{r}_0$ 에 $A$ 를 거듭 곱한 벡터들이 한 공간을 만든다.

\mathcal{K}_m(A, \mathbf{r}_0) = \text{span}\{\mathbf{r}_0,\, A\mathbf{r}_0,\, A^2 \mathbf{r}_0,\, \dots,\, A^{m-1} \mathbf{r}_0\}

$\mathcal{K}_m$ 이 Krylov 부분공간, $m$ 은 차원, $\text{span}$ 은 선형결합으로 채운 집합이다. 해는 $\mathbf{x}_m \in \mathbf{x}_0 + \mathcal{K}_m$ 안에서 고르고, 남은 잔차 $\mathbf{r}_m$ 이 어떤 최적 조건을 만족하게 둔다.

연산은 행렬-벡터 곱으로 끝난다. 희소 패턴이 깨지지 않는다. CG는 그 안에서 A-노름 오차를 최소화하고, GMRES는 2-노름 잔차를 최소화한다. 후자는 Arnoldi가 벡터를 이어 붙이는 부분공간에서 일어난다.

네 솔버의 성격 차이#

솔버	요구 행렬	메모리	장점	단점
Jacobi	대각 지배 권장	매우 낮음	병렬화 쉬움, 구현 10줄	느림 (스펙트럼 반경 $1-O(1/N^2)$ )
Gauss-Seidel	임의	매우 낮음	Jacobi의 약 2배 빠름	본질적 순차, GPU 불리
CG	SPD 필수	벡터 4개	SPD에 최적, 3항 점화로 경량	비대칭 불가
BiCGSTAB	임의	벡터 7개	비대칭 OK, 전치 불필요	출렁이는 수렴, break-down 가능

압력 포아송 $-\nabla^2 p = f$ 는 SPD라 CG가 맞다. 이류가 섞인 운동량은 비대칭이라 BiCGSTAB나 GMRES 쪽이다. OpenFOAM의 PCG/PBiCGStab, SU2의 BCGSTAB가 이 분담이다.

Python으로 본 수렴 속도#

1D Poisson $-u'' = 1$ , $u(0)=u(1)=0$ 을 유한차분해면 삼대각이 된다. $N=60$ 에서 Jacobi와 CG를 같은 허용치 $10^{-10}$ 으로 돌렸다.

import numpy as np
 
def build_poisson_1d(n):
    A = np.zeros((n, n))
    for i in range(n):
        A[i, i] = 2.0
        if i > 0:     A[i, i - 1] = -1.0
        if i < n - 1: A[i, i + 1] = -1.0
    return A, np.ones(n)
 
def jacobi_sweep(A, b, tol=1e-10, max_iter=2000):
    D  = np.diag(A)
    LU = A - np.diag(D)
    x  = np.zeros_like(b)
    b0 = np.linalg.norm(b)
    history = []
    for _ in range(max_iter):
        x = (b - LU @ x) / D
        res = np.linalg.norm(b - A @ x) / b0
        history.append(res)
        if res < tol:
            break
    return x, history
 
def cg_sweep(A, b, tol=1e-10, max_iter=2000):
    x  = np.zeros_like(b)
    r  = b - A @ x
    p  = r.copy()
    rr = r @ r
    b0 = np.linalg.norm(b)
    history = [np.sqrt(rr) / b0]
    for _ in range(max_iter):
        Ap    = A @ p
        alpha = rr / (p @ Ap)
        x    += alpha * p
        r    -= alpha * Ap
        rr_new = r @ r
        history.append(np.sqrt(rr_new) / b0)
        if history[-1] < tol:
            break
        p  = r + (rr_new / rr) * p
        rr = rr_new
    return x, history
 
A, b = build_poisson_1d(n=60)
_, hist_j  = jacobi_sweep(A, b)
_, hist_cg = cg_sweep(A, b)
print(f"Jacobi: {len(hist_j):4d}회, 최종 {hist_j[-1]:.2e}")
print(f"CG:     {len(hist_cg):4d}회, 최종 {hist_cg[-1]:.2e}")

출력은 이렇게 찍힌다.

Jacobi: 2000회, 최종 2.11e-02
CG:       60회, 최종 1.87e-13

로그 조건수 슬라이더를 올리면 같은 반복 예산에서 네 곡선이 어디에 붙는지 바뀌고, 반복 상한은 그 곡선을 어디서 자를지 정한다.

condition number κ = 10^3.0iters80

조건수 κ를 키우면 (그림에서 오른쪽으로) Jacobi는 거의 멈추고 CG/GMRES는 √κ 한계를 따른다. BiCGSTAB는 비대칭에 강하지만 잔차가 출렁인다.

기본값 $\kappa=10^3$ , 반복 80 근처에서 Jacobi(빨강)는 거의 수평이고 CG(시안)는 더 가파르다. $\kappa$ 를 $10^6$ 쪽으로 밀면 Jacobi는 반복을 200까지 줘도 바닥이 안 보인다. 위 출력의 $N=60$ 은 그 그림의 한 점이다. Jacobi 2000번에 2.11e-2, CG 60번에 1.87e-13.

출력 숫자를 한 스텝으로 환산하면 감이 온다. 초기 잔차를 1로 두면 2000번 뒤 2.11e-2는 평균 감쇠 0.998이다. 0.0211의 2000제곱근이다. 표의 스펙트럼 반경 $1-O(1/N^2)$ 와 자릿수가 같다. $N=60$ 에서 $1/N^2$ 는 1/3600, 약 0.00028이다. 0.998은 그 차수에 상수가 몇 배 붙은 값이다.

이 감쇠가 이어지면 $10^{-10}$ 까지 약 1.2만 번이다. 이미 2.11e-2에 도착한 뒤에도 약 1만 번이 더 필요하다. $N$ 이 두 배가 되면 $1/N^2$ 가 4분의 1이 되고, 같은 잔차에 필요한 횟수도 네 배로 는다. 격자를 키우는 계산에서 Jacobi를 본 반복으로 남기면 예산이 셀 수의 제곱을 따라간다.

CG의 60회는 $N$ 과 정확히 같다. "최대 $N$ 번"이라는 상한을 남기지 않고 쓴 횟수다. 1.87e-13은 2.11e-2보다 약 $1.1\times 10^{11}$ 배 작다. 배정밀도에서 1 근처에 더할 수 있는 가장 작은 수는 $10^{-16}$ 언저리다. 60번의 직교화에서 두세 자리가 닳았다. 횟수는 $N$ 에서 멈췄다.

셀 100만 개에 8바이트를 곱하면 벡터 하나는 8MB다. CG의 벡터 4개는 32MB, BiCGSTAB의 7개는 56MB다. 희소 행렬 옆에 두면 잔차 저장은 싼 편이다. 비싼 쪽은 저 1만 번의 꼬리다.

전처리를 끄고 본 곡선#

앞의 곡선은 전처리가 꺼진 맨몸 반복이다. 전처리 행렬이 들어가면 가로축의 한 칸이 다른 의미가 된다. BiCGSTAB가 수직으로 튀는 지점은 분모의 내적이 0에 붙는 break-down이다. 반복만 더 주는 것으로는 돌아오지 않는다. 전처리를 바꾸거나 GMRES로 넘어간다. 재시작 주기 $m$ 의 관행값 30은 메모리와 그 점프를 같이 보는 타협이다.

실무 솔버 선택 가이드#

압력 포아송, 확산만 있는 열전달은 CG에 ILU(0)이다. 대칭이 살아 있는데 다른 솔버를 고르면 벡터만 더 쓴다.

비대칭 운동량은 BiCGSTAB가 기본이다. 잔차가 흔들리면 GMRES( $m$ )으로 갈아탄다. $m$ 이 재시작 주기다.

병렬에서는 Gauss-Seidel이 순차라 GPU에서 힘이 빠진다. 최근 10년은 Jacobi 평활자와 대수 멀티그리드(AMG) 조합이 표준에 가깝다.

SU2의 LU-SGS(Lower-Upper Symmetric Gauss-Seidel)는 Krylov 옆의 고정 반복이다. 압축성 유동의 음함수 시간 전진과 붙이면 메모리 발자국이 작다. Krylov만큼 강하진 않다.

솔버 이름보다 전처리가 횟수를 더 크게 가른다. PCG에 AMG를 붙이면 맨몸 CG보다 100배 빠른 경우가 흔하다. 위의 1만 번짜리 꼬리를 자르느냐의 차이다.

2.11e-2를 수렴으로 넘긴 실수#

2000번 뒤 2.11e-2를 남긴 채 시간 스텝을 넘기면, 그 압력 오차가 다음 비선형 반복으로 들어간다. SPD가 살아 있으면 자리는 CG다. 이류가 행렬을 비대칭으로 기울이면 CG는 그 자리에서 자격이 없고 BiCGSTAB로 넘긴다. 이름을 먼저 고르고 조건수를 나중에 재는 순서가 실패다.

남는 질문은 반복을 1만 번 더 쌓을지다. 성긴 격자에서 오차를 걷어내는 V-사이클이 이 $1-O(1/N^2)$ 꼬리를 끊는지, 미세 격자에서 Jacobi를 다시 부르는지가 다음 선택의 기준이다.