[论文评述] 手推雅可比矩阵推到心累？——算子重载自动微分（ADOO）

手推 HLLC 通量雅可比,三张纸丢掉了。乘积法则漏了一项,代码无声地发散。通量没错。错的是导数。隐式(implicit)CFD 有一半的活,是把空间离散的雅可比矩阵填准。Fraysse 等人(2019)的 ADOO(算子重载自动微分)把这份导数交给代码。Godunov 那种内部带求根的格式,也不再写手推公式。Newton 从这里赚多少,收敛曲线自己会说。

把有限差分步长 h 调小,误差曲线先降一次,随后重新爬升。拖动 x,对偶路径仍贴在底上。

谷底大约在 $h\sim 10^{-8}$。离开那个口袋,只有橙色曲线会再抬起来。

论文信息

• 标题: Automatic Differentiation using Operator Overloading (ADOO) for implicit resolution of hyperbolic single phase and two-phase flow models
• 作者: G. Fraysse 等
• 年份: 2019
• 关键词: automatic differentiation, implicit, two-phase, finite volume, unstructured meshes

一句话总结：通量计算的代码保持不动，只换数据类型，就取出精确的雅可比。

雅可比为什么是个麻烦

隐式时间推进在每一步用Newton迭代求解非线性残差 $\mathbf{R}(\mathbf{Q}) = 0$。

这里 $\mathbf{J}$ 是残差的雅可比（对各守恒变量偏导构成的矩阵），$\delta\mathbf{Q}$ 是更新量。要让Newton二次收敛，$\mathbf{J}$ 必须精确。

有三条路，每条都有软肋。手推（analytic）精确，但格式一改就得重推。遇到AUSM+这类分支繁多、或Godunov这类含迭代的通量，推导本身就是噩梦。有限差分（finite difference）能复用代码，但步长 $h$ 得在舍入误差和截断误差之间进退两难。

Newton-Krylov-matrix-free只对雅可比-向量乘积做有限差分，因此从不组装矩阵。但要停住没有预处理矩阵就会转上千次的 Krylov，终究还是需要矩阵的真实元素。ADOO用计算导数而非近似导数斩断这个两难。

对偶数：给数值挂上导数

想法很朴素。把实数 $x$ 换成一个含两个分量 $(v, dv)$ 的对象：$v$ 是值，$dv$ 是它在该点的导数。这就是对偶数。

算术规则不过是乘积法则和链式法则的照搬。

每条规则左边是值，右边更新导数分量。以 $f(x) = \sin(x^2) + x^2$ 为例。手推得 $f'(x) = \cos(x^2)\cdot 2x + 2x$。用对偶数，把 $x$ 的导数分量设为 $1$（因为 $dx/dx = 1$），照常运行代码，最后那个对象的 $dv$ 就是 $f'(x)$。

from dataclasses import dataclass
import math
 
@dataclass
class Dual:
    v: float   # 值
    d: float   # 导数分量
 
    # 算子重载 —— 每条规则即乘积法则/链式法则（论文 III.2）
    def __add__(self, o):
        o = o if isinstance(o, Dual) else Dual(o, 0.0)
        return Dual(self.v + o.v, self.d + o.d)
 
    def __mul__(self, o):
        o = o if isinstance(o, Dual) else Dual(o, 0.0)
        return Dual(self.v * o.v, self.v * o.d + self.d * o.v)
 
    def __truediv__(self, o):
        o = o if isinstance(o, Dual) else Dual(o, 0.0)
        return Dual(self.v / o.v, (self.d * o.v - self.v * o.d) / (o.v * o.v))
 
def sin_d(a: Dual) -> Dual:
    return Dual(math.sin(a.v), math.cos(a.v) * a.d)
 
# 在 x=1.3 处对 f(x) = sin(x^2) + x^2 求导
x = Dual(1.3, 1.0)          # 种子: dx/dx = 1
f = sin_d(x * x) + x * x    # 原式原封不动
print(f.v, f.d)             # 值, f'(1.3)
print(math.cos(1.3**2) * 2 * 1.3 + 2 * 1.3)  # 与手推对照

输出的 f.d 与手推导数吻合到舍入量级。这里手推的公式只是校验用，从未进入实际计算——这正是关键。

为什么不能信任有限差分

ADOO的真正价值，要和有限差分并排才看得清。中心差分 $[f(x+h)-f(x-h)]/2h$ 的误差是两股力的角力。$h$ 大时截断误差（Taylor展开余项）主导；$h$ 太小时舍入误差（两个相近数相减丢失有效位）爆炸。于是误差对 $h$ 画出一条U形。最优点大致在 $h \sim \sqrt{\epsilon} \approx 10^{-8}$ 附近，可即便如此也因问题而异。

上面的橙色曲线就是那条 U，青色虚线是没有步长的底。AD 根本没有步长，也就没有可调的步长。

Euler通量雅可比：手推 vs AD#

现在离开标量。一维可压缩Euler有守恒变量 $\mathbf{Q} = [\rho, \rho u, \rho E]^\top$ 和通量 $\mathbf{F}(\mathbf{Q})$。理想气体下通量雅可比 $\partial\mathbf{F}/\partial\mathbf{Q}$ 是那个著名的 $3\times3$ 矩阵，但手推时 $\gamma$ 与动能项纠缠不清。把对偶数扩成向量（让每个 $dv$ 变成三分量数组），一次运行就把整个矩阵填满。

import numpy as np
 
class DualVec:
    """值 + 对3个自变量的梯度（种子向量）。"""
    def __init__(self, v, grad):
        self.v = v
        self.g = np.asarray(grad, float)  # ∂(this)/∂Q, 长度3
    def __add__(s, o):  return DualVec(s.v + o.v, s.g + o.g)
    def __sub__(s, o):  return DualVec(s.v - o.v, s.g - o.g)
    def __mul__(s, o):
        if isinstance(o, DualVec):
            return DualVec(s.v * o.v, s.v * o.g + s.g * o.v)  # 乘积法则
        return DualVec(s.v * o, s.g * o)
    def __truediv__(s, o):
        return DualVec(s.v / o.v, (s.g * o.v - s.v * o.g) / (o.v * o.v))
 
def euler_flux_jacobian(Q, gamma=1.4):
    # 为每个守恒变量种子: rho -> (Q0, e0), 等
    rho  = DualVec(Q[0], [1, 0, 0])
    rhou = DualVec(Q[1], [0, 1, 0])
    rhoE = DualVec(Q[2], [0, 0, 1])
 
    u = rhou / rho                                  # 速度
    kinetic = rhou * u * 0.5                         # ½ρu²
    p = (rhoE - kinetic) * (gamma - 1.0)             # 压力（理想气体）
 
    F0 = rhou                                        # ρu
    F1 = rhou * u + p                                # ρu² + p
    F2 = (rhoE + p) * u                              # (ρE + p)u
 
    # 每个F分量的 .g 就是该行的雅可比（对应论文 式9~10）
    return np.array([F0.g, F1.g, F2.g])
 
Q = np.array([1.2, 0.6, 3.0])   # ρ, ρu, ρE
J_ad = euler_flux_jacobian(Q)
 
# 与手推精确雅可比对照
rho, mom, E = Q
u = mom / rho; g = 1.4
J_ref = np.array([
    [0, 1, 0],
    [0.5*(g-3)*u*u, (3-g)*u, g-1],
    [((g-1)*u**3 - g*u*E/rho), (g*E/rho - 1.5*(g-1)*u*u), g*u],
])
print("最大误差:", np.abs(J_ad - J_ref).max())   # ~1e-15

euler_flux_jacobian 走的是黎曼求解器实际用的通量代码里同一套算术。我们从没推过雅可比，精度却是机器精度。想换成AUSM+？只替换通量函数，雅可比自动跟上。

把手推公式代进代码里的状态，矩阵就收成数字。$Q=(1.2,\ 0.6,\ 3.0)$、$\gamma=1.4$ 时，$u=0.6/1.2=0.5$。动量行是 $0.5(\gamma-3)u^2$、$(3-\gamma)u$、$\gamma-1$，也就是 $(-0.2,\ 0.8,\ 0.4)$。能量行更长。$(\gamma-1)u^3=0.05$，$\gamma u E/\rho=1.75$，第一项便是 $0.05-1.75=-1.70$。从 $\gamma E/\rho=3.5$ 减去 $1.5(\gamma-1)u^2=0.15$，中间项是 $3.35$。最后一项是 $\gamma u=0.7$。质量行仍是 $(0,\ 1,\ 0)$。代码注释把这张矩阵和 AD 的最大差写成 $\sim 10^{-15}$。标量核对停在同一处。$x=1.3$ 时 $x^2=1.69$、$2x=2.6$，手推导数是 $\cos(1.69)\cdot 2.6+2.6$，f.d 只在舍入范围内和它分开。

Newton收敛：精确带来的回报#

精确雅可比为何重要，收敛曲线自会说话。近似雅可比把Newton从二次收敛拉低到一次。迭代次数攀升，大CFL下干脆发散。论文报告，凭借精确雅可比，即便CFL 20也能在不到10次迭代内把残差压到 $10^{-6}$——这是二次收敛的签名。

在下面把Jacobian误差滑块从0往上推。

误差为0%（AD精确雅可比）时，残差每迭代一次位数近乎翻倍下跌，直线俯冲——二次收敛的签名。只给20%误差，曲线就躺成一条平缓直线（一次），达到同样精度要多得多的迭代。把Krylov迭代也算进去，这个差距就是墙钟时间。

批判性思考：ADOO的阴影#

ADOO并非免费。算子重载为每个标量创建对象、乘数组。Fraysse的论文也承认，前向模式的开销与自变量数目成正比——在块很大的多相（multiphase）系统里，这些数组会变重。这正是源码变换（ADSCT，如Tapenade）可借编译期优化更快的原因。

复现时我还撞到一个实务问题。微分Godunov这类含求根迭代的格式时，不仅值要收敛，导数分量也要收敛。导数通常比值收敛得慢，若停止条件只看值，雅可比就会被污染。论文指出的这一句，替我省下好几天。

从OpenFOAM的角度看，这套方法并不陌生。有人试过用对偶标量类型组装 blockLduMatrix，SU2也早已用源码变换AD构造伴随（adjoint）。要点一致——由人手推雅可比的时代正在落幕。

只看函数值就停掉求根，雅可比会被弄脏

通量内部有求根时，Godunov 是常见的那种，值和导数的收敛速度并不一样。导数更慢。停止条件只看值，通量看着对，矩阵元素已经脏了。论文里的精确雅可比，在 CFL 20 下不到 10 次迭代就能把残差压到 $10^{-6}$。混进 20% 的误差，二次就躺成直线。原本每次迭代位数大致翻倍的曲线，变成斜率差不多不变的一条。

前向模式的开销与自变量个数成正比。多相块一大，对偶数组成的数组就沉，源码变换可能更便宜。分支很多的通量，比如用幽灵流体填上界面的那一类，手推会先死掉。停止判据不扩展到导数分量，这份雅可比就会带着错误走进 Newton。