在非结构网格上拯救梯度的权重 — 移动最小二乘(MLS)近似

近点和远点用同一个权重,梯度会跳到 −4.1;再把 h 缩小,A 就往奇异一侧靠。

非结构网格上梯度崩坏的那天

下午两点,在非结构四面体网格上,LSQ(最小二乘)梯度在某一个网格单元突然跳到负值。相邻单元为 0.3、下一个为 0.5,而中间那个却是 -4.1。回到网格仔细看,那一格的七个邻居距离参差不齐。一个近邻点与六个远邻点以相同的权重进入正规方程。

解决方案在 1981 年就已经存在,只是出自另一个学科的另一个名字。

拨动 h、多项式阶数或样本数 N,导数曲线的 RMSE 会跟着动。

Derivative recovery from 25 scattered samples of u(x) = sin(2πx)

kernel radius h: 0.070polynomial order: 2sample count N: 25

RMSE on derivative (x ∈ [0.05, 0.95]): 0.222

把 $h$ 缩到 0.025,某些评估点的活动节点过少, $\mathbf{A}$ 接近奇异,RMSE 爆炸。反过来把 $h$ 调到 0.20,函数的振幅被平均掉,导数变钝。两个极端之间有一段狭窄的 sweet spot。

1981 年,从曲面拟合借来的式子#

Lancaster 和 Salkauskas 在一篇曲面重建论文中做了一个约定:"在每个评估点,只用该点附近的数据重新拟合多项式。"他们把它命名为 moving least-squares(MLS)。这本是图形学里用来构造光滑曲面的工具。

二十年后,Belytschko 的 EFG(Element-Free Galerkin)把同一个公式借走,作为无网格结构分析的骨架。又过十年,非结构 FVM 把它取来做梯度重构。SPH 的一致性修正导数也站在同一位置。出发点不是网格而是离散点,这一事实让每个领域都召唤出同一个公式。

权重函数把距离转换为影响

MLS 的出发式很简单。在评估点 $\mathbf{x}^*$ 附近,用多项式基底拟合函数。

u^h(\mathbf{x}) = \mathbf{p}(\mathbf{x})^\top \mathbf{a}(\mathbf{x}^*)

$\mathbf{p}$ 是多项式基向量(2D 一阶时为 $[1, \xi, \eta]$ ,二阶时为 $[1, \xi, \eta, \xi^2/2, \xi\eta, \eta^2/2]$ ), $\mathbf{a}$ 是系数向量。

关键在于 $\mathbf{a}$ 是 评估点 $\mathbf{x}^*$ 的函数。评估点一动,系数就被重新求解。这就是 "moving" 最小二乘的名字由来。

系数通过最小化加权残差平方和 $J$ 确定。

J(\mathbf{a}) = \sum_i w(\|\mathbf{x}^* - \mathbf{x}_i\|)\, \bigl(\mathbf{p}(\mathbf{x}_i)^\top \mathbf{a} - u_i\bigr)^2

权重函数 $w$ 起决定性作用。靠近评估点的点权重大,远处的点权重几乎为零。高斯函数 $w = \exp(-(r/h)^2)$ 、三次样条、四次多项式 — 它们有一个共同点:距离超过 $h$ 后衰减到零。这个 $h$ 就是 kernel radius(支撑半径)。

正规方程与形状函数的登场

由 $\partial J / \partial \mathbf{a} = 0$ 得到正规方程:

\mathbf{A}(\mathbf{x}^*)\,\mathbf{a} = \mathbf{B}(\mathbf{x}^*)\,\mathbf{u}

其中

\mathbf{A} = \sum_i w_i\, \mathbf{p}(\mathbf{x}_i)\, \mathbf{p}(\mathbf{x}_i)^\top, \quad \mathbf{B}_i = w_i\, \mathbf{p}(\mathbf{x}_i)

$\mathbf{u}$ 为节点值向量, $w_i = w(\|\mathbf{x}^* - \mathbf{x}_i\|)$ 。把系数代回原式:

u^h(\mathbf{x}^*) = \mathbf{p}(\mathbf{x}^*)^\top \mathbf{A}^{-1} \mathbf{B}\, \mathbf{u} = \sum_i \Phi_i(\mathbf{x}^*)\, u_i

$\Phi_i(\mathbf{x}^*)$ 就是 形状函数(shape function)。但这个形状函数带来一个虽小却麻烦的性质 — $\Phi_i(\mathbf{x}_j) \neq \delta_{ij}$ 。也就是说,在节点上插值并不精确。这在 EFG 中使得 Dirichlet 边界条件的施加变得棘手。在 CFD 梯度重构中并非大问题:我们要的是光滑的导数,而不是节点值的精确插值。

Python — 用 25 个点重新对 sin 函数求导#

将基底取为 Taylor 形式 $p_k(\xi) = \xi^k/k!$ ,系数本身就成为导数 — $a_k \approx u^{(k)}(\mathbf{x}^*)$ 。

import numpy as np
from math import factorial
 
def mls_taylor_1d(xs, us, x_star, h, order=2):
    """在评估点 x_star 处使用 Taylor 基底的加权最小二乘。
    返回的 a[k] 是 u^(k)(x_star) 的近似。
    """
    m = order + 1
    A = np.zeros((m, m))
    B = np.zeros(m)
    for x_i, u_i in zip(xs, us):
        dx = x_i - x_star
        w = np.exp(-(dx / h) ** 2)
        if w < 1e-10:
            continue
        p = np.array([(dx ** k) / factorial(k) for k in range(m)])
        A += w * np.outer(p, p)
        B += w * p * u_i
    return np.linalg.solve(A, B)
 
# 验证 — 恢复 sin(2 pi x) 的导数
np.random.seed(7)
xs = np.sort(np.random.uniform(0, 1, 25))
us = np.sin(2 * np.pi * xs)
 
x_eval = np.linspace(0.05, 0.95, 200)
du_hat = np.array([mls_taylor_1d(xs, us, xq, h=0.07, order=2)[1] for xq in x_eval])
du_true = 2 * np.pi * np.cos(2 * np.pi * x_eval)
 
rmse = np.sqrt(np.mean((du_hat - du_true) ** 2))
print(f"derivative RMSE: {rmse:.3e}")
# derivative RMSE: 2.1e-01   (h=0.07, order=2)

25 个随机点下,导数 RMSE 约 0.21。相对于振幅 $2\pi \approx 6.28$ ,相对误差大约 3%。把 order 提高到 3 并将 $h$ 设为 0.10 会进一步降低误差。但 $h$ 过大会带来另一种问题。

高斯 $w=\exp(-(r/h)^{2})$ 里,权重掉到一半的位置是 $(r/h)^{2}=\ln 2$ 。 $h=0.07$ 时这个半径大约 0.058。区间 0 到 1 上放 25 个点,平均间距是 $1/24\approx 0.042$ ,半权重半径大约一格邻居。代码在 $w<10^{-10}$ 时丢掉点,距离大约 $4.8h$ 。 $h=0.07$ 时是 0.34, $h=0.025$ 时是 0.12。

order 为 2 时未知数是 3 个。 $h=0.025$ 的窗口左右加起来约 0.24,按间距 0.042 只能罩住五六个点。种子 7 那样间距参差,就会出现低于平均的评估点,那里 $\mathbf{A}$ 的条件数会跳。导数 RMSE 0.21 除以振幅 $2\pi\approx 6.28$ ,得到 0.033,大约 3%。

把 $h$ 取成平均间距的 1.5 倍, $0.042\times 1.5\approx 0.063$ ,和代码里的 0.07 几乎重合。推到 3 倍、约 0.13,半权重半径越过两格邻居,高阶导数变钝。单元间距各不相同的非结构网格上,一个全局 $h$ 没法同时套进这个窗口。

去掉权重的最小二乘在歪斜单元上怎么裂开,写在Green–Gauss 与加权最小二乘在非结构网格上分道的地方。那天的 −4.1,就是近点和远点用同一个系数进了正规方程。

核半径 h — 大小两端各自损坏什么#

在梯度重构中, $h$ 在两件事之间做权衡。

小 $h$ : $\mathbf{A}$ 的条件数失控。1〜2 个近邻点承担所有权重,用少于 $m$ 个活动点去拟合 $m$ 阶多项式,结果对噪声非常敏感。
大 $h$ : 远处点的信息混入并被平均化。零阶矩还能保留,高阶导数则被削弱。

实务上, $h$ 取 局部平均节点间距的 1.5〜3 倍,使支撑域内包含约 1.5〜2 倍 $m$ 的活动节点。在非结构网格上必须按单元自适应调整 $h$ — 全局 $h$ 几乎总是失败。

扩展到二维。28 个散布节点上给出 $u = \sin(2\pi x)\cos(2\pi y)$ 。拖动黄色评估点,观察权重的变化。

u(x, y) = sin(2πx) cos(2πy) sampled at 28 scattered nodes

kernel radius h: 0.180polynomial order: 2

exact u(x*): -0.000

MLS u^h(x*): 0.039

error: 0.039

active nodes inside kernel: 28

把评估点拖到左上角空白处,活动节点降到 $m$ 以下, $\mathbf{A}$ 变奇异。再拖回节点密集区,误差又变小。这正是 MLS 在非结构网格上的结构性弱点,也是实际中必须使用自适应 $h$ 的原因。

在 SPH、EFG、FVM 中留下的痕迹#

MLS 今天以不同的名字活在三个地方。

SPH 一致性修正: 标准 SPH 微分在非均匀粒子分布下连一阶精度都会丢失。corrected SPH 与 reproducing-kernel SPH 通过加权最小二乘修正恢复一阶一致性 — 本质上就是一阶 MLS。
EFG / MLPG 无网格结构分析: 形状函数就是 MLS。节点处插值不精确使 Dirichlet 边界条件的施加变得棘手,Lagrange 乘子法或惩罚法是标准对策。
非结构 FVM 梯度重构: Compact Least-Squares Reconstruction 在单元中心上求解正规方程,不仅取出梯度还提取曲率。OpenFOAM 的 leastSquares 选项也是一阶多项式 MLS 的变体。

三个领域解的是同一组正规方程。变化的只有权重函数和多项式基底的维数。

回到梯度崩坏的那天

回到开头的那个单元。LSQ 梯度之所以跳到 -4.1,原因很简单:两个近邻点与五个远邻点以相同权重进入了正规方程。引入距离权重之后,-4.1 变成 0.4,发散消失了。

权重函数构造的局部性、每个评估点重新求解的正规方程、形状函数的非可分性 — 这三件事是 MLS 从图形学里借给我们的骨架。今天所有无网格积分与一半的非结构积分,都站在这副骨架之上。

用同一权重解正规方程时梯度跳到 −4.1 的那个单元#

两个近邻和五个远邻以同一个系数进去,那天的梯度就跳到了 −4.1。加上距离权重之后才回到 0.4。把 $h$ 取得远小于节点间距,这次换成 $\mathbf{A}$ 往奇异一侧靠。符号可以对,大小还会再晃。

梯度救回来之后,下一道判断是这个梯度允许进单元到什么程度。非结构网格上二阶重构造出假极值的那段,写在Barth–Jespersen 限制器切开的位置。权重把梯度还回来,若不再加上限制器就原样送进去,那个单元还会再长出极值。