CFL을 1로 두면 만 번, 100으로 두면 폭발한다 — Switched Evolution Relaxation

정상상태 해석에서 CFL 수를 1로 두면 잔차가 만 번을 돌아도 안 떨어진다. 100으로 올리면 첫 반복에서 해가 폭발한다. 둘 다 흔한 풍경이다. 이 글은 그 사이를 자동으로 오가는 기법인 SER(Switched Evolution Relaxation)을 다룬다. 잔차가 줄어드는 만큼 CFL을 스스로 키워, 초반엔 안전하게 후반엔 Newton처럼 빠르게 착지하는 방법을 작은 정상 Burgers 솔버로 직접 만져본다.

여기서 다룰 CFL 수(시간 간격을 격자·파속으로 무차원화한 양)는 정상상태 솔버에서 더 이상 정확도가 아니라 수렴 속도와 안정성을 동시에 쥔 손잡이다.

정상상태는 가짜 시간으로 쫓는다#

정상상태 방정식은 $R(U) = 0$ 한 줄이다. $R$은 대류·확산·소스를 모은 공간 잔차다. 이걸 곧장 풀려면 비선형 대수계를 풀어야 한다. 대신 대부분의 CFD 코드는 가짜 시간 항을 붙인다.

여기서 $\tau$는 물리적 의미가 없는 의사시간(pseudo-time)이다. $\tau \to \infty$에서 $\partial U / \partial \tau \to 0$이 되면 $R(U)=0$, 즉 우리가 원하던 정상해다. 시간 정확도는 버리고 정상해만 노린다. 이를 의사시간 연속법(pseudo-transient continuation)이라 부른다.

가짜 시간을 음함수(backward Euler)로 전진시키고 잔차를 선형화하면 매 반복은 이렇게 된다.

여기서 $\partial R/\partial U$는 플럭스 Jacobian, $\Delta U = U^{n+1}-U^n$이다. 매 반복 좌변의 행렬을 한 번 풀어 $\Delta U$를 얻는다.

작게 시작해 Newton으로 착지하는 길#

이 한 줄에 SER의 모든 직관이 들어 있다. $\Delta\tau$를 양 극단으로 보내보자.

$\Delta\tau \to 0$이면 $I/\Delta\tau$가 행렬을 압도해 $\Delta U \approx -\Delta\tau\, R$이 된다. 이건 명시적 전진과 똑같다. 매우 안정적이지만 거북이 걸음이다.

$\Delta\tau \to \infty$이면 $I/\Delta\tau$가 사라져 $(\partial R/\partial U)\,\Delta U = -R$만 남는다. 이건 정확히 Newton 반복이다. 해 근처에서 2차 수렴(오차 제곱으로 줄어듦)을 한다. 단, 초기 추측이 나쁘면 발산한다.

그래서 이상적인 전략은 분명하다. 초반엔 $\Delta\tau$를 작게 두어 비선형성을 견디고, 잔차가 줄어 해 근처에 들어서면 $\Delta\tau$를 키워 Newton의 2차 수렴을 빨아들인다. 문제는 "언제 얼마나 키울 것인가"이다. SER은 그 판단을 잔차에게 맡긴다.

SER — 잔차의 역수로 CFL을 키운다#

Mulder와 van Leer가 1985년에 제안한 규칙은 놀랄 만큼 단순하다. 현재 CFL을 초기 잔차와 현재 잔차의 비로 키운다.

$\mathrm{CFL}_0$는 시작 CFL, $\lVert R_n \rVert$은 $n$번째 잔차의 L2 노름, $p$는 보통 1 안팎의 지수, $\mathrm{CFL}_{\max}$는 폭주를 막는 상한이다. 잔차가 10배 줄면 ($p=1$일 때) CFL도 10배 커진다. 의사시간 간격은 이 CFL과 국소 파속으로 환산한다.

분모는 대류 파속 $|u_i|$와 확산 파속 $2\nu/\Delta x$의 합이다. 셀마다 다르게 줄 수 있어 국소 시간전진(local time-stepping)과 자연스럽게 결합된다.

아래 시뮬레이션에서 직접 파라미터를 조작해보자. 고정 CFL과 SER을 같은 정상 Burgers 문제에 나란히 돌린 잔차 이력이다.

fixed CFL을 1로 내리면 호박색 곡선이 좀처럼 바닥에 닿지 못한다. 20까지 올리면 발산(diverged)으로 튄다. 반면 SER(청록색)은 $\mathrm{CFL}_0=1$로 안전하게 출발하고도 후반에 급격히 꺾여 내려간다. 지수 $p$를 1에서 1.5로 올리면 더 공격적으로 빨라지지만, 너무 키우면 SER조차 초반에 흔들린다.

Python — 정상 Burgers를 두 전략으로 푼다#

점성 Burgers의 정상해 $u u_x = \nu u_{xx}$를 푼다. 경계조건은 $u(0)=1$, $u(1)=-1$로, 가운데에 정지 충격층이 생긴다. 의사시간 음함수 전진 한 스텝은 삼중대각 행렬을 Thomas 알고리즘으로 푼다.

import numpy as np
 
def steady_residual(u, nu, dx):
    """R(u) = u u_x - nu u_xx, 경계 u(0)=1, u(1)=-1."""
    um = np.concatenate(([1.0], u[:-1]))   # 왼쪽 이웃
    up = np.concatenate((u[1:], [-1.0]))   # 오른쪽 이웃
    conv = u * (up - um) / (2 * dx)
    diff = nu * (up - 2 * u + um) / dx**2
    return conv - diff
 
def thomas_solve(a, b, c, d):
    """삼중대각계 (하/주/상 대각 a,b,c, 우변 d) 풀이."""
    n = len(d)
    b, d = b.copy(), d.copy()
    for i in range(1, n):
        m = a[i] / b[i - 1]
        b[i] -= m * c[i - 1]
        d[i] -= m * d[i - 1]
    x = np.empty(n)
    x[-1] = d[-1] / b[-1]
    for i in range(n - 2, -1, -1):
        x[i] = (d[i] - c[i] * x[i + 1]) / b[i]
    return x
 
def ptc_update(u, nu, dx, dt):
    """의사시간 음함수 한 스텝: (I/dt + dR/du) du = -R."""
    n = len(u)
    um = np.concatenate(([1.0], u[:-1]))
    up = np.concatenate((u[1:], [-1.0]))
    a = -u / (2 * dx) - nu / dx**2                       # 하대각
    b = (up - um) / (2 * dx) + 2 * nu / dx**2 + 1.0 / dt  # 주대각
    c = u / (2 * dx) - nu / dx**2                        # 상대각
    du = thomas_solve(a, b, c, -steady_residual(u, nu, dx))
    return u + du
 
def ser_schedule(r0, rn, cfl0, p, cfl_max):
    return min(cfl0 * (r0 / rn) ** p, cfl_max)
 
def converge(nu=0.02, ser=True, cfl0=1.0, p=1.0, cfl_max=1e5,
             N=80, tol=1e-9, max_iter=600):
    dx = 1.0 / (N + 1)
    x = np.linspace(dx, 1 - dx, N)
    u = 1 - 2 * x                          # 선형 초기 추측
    r0 = np.linalg.norm(steady_residual(u, nu, dx)) / np.sqrt(N)
    rn = r0
    for it in range(1, max_iter + 1):
        cfl = ser_schedule(r0, rn, cfl0, p, cfl_max) if ser else cfl0
        dt = cfl * dx / (np.abs(u).max() + 2 * nu / dx)
        u = ptc_update(u, nu, dx, dt)
        rn = np.linalg.norm(steady_residual(u, nu, dx)) / np.sqrt(N)
        if rn / r0 < tol:
            return u, it
    return u, max_iter
 
for tag, kw in [("fixed CFL=3", dict(ser=False, cfl0=3.0)),
                ("SER  p=1.0", dict(ser=True, cfl0=1.0, p=1.0))]:
    _, iters = converge(**kw)
    print(f"{tag:14s} -> {iters:4d} iters")

출력은 다음과 같다.

fixed CFL=3    ->  588 iters
SER  p=1.0     ->   34 iters

같은 정확도의 정상해에 도달하는 데 고정 CFL은 수백 번, SER은 수십 번이면 충분하다. 차이는 후반부에서 벌어진다. SER이 CFL을 수천까지 끌어올려 사실상 Newton 반복으로 변신하기 때문이다.

고정 CFL과 SER를 나란히 세우면#

SER이 어떻게 CFL을 끌어올리는지 스케줄 자체를 보자. 아래 그래프에서 직접 조작해보자. 청록색이 CFL 스케줄, 분홍색이 잔차다.

초반 잔차가 평평한 구간에선 CFL도 거의 $\mathrm{CFL}_0$에 머문다. 잔차가 무너지기 시작하는 순간 CFL이 거의 수직으로 치솟는다. 이게 SER의 본질이다. 해가 좋아질수록 더 큰 발걸음을 허락한다. cap을 20으로 낮추면 CFL이 천장에 막혀 후반 Newton 가속을 잃고, 다시 반복 수가 늘어난다.

현장에서 SER를 켤 때 빠지는 함정#

첫째, $\mathrm{CFL}_{\max}$는 반드시 둔다. 무한대로 풀어두면 잔차가 잠깐 튀는 순간 CFL이 폭주해 해를 날린다.

둘째, 잔차가 비단조(non-monotone)로 출렁일 때 $R_0/R_n$이 1보다 작아져 CFL이 거꾸로 줄 수 있다. 이건 버그가 아니라 안전장치다. 잔차가 다시 나빠지면 발걸음을 줄이는 게 옳다.

셋째, $p$는 공짜 가속이 아니다. $p>1$은 잔차가 조금만 줄어도 CFL을 크게 키워 위험을 키운다. 견고한 시작은 $p=1$, $\mathrm{CFL}_0 \approx 1$이다.

넷째, Jacobian이 부정확하면(근사 Jacobian, 또는 행렬프리) 큰 CFL에서 Newton의 2차 수렴이 1차로 떨어진다. 이때는 $\mathrm{CFL}_{\max}$를 낮춰 안정성을 사는 편이 낫다.

기억할 점#

• 정상상태 솔버에서 CFL은 정확도가 아니라 수렴 속도·안정성의 손잡이다. 작으면 거북이, 크면 폭발.
• SER은 $\mathrm{CFL}_n = \min(\mathrm{CFL}_0 (\lVert R_0\rVert/\lVert R_n\rVert)^p,\ \mathrm{CFL}_{\max})$로, 잔차가 줄면 자동으로 발걸음을 키운다.
• 의사시간 음함수 반복은 $\Delta\tau\to\infty$에서 Newton이 된다. SER은 그 한계로 부드럽게 데려가는 스케줄이다.