積分点を一つ減らしたら解が発散した — DGの求積下限とテイラー基底

DGの体積項は次数 $2p-1$ の多項式です。ガウス $n$ 点は $2n-1$ まで厳密なので下限は $n = p$ であり、それを下回ると精度が落ちるのではなくスキームそのものが壊れます。

積分点を一つ抜いたら答えが丸ごと消えた

不連続ガレルキン（DG, Discontinuous Galerkin — セルごとに独立した多項式を置き、面のフラックスで繋ぐ高次手法）コードのセル積分を、ガウス3点から2点へ減らしたことがあります。計算上はセルあたりの積分コストが3分の1減ります。走らせてみると、L2誤差は小数点以下13桁まで同じでした。欲が出て1点まで下げました。今度は精度が一次落ちたのではなく、解が一周もせずに発散しました。

境界線は格子幅にもCFL数にもありませんでした。被積分関数の多項式次数にあったのです。この記事では、その境界線がどこにあり、なぜそこなのか、そして任意格子でその線を守るには基底関数をどう取るべきかを扱います。根拠は1次元DG-P2ソルバと質量行列の条件数計算の二つです。

DG order p を上げると緑のまま残るガウス点の下限も上がり、Gauss points n をその下へ下げるとバッジが赤へ変わります。

DG order p2Gauss points n2test mode i2

u_h shape0.60true 0.0000 · quad 0.0000

Set p = 2 and drag Gauss points from 3 down to 2: the badge stays green and the error bar stays empty, because the integrand only has degree 3. Drag to 1 and it turns red — no shape of u_h will bring it back. The minimum is n = p, and it is a property of the integrand’s degree, not of how fine the mesh is.

$n \ge p$ の間はバッジが緑で、誤差バーは空のままです。u_h shape を揺らしても同じです。 $n$ をあと一段下げた瞬間、バーに誤差が入ります。

Q1. DGは有限要素なのか、有限体積なのか#

両方です。セル一つだけを見れば有限要素、セル境界だけを見れば有限体積です。

保存形方程式に試験関数 $\phi_i$ を掛け、セル $\Omega_e$ 上で積分して部分積分すると次の形になります。

\int_{\Omega_e} \frac{\partial U_h}{\partial t}\,\phi_i \, d\Omega - \int_{\Omega_e} \mathbf{F}(U_h)\cdot\nabla\phi_i \, d\Omega + \oint_{\partial\Omega_e} \phi_i\, \hat{\mathbf{F}}(U_h^-, U_h^+)\cdot\mathbf{n}\, dS = 0

$U_h$ はセル内部の近似解、 $\mathbf{F}$ は対流フラックス、 $\hat{\mathbf{F}}$ は両側のトレース $U_h^-, U_h^+$ から作る数値フラックス、 $\mathbf{n}$ は面法線です。粘性フラックスとソース項はそれぞれ項が一つずつ増えますが、構造は同じです。

重要なのは、式が二つの塊に分かれる点です。体積積分はセル内部で閉じ、隣と通信するのは面積分だけです。そこに入るのはリーマンソルバが作った単一値のフラックスです。有限体積法がセル平均一つで行うことを、DGは多項式係数の複数個で行っているにすぎません。セル内に多項式を置き、面だけで数値フラックスを繋ぐ骨格は非構造格子の不連続ガレルキン離散化と同じです。フラックス差分を失えば、衝撃波速度も一緒に外れます。

近似解を基底関数の線形結合で書けば

U_h(\mathbf{x}, t) = \sum_{k=1}^{K} U_k(t)\, b_k(\mathbf{x})

となり、時間項は質量行列 $M_{ik} = \int_{\Omega_e} b_i b_k \, d\Omega$ になります。セルごとに $K \times K$ の小さな行列が一つ。隣と混ざらないので、セル単位で逆行列を先に求めて保存できます。DGが並列化に強い理由の大きな部分がこれです。

Q2. 積分は何次まで厳密であるべきか#

体積積分の被積分関数の次数を数えれば答えが出ます。

$p$ 次の多項式空間を使うとします。 $U_h$ は次数 $p$ 、試験関数 $\phi_i$ も最大 $p$ 次なので $\nabla\phi_i$ は $p-1$ 次です。線形フラックスなら積の次数はこうなります。

\deg\left(\mathbf{F}(U_h)\cdot\nabla\phi_i\right) = p + (p-1) = 2p - 1

ガウス・ルジャンドル $n$ 点公式は次数 $2n-1$ まで厳密です。二つを繋ぐと下限が現れます。

2n - 1 \ge 2p - 1 \quad \Longrightarrow \quad n \ge p

原資料にある「最低 $2K-1$ オーダーの積分を行わないと次数が低下する」という一文がこれです。格子をいくら細かくしてもこの不等式は動きません。多項式の次数はセル幅と無関係だからです。

注意点が二つあります。質量行列の被積分関数は $b_i b_k$ なので次数が $2p$ で、下限は $n \ge p+1$ と一段高くなります。さらにフラックスが非線形なら $\mathbf{F}(U_h)$ はそもそも多項式ではありません。CockburnとShuが体積 $2p$ 次・面 $2p+1$ 次を推奨したのはそのためです。実務では曲面要素のヤコビアンまで掛かるので、さらに余裕を取ります。

Q3. 1点まで下げると何が壊れるのか#

測ってみるのが早いです。周期領域 $[0, 2\pi]$ で $u_t + u_x = 0$ をDG-P2で解きます。基底はルジャンドル、時間進行はSSP-RK3、面フラックスは風上です。質量行列は解析的に入れて、体積積分の点数だけを変数として残します。

from math import pi, sin, exp, log, sqrt, ceil
 
GAUSS = {                                   # Gauss-Legendre on [-1,1]: exact to degree 2n-1
    1: ([0.0], [2.0]),
    2: ([-0.5773502691896257, 0.5773502691896257], [1.0, 1.0]),
    3: ([-0.7745966692414834, 0.0, 0.7745966692414834], [5/9, 8/9, 5/9]),
    6: ([-0.9324695142031521, -0.6612093864662645, -0.2386191860831969,
          0.2386191860831969,  0.6612093864662645,  0.9324695142031521],
        [0.1713244923791704, 0.3607615730481386, 0.4679139345726910,
         0.4679139345726910, 0.3607615730481386, 0.1713244923791704]),
}
PHI  = [lambda s: 1.0, lambda s: s,   lambda s: 1.5*s*s - 0.5]   # ルジャンドルモード, p = 2
DPHI = [lambda s: 0.0, lambda s: 1.0, lambda s: 3.0*s]
K = 3
 
def dg_rhs(U, h, nq):
    """u_t + u_x = 0 のDG半離散残差。面は風上フラックス。"""
    xq, wq = GAUSS[nq]
    N = len(U)
    uR = [sum(U[j][i]*PHI[i](1.0) for i in range(K)) for j in range(N)]   # 右トレース
    R  = []
    for j in range(N):
        fR = uR[j]                       # a = 1 > 0 なので面は左セルの値を取る
        fL = uR[j-1]
        row = []
        for i in range(K):
            vol = 0.0
            for xk, wk in zip(xq, wq):
                uh = sum(U[j][m]*PHI[m](xk) for m in range(K))
                vol += wk*DPHI[i](xk)*uh
            surf = PHI[i](1.0)*fR - PHI[i](-1.0)*fL
            row.append((vol - surf)*(2*i+1)/h)                # M_ii = h/(2i+1)
        R.append(row)
    return R
 
def run_dg(N, nq, T=1.0, cfl=0.05):
    h  = 2*pi/N
    xc = [h*(j + 0.5) for j in range(N)]
    xg, wg = GAUSS[6]
    u0 = lambda x: exp(sin(x))
    U  = [[(2*i+1)/2*sum(w*PHI[i](s)*u0(xc[j] + h/2*s) for s, w in zip(xg, wg))
           for i in range(K)] for j in range(N)]
    nt = int(ceil(T/(cfl*h/5))); dt = T/nt
    for _ in range(nt):                                        # SSP-RK3
        R0 = dg_rhs(U, h, nq)
        U1 = [[U[j][i] + dt*R0[j][i] for i in range(K)] for j in range(N)]
        R1 = dg_rhs(U1, h, nq)
        U2 = [[0.75*U[j][i] + 0.25*(U1[j][i] + dt*R1[j][i]) for i in range(K)] for j in range(N)]
        R2 = dg_rhs(U2, h, nq)
        U  = [[(U[j][i] + 2*(U2[j][i] + dt*R2[j][i]))/3 for i in range(K)] for j in range(N)]
    e2 = 0.0
    for j in range(N):
        for s, w in zip(xg, wg):
            uh = sum(U[j][i]*PHI[i](s) for i in range(K))
            e2 += w*(uh - u0(xc[j] + h/2*s - T))**2*h/2
    return sqrt(e2)
 
print("nq  exact-to-deg |   N=10       N=20       N=40    | order")
for nq in (1, 2, 3):
    e = [run_dg(N, nq) for N in (10, 20, 40)]
    print(f" {nq}       {2*nq-1}      | {e[0]:.3e}  {e[1]:.3e}  {e[2]:.3e} |  {log(e[1]/e[2], 2):.2f}")

nq  exact-to-deg |   N=10       N=20       N=40    | order
 1       1      | 1.170e+01  1.302e+01  1.186e+01 |  0.13
 2       3      | 5.989e-03  7.369e-04  9.211e-05 |  3.00
 3       5      | 5.989e-03  7.369e-04  9.211e-05 |  3.00

三行を順に読んでみましょう。 $n=2$ と $n=3$ は三つの格子すべてで表示桁まで同一です。実際には13桁目の有効数字で分かれ、その差は丸め誤差です。被積分関数が次数3なので、2点公式がすでに厳密な値を返しているからです。点を増やしても得るものはありません。

$n=1$ の行は性格が違います。誤差は $10^1$ の規模で、格子を4倍細かくしても縮みません。収束次数0.13は「一次に落ちた」ではなく「収束しない」という意味です。不足求積（under-integration）は毎ステップ誤った体積項を食わせ、その誤差が時間方向に増幅されます。次のシミュレーションがその過程をそのまま見せてくれます。

Gauss points n2cells N16play rate0.50

t = 0.00 · L2 0.00e+0

Watch the face jumps first: they are tiny while the rule is consistent, and they are what the upwind flux has to reconcile. Now drag Gauss points from 3 to 2 — nothing moves, the L2 readout does not budge. Drag to 1 and the parabolas tear apart within a fraction of a revolution. Adding cells only makes it happen sooner.

まず Gauss points を3から2へ下げてもL2誤差の表示が動かないことを確認してください。次に1へ下げると、セルごとの放物線が一周もせずに裂けます。cells N を増やすとより早く壊れます。

出力の三行を次数で分けると、下限が表と一致します。P2の体積被積分は $2p-1=3$ 次です。ガウス $n$ 点の厳密次数は $2n-1$ なので、 $n=2$ は3次、 $n=3$ は5次まで正確です。3次をすでに積分できる2点公式は、3点と同じ積分値を返します。だから N=10, 20, 40 の L2 が $5.989\times 10^{-3}$ 、 $7.369\times 10^{-4}$ 、 $9.211\times 10^{-5}$ で重なります。次数は $\log_2(7.369\times 10^{-4}/9.211\times 10^{-5})$ です。比は約 8.00 で $\log_2 8 = 3$ なので、出力の 3.00 がその割り算です。

$n=1$ の厳密次数は 1 で、 $2p-1=3$ に届きません。誤差は $1.170\times 10^{1}$ 、 $1.302\times 10^{1}$ 、 $1.186\times 10^{1}$ です。N=20 から N=40 への比は $1.302/1.186 \approx 1.098$ で、 $\log_2 1.098 \approx 0.135$ 、出力の 0.13 と一致します。格子を4倍にしても誤差が 10 付近に残る理由がこの比です。13桁まで同じという話は $n=2$ と $n=3$ の丸めで、 $n=1$ とは桁そのものが違います。

条件数の列も同じ読み方です。 $h$ が 1、 $10^{-1}$ 、 $10^{-2}$ 、 $10^{-3}$ のとき、正規化前の値は $7.225\times 10^{2}$ 、 $7.200\times 10^{6}$ 、 $7.200\times 10^{10}$ 、 $7.200\times 10^{14}$ です。長さが10分の1になるたびに条件数は約 $10^{4}$ 倍です。 $7.200\times 10^{6}/7.225\times 10^{2} \approx 9.97\times 10^{3}$ なので、最初の段だけは 7.225 のせいでぴったり $10^{4}$ ではありません。その下の段は $10^{4}$ です。 $p=2$ では指数が $2p=4$ であることと合います。正規化の列は $7.225\times 10^{2}$ のままです。

Q4. なぜテイラー基底なのか#

ここまでは1次元だったので楽でした。実際の格子には四面体・六面体・プリズム・ピラミッド・多面体が混在します。標準的な有限要素は形状ごとに参照要素へ写像し、その上で形状関数を定義します。ヤコビアン作業が形状ごとに一式必要になるということです。多面体には参照要素がそもそも存在しません。

Luoらが提案したテイラー基底は写像を飛ばします。セル中心 $\mathbf{x}_c$ でそのままテイラー展開するだけです。

U_h = \bar{U} + \left.\frac{\partial U}{\partial x}\right|_c (x - x_c) + \left.\frac{\partial U}{\partial y}\right|_c (y - y_c) + \left.\frac{\partial^2 U}{\partial x^2}\right|_c \frac{(x - x_c)^2}{2} + \cdots

各項からその項自身のセル平均を引いておくと、先頭係数 $\bar{U}$ がちょうどセル平均になります。この性質が実務では大きいのです。 $p=0$ とすればDGは有限体積法と完全に一致し、有限体積用のリミターをそのまま載せられます。有限体積のリミターがDGコードで再利用される通り道がこれです。セル形状を問わないので、混合格子でもコードが一式で済みます。

代わりに対価が一つあります。 $(x-x_c)^k$ をそのまま使うと質量行列の成分が $h^{k+l+1}$ でスケールします。条件数がセル幅とともに暴走します。境界層格子のように $h$ が $10^{-3}$ 程度になるとどうなるか測ってみましょう。

from math import factorial, sqrt
 
def taylor_mass(h, K, scale):
    """幅 h のセル上でのテイラー基底 b_k = ((x-xc)/scale)^k / k! の質量行列。"""
    M = [[0.0]*K for _ in range(K)]
    for i in range(K):
        for j in range(K):
            n = i + j
            if n % 2:                                   # 重心まわりの奇数モーメントは0
                continue
            M[i][j] = (h/scale)**n * h / (2**n * (n+1) * factorial(i) * factorial(j))
    return M
 
def jacobi_eig(A, sweeps=60):
    """対称行列の固有値 — 巡回ヤコビ回転。"""
    K = len(A); A = [row[:] for row in A]
    for _ in range(sweeps):
        for p in range(K-1):
            for q in range(p+1, K):
                if abs(A[p][q]) < 1e-300:
                    continue
                th = 0.5*(A[q][q]-A[p][p])/A[p][q]
                t  = (1 if th >= 0 else -1)/(abs(th)+sqrt(th*th+1))
                c  = 1/sqrt(t*t+1); s = t*c
                for k in range(K):
                    akp, akq = A[k][p], A[k][q]
                    A[k][p], A[k][q] = c*akp - s*akq, s*akp + c*akq
                for k in range(K):
                    apk, aqk = A[p][k], A[q][k]
                    A[p][k], A[q][k] = c*apk - s*aqk, s*apk + c*aqk
    return [A[k][k] for k in range(K)]
 
print(" h        raw Taylor      normalized")
for h in (1.0, 1e-1, 1e-2, 1e-3):
    out = []
    for scale in (1.0, h):
        ev = [abs(v) for v in jacobi_eig(taylor_mass(h, 3, scale))]
        out.append(max(ev)/min(ev))
    print(f" {h:<8.0e} {out[0]:.3e}       {out[1]:.3e}")

 h        raw Taylor      normalized
 1e+00    7.225e+02       7.225e+02
 1e-01    7.200e+06       7.225e+02
 1e-02    7.200e+10       7.225e+02
 1e-03    7.200e+14       7.225e+02

$h$ が10分の1になるたびに条件数が $10^4$ 倍になります。 $p=2$ では指数が $2p$ です。 $h = 10^{-3}$ なら $7.2 \times 10^{14}$ に達し、倍精度の $10^{16}$ の余裕をほぼ使い切ります。セル幅で正規化した右列は $h$ に関係なく722を保ちます。セル内で $\Delta x$ で割る一行がその差を作ります。 $p=3$ では指数が6になるので、正規化なしでは実用格子では使えません。

Q5. あらかじめ表に入れておくものは何か#

DGコードの初期化段階は、実質的に表を作る作業です。順序はこうなります。

セルを形状別に分類する — 四面体／六面体／プリズム／ピラミッド／多面体。
面を形状別に分類する — 三角形／四角形／多角形。
形状ごとに必要な次数のガウス求積則を用意する。
各ガウス点で基底関数の値とその勾配を計算して保存する。

3次元における $p$ 次完全多項式空間の自由度は $\binom{p+3}{3}$ です。

$p$	0	1	2	3	4
セルあたりモード数 $K$	1	4	10	20	35

原資料の (1,4,10,20,35) がこの行で、*3 が付いている方は各モードの勾配3成分です。3次元圧縮性解析では保存変数が5個なので、 $p=2$ の六面体格子では状態ベクトルだけでセルあたり $5 \times 10 \times 8 = 400$ バイトです。さらにガウス点ごとの基底値が加わります。 $p=2$ の六面体で体積求積を $3^3 = 27$ 点にすれば、セルあたり $27 \times 10$ 個の実数が追加で必要になります。

この表をセルごとに個別に持つ必要はありません。参照座標での基底値は形状が同じなら同一だからです。形状ごとに一式だけ作り、セルにはヤコビアンと重心・寸法だけを持たせれば足ります。多面体だけが例外的に自分の表を持ちます。

P1からP2へ上げるとき、どこで代金を払うのか#

$p$ を1から2へ上げると、3次元のモード数は4から10へ、メモリは2.5倍です。体積点の下限も $n \ge p$ に従い、テンソル積なら $n^3$ なので点は8倍になります。安定CFLはおよそ $1/(2p+1)$ で、 $p=1$ の $1/3$ から $p=2$ の $1/5$ へ減り、時間刻みは5分の3になります。テイラー基底なら $\Delta x^k$ の指数が大きくなり、 $h=10^{-3}$ の条件数 $7.2\times 10^{14}$ を避けるには正規化が欠かせません。

滑らかな解が広ければ誤差は $h^{p+1}$ で減るので、 $p$ を上げるほうが格子を細かくするより安くなります。衝撃波が支配するなら、リミターがその利得を削ります。セルを事後に疑う側はMOOD subcellリミターがDGオイラーでやる仕事と同じ場所です。点数を $n < p$ までケチる取引だけは残りません。その線の下では、スキームが別の方程式を解きます。