積分点を一つ減らしたら解が発散した — 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 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.

DG order p と Gauss points n を別々に動かしてみてください。 $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は多項式係数の複数個で行っているにすぎません。ですから保存形と原始形が分かれる場所で見た話がそのまま当てはまります。フラックス差分の構造を失えば、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 を増やすとより早く壊れます。

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は有限体積法と完全に一致し、有限体積用のリミターをそのまま載せられます。Barth–Jespersen・Venkatakrishnanリミターが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$ に従って一緒に上がります。3次元のテンソル積なら $n^3$ なので点数は8倍です。第二に、陽的時間進行の安定CFLがおよそ $1/(2p+1)$ で減り、時間刻みが5分の3に短くなります。第三に、テイラー基底を使うなら正規化定数 $\Delta x^k$ の指数が大きくなり、条件数の管理が必須になります。

三つの費用を払う価値があるかは問題が決めます。滑らかな解が広く分布する問題なら、 $p$ を上げるほうが格子を細かくするより安上がりです。誤差が $h^{p+1}$ で減るからです。衝撃波が支配する問題なら、リミターが $p$ の利得の大半を削ります。いずれにせよ、求積点を節約しようとして $n < p$ まで下げることだけは得になりません。その線の下では精度が少し悪くなるのではなく、スキームが別の方程式を解きます。