格子ペクレ数を2に上げたら解が-0.33まで落ちた — 最小化すべきエネルギーが消えた場所

一本の棒では四つの方法が同じ答えを出した

有限要素法の講義ノート第1章は、同じ棒の問題を四回解きます。直接剛性法、最小ポテンシャルエネルギー、重み付き残差法、そしてガラーキン法です。四回とも同じ5×5の剛性行列が出てきます。ノートにはこう書かれています。どの方法を使っても結果に大きな違いはなく、正確な値が得られる、と。

しかしこの一文には条件が付いています。その条件は棒の問題の中に隠れていて表に出てきません。流れの方程式はその条件を満たしません。条件が壊れたときにどの保証が消え、安定化が何を返して何を返せないのかを数字で測ります。advection a をゼロから上げると、ピンクの点の一周の仕事は0.000からちょうど $2\pi a$ に変わります。

灰色の等エネルギー線を横切らなかった黄色の球は、a が大きくなると螺旋を描いて逆に登ります。diffusion s を端まで動かしても一周の仕事はそのままです。

最小化するには行列が対称でなければならない

最小ポテンシャルエネルギー法は、全ポテンシャルエネルギーを節点未知数で微分してゼロと置きます。離散化した形はこうなります。

$\boldsymbol{\phi}$は節点未知数ベクトル、$\mathbf{K}$は剛性行列、$\mathbf{f}$は荷重ベクトルです。ここで$i$番目の成分で微分します。

$K_{ij} + K_{ji}$が出てきたことに注目してください。この式が$\mathbf{K}\boldsymbol{\phi} = \mathbf{f}$になるためには$K_{ij} = K_{ji}$である必要があります。対称でなければ、最小化が実際に解いているのは対称化された行列$\tfrac{1}{2}(\mathbf{K}+\mathbf{K}^{T})$であり、それは元の方程式ではありません。

同じ場所からもっと根本的な話も出てきます。残差場$\mathbf{r}(\boldsymbol{\phi}) = \mathbf{f} - \mathbf{K}\boldsymbol{\phi}$が何らかのスカラー関数の勾配であるためには、そのヤコビアンが対称でなければなりません。対称でなければポテンシャル関数はそもそも存在しません。存在を確かめる一番速い方法は閉じた経路を一周することです。勾配場なら一周した仕事は必ずゼロになります。上で a を上げたとき球が等エネルギー線を逆に登ったのが、そのゼロが壊れた場所です。

移流項が壊す量はちょうど$u$です#

1次元移流拡散方程式$u\,\phi' = \epsilon\,\phi''$の弱形式を書くと、二つの項の性格の違いが見えてきます。

$w$は試験関数、$\epsilon$は拡散係数、$u$は移流速度です。第1項は$w$と$\phi$を入れ替えても変わりません — 対称です。後ろの項は部分積分すると符号が反転します。境界で消える試験関数に対して$\int w\,u\,\phi'\,dx = -\int \phi\,u\,w'\,dx$ですから、移流項は純粋に反対称な成分です。

線形要素で組み立てると、この構造がそのまま係数に残ります。要素サイズを$h$とすると、内部節点の行はこうなります。

拡散は両隣に同じ値を入れ、移流は片側に$+u/2$、反対側に$-u/2$を入れます。だから対称からのずれはちょうど一つの数になります。

格子をいくら細かくしてもこの値は$u$のままです。$h$が入っていないからです。流れがある限り、最小ポテンシャルエネルギー原理は戻ってきません。剛性を6×6のVoigt行列に畳めたのは、行列が対称だったときだけです。移流が段差の近くで2次スキームを揺らすときの処方は流束制限子のTVD移流の側にあります。

Pythonで同じ格子に三つの手法を立てた#

要素10個、$u = 1$、境界条件$\phi(0)=0$と$\phi(1)=1$。人工拡散係数$\beta$一つで三つのスキームを作ります。$\epsilon_{\text{eff}} = \epsilon + \beta\,u\,h/2$として、$\beta = 0$なら純粋ガラーキン、$\beta = 1$なら完全風上差分、$\beta = \coth(Pe_h) - 1/Pe_h$ならSUPGです。

import numpy as np
 
def assemble_ad(n, vel, eps, beta):
    """線形要素上の1次元移流拡散。beta は人工拡散係数。"""
    h = 1.0 / n
    eps_eff = eps + beta * vel * h / 2.0
    kd = (eps_eff / h) * np.array([[1.0, -1.0], [-1.0, 1.0]])   # 拡散: 対称
    ka = (vel / 2.0) * np.array([[-1.0, 1.0], [-1.0, 1.0]])     # 移流: 反対称
    K = np.zeros((n + 1, n + 1))
    for e in range(n):
        K[e:e + 2, e:e + 2] += kd + ka
    return K
 
def solve_bvp(K):
    n = K.shape[0] - 1
    A, b = K.copy(), np.zeros(n + 1)
    A[0, :], A[0, 0], b[0] = 0.0, 1.0, 0.0
    A[n, :], A[n, n], b[n] = 0.0, 1.0, 1.0
    return np.linalg.solve(A, b)
 
def exact_ad(x, pe):
    return (np.exp(pe * (x - 1.0)) - np.exp(-pe)) / (1.0 - np.exp(-pe))
 
def skew_ratio(K):
    return np.linalg.norm(K - K.T) / np.linalg.norm(K + K.T)
 
def loop_work(K, m=20000):
    """自由度空間の単位円を一周したときの残差場 -K.phi の仕事。"""
    t = np.linspace(0.0, 2.0 * np.pi, m, endpoint=False)
    path = np.stack([np.cos(t), np.sin(t)])          # 位置
    tang = np.stack([-np.sin(t), np.cos(t)])         # dl / dt
    return float(np.sum(np.sum(-(K @ path) * tang, axis=0)) * (2.0 * np.pi / m))
 
n, vel = 10, 1.0
x = np.linspace(0.0, 1.0, n + 1)
print(f"{'Pe_h':>5} {'scheme':>9} {'max|err|%':>10} {'min phi':>9} {'skew/sym':>9} {'loop W':>8}")
for pe_h in [0.5, 1.0, 2.0, 5.0]:
    eps = vel / (2.0 * n * pe_h)
    ex = exact_ad(x, vel / eps)
    for name, beta in [("Galerkin", 0.0), ("upwind", 1.0),
                       ("SUPG", 1.0 / np.tanh(pe_h) - 1.0 / pe_h)]:
        K = assemble_ad(n, vel, eps, beta)
        phi = solve_bvp(K)
        print(f"{pe_h:5.1f} {name:>9} {100 * np.max(np.abs(phi - ex)):10.2f} "
              f"{phi.min():9.4f} {skew_ratio(K):9.4f} {loop_work(K[4:6, 4:6]):8.4f}")
 
K0 = assemble_ad(n, 0.0, 0.1, 0.0)
print(f"\nvel = 0 : skew/sym = {skew_ratio(K0):.2e},  loop W = {loop_work(K0[4:6, 4:6]):.2e}")
print(f"pi * u  = {np.pi * vel:.4f}")

 Pe_h    scheme  max|err|%   min phi  skew/sym   loop W
  0.5  Galerkin       3.45    0.0000    0.2924   3.1416
  0.5    upwind      13.17    0.0000    0.1954   3.1416
  0.5      SUPG       0.00   -0.0000    0.2704   3.1416
  1.0  Galerkin      13.53    0.0000    0.5774   3.1416
  1.0    upwind      19.80    0.0000    0.2924   3.1416
  1.0      SUPG       0.00   -0.0000    0.4428   3.1416
  2.0  Galerkin      35.17   -0.3334    1.1010   3.1416
  2.0    upwind      18.17   -0.0000    0.3885   3.1416
  2.0      SUPG       0.00    0.0000    0.5572   3.1416
  5.0  Galerkin      69.61   -0.6961    2.1517   3.1416
  5.0    upwind       9.09    0.0000    0.4836   3.1416
  5.0      SUPG       0.00    0.0000    0.5773   3.1416
 
vel = 0 : skew/sym = 0.00e+00,  loop W = -4.29e-16
pi * u  = 3.1416

Pe_h スライダーを1の向こうへ押すと、右の a_E が赤く越えて節点値が0の下へ沈みます。

beta をゼロのまま Pe_h を2にすると、表の min phi のように -0.3334 付近が見えます。SUPG を押すと誤差は0 %に張りつき、ピンクの反対称の行は3.1416から動きません。

境界値は0と1に固定されているのに、$Pe_h = 2$のガラーキン解は$-0.3334$まで下がります。$Pe_h = 5$では$-0.6961$です。物理的に出てはいけない値です。

$Pe_h = 1$で係数が一つ符号を変える#

格子ペクレ数を$Pe_h = u h / (2\epsilon)$と定義すると、先ほどの行の右側の係数はこう書き直せます。

$Pe_h > 1$になると非対角成分が正になります。この瞬間に行列はM行列であることをやめ、離散最大値原理 — 内部の解が境界値の作る区間を出ないという保証 — も一緒に消えます。表でmin phiが$Pe_h = 1.0$まではゼロで、2.0で負になる理由がこれです。上で a_E が赤く越えた瞬間が、この符号の変化です。

安定化は対称性を取り戻さない

ここでよくある読み違いがあります。風上差分やSUPGが「安定性を回復させる」という言い方は、対称性と最小化原理を取り戻すという意味ではありません。表のloop W列を見てください。三つのスキーム、四つの$Pe_h$ — 全12行がすべて3.1416です。この値は$\pi u$であり、$\epsilon$も$\beta$も入っていません。

理由は単純です。$\beta$が大きくするのは$\epsilon_{\text{eff}}$だけで、それは行列の対称部分です。反対称部分の$u$には手を触れません。skew/sym比が$Pe_h = 5$で2.1517から0.4836に減ったのも、反対称が小さくなったのではなく、分母の対称部分が大きくなったからです。

$Pe_h = 0.5$ のガラーキンは max|err|% が3.45、min phi が0.0000、skew/sym が0.2924です。$Pe_h = 1.0$ では誤差が13.53まで上がっても最小値はまだ0.0000です。$K_{i,i+1} = (\epsilon/h)(Pe_h - 1)$ で $Pe_h = 1$ なら右側の係数はゼロで、1を越えた瞬間に正になります。$Pe_h = 2$ のガラーキンの最小値は -0.3334、誤差は35.17 %、skew/sym は1.1010です。同じ行の風上差分は最小値 -0.0000、誤差18.17 %、skew/sym 0.3885です。SUPGは誤差0.00、skew/sym 0.5572です。三行の loop W はすべて3.1416です。$u = 1$ のとき $\pi u$ が3.1416なので、$\beta$ が0でも1でも $\coth(Pe_h) - 1/Pe_h$ でも一周の仕事は同じです。$Pe_h = 5$ のガラーキンは最小値 -0.6961、誤差69.61 %、skew/sym 2.1517です。風上差分は誤差を9.09 %へ、skew/sym を0.4836へ下げます。小さくなったのは分母の対称部分です。速度ゼロの対照は skew/sym 0.00e+00、loop W -4.29e-16です。要素10個、$u = 1$ なら $h = 0.1$ で、$Pe_h = uh/(2\epsilon)$ から $Pe_h = 2$ の $\epsilon$ は0.025です。その $\epsilon$ を $\beta$ で膨らませても $K_{i,i+1} - K_{i+1,i} = u$ は1のままです。

つまり安定化が実際に買っているのは、より弱い保証です。最小化原理(エネルギーノルムでの最良近似)を失った代わりに、M行列の性質と離散最大値原理を手に入れます。代金は精度で払います。$Pe_h = 2$で風上差分は振動を消しましたが、誤差は18.17 %です。一方SUPGは$\beta$を必要な分だけ入れるので節点で正確です。

この値は$Pe_h \to 0$でゼロへ、$Pe_h \to \infty$で1へ向かいます。拡散が支配するなら安定化を切り、移流が支配するなら完全風上差分へ行くという意味です。ただし節点での正確さは1次元定数係数問題の特権です。2次元では流線方向にだけ人工拡散を入れる本来のSUPGの形が必要で、そのときも「正確」は消えます。

有限体積法が同じ場所で使う別の名前

この計算は有限要素法の言葉で行いましたが、結論は離散化手法に依りません。有限体積法でも中心差分の対流項はまったく同じステンシルを作り、同じ$Pe_h = 1$で係数が符号を変えます。そこで登場するのが1次風上差分で、それが持ち込む数値拡散は$\beta = 1$のときの$u h / 2$と同じ値です。

二つの世界が分かれるのは保存性です。有限体積法の風上差分は面フラックスを変える方式なので、総量保存が保たれます。保存形と原始形が衝撃波速度で分かれたのと同じ種類の分岐です。有限要素法の人工拡散は剛性行列に項を足すので、何を保存しているかを別に見る必要があります。

重み付き残差の系列の中で別の選択肢を選ぶと、また別の代価を払います。最小二乗法は残差ノルムを最小化するので、常に対称正定値行列を作ります — 最小化原理が戻ってきます。その代わり条件数が二乗され、線形要素では要素内部の2階微分がゼロなので拡散項がまるごと消えます。基底と積分則が定式化の性質を静かに変える場所です。

$Pe_h$ が2なら、細分化の前に有界性を買う#

講義ノートの「どの方法でも同じ」という文は、自己随伴作用素の内側でだけ当たります。拡散、弾性、ポテンシャル流はその中です。移流が入ると $K_{i,i+1} - K_{i+1,i} = u$ で対称が壊れ、loop W は $\epsilon$ をどう変えても3.1416、つまり $\pi u$ に固定されます。そのときCGが約束していたエネルギーノルムの最良近似はありません。非対称な行列はKrylov側、JacobiからBiCGSTABまでへ渡すほうが正直です。

min phi が負ならセルペクレ数は1を超えています。$Pe_h = 2$ でガラーキンは -0.3334、誤差35.17 %です。風上差分は振動を消しますが誤差18.17 %を残します。SUPGの節点誤差は0.00なのに、skew/sym は1.1010から0.5572へ下がっただけで、反対称そのものはそのままです。3次元で $Pe_h$ を5から1へ下げるとセルは125倍です。その倍数を持てないときに買うのは精度ではなく有界性です。