歪んだ格子で勾配がずれるとき — Green–Gaussと重み付き最小二乗による再構成

線形場は乱れたステンシルでも小数6桁まで一致し、曲面は遠い近傍をn=1で抑えなければ誤差が残る。

残差が妙に寝てきて、そのまま発散しました。ログを追うと、壁近くのあるセルで温度勾配が実際とほぼ直角にずれていました。格子を見ると、そのセルだけがやけに扁平でした。値は正常です。勾配を作る手法が、格子の形につまずいたのです。

有限体積法（FVM）でセル勾配は、拡散フラックス、2次再構成、リミッタの入力になります。セル中心値だけでは足りません。近傍の情報を集めて $\nabla\phi_C$ を推定する必要があります。比べるのは Green–Gauss と最小二乗の二つです。歪んだ格子でどちらが崩れるか、距離の重みが誤差をどれだけ戻すかを、出力の数字で追います。

anisotropy を 0.2 付近まで下げるとステンシルが扁平になります。weight exponent を 0 から 1 に上げると、オレンジの矢印が緑の真値へ戻ります。

angle error 16.8° · |∇φ| 0.64 vs 1.08 true — the reconstructed gradient is skewed

Stencil anisotropy: 0.35Field curvature κ: 1.40Weight exponent n: 0.0Neighbor clustering: 0.60

curvature κ も一緒に上げると、二つの矢印の間隔が開きます。扁平なテンプレートで、重みが向きを戻すかどうかだけを見ます。

Green–Gauss: 面で囲んで平均する#

発散定理をセルにそのまま適用します。体積積分を面積分に変えます。

\nabla \phi_C = \frac{1}{V_C} \sum_{f} \phi_f \, \mathbf{n}_f \, A_f

$V_C$ はセル体積、 $\mathbf{n}_f$ は面の外向き法線、 $A_f$ は面の面積、 $\phi_f$ は面上の値です。

面の値 $\phi_f$ は直接はわかりません。セル中心値から補間します。もっとも一般的な線形補間は次のとおりです。

\phi_f = g_f\, \phi_C + (1 - g_f)\, \phi_N

$g_f$ は面から近傍中心 $N$ までの距離比で決めた重みです。計算が安く、直感的です。構造格子では2次精度が出ます。

歪んだ格子がGreen–Gaussを揺さぶる#

問題は $g_f$ に隠れています。線形補間は、面中心と2つのセル中心が一直線上にあると仮定します。非構造格子では、この仮定がほぼ必ず崩れます。

面中心が2つのセル中心を結ぶ線から外れていると（ $\text{skewness}$ ）、補間された $\phi_f$ が偏ります。セルが扁平だったり、近傍の大きさが大きく異なったりすると、偏りはさらに大きくなります。この誤差は反復のたびに積み重なります。結局、勾配が見当違いの方向を指します。

補正項を付けると緩和されます。しかし補正には、また勾配が必要です。循環が生じます。反復補正は高価で、歪みが激しいと収束が遅くなります。

最小二乗: 近傍に平面をあてはめる

発想を変えます。面を経由しません。セル $C$ 周辺の値の分布に、直接1枚の平面をあてはめます。

近傍 $k$ がオフセット $\mathbf{d}_k = \mathbf{x}_k - \mathbf{x}_C$ にあるとき、線形近似は $\phi_k \approx \phi_C + \nabla\phi_C \cdot \mathbf{d}_k$ です。すべての近傍について、この誤差の二乗和を最小化します。

\min_{\nabla\phi_C} \sum_k w_k \left( \phi_k - \phi_C - \nabla\phi_C \cdot \mathbf{d}_k \right)^2

$w_k$ は近傍ごとの重みです。微分して0と置くと正規方程式が得られます。

\left( \sum_k w_k\, \mathbf{d}_k \mathbf{d}_k^{\mathsf T} \right) \nabla\phi_C = \sum_k w_k \left( \phi_k - \phi_C \right) \mathbf{d}_k

2次元なら $2\times2$ のシステムです。書き下すと次のようになります。

\begin{pmatrix} \sum w\, d_x^2 & \sum w\, d_x d_y \\[2pt] \sum w\, d_x d_y & \sum w\, d_y^2 \end{pmatrix} \begin{pmatrix} \partial_x\phi \\[2pt] \partial_y\phi \end{pmatrix} = \begin{pmatrix} \sum w\, \Delta\phi\, d_x \\[2pt] \sum w\, \Delta\phi\, d_y \end{pmatrix}

$\Delta\phi = \phi_k - \phi_C$ です。左辺の行列は近傍の配置だけで決まります。セルごとに一度計算して保存できます。利点は大きいです。面補間の仮定がありません。歪みに鈍感です。近傍が3個以上あれば必ず解けます。

面補間が崩れる地点は、歪んだセルでの Green–Gauss と最小二乗と同じ場所です。近傍が3個でも正規方程式は解けますが、その3点が一直線なら横方向の答は空になります。

距離で重み付けする

重み $w_k$ が精度を左右します。よくある選択は距離の逆数です。

w_k = \frac{1}{\lvert \mathbf{d}_k \rvert^{\,n}}

$n=0$ ならすべての近傍を同じに見ます。遠くにある近傍は値の差が大きく、その分だけ曲率誤差も大きくなります。その近傍が平面のあてはめを引っ張ります。 $n$ を1に上げると、遠い近傍の発言権が減ります。局所線形の仮定がよく当てはまる近い近傍に重みが乗ります。非粘性の問題では通常 $n=1$ を使います。遠い近傍の票を削り、曲率誤差を落とす選択です。

Python — 線形場は正確に、曲面は重み付けで#

最小二乗の力を2回に分けて確認します。まず線形場 $\phi = 1.3x - 0.7y + 2$ の勾配を、でたらめに散らばったステンシルで再構成します。理論上、誤差は0のはずです。次に曲面 $\phi = e^{0.7x}\cos(0.9y)$ で、重みなしの側と $n=1$ の重み付けを競わせます。

import numpy as np
 
def weighted_lsq_gradient(offsets, dphi, n):
    # offsets: (K,2) 近傍オフセット d_k; dphi: (K,) 差分 phi_k - phi_C
    w = 1.0 / np.maximum(np.linalg.norm(offsets, axis=1), 1e-12) ** n
    M = (offsets[:, :, None] * offsets[:, None, :] * w[:, None, None]).sum(axis=0)
    b = (offsets * (dphi * w)[:, None]).sum(axis=0)
    return np.linalg.solve(M, b)
 
xc = np.array([0.4, -0.2])
 
# (A) 線形場 — どんなステンシルでも最小二乗は正確
lin = lambda x, y: 1.30 * x - 0.70 * y + 2.0
scatter = np.array([[0.45, 0.30], [-0.40, 0.12], [0.05, -0.48], [-0.33, -0.29],
                    [0.90, 0.15], [-0.11, 0.44], [0.28, -0.09]])   # 不規則なステンシル
dlin = np.array([lin(*(xc + d)) - lin(*xc) for d in scatter])
gA = weighted_lsq_gradient(scatter, dlin, 1.0)
print(f"linear field   LSQ = ({gA[0]:+.6f}, {gA[1]:+.6f})   exact = (+1.300000, -0.700000)")
 
# (B) 曲面 — 距離による重み付けが再構成誤差を削る
phi  = lambda x, y: np.exp(0.7*x) * np.cos(0.9*y)
grad = lambda x, y: np.array([0.7*np.exp(0.7*x)*np.cos(0.9*y),
                             -0.9*np.exp(0.7*x)*np.sin(0.9*y)])
near = np.array([[0.10, 0.05], [-0.08, 0.09], [0.06, -0.10], [-0.09, -0.06], [0.11, 0.10]])
offs = np.vstack([near, [0.90, 0.70]])          # 遠くに離れた近傍を1つ
dphi = np.array([phi(*(xc + d)) - phi(*xc) for d in offs])
gt   = grad(*xc)
for n in (0.0, 1.0):
    g = weighted_lsq_gradient(offs, dphi, n)
    err = np.linalg.norm(g - gt)
    print(f"curved field   n={n:.0f}  LSQ = ({g[0]:+.3f}, {g[1]:+.3f})   |error| = {err:.4f}")
print(f"curved field   true = ({gt[0]:+.3f}, {gt[1]:+.3f})")

出力は次のようになります。

linear field   LSQ = (+1.300000, -0.700000)   exact = (+1.300000, -0.700000)
curved field   n=0  LSQ = (+0.885, +0.200)   |error| = 0.0293
curved field   n=1  LSQ = (+0.891, +0.205)   |error| = 0.0222
curved field   true = (+0.911, +0.213)

線形場では、ステンシルがどれだけ歪んでも小数点以下6桁まで一致します。出力の (+1.300000, -0.700000) がその値です。線形関数は近傍の位置に関係なく、差分と勾配の内積がぴったり合います。不規則なステンシルに (0.90, 0.15) が入っていても誤差は出ません。歪みはこの場では解を押しません。

曲面は違います。真値は (+0.911, +0.213) です。n=0 の (+0.885, +0.200) は x が 0.026、y が 0.013 足りません。n=1 の (+0.891, +0.205) では不足が 0.020 と 0.008 になります。ノルムは 0.0293 から 0.0222 へ下がります。差 0.0071 は約 24% で、4分の1という言い方と揃います。遠い近傍 (0.90, 0.70) の成分は、近い五点の 0.1 前後よりずっと大きいです。n=0 ではその点が同じ票で平面を引っ張ります。n=1 では距離が票を削ります。

近傍が一列に並んだら逆行列を信じない

最小二乗が歪みに強いからといって万能ではありません。左辺の行列 $M = \sum w_k \mathbf{d}_k \mathbf{d}_k^{\mathsf T}$ が病的になることがあります。近傍が一直線上に並ぶと、 $M$ は特異に近づきます。ステンシルの方向には勾配がよく定まりますが、横方向にはほぼ自由になります。条件数が爆発します。

Collinearity を 0.9 付近まで上げると、近傍が一列に集まり、cond(M) が数十から無限大へ跳ねます。

cond(M) = 1.4 (λ_max 3.14 / λ_min 2.271) — well-conditioned

Collinearity: 0.20Neighbor count: 6

楕円が長く薄く伸びる軸が、勾配を定められない軸です。壁に張り付いた境界層セルがこの形です。対角はオフセット二乗の和なので、近傍が一方向にしか長いと反対の対角が空きます。境界セルでは近傍が流体側に寄ります。壁の反対側の近傍や面中心をステンシルに足すと、空いた軸が少し埋まり、条件数が下がります。格子が固定なら $M$ は速度場が変わってもそのままです。毎ステップ組み直さず、逆行列か QR を保存します。距離の指数 $n$ は問題ごとに違います。粘性項が大きいときは $n=2$ の方がよい場合もあります。

移動最小二乗は、この重みを近傍ごとに組み直す側へもう一歩進みます。Green–Gauss は安い一方でスキューに弱く、重み付き最小二乗はスキューに強いものの、ステンシルが一列なら行列が先に病みます。手法を選ぶ前に、そのセルの近傍が平面に広がっているかを見ます。