같은 코드가 0.79차와 2.00차를 동시에 냈다 — 검증 케이스가 실제로 재는 것

불연속이 든 해에서 잰 수렴차수는 스킴의 차수가 아니다. 차수는 매끄러운 해에서만 재고, Sod 충격관은 충격파 위치와 단조성을 본다.

2차 스킴이 0.79차를 내놓았다#

압축성 코드 검증 목록에는 세 가지가 적혀 있었다. Sod 충격관, 등엔트로피 와동 이송, 충격파–와동 상호작용. 순서대로 Sod부터 걸었다.

격자를 100에서 800까지 네 단계로 조이고 밀도의 $L_1$ 오차를 쟀다. 수렴차수가 0.79로 나왔다. MUSCL 재구성에 HLLC 플럭스, 시간은 2단 SSP Runge–Kutta. 설계 차수는 2다.

이틀 동안 재구성과 경계조건을 뒤졌다. 버그는 없었다. 0.79는 정상이다.

이 글은 왜 정상인지를 숫자로 보여준다. 같은 코드가 다른 테스트에서 2.001을 내놓는 것까지 확인하고, 세 검증 케이스가 각각 무엇을 잴 수 있는지 정리한다.

수렴차수는 스킴이 아니라 해의 성질이다#

차수 증명은 테일러 전개에서 나온다. 전개를 쓰려면 해가 그만큼 미분 가능해야 한다. 불연속에서 $\rho'$ 는 존재하지 않는다. 전제가 깨지면 결론도 따라 깨진다.

그래서 "2차 스킴"은 "매끄러운 해에서 2차"의 줄임말이다. 괄호 안을 빼면 문장이 틀린다. N을 올리면 번진 폭만 줄고, show Linf를 누르면 같은 해가 다른 기울기로 읽힌다.

N = 100 smear = 3 cells t = 0.200 
at N = 100: L1 0.00e+0  ·  Linf 0.00e+0  ·  max error sits at x = 0.000

번진 구간의 폭은 격자와 함께 좁아진다. 불연속을 건너는 높이는 그대로다. $L_1$ 기울기와 $L_\infty$ 기울기가 한 화면에서 갈라진다. 수치해는 정확해를 격자 폭에 비례해 번지게 한 모형이다.

번진 불연속이 $L_1$ 에 남기는 면적#

충격포착 스킴은 불연속을 몇 개 셀에 걸쳐 번지게 푼다. 번지는 셀 수 $\alpha$ 는 격자를 조여도 거의 일정하다. 물리적 번짐 폭은 $\alpha h$ 다.

\|e\|_{1} = \int_0^1 |\rho_h - \rho|\,\mathrm{d}x \approx |\Delta\rho|\,\alpha h

$\Delta\rho$ 는 불연속 양쪽의 밀도 차, $\alpha$ 는 번짐이 걸친 셀 수, $h$ 는 격자 간격이다. 높이는 고정이고 폭만 줄어든다. 면적은 $h$ 에 비례한다. 1차다.

Sod 해에는 불연속이 둘 있다. 충격파와 접촉면이다. 둘은 서로 다르게 번진다.

충격파에는 스스로 날을 세우는 압축 기구가 있다. 번짐 폭이 셀 몇 개로 묶여 $O(h)$ 를 낸다. 접촉면에는 그 기구가 없다. 수치 확산이 시간에 걸쳐 쌓여 폭이 $h$ 보다 느리게 줄어든다. 모노톤 스킴에서 알려진 추정은 $O(h^{2/3})$ 수준이다.

측정값 0.79는 그 둘의 혼합이다. 2/3과 1 사이에 앉는다.

$L_\infty$ 가 0.083에서 멈춰 선 이유#

같은 계산을 최대 노름으로 다시 읽으면 숫자가 아예 움직이지 않는다.

\|e\|_{\infty} = \max_i |\rho_{h,i} - \rho(x_i)| \approx \beta\,|\Delta\rho|

번진 구간 안에는 정확해와 수치해가 $|\Delta\rho|$ 의 일정 비율만큼 벌어진 셀이 반드시 하나 있다. 그 비율 $\beta$ 는 격자와 무관하다.

측정값은 0.0678 → 0.0830 → 0.0834 → 0.0857이었다. 출력에 적힌 값은 $6.7837\times 10^{-2}$ , $8.3037\times 10^{-2}$ , $8.3428\times 10^{-2}$ , $8.5667\times 10^{-2}$ 이고, 차수는 −0.292, −0.007, −0.038이다. N=200에서 이미 0.083 근처에 붙어, 800까지 조여도 0.0857을 넘기지 않는다. 격자를 조여도 폭이 셀 몇 개에 고정되는 일은 충격파를 일부러 번지게 둔 인공점성과 같은 자리다. 거기서는 점성 항이 폭을 고정했고, 여기서는 노름이 높이를 고정한다.

$L_1$ 쪽 숫자는 다르게 움직인다. N=100에서 $6.5089\times 10^{-3}$ , 200에서 $3.7579\times 10^{-3}$ 에 차수 0.792, 400에서 $2.1461\times 10^{-3}$ 에 0.808, 800에서 $1.2420\times 10^{-3}$ 에 0.789다. 0.79 근처에서 흔들리되 1을 넘지 않는다. 매끄러운 파에서 제한자를 끈 $L_1$ 은 같은 네 격자에서 $2.0251\times 10^{-3}$ , $5.0589\times 10^{-4}$ (2.001), $1.2630\times 10^{-4}$ (2.002), $3.1562\times 10^{-5}$ (2.001)이다. minmod를 켜면 $5.6574\times 10^{-3}$ , $2.0125\times 10^{-3}$ (1.491), $5.7467\times 10^{-4}$ (1.808), $1.5736\times 10^{-4}$ (1.869)이다. N=100에서 두 오차의 비는 이미 말한 2.8배다. N=800에서도 끈 쪽 $3.1562\times 10^{-5}$ 와 켠 쪽 $1.5736\times 10^{-4}$ 가 나란히 남는다. 차수 1.869는 2로 가는 도중이지, 설계 차수가 바뀐 것은 아니다.

같은 코드, 같은 스킴, 매끄러운 해#

압축성 오일러 방정식에도 쓰기 쉬운 매끈한 정확해가 있다. 균일한 압력과 균일한 속도 위에 밀도만 흔든 장이다. 접촉파 하나뿐이라 그 모양이 그대로 실려 간다.

\rho(x,t) = \rho_0(x - u_0 t), \qquad u \equiv u_0, \qquad p \equiv p_0

$\rho_0$ 는 초기 밀도 분포, $u_0$ 는 균일 속도, $p_0$ 는 균일 압력이다. 주기 경계에서 한 바퀴 돌린 뒤 초기장과 비교하면 된다. 2차원 등엔트로피 와동이 맡는 역할과 같다.

가우스 봉우리를 $u_0 = 1$ 로 한 바퀴 이송했다. 제한자를 끈 중앙 기울기에서 2.001, 2.002, 2.001이 나왔다. 설계 차수가 그대로 보인다. 봉우리만 사각파 쪽으로 세우면 스킴 줄은 그대로인데 차수가 1 쪽으로 내려앉는다. 그 슬라이더는 출력 표 바로 뒤에 있다.

minmod이 2.001을 1.491로 깎는 자리#

실전 코드에는 제한자가 켜져 있다. minmod은 두 차분의 부호가 다르면 기울기를 0으로 만든다.

\sigma_i = \operatorname{minmod}(\rho_i - \rho_{i-1},\; \rho_{i+1} - \rho_i)

$\sigma_i$ 는 셀 $i$ 의 재구성 기울기다. 봉우리 꼭대기에서는 왼쪽 차분이 양수이고 오른쪽 차분이 음수다. 조건에 걸려 $\sigma_i = 0$ 이 된다. 그 셀만 1차 재구성으로 떨어진다.

측정 차수는 1.491 → 1.808 → 1.869였다. 2로 돌아가는 중이다. 차수가 바뀐 것이 아니라 극값 한 점의 1차 오차가 아직 남아 있는 것이다. 같은 격자에서 제한자를 끄면 오차가 2.8배 작다.

실무 함의는 하나다. 차수 검증은 제한자를 끄거나 극값 보존 제한자로 하고, 생산 계산은 제한자를 켠다. 두 조건을 한 표에 섞어 적으면 어느 쪽도 못 믿는다.

코드 한 벌로 네 숫자를 같이 뽑는다#

같은 솔버에 Sod와 매끄러운 파를 차례로 걸었다. 기준해는 리만 문제의 정확해와 평행이동한 가우스다.

import numpy as np
G = 1.4
 
def riemann_star_pressure(rl, ul, pl, rr, ur, pr):
    al, ar = np.sqrt(G*pl/rl), np.sqrt(G*pr/rr)
    def fk(p, r, u, a, pk):
        if p > pk:
            A, B = 2.0/((G+1)*r), (G-1)/(G+1)*pk
            f = (p-pk)*np.sqrt(A/(p+B))
            df = np.sqrt(A/(B+p))*(1.0-0.5*(p-pk)/(B+p))
        else:
            f = 2*a/(G-1)*((p/pk)**((G-1)/(2*G))-1.0)
            df = 1.0/(r*a)*(p/pk)**(-(G+1)/(2*G))
        return f, df
    p = 0.5*(pl+pr)
    for _ in range(80):
        f1, d1 = fk(p, rl, ul, al, pl)
        f2, d2 = fk(p, rr, ur, ar, pr)
        dp = (f1+f2+(ur-ul))/(d1+d2)
        p = max(p-dp, 1e-10)
        if abs(dp) < 1e-13: break
    f1, _ = fk(p, rl, ul, al, pl)
    f2, _ = fk(p, rr, ur, ar, pr)
    return p, 0.5*(ul+ur)+0.5*(f2-f1)
 
def sample_exact_riemann(x, t, rl, ul, pl, rr, ur, pr):
    ps, us = riemann_star_pressure(rl, ul, pl, rr, ur, pr)
    al, ar = np.sqrt(G*pl/rl), np.sqrt(G*pr/rr)
    s = x/t
    rho = np.empty_like(s)
    for i, si in enumerate(s):
        if si <= us:                                  # 접촉면 왼쪽
            if ps > pl:                               # 왼쪽 충격파
                rsl = rl*((ps/pl+(G-1)/(G+1))/((G-1)/(G+1)*ps/pl+1))
                ss = ul-al*np.sqrt((G+1)/(2*G)*ps/pl+(G-1)/(2*G))
                rho[i] = rl if si < ss else rsl
            else:                                     # 왼쪽 팽창파
                asl = al*(ps/pl)**((G-1)/(2*G))
                if si < ul-al: rho[i] = rl
                elif si > us-asl: rho[i] = rl*(ps/pl)**(1/G)
                else:
                    c = 2/(G+1)+(G-1)/((G+1)*al)*(ul-si)
                    rho[i] = rl*c**(2/(G-1))
        else:
            if ps > pr:
                rsr = rr*((ps/pr+(G-1)/(G+1))/((G-1)/(G+1)*ps/pr+1))
                ss = ur+ar*np.sqrt((G+1)/(2*G)*ps/pr+(G-1)/(2*G))
                rho[i] = rr if si > ss else rsr
            else:
                asr = ar*(ps/pr)**((G-1)/(2*G))
                if si > ur+ar: rho[i] = rr
                elif si < us+asr: rho[i] = rr*(ps/pr)**(1/G)
                else:
                    c = 2/(G+1)-(G-1)/((G+1)*ar)*(ur-si)
                    rho[i] = rr*c**(2/(G-1))
    return rho
 
def to_prim(U):
    r = U[0]; u = U[1]/r; p = (G-1)*(U[2]-0.5*r*u*u)
    return r, u, p
 
def hllc_flux(rL, uL, pL, rR, uR, pR):
    aL, aR = np.sqrt(G*pL/rL), np.sqrt(G*pR/rR)
    EL = pL/(G-1)+0.5*rL*uL**2; ER = pR/(G-1)+0.5*rR*uR**2
    sL = np.minimum(uL-aL, uR-aR); sR = np.maximum(uL+aL, uR+aR)
    sM = (pR-pL+rL*uL*(sL-uL)-rR*uR*(sR-uR))/(rL*(sL-uL)-rR*(sR-uR))
    FL = np.array([rL*uL, rL*uL**2+pL, uL*(EL+pL)])
    FR = np.array([rR*uR, rR*uR**2+pR, uR*(ER+pR)])
    UL = np.array([rL, rL*uL, EL]); UR = np.array([rR, rR*uR, ER])
    UsL = rL*(sL-uL)/(sL-sM)*np.array([np.ones_like(rL), sM, EL/rL+(sM-uL)*(sM-pL/(rL*(sL-uL)))])
    UsR = rR*(sR-uR)/(sR-sM)*np.array([np.ones_like(rR), sM, ER/rR+(sM-uR)*(sM-pR/(rR*(sR-uR)))])
    F = np.where(sL >= 0, FL, np.where(sM >= 0, FL+sL*(UsL-UL),
        np.where(sR >= 0, FR+sR*(UsR-UR), FR)))
    return F
 
def minmod_slope(q, limited):
    dL = q[1:-1]-q[:-2]; dR = q[2:]-q[1:-1]
    if not limited:
        return 0.5*(dL+dR)
    return np.where(dL*dR <= 0, 0.0, np.where(np.abs(dL) < np.abs(dR), dL, dR))
 
def residual(U, h, periodic, limited):
    V = np.vstack(to_prim(U))                      # 3 x (N+4)
    sl = np.zeros_like(V)
    for k in range(3):
        sl[k, 1:-1] = minmod_slope(V[k], limited)
    VL = V[:, 1:-1]+0.5*sl[:, 1:-1]                # 셀 i의 오른쪽 면
    VR = V[:, 2:]-0.5*sl[:, 2:]                    # 셀 i+1의 왼쪽 면
    F = hllc_flux(VL[0], VL[1], VL[2], VR[0], VR[1], VR[2])
    R = np.zeros_like(U)
    R[:, 2:-2] = -(F[:, 1:-1]-F[:, :-2])/h
    return R
 
def fill_ghost(U, periodic):
    if periodic:
        U[:, :2] = U[:, -4:-2]; U[:, -2:] = U[:, 2:4]
    else:
        U[:, 0] = U[:, 2]; U[:, 1] = U[:, 2]
        U[:, -1] = U[:, -3]; U[:, -2] = U[:, -3]
 
def evolve(U, h, tend, periodic, limited, cfl=0.4):
    t = 0.0
    while t < tend-1e-14:
        fill_ghost(U, periodic)
        r, u, p = to_prim(U[:, 2:-2])
        dt = cfl*h/np.max(np.abs(u)+np.sqrt(G*p/r))
        dt = min(dt, tend-t)
        U1 = U+dt*residual(U, h, periodic, limited); fill_ghost(U1, periodic)
        U = 0.5*(U+U1+dt*residual(U1, h, periodic, limited))
        t += dt
    fill_ghost(U, periodic)
    return U
 
def order_table(ns, errs):
    return [None]+[np.log(errs[i-1]/errs[i])/np.log(ns[i]/ns[i-1]) for i in range(1, len(ns))]
 
def run_sod(ns, tend=0.2):
    e1, ei = [], []
    for N in ns:
        h = 1.0/N
        x = (np.arange(N)+0.5)*h
        r = np.where(x < 0.5, 1.0, 0.125); u = np.zeros(N)
        p = np.where(x < 0.5, 1.0, 0.1)
        U = np.zeros((3, N+4))
        U[:, 2:-2] = np.vstack([r, r*u, p/(G-1)+0.5*r*u*u])
        U = evolve(U, h, tend, False, True)
        re = sample_exact_riemann(x-0.5, tend, 1.0, 0.0, 1.0, 0.125, 0.0, 0.1)
        d = np.abs(U[0, 2:-2]-re)
        e1.append(np.sum(d)*h); ei.append(np.max(d))
    return e1, ei
 
def run_wave(ns, limited, tend=1.0):
    e1 = []
    for N in ns:
        h = 1.0/N
        x = (np.arange(N)+0.5)*h
        prof = lambda z: 1.0+0.5*np.exp(-60.0*(np.mod(z-0.5, 1.0)-0.5)**2)
        r = prof(x); u = np.ones(N); p = np.ones(N)
        U = np.zeros((3, N+4))
        U[:, 2:-2] = np.vstack([r, r*u, p/(G-1)+0.5*r*u*u])
        U = evolve(U, h, tend, True, limited)
        e1.append(np.sum(np.abs(U[0, 2:-2]-prof(x-tend)))*h)
    return e1
 
 
def rankine_hugoniot(m1):
    """정지 충격파의 랭킨-위고니오 비 (충격파-와동 설정)"""
    rr = ((G+1)*m1**2)/((G-1)*m1**2+2)
    pr = (2*G*m1**2-(G-1))/(G+1)
    m2 = np.sqrt((1+0.5*(G-1)*m1**2)/(G*m1**2-0.5*(G-1)))
    return rr, pr, m2
 
NS = [100, 200, 400, 800]
s1, si = run_sod(NS)
wl, wc = run_wave(NS, True), run_wave(NS, False)
o_s1, o_si, o_wl, o_wc = (order_table(NS, e) for e in (s1, si, wl, wc))
fmt = lambda o: "    -  " if o is None else f"{o:7.3f}"
print("Sod shock tube (two discontinuities), t = 0.2")
print("    N      L1 error   order     Linf error  order")
for i, N in enumerate(NS):
    print(f" {N:4d}   {s1[i]:.4e} {fmt(o_s1[i])}   {si[i]:.4e} {fmt(o_si[i])}")
print()
print("smooth density wave (exact = translated Gaussian), t = 1.0")
print("    N      L1 minmod  order     L1 unlim.   order")
for i, N in enumerate(NS):
    print(f" {N:4d}   {wl[i]:.4e} {fmt(o_wl[i])}   {wc[i]:.4e} {fmt(o_wc[i])}")
ps, us = riemann_star_pressure(1.0, 0.0, 1.0, 0.125, 0.0, 0.1)
print(f"\nSod star state   p* = {ps:.6f}  u* = {us:.6f}")
rr, pr, m2 = rankine_hugoniot(1.2)
print(f"RH at M1 = 1.2   rho2/rho1 = {rr:.4f}  p2/p1 = {pr:.4f}  M2 = {m2:.4f}")
amp = (G-1)*25.0/(8*G*np.pi**2)
for r in (2.5, 5.0):
    print(f"vortex tail r = {r:.1f}  exp(1-r^2) = {np.exp(1-r*r):.3e}  dT(beta=5) = {amp*np.exp(1-r*r):.3e}")

출력:

Sod shock tube (two discontinuities), t = 0.2
    N      L1 error   order     Linf error  order
  100   6.5089e-03     -     6.7837e-02     -  
  200   3.7579e-03   0.792   8.3037e-02  -0.292
  400   2.1461e-03   0.808   8.3428e-02  -0.007
  800   1.2420e-03   0.789   8.5667e-02  -0.038
 
smooth density wave (exact = translated Gaussian), t = 1.0
    N      L1 minmod  order     L1 unlim.   order
  100   5.6574e-03     -     2.0251e-03     -  
  200   2.0125e-03   1.491   5.0589e-04   2.001
  400   5.7467e-04   1.808   1.2630e-04   2.002
  800   1.5736e-04   1.869   3.1562e-05   2.001
 
Sod star state   p* = 0.303130  u* = 0.927453
RH at M1 = 1.2   rho2/rho1 = 1.3416  p2/p1 = 1.5133  M2 = 0.8422
vortex tail r = 2.5  exp(1-r^2) = 5.248e-03  dT(beta=5) = 4.747e-04
vortex tail r = 5.0  exp(1-r^2) = 3.775e-11  dT(beta=5) = 3.415e-12

sharpness $m$ 을 올리면 봉우리가 가우스에서 사각파로 바뀌고, 측정 차수가 2 근처에서 1 쪽으로 내려앉는다.

sharpness m = 2 N = 200 smear = 3 cells 
m = 2 gives a measured L1 order of 2.00 — the scheme never changed.

제한자를 끈 줄의 2.001이 $m$ 과 함께 1 쪽으로 미끄러진다. show Linf로 바꾸면 같은 계산이 2에서 0까지 떨어진다. 스킴 줄은 그대로다.

$p^{*} = 0.303130$ , $u^{*} = 0.927453$ 은 교과서 Sod 값이다. 충격파 속도는 1.7522가 된다. 이 숫자가 맞는지가 Sod 충격관이 실제로 보는 것이다.

등엔트로피 와동이 맡는 자리#

2차원에서 차수를 재는 표준 케이스가 등엔트로피 와동 이송이다. $[0,10]^2$ 주기 영역에 균일 유동을 깔고, 그 위에 와동 교란을 더한다.

\delta u = -\frac{\beta}{2\pi}(y - y_c)\,e^{(1-r^2)/2}, \qquad \delta v = \frac{\beta}{2\pi}(x - x_c)\,e^{(1-r^2)/2}

$\beta$ 는 와동 강도, $r$ 은 와동 중심까지의 거리, $(x_c, y_c)$ 는 중심 좌표다. 온도는 등엔트로피 관계로 묶인다.

T = 1 - \frac{(\gamma-1)\beta^2}{8\gamma\pi^2}e^{1-r^2}, \qquad \rho = T^{1/(\gamma-1)}, \qquad p = \rho^{\gamma}

이 장은 자유류 속도로 평행이동하는 정확해다. 매끄럽고 불연속이 없다. 그래서 차수를 재도 된다.

함정은 꼬리다. 교란은 $e^{1-r^2}$ 로 줄지만 0이 되지는 않는다. 영역 경계에서 잘라 버리면 "정확해"가 더 이상 정확하지 않다.

$[0,10]^2$ 에서 중심까지 거리는 5다. $e^{1-r^2}$ 가 $3.8 \times 10^{-11}$ , $\beta = 5$ 일 때 온도 교란이 $3.4 \times 10^{-12}$ 다. 2차 스킴에는 충분히 안전하다. 영역을 $[0,5]^2$ 로 줄이면 $r = 2.5$ 에서 각각 $5.2 \times 10^{-3}$ 과 $4.7 \times 10^{-4}$ 가 된다. 고차 스킴은 격자를 조이다 이 바닥에 먼저 닿는다. 차수가 거기서 무너진다.

충격파–와동은 차수를 재지 않는다#

세 번째 케이스는 성격이 다르다. 왼쪽에 $M = 1.2$ 정지 충격파를 세우고 오른쪽 상태를 랭킨–위고니오 관계로 정한 뒤, 그 충격파로 등엔트로피 와동을 통과시킨다.

이 문제에는 닫힌 해가 없다. 차수를 뽑을 수 없다.

확인할 숫자는 충격파 전후의 비다. $M_1 = 1.2$ , $\gamma = 1.4$ 에서 $\rho_2/\rho_1 = 1.3416$ , $p_2/p_1 = 1.5133$ , $M_2 = 0.8422$ 다. 코드가 이 숫자를 내놓지 못하면 충격파 뒤의 그림은 전부 의미가 없다. 주기 방향의 질량과 운동량 적분은 기계 정밀도 수준에 붙어 있어야 한다. 음향파 진폭이나 와도 적분 같은 범함수는 세 격자에서 재고, 가장 고운 격자를 기준해로 삼는다. 그건 차수가 아니라 격자 수렴성이다.

리포트에 숫자 하나를 적을 때 고르는 것#

"2차 정확도를 확인했다"는 한 문장에는 선택이 셋 숨어 있다. 어떤 해, 어떤 노름, 제한자를 켰는지다.

케이스	정확해	매끄러움	잴 수 있는 것
Sod 충격관	있음 (리만)	불연속	충격파 위치·속도, 단조성, 양의성
등엔트로피 와동	있음 (평행이동)	매끄러움	설계 수렴차수
충격파–와동	없음	혼합	보존량, RH 정합, 격자 수렴성

Sod에서 0.79가 나왔다면 고칠 것은 스킴이 아니라 리포트 문장이다. 와동에서 2.00이 나왔다면 제한자 상태를 함께 적어야 한다. 매끄러운 해석해 위에서 차수를 재는 일은 포아즈이유 관의 해석해와 맞대는 검증이 하는 일과 같다. 어디서 잴지 고르는 일이 그 절반이다.