[论文评论] 把一次 Newton 拆成两次,自由水面就能解 — Casulli–Zanolli (2012)

$V(\eta) + T\eta = \mathbf{b}$ 。只有一行的方程。第一次看见它时,以为一次 Newton 就能收尾。直到撞上井缘——水深跨过零的胞——Jacobian 在那里悄无声息地变零。同一个胞,这一时刻浸在水里,下一时刻又彻底干掉(wetting/drying)。所有基于 Newton 的自由水面代码都在这一关栽过跟头。Casulli 与 Zanolli 在 2012 年所做的事,是把 Newton 本身切成两半。外层迭代自下而上推进,内层迭代自上而下压回——两侧夹住,单调收敛被严格保证。

论文信息#

作者: Vincenzo Casulli, Paola Zanolli
期刊: Journal of Computational and Applied Mathematics, vol. 236 (2012), pp. 3937–3947
DOI: 10.1016/j.cam.2012.04.025
关键词: mildly nonlinear systems / wetting and drying / free-surface hydrodynamics / Richards equation / confined–unconfined aquifers

Q1. 自由水面上 Newton 为何频频死掉#

方程看起来温和。

\mathbf{V}(\eta) + T\eta = \mathbf{b}

其中 $\eta \in \mathbb{R}^N$ 是各胞的 piezometric head, $T$ 是对称半正定稀疏矩阵(离散扩散算子), $\mathbf{V}$ 是对角非线性函数——每胞体积 $V_i(\eta_i) = \int_{-\infty}^{\eta_i} a_i(z)\,dz$ 。 $a_i$ 只需是 有界变差(允许不连续)。三类典型的自由水面例:

wetting/drying: $V(\eta) = \max(0, \eta)$ 。干胞处 $V'=0$ 。
confined–unconfined aquifer: $V(\eta) = \max(0, \min(c-(-h), \eta-(-h)))$ 。两个折点。
Richards 方程(mixed form): $V(\eta)$ 取决于土壤,光滑但近乎平坦的区间宽。

三者的共性是: $V'(\eta)$ 在整段区间上为零。标准 Newton 用 $J = V'(\eta^{(k)}) + T$ 的逆;当 $V'=0$ 且 $T$ 的相应行弱耦合(或落入 $T$ 的零空间)时, $J$ 退化。初值若不够接近真解便发散——这正是自由水面代码三十年来的顽疾。

Q2. 把 $\mathbf{V}(\eta)$ 拆成两段 — Jordan 分解#

第一招是函数分解。 $a_i$ 有界变差意味着可写成两个非减函数之差(Jordan 分解)。

a_i(\eta) = p_i(\eta) - q_i(\eta), \quad p_i, q_i \ge 0,\ \text{非减}

$V_i$ 自动随之分解。

V_i(\eta) = V_{1,i}(\eta) - V_{2,i}(\eta), \quad V_{j,i}(\eta) = \int_{-\infty}^\eta p_i\ \text{or}\ q_i

$V_1, V_2$ 二者均非减且有界。wetting/drying 例中 $p(\eta)=1\,[\eta \ge 0]$ , $q\equiv 0$ ;confined–unconfined 例中 $p(\eta)=1\,[\eta \ge -h]$ , $q(\eta)=1\,[\eta > c]$ 。一旦分开,每一段都单调——这一句话是后续所有内容的起点。

Q3. 一内一外,两次线性化#

我们求解

V_1(\eta) - V_2(\eta) + T\eta = \mathbf{b}

外层(outer)迭代将 $V_2$ 在 $\eta^{n-1}$ 处线性化。

V_1(\eta^n) + (T - Q^{n-1})\eta^n = \mathbf{b} + V_2^{n-1} - Q^{n-1}\eta^{n-1}

其中 $Q^{n-1} = \mathrm{diag}(q_i(\eta_i^{n-1}))$ 。起点为 $\eta^0 \le \ell$ (下界栏)。这一式仍是非线性的,故内层(inner)迭代将 $V_1$ 线性化:

(T + P^{n,m-1} - Q^{n-1})\eta^{n,m} = P^{n,m-1}\eta^{n,m-1} - V_1^{n,m-1} + \mathbf{d}^{n-1}

起点为 $\eta^{n,0} \ge u$ (上界栏)。设计 $\ell, u$ 使 $T+P-Q$ 成为 Stieltjes 矩阵(对称 M-矩阵)是关键。回报是:

内层从上方单调递减( $\eta^{n,m+1} \le \eta^{n,m}$ )
外层自下方单调递增( $\eta^{n+1} \ge \eta^n$ )
二者夹住唯一解

收敛性由不等式证明,而非浮点逻辑。Newton "需接近真解" 的弱点彻底消失。在下方亲手拨动:

T = 0.50b = 0.90

step 1/4 · outer n=0 · eta=-0.4000 · |r|=1.10e+0

Pink = outer iterate (rises from below); green = inner iterate (falls from above). They squeeze toward the true root of V(eta)+T*eta=b with V(eta)=max(0,min(1,eta)).

把 T 降到 0.05, $T+P-Q$ 接近奇异,但迭代仍在 5 步内结束。把 b 推过 1.5,根落在 $V$ 的平坦尾部(Newton 会永久振荡的地方)——内层一步停住。

Q4. 在 1D 井里放一台泵 — Python 实现#

把论文 §6 的二维 confined–unconfined aquifer 降到一维,从头写一遍。试验问题是抛物面底/顶+中心泵:随天数推进,phreatic 区域先扩张,继而干区出现。CFD 经典玩具问题(Sod、Karman、blast)在此不适用,直接采用论文本身的设置。

import numpy as np
 
# 每胞体积: V_i(eta) = a0 * clamp(eta, -h_i, c_i) (按 -h_i 平移)
def cell_volume(eta, h, c, a0=0.3):
    return a0 * (np.clip(eta, -h, c) - (-h))
 
# Jordan 分解给出的对角 Jacobian
def p_diag(eta, h, a0=0.3):
    # V1 = a0 * max(0, eta - (-h)) 的导数
    return np.where(eta >= -h, a0, 0.0)
 
def q_diag(eta, c, a0=0.3):
    # V2 = a0 * max(0, eta - c) 的导数
    return np.where(eta > c, a0, 0.0)
 
def nested_newton_step(eta_prev, h, c, T, dt, phi, tol=1e-6):
    """单时间步的单调嵌套 Newton 求解。
 
    求解 V(eta_new) + T eta_new = V(eta_prev) + dt*phi,
    其中 T 来源于 -dt*Laplacian, 对称正定。
    """
    N = len(eta_prev)
    V_prev = cell_volume(eta_prev, h, c)
    rhs0 = V_prev + dt * phi
 
    # 外层初值: 下界栏之下 (所有胞干)
    eta = np.full(N, -h.max() - 1.0)
    for n_out in range(20):
        Q = np.diag(q_diag(eta, c))
        V2 = np.where(eta > c, 0.3 * (eta - c), 0.0)
        d = rhs0 + V2 - Q @ eta
 
        # 内层初值: 上界栏之上 (所有胞加压)
        eta_in = np.full(N, c.max() + 1.0)
        for m_in in range(30):
            P = np.diag(p_diag(eta_in, h))
            V1 = np.where(eta_in > -h, 0.3 * (eta_in - (-h)), 0.0)
            f = P @ eta_in - V1 + d
            A = T + P - Q                     # Stieltjes (对称 M-矩阵)
            eta_in = np.linalg.solve(A, f)    # 一维: 稠密; 二维用 PCG
            r_in = V1 + (T - Q) @ eta_in - d
            if np.max(np.abs(r_in)) < tol:
                break
        eta = eta_in
        V_cur = cell_volume(eta, h, c)
        r_out = V_cur + T @ eta - rhs0
        if np.max(np.abs(r_out)) < tol:
            return eta, n_out + 1, m_in + 1
    return eta, 20, m_in + 1
 
# 一维井 (抛物面剖面)
N = 41
L = 1000.0
x = np.linspace(0, L, N)
h_arr = 10 * np.maximum(0, 1 - ((x - L/2) / (L/2))**2)
c_arr = h_arr.copy()
 
# 扩散矩阵 T = dt*kappa/dx^2 * Laplacian (零通量边界)
dx = L / (N - 1)
dt = 86400.0     # 1 日
kappa = 1.0
alpha = dt * kappa / dx**2
T = -alpha * np.eye(N, k=-1) + 2*alpha*np.eye(N) - alpha*np.eye(N, k=1)
T[0, 0] = alpha; T[0, 1] = -alpha
T[-1, -1] = alpha; T[-1, -2] = -alpha
 
# 初始: 顶部全域加压
eta = c_arr.copy()
phi = np.zeros(N); phi[N//2] = -8e-4   # 中心泵 (汇)
 
for day in range(10):
    eta, n_out, m_in = nested_newton_step(eta, h_arr, c_arr, T, dt, phi)
    pre = int(np.sum(eta >= c_arr - 1e-6))
    dry = int(np.sum(eta <= -h_arr + 1e-6))
    print(f"day {day+1}: pre={pre} dry={dry} outer={n_out} inner_last={m_in}")

典型运行结果:

day 1: pre=39 dry=0 outer=5 inner_last=5
day 2: pre=37 dry=0 outer=4 inner_last=4
...
day 9: pre=0 dry=8 outer=1 inner_last=5
day 10: pre=0 dry=12 outer=1 inner_last=5

**外层平均 3 次,内层约 5 次。**与论文 Table 1 在定性上一致。最有意思的是加压向 phreatic 过渡期——内外两层常常仅一次便停(论文 Remark 6:若 $u \le \tilde\eta^{1,1} \le \ell$ ,根已经被命中)。

Q5. 把二维剖面亲自演化#

同一算法,在一维井剖面上按日推进。泵强为滑条,按 "+1 day" 推进一日。

pump flux phi = 0.0008 (m/s per cell)

day = 0 · pressurized cells = 58 · phreatic = 0 · dry = 2 · total water V = 3998.9 m

Cyan = pressurized (confined) region; teal = phreatic; black = dry. Yellow line is piezometric head eta(x). Pump near the center steadily drains the aquifer; the wet/dry boundary advances inward as days pass.

把泵推到 0.0025 以上跑十日,中心会出现一片向外扩张的黑色干区。求解器并不关心湿胞集合每日都在改变。若是单一 Newton,湿/干模板一旦变就要重建 Jacobian。嵌套循环的代码原样不动,继续解算。

不能做到什么 — 一段批判#

并非万能。几点警示:

$P-Q = 0$ 的胞。干胞处 $T$ 对应行也归零时,矩阵会失去 T2(不可约)条件。论文建议把干胞从系统中剔除,但有效胞集合的簿记并非小事。
并不快。每一时间步外层×内层平均 15–25 次线性求解。当带 CG 的 implicit Euler 能收敛时,后者往往更快。代价换来的是零发散事故。
不能直接套到可压缩 Euler。本方法仅限对角的 $V$ 且 $a$ 有界变差的轻度非线性。Riemann 求解器中的非线性无法这样分解。
内层求解仍要选择。每一内层系统是 SPD,但如何预条件是另一门学问。论文推荐 PCG;大规模网格上通常优先多重网格。

复现可能性评分#

数学(§2–§5):自洽。一支笔即可重走所有不等式。A.
算法盒(1, 2):伪代码完备。A.
数值实验(§6):网格、时间步、所有系数皆明示。30 行 Python 完成定性复现。A.
工程落地:OpenFOAM 缺一线 wetting/drying 模块,但 meshFlow/overFlowSolver 系可植入。海洋与河口模型(SCHISM、TELEMAC)早有类似思想。

接下来值得读的#

Brugnano & Casulli (2008/2009) — 分段线性系统的有限步终止定理,本文直接前作。
Casulli (2009) — 同作者的高分辨率 wetting/drying 算法。
Casulli & Zanolli (2010) — Richards 方程 mixed form 上的嵌套 Newton(光滑 $a(\eta)$ 情形)。