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

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

把 T 降到 0.05,把 b 推过 1.5。标量嵌套 Newton 从上方和下方同时收紧。

T 等于 0.05 时 $T+P-Q$ 近乎奇异,迭代仍在 5 步内停下。b 超过 1.5 时根落在 $V$ 的平坦尾部,也就是 $\eta > 1$,内层一步即停。

论文信息

• 作者: 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 为何频频死掉#

方程看起来温和。

其中 $\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 分解)。

$V_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. 一内一外,两次线性化#

我们求解

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

其中 $Q^{n-1} = \mathrm{diag}(q_i(\eta_i^{n-1}))$。起点为 $\eta^0 \le \ell$(下界栏)。这一式仍是非线性的,故内层(inner)迭代将 $V_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 "需接近真解" 的弱点彻底消失。上方滑块上的 5 步和 1 步,就是这条不等式的次数。内层矩阵保持 Stieltjes 时,求解和 V 循环先清掉长波的地方是同一类 SPD 问题。

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

把泵强和 +1 day 往上推,上面输出里的 pre 与 dry 会在剖面上摊开。

超过 0.0025 再跑十日,中心会出现一片黑色干区。有效胞数每天都变,嵌套循环并不为此重建雅可比。

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

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

只按胞数重读这十天。第 1 天加压胞是 41 个里的 39 个,干胞 0,外层 5,最后一次内层 5。第 2 天加压降到 37,外层 4。打印跳过的日子不作估计。第 9 天加压 0,干胞 8,外层 1,内层 5。第 10 天只有干胞增到 12,外层仍是 1。泵在中心胞, $\phi=-8\times 10^{-4}$,时间步 86400,长度 1000 分成 41 个胞。$a_0=0.3$ 只在湿区间留在 $P$ 里。

水头贴着顶板 $c$ 时 $q$ 打开,$Q$ 等于 $a_0$。落到 $c$ 以下则 $Q=0$。落到井底 $-h$ 以下则 $P$ 也是 0。第 9 天外层从 5 降到 1,是贴近 Remark 6 的信号。内层仍是 5 次。干胞让 $P-Q=0$ 的行变多,内层线性求解的次数并不因此减少。

一步里外层平均 3 次、内层约 5 次,线性系就是 15 到 25 个。矩阵每次都是 Stieltjes,也就是对称 M-矩阵。把发散次数钉在 0 的代价是,不剔除活动胞的话零行会留下。第 10 天的 12 个干胞是 41 个里的一部分。剔除则集合每天变,留下则奇异行可能变多。标量上 $T=0.05$ 时迭代仍在 5 步内结束,b 超过 1.5 时内层一步停下。和水井第 9 天外层 1 次是同一种停止。

不能做到什么 — 一段批判

并非万能。几点警示:

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

复现可能性评分

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

何时把干胞从系统里拿掉

$P-Q=0$ 的行变多,T2 就破。论文建议拿掉干胞。第 10 天干胞 12 个,一共 41 个。拿掉则活动集合每天变,留下则零行还在。嵌套 Newton 两边都保住不等式收敛。变的是线性系的尺寸。

Brugnano 与 Casulli(2008/2009)的有限步终止,说的是分段线性。$a(\eta)$ 光滑的 Richards 归 Casulli 与 Zanolli(2010)。波做出来的自由水面在 水深把开尔文楔撑开的计算 里,和这口井的干胞只共享一个词。