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

Krylov 부분공간 4대 솔버의 성격과 수렴 속도를 코드·시뮬레이션으로 비교한다

1952년 Cornelius Lanczos와 Magnus Hestenes, Eduard Stiefel은 한 학회에서 우연히 같은 결과를 독립적으로 들고 나타났다. 대칭 양정부호(SPD) 행렬을 $N$ 번 반복만에 푸는 알고리즘이었다. 지금 우리가 켤레기울기법(Conjugate Gradient, CG) 이라 부르는 그 방법이다. 흥미로운 건 70년이 지난 오늘도 CFD 솔버의 심장에는 이 아이디어가 거의 그대로 살아 있다는 점이다 — Krylov 부분공간에서 해를 찾는다. 이 포스트는 CFD에서 가장 자주 쓰이는 네 솔버(Jacobi, Gauss-Seidel, CG, BiCGSTAB)의 수렴 속도를 Python과 인터랙티브 시뮬레이션으로 비교한다. 끝까지 읽으면 "이 문제엔 어느 솔버가 맞는가"를 실무에서 바로 판단할 수 있다.

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

Navier–Stokes를 이산화하면 매 스텝마다 $A\mathbf{x}=\mathbf{b}$ 형태의 큰 선형 시스템이 떨어진다. 여기서 $A$ 는 계수 행렬, $\mathbf{x}$ 는 셀 중심 변수 벡터, $\mathbf{b}$ 는 우변 소스다. 3D 격자 100만 셀이면 $A$ 는 $10^6 \times 10^6$ 크기의 희소 행렬이 된다. LU 분해(행렬을 하삼각·상삼각 곱으로 쪼개는 직접법)는 이 크기에서 메모리와 CPU를 함께 태운다. OpenFOAM, SU2, Gerris 같은 오픈소스 CFD 코드가 예외 없이 반복법(iterative method) 을 쓰는 이유다.

반복법의 기본 아이디어는 단순하다. 초기 추측 $\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}$ 은 선형결합으로 채운 공간이다. Krylov 계열 솔버는 해를 $\mathbf{x}_m \in \mathbf{x}_0 + \mathcal{K}_m$ 안에서 찾되, 남은 잔차 $\mathbf{r}_m$ 이 어떤 최적성 조건을 만족하도록 결정한다.

해를 부분공간에 가두면 이점 두 가지가 있다. 첫째, 행렬-벡터 곱 외에 다른 연산이 없다 — 희소 구조를 그대로 살린다. 둘째, 이론적으로 해당 부분공간 안에서 최적인 해를 뽑을 수 있다. CG는 A-노름 오차를, GMRES는 2-노름 잔차를 Krylov 공간 안에서 정확히 최소화한다.

네 솔버의 성격 차이#

솔버	요구 행렬	메모리	장점	단점
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$ 을 유한차분으로 이산화하면 삼대각 행렬이 떨어진다. Jacobi와 CG를 직접 비교해 보자.

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

Jacobi가 2000번을 돌아도 $10^{-2}$ 를 못 넘는 동안, CG는 $N=60$ 번 만에 기계 정밀도에 도달한다. 이론적으로 CG는 최대 $N$ 번에 정확해에 닿는다 (부동소수점 오차 무시 시). SPD 앞에서 CG는 거의 초능력이다.

직접 수렴 경쟁을 보자#

아래 시뮬레이션에서 직접 조작해보자. 슬라이더로 격자 크기 $N$ 과 최대 반복 횟수를 바꾸면 네 솔버의 잔차 곡선이 로그 스케일로 다시 그려진다.

격자 N40최대 반복200

$N$ 을 80까지 키우면 Jacobi와 Gauss-Seidel의 곡선이 거의 수평에 가까워진다. 반면 CG는 계단식으로 떨어져 $N$ 근처에서 급격히 수렴한다. BiCGSTAB는 출렁이며 내려오는데, 이 "들쭉날쭉함"이 break-down(분자 또는 분모가 0으로 폭발하는 수치 파괴) 위험의 시각적 신호다. 실전에서 BiCGSTAB가 갑자기 발산하면 전처리자를 강화하거나 GMRES로 갈아타는 게 정석이다.

실무 솔버 선택 가이드#

압력 포아송, 확산만 있는 열전달: CG + ILU(0) 전처리자. 대칭이 보장되면 다른 선택지는 낭비다.

비대칭 운동량, 이류-확산 혼합: BiCGSTAB가 기본. 수렴이 흔들리거나 막히면 GMRES( $m$ )으로 갈아탄다. $m$ 은 재시작 주기이며 메모리-속도 트레이드오프를 결정한다. 업계 관행은 $m=30$ 근처다.

병렬 CFD: Gauss-Seidel은 순차 특성 때문에 GPU에서 힘을 못 쓴다. 최근 10년 표준은 Jacobi 기반 평활자(smoother)와 대수 멀티그리드(AMG)의 조합이다.

SU2의 LU-SGS: SU2는 Krylov 외에 LU-SGS(Lower-Upper Symmetric Gauss-Seidel)라는 고정 반복 방법도 탑재한다. 압축성 유동에서 묵시적(implicit) 시간 전진과 결합할 때 메모리 발자국이 훨씬 작다. 대신 Krylov만큼 강건하지는 않다.

Tip: 실무에선 "솔버 선택"보다 "전처리자 선택"이 수렴 속도를 훨씬 크게 바꾼다. CG를 맨몸으로 쓰기보다 PCG(Preconditioned CG) + AMG 조합이 100배 빠른 경우가 흔하다.

핵심 3줄 요약#

Jacobi와 Gauss-Seidel은 단순하지만 CFD 실제 문제에서 수렴률이 $1-O(1/N^2)$ 로 형편없다. 문제가 커질수록 반복 횟수가 제곱으로 는다.
Krylov 계열(CG, BiCGSTAB, GMRES)은 행렬-벡터 곱만으로 $N$ 번 안에 기계 정밀도에 도달한다. SPD는 CG, 비대칭은 BiCGSTAB 또는 GMRES가 실무 기본값이다.
솔버보다 전처리자가 중요하다. ILU(0), AMG, Jacobi 스무더 — 이 셋의 조합이 10배에서 100배 속도를 가른다.