同じコードが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$ は一様圧力です。周期境界で一周させて初期場と比べるだけです。二次元の等エントロピー渦が担う役割と同じです。

ガウス型の山を $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衝撃波管が本当に見ているところです。

等エントロピー渦が担う役割

二次元で次数を測る標準ケースが等エントロピー渦の移流です。 $[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が出たなら、制限子の設定を並べて書くべきです。滑らかな厳密解の上で次数を測る作業は、ポアズイユ管の厳密解と突き合わせる検証と同じ仕事です。どこで測るかを選ぶことがその半分です。