정상 상태용 예조건을 이중 시간 루프에 넣었더니 내부 반복이 19배 늘었다 — 저마하 음향 펄스 실험

빨라지라고 켠 옵션이 19배 느리게 만들었다#

저마하 압축성 코드에서 예조건(preconditioning)은 "켜면 빨라지는" 옵션으로 통한다. 그래서 비정상 계산에도 그대로 켜고 싶어진다. 실제로 켜 보면 두 가지 일이 생긴다. 실시간에 넣으면 압력 펄스가 제자리에 멈춘다. 이중 시간(dual time) 루프 안에 넣으면 답은 돌아오지만, 한 스텝의 내부 반복이 16회에서 305회로 늘어난다.

이 글은 2018년에 정리한 코드 구조 슬라이드에서 출발한다. 슬라이드에는 "예조건은 정상 상태에만 쓴다. 비정상에서는 물리적으로 맞지 않다"는 한 줄이 있다. 바로 다음 장에는 "Implicit + Unsteady + dual time + Preconditioning"이 최종 형태로 적혀 있다. 이 두 줄이 어떻게 함께 성립하는지, 그리고 그 조합이 언제 이득이고 언제 손해인지를 1차원 실험으로 재 본다.

지배방정식 여덟 갈래#

슬라이드 첫 장은 경우의 수를 센다. 예조건 유무, 정상·비정상, 양해·음해. 2×2×2로 여덟 갈래다. 코드는 이 여덟 갈래를 하나의 식으로 묶으려 한다.

$\tau$는 가짜 시간(pseudo time, 수렴용으로만 쓰는 반복 변수), $t$는 실제 시간이다. $\Gamma$는 예조건 행렬, 둘째 항은 2차 후방차분(BDF2)이다. 정상 문제는 둘째 항을 지우면 되고, 예조건이 없으면 $\Gamma=I$를 넣으면 된다. 가짜 시간에서 수렴하면 첫 항이 0이 된다. 그러면 $\Gamma$가 무엇이든 남는 식은 BDF2로 이산화한 원래 방정식이다. 이것이 슬라이드의 두 줄이 동시에 참인 이유다.

Γ가 바꾸는 것은 파동 속도다#

선형화한 1차원 오일러 방정식을 $Q=(p,u)$로 쓰자. 밀도와 음속을 1로 두면 $A=\begin{pmatrix}u_0&1\\1&u_0\end{pmatrix}$이고 고유값은 $u_0\pm c$다. Weiss–Smith 형태의 예조건은 압력의 시간 미분 계수만 바꾼다. 여기서는 $\Gamma=\mathrm{diag}(\delta,1)$로 단순화한다.

$\delta=1$이면 원래의 $u_0\pm c$가 나온다. 정상 상태용 선택 $\delta=M^2$를 넣으면 근호 안이 $\approx 5u_0^2$가 되어 $\lambda\approx u_0(1\pm\sqrt5)/2$가 된다. 가장 빠른 파동과 유속의 비가 $1/M$에서 약 1.6으로 줄어든다. 정상 문제에서 수렴이 빨라지는 이유가 이 비율이다.

아래 시뮬레이션에서 $M$과 $\delta$를 직접 조작해보자.

$\delta$를 1에서 $M^2$로 내리면 빨간 펄스가 둘로 갈라지지 않고 거의 제자리에 머문다. 아래 그래프에서 두 고유값이 $u\pm c$ 점선에서 $u$ 근처로 내려오는 것이 그 이유다.

이 식을 실제 시간으로 적분하면 음파가 음속이 아니라 유속 수준으로 움직인다. 슬라이드가 "비정상에서는 물리적으로 맞지 않다"고 쓴 것이 이것이다.

이중 시간: Γ를 가짜 시간 쪽으로 옮긴다#

해법은 $\Gamma$를 실제 시간 미분에서 떼어 가짜 시간 미분에 붙이는 것이다. 각 실제 시간 스텝마다 내부 반복을 돌린다. 내부 반복은 점-음해(point-implicit) 갱신으로 실제 시간 항만 음해로 처리한다.

$R$은 위 식의 BDF2 잔차 전체, $\Delta\tau$는 예조건 파동 속도로 정한 가짜 시간 간격이다. 공간 플럭스의 수치 소산은 $\Gamma^{-1}|\Gamma A|$를 쓴다. 예조건 코드는 소산 행렬도 함께 바꾸기 때문이다. 수렴 판정은 잔차의 최대 노름이 처음의 $10^{-6}$배로 떨어질 때다. 판정 기준이 반복 수를 어떻게 속이는지는 Jacobi 증분 판정 글에서 다뤘다.

코드: 저마하 음향 펄스, 세 가지 Γ#

$M=0.05$, 격자 100칸의 주기 영역에 가우스 압력 펄스를 둔다. 정답은 $u_0\pm c$로 갈라지는 두 음파다. [A]는 $\Gamma$를 실제 시간에 넣은 경우다. [B]는 BDF2 이중 시간에서 $\Gamma$ 없음, $\delta=M^2$, 그리고 $\Delta t$로 제한한 $U_r$의 세 경우다. 음향 CFL(음속 기준 쿠랑 수)을 1, 4, 20으로 바꾼다.

import numpy as np
 
N, M, T = 100, 0.05, 0.3          # cells, Mach number, final time (c = rho = 1)
dx = 1.0 / N
x = (np.arange(N) + 0.5) * dx
 
def pulse(s):                      # periodic Gaussian pressure pulse
    return np.exp(-(((s - 0.5 + 0.5) % 1 - 0.5) / 0.08) ** 2)
 
def exact_p(t):                    # splits into two acoustic waves u0 +/- c
    return 0.5 * (pulse(x - (M + 1) * t) + pulse(x - (M - 1) * t))
 
def precond_mats(delta):
    A = np.array([[M, 1.0], [1.0, M]])          # linearized Euler, q = (p, u)
    G = np.diag([delta, 1.0])                   # Gamma: scales dp/dt only
    lam, R = np.linalg.eig(G @ A)
    D = np.linalg.inv(G) @ (R @ np.diag(abs(lam)) @ np.linalg.inv(R)).real
    return A, G, D, abs(lam).max()              # D = Gamma^-1 |Gamma A|
 
def flux_res(q, A, D):             # upwind finite-volume residual dF/dx
    qr = np.roll(q, -1, 1)
    F = A @ (q + qr) / 2 - D @ (qr - q) / 2
    return (F - np.roll(F, 1, 1)) / dx
 
def real_time_precond(delta):      # Gamma put directly into physical time (RK2)
    A, G, D, lam = precond_mats(delta)
    ns = int(np.ceil(T / (0.5 * dx / lam))); dt = T / ns
    q = np.vstack([pulse(x), 0 * x])
    for _ in range(ns):
        k1 = -G @ flux_res(q, A, D)
        k2 = -G @ flux_res(q + dt * k1, A, D)
        q = q + dt / 2 * (k1 + k2)
    return ns, np.abs(q[0] - exact_p(T)).max(), x[np.argmax(q[0])]
 
def dual_time_run(cfl, delta, tol=1e-6):
    A, G, D, lam = precond_mats(delta)
    ns = int(round(T / (cfl * dx / (1 + M)))); dt = T / ns
    dtau = 0.9 * dx / lam                       # pseudo-time step from Gamma A
    g = np.array([[delta], [1.0]])
    q = np.vstack([pulse(x), 0 * x]); qm = q.copy(); its = 0
    for n in range(ns):
        qn = q.copy(); a = 1.5 / dt if n else 1.0 / dt
        for k in range(20000):                  # inner pseudo-time loop
            tt = (3 * q - 4 * qn + qm) / (2 * dt) if n else (q - qn) / dt
            r = tt + flux_res(q, A, D)          # BDF2 unsteady residual
            rn = np.abs(r).max()
            if k == 0: r0 = rn
            if rn < tol * r0: break
            q = q - r / (1 / (g * dtau) + a)    # point-implicit in real-time term
        its += k; qm = qn
    return its / ns, np.abs(q[0] - exact_p(T)).max()
 
print(f"M={M}, N={N}, T={T}, exact peaks at x={(0.5+(M+1)*T)%1:.3f}, {(0.5+(M-1)*T)%1:.3f}")
print("[A] Gamma in physical time (explicit RK2)")
for name, d in [("no Gamma", 1.0), ("Gamma d=M^2", M * M)]:
    ns, err, xp = real_time_precond(d)
    print(f"  {name:12s} steps={ns:3d}  max|p-p_exact|={err:.3f}  peak x={xp:.3f}")
print("[B] dual time BDF2, mean inner iterations per step / max|p-p_exact|")
print("  CFL_ac   no Gamma         Gamma d=M^2      Ur=max(u,L/(pi dt)) -> d")
for cfl in [1, 4, 20]:
    dt = T / int(round(T / (cfl * dx / (1 + M))))
    d_un = min(1.0, max(M, 1.0 / (np.pi * dt)) ** 2)   # L = domain length
    row = [dual_time_run(cfl, d) for d in (1.0, M * M, d_un)]
    print(f"  {cfl:5d}" + "".join(f"   {i:6.1f} / {e:.3f}" for i, e in row) + f"   d={d_un:.2f}")

M=0.05, N=100, T=0.3, exact peaks at x=0.815, 0.215
[A] Gamma in physical time (explicit RK2)
  no Gamma     steps= 63  max|p-p_exact|=0.145  peak x=0.215
  Gamma d=M^2  steps=  5  max|p-p_exact|=0.947  peak x=0.495
[B] dual time BDF2, mean inner iterations per step / max|p-p_exact|
  CFL_ac   no Gamma         Gamma d=M^2      Ur=max(u,L/(pi dt)) -> d
      1     16.4 / 0.149    304.7 / 0.317     16.4 / 0.149   d=1.00
      4     50.1 / 0.221   1093.4 / 0.314     50.1 / 0.221   d=1.00
     20    199.5 / 0.344   2578.0 / 0.344    199.5 / 0.344   d=1.00

위 코드를 브라우저에서 직접 돌려보자.

기본값에서 두 경우가 같은 반복 예산으로 경주한다. 노란 쪽은 16.4회/스텝, 빨간 쪽은 304.7회/스텝 근처에서 끝난다. $\delta$를 1 쪽으로 올리거나 CFL_ac를 키우면 [B] 표의 두 열이 어떻게 좁혀지는지 볼 수 있다.

[A]: 다섯 스텝 만에 끝났지만 펄스는 0.495에 있다#

실시간에 $\Gamma$를 넣으면 가장 빠른 파동 속도가 0.08로 줄어든다. 그래서 63스텝이 5스텝이 된다. 그러나 정답의 두 봉우리는 0.215와 0.815에 있다. 계산된 봉우리는 0.495, 출발점 근처다. 오차 0.947은 펄스 높이와 거의 같다. 빠르게 계산한 것은 다른 방정식의 해였다.

[B]: 내부 반복 동안 압력이 거의 움직이지 않았다#

이중 시간은 이 물리 오류를 지운다. 그 대가는 반복 수였다. 원인은 점-음해 갱신의 압력 계수 $1/(\delta\Delta\tau)$에 있다. $\delta=M^2$이면 $\Delta\tau$는 약 13배 커진다. 하지만 압력이 실제로 받는 보폭 $\delta\Delta\tau$는 $0.0025\times13\approx0.032$배, 즉 31분의 1로 줄어든다.

정상 문제에서는 이것이 맞는 거래다. 느린 대류파가 영역을 빠져나가는 시간이 수렴을 좌우하고, $\Gamma$는 그 파동을 음파만큼 빠르게 만든다. 이 실험의 내용물은 음파다. 게다가 실제 시간 항 $3/(2\Delta t)$가 이미 모든 모드를 고르게 감쇠시킨다. 고칠 강성(stiffness, 고유값 크기 차이로 생기는 수렴 지연)이 없는데 압력만 31배 느리게 만들었다. 그 결과가 304.7회다. CFL_ac를 20으로 키워도 2578회 대 199.5회로 격차가 남는다.

Ur에 Δt를 넣으면 Γ는 스스로 꺼진다#

Weiss와 Smith는 비정상용 기준 속도를 $U_r=\max(u, L/(\pi\Delta t))$로 제안했다. $L$은 영역 길이, $\Delta t$는 실제 시간 간격이다. 시간 간격이 짧아 음파가 해상될수록 $U_r$이 음속을 넘고 $\delta=\min(1,U_r^2/c^2)$는 1이 된다. 표의 마지막 열은 세 CFL 모두에서 $\delta=1.00$, 반복 수는 $\Gamma$ 없음과 같다.

즉 제대로 만든 비정상 예조건은 음향이 중요한 문제에서 꺼진다. 켜져 있어야 하는 것은 $\Delta t$가 $L/u$ 수준으로 커서 음파를 아예 해상하지 않는 문제다. 그때는 $\Delta t$가 음파 주기의 수십 배이고, 내부 반복이 다시 정상 문제처럼 대류파에 묶인다. 이 두 영역의 차이는 저마하 점근 보존 글이 다룬 "음파를 풀 것인가 걸러낼 것인가"와 같은 질문이다.

답이 0.149와 0.317로 갈린 이유#

이중 시간으로 수렴했는데 $\delta=M^2$의 오차는 0.317, $\Gamma$ 없음은 0.149다. 시간 미분의 $\Gamma$는 수렴하면 사라진다. 소산 행렬 $\Gamma^{-1}|\Gamma A|$는 사라지지 않는다. $\delta$가 작으면 이 행렬의 압력 성분이 커져 음파를 더 깎는다.

저마하 정상 유동에서는 이 소산이 오히려 정확도를 살린다. 음파가 주인공인 비정상 계산에서는 손해다. 이중 시간이 지켜주는 것은 시간 미분뿐이다. 플럭스에 들어간 $\Gamma$는 별도로 판단해야 한다.

서브루틴 다섯 개의 순서#

슬라이드 8장은 최종 형태(음해 + 비정상 + 이중 시간 + 예조건)의 호출 순서를 적는다. 위 Python 내부 루프와 짝지으면 이렇다.

1. RHS.f90 → cflux.f90, vflux.f90: 대류·점성 플럭스. 코드의 flux_res.
2. hsource.f90: 실제 시간 잔차. 코드의 tt(BDF2 항).
3. hsource.f90: 소스 항. 이 실험에는 없다.
4. 1~3을 합친 RHS. 코드의 r.
5. LHS.f90: 좌변 역행렬. 원래 코드는 LU-SGS나 Gauss–Seidel로 풀고, 위 코드는 대각 점-음해만 쓴다.

슬라이드는 5단계의 음해 처리가 "잔차 평활(residual smoothing)과 같은 개념"이라고 적는다. 정확하다. 내부 반복의 LHS는 수렴한 답을 바꾸지 않는다. 바꾸는 것은 수렴까지의 경로와 반복 수뿐이다. $\Gamma$도 시간 미분에 있는 한 같은 부류다.

비정상 코드에서 Γ를 켜기 전에 계산할 비율 하나#

먼저 $c\Delta t/L$을 계산한다. 이 값이 1보다 크게 작으면 음파를 해상하는 계산이다. 이때는 $U_r$에 $L/(\pi\Delta t)$ 하한을 넣어 $\Gamma$가 꺼지게 둔다. 이 값이 $1/M$ 수준이면 음파를 걸러내는 계산이고, 그때만 $\Gamma$가 내부 반복을 줄인다.

다음으로 소산 행렬을 따로 본다. 시간 미분의 $\Gamma$와 플럭스의 $\Gamma$는 서로 다른 스위치로 두는 편이 안전하다. 마지막으로 한 스텝의 내부 반복 수를 $\Gamma$ 켜고 끄고 한 번씩 찍어 본다. 이 실험처럼 19배 차이가 나면 슬라이드의 한 줄을 다시 읽을 때다.