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

非構造格子でセル勾配を求める2つの手法と、歪みに強いのはどちらか

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

有限体積法（FVM）でセル勾配は、拡散フラックス、2次再構成、リミッタの入力になります。セル中心値だけでは足りません。近傍の情報を集めて $\nabla\phi_C$ を推定する必要があります。この記事では、その推定を行う2つの手法 — Green–Gaussと最小二乗 — を比較します。なぜ歪んだ格子で一方が崩れるのか、そして距離による重み付けが何を直すのかを、コードとともに見ていきます。

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個以上あれば必ず解けます。

距離で重み付けする

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

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

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

anisotropyを0.2付近まで下げてステンシルを扁平にし、curvature κを上げると、オレンジ色の再構成勾配が緑色の真値から離れていきます。ここでweight exponent nを0から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つあるだけで誤差が大きくなります。 $n=1$ の距離重み付けがその近傍を抑え、誤差を4分の1ほど減らします。

コードを書くときに陥ってはならない罠

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

下のシミュレーションで直接操作してみましょう。

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

Collinearity: 0.20Neighbor count: 6

Collinearityを0.9付近まで上げると、近傍が一列に集まり、cond(M)が数十から無限大へ跳ね上がります。楕円が長く薄く伸びる方向こそ、勾配が不正確になる方向です。壁にぴったり張り付いた境界層セルでよく起きる状況です。

現場で押さえるべき3つ。第一に、境界セルでは近傍が片側に偏ります。壁の反対側の近傍や面中心をステンシルに追加して、条件数を下げます。第二に、 $M$ を毎ステップ組み直さず、格子が固定なら逆行列（またはQR分解）を事前に保存します。第三に、距離重み $n$ は問題に合わせて調整します。粘性項が支配的なら $n=2$ の方がよいこともあります。

ひと言で残します。Green–Gaussは安いが歪みに弱く、重み付き最小二乗は歪みに強いがステンシルが一列に並ぶと病的になる — 格子を先に見て手法を選びましょう。