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

정상상태 잔차는 CFL에 예민하다. 1로 두면 만 번을 돌아도 잘 안 떨어진다. 100으로 올리면 첫 반복에서 해가 터진다. 그 사이를 잔차가 고르는 방법이 SER(Switched Evolution Relaxation)이다. 줄면 키우고, 튀면 다시 줄인다.

fixed CFL 슬라이더를 1에서 20으로 밀면 호박색 잔차가 diverged로 튀고, SER은 시작 CFL과 지수만 바꾼다.

$\mathrm{CFL}_0=1$에 $p=1$이면 청록 곡선이 후반에 꺾여 내려간다. $p$를 1.5로 올리면 더 가파르지만, 초반 잔차가 출렁이는 순간 청록도 같이 흔들린다.

여기서 다룰 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)과 자연스럽게 결합된다.

격자 수 N=80이면 간격은 1/81이다. ν=0.02에서 확산 파속 $2\nu/\Delta x$는 3.24까지 올라간다. 초기 추측의 최댓값은 1을 넘지 못한다. 분모는 대류 $|u_i|$보다 확산이 크다. CFL을 3으로 고정하면 $\Delta\tau$는 간격 1/81을 4 남짓으로 나눈 뒤 3을 곱한 값이라, 0.01을 밑돈다. 같은 분모에 상한 1e5를 넣으면 $\Delta\tau$는 수백이 된다. 좌변의 $I/\Delta\tau$는 그때 거의 사라진다. 코드는 칸마다 $|u_i|$를 쓰지 않는다. 속도 절댓값의 최댓값을 모든 칸이 나눠 갖는다. 충격층 한가운데, 속도가 0에 가까운 칸도 가장 빠른 칸의 파속을 빌린다. 본문의 국소 시간전진보다 한 박자 보수적이다.

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

점성 Burgers의 정상해 $u u_x = \nu u_{xx}$를 푼다. 경계조건은 $u(0)=1$, $u(1)=-1$로, 가운데에 정지 충격층이 생긴다. 의사시간 음함수 전진 한 스텝은 삼중대각 행렬을 Thomas 알고리듬으로 푼다. 1D 확산에서 그 역추적이 무조건 안정을 만드는 자리와, 여기서 좌변에 시간항을 더하는 자리는 같다.

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

cap 슬라이더를 내리면 출력의 34걸음이 어디서 다시 길어지는지, 청록 CFL과 분홍 잔차가 같이 움직인다.

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

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

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

588을 34로 나누면 약 17배다. 공차는 현재 잔차를 초기 잔차로 나눈 값이 1e-9보다 작아지는 지점이다. $p=1$, 시작 CFL이 1, 상한이 1e5이면 잔차 비가 1e5가 되는 순간 CFL이 천장에 붙는다. 그때 잔차 비는 1e-5다. 공차 1e-9까지는 네 자릿수가 더 남는다. 그 구간은 CFL이 1e5에 고정된 채 행렬이 Newton 쪽에 붙는 걸음이다. 고정 CFL 3에는 그 구간이 없어서 588걸음을 다 쓴다.

스케줄은 한 줄로 적힌다. $\mathrm{CFL}_n = \min(\mathrm{CFL}_0 (\lVert R_0\rVert/\lVert R_n\rVert)^p,\ \mathrm{CFL}_{\max})$. 잔차가 줄면 발걸음이 커지고, 상한에서 멈춘다.

현장에서 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}$를 낮춰 안정성을 사는 편이 낫다.

상한을 비우면 후반이 Newton 전에 죽는다#

$\mathrm{CFL}_{\max}$를 비우면 잔차가 한 번 튀는 순간 CFL이 잔차 비를 따라 끝없이 커진다. 의사시간 음함수 반복은 $\Delta\tau\to\infty$에서 Newton이 된다. 그 2차는 해 근처에서만 나온다. 초반에 그 한계로 뛰어들면 행렬이 해를 밀어 올린다. 코드의 상한 1e5는 잔차가 십만 배 줄어든 뒤에야 그 영역에 들어가게 막는 문이다.

행렬을 매번 정확히 만들 수 없을 때도 같은 문이 필요하다. 큰 CFL의 2차가 1차로 떨어지면, 상한을 더 낮추는 쪽이 남는다. 스윕으로 그 행렬을 근사하는 길은 양함수가 터진 뒤의 LU-SGS에 따로 있다. 거기서는 행렬을 조립하지 않고도 음함수 걸음을 유지한다.