[論文レビュー] 嵐が来る前にシミュレーションを終わらせる — Casulli

対角優位の余裕が1のまま残るときだけ、波速は安定条件から外れます。

1990年6月、ハンブルクの BAW は Jade–Weser 河口の潮位13日分を合わせ直す必要がありました。三角形は15万個です。時間刻みは5分です。実時間より速く終わらないと、次の予報に計算時間が残りません。

陽に進めると自由表面波の波速 $c=\sqrt{gh}$ が CFL を握ります。水深28 m で $c \approx 16.6$ m/s です。最短辺が $\sqrt{10}$ m なら許される刻みは0.2秒です。13日は1,123,200秒です。0.2秒で切ると5,616,000回で、論文の560万ステップと一致します。5分なら3,744回です。比はちょうど1500倍です。300秒を0.2秒で割った値なので、期間を替えても倍率は動きません。

その倍率を持たせるのが Casulli–Walters(2000) の半陰解法です。自由表面の勾配だけを新しい時刻に置き、移流は古い時刻に残します。

CFL スライダーを1より上げると赤い陽解法の曲線が壊れ、θ を上げるとシアンの曲線の振幅だけが落ちます。

CFL ratio (Δt / Δt_explicit): 2.00θ (implicitness): 0.60

CFL > 1 ⇒ explicit blows up; semi-implicit stays bounded for any CFL when θ ≥ 1/2.

CFL=6 でもシアンの峰は残ります。θ を1まで押すとその峰はすぐ消え、振幅が残る実務の帯は0.55から0.6です。

論文メタデータ

著者: Vincenzo Casulli (Univ. of Trento), Roy A. Walters (USGS)
誌名: International Journal for Numerical Methods in Fluids, 32, 331–348 (2000)
題目: "An unstructured grid, three-dimensional model based on the shallow water equations"
キーワード: semi-implicit, shallow water, unstructured orthogonal grid, wetting/drying

なぜ陽解法は河口で行き詰まるのか

水深10 m なら $c \approx 10$ m/s、100 m なら $c \approx 31$ m/s です。流速 $u$ はたいてい1 m/s を下回ります。潮が区間を渡る時間は $L/u$ なのに、陽的な刻みは $\Delta x/c$ に縛られます。 $c/u \approx 10\sim30$ なら、水が一度過ぎるあいだにソルバは時間をそれより細かく切ります。

非構造格子はその制限を最短辺へ寄せます。Big Lost River は1.2 km × 1 km に三角形12,622個で、最短辺は約5 mです。陽的な限界は0.5秒ほどです。5 m を0.5秒で割ると10 m/s です。深い澪の16 m/s が時計を止めているのではありません。浅い場所の短い辺が止めます。

遅い項は陽で、速い項は陰です。自由表面の圧力勾配だけを陰に回せば、 $\sqrt{gh}$ は安定条件から外れます。残る制限は流速の側です。

直交非構造格子 — Voronoi–Delaunay の対#

非構造格子に有限差分をそのまま載せると、二箇所が壊れます。セル間隔が形ごとに変わり、勾配が片側へ寄ります。曲線座標で避けると、変換項が精度と安定性を一緒に削ります。

Casulli と Walters は格子の側を制限しました。直交非構造格子です。隣り合う中心を結ぶ線分が、共有辺と垂直です。Voronoi 多角形とその双対である Delaunay 三角分割は、鋭角三角形だけのときこの条件を満たします。鈍角が一つ入ると外心は三角形の外へ出ます。中心を結ぶ線は辺と垂直のままでも、セルの内部を通りません。外接円を空に保つ点挿入が、その鋭角をどれだけ長く保つかが、この差分が成り立つ時間です。

直交が渡すものは差分一行です。辺の法線速度 $u_j$ だけを置けば、 $\partial \eta / \partial n$ は $(\eta_{i_2} - \eta_{i_1})/\delta_j$ になります。曲率の補正はありません。

半陰解法分割: θ 法で 2 つの時間スケールを切り離す#

θ 法は陰的さの係数 $\theta \in [1/2, 1]$ で時間を分けます。運動量では自由表面勾配を $\theta \eta^{n+1} + (1-\theta)\eta^n$ とし、自由表面式では速度を $\theta u^{n+1} + (1-\theta)u^n$ とします。

離散化した運動量です。鉛直粘性と風応力が入っています。

\mathbf{A}_j^n \mathbf{U}_j^{n+1} = \mathbf{G}_j^n - \theta g \frac{\Delta t}{\delta_j}\!\left[\eta_{i(j,2)}^{n+1} - \eta_{i(j,1)}^{n+1}\right]\Delta\mathbf{Z}_j^n

$\mathbf{U}_j^{n+1}$ は辺 $j$ の法線速度、 $\mathbf{A}_j^n$ は鉛直粘性・風・底面摩擦を集めた三重対角、 $\mathbf{G}_j^n$ は移流・コリオリ・水平摩擦です。辺ごとに独立なので、 $N_s$ 個の $N_z \times N_z$ 三重対角に分かれます。

自由表面式は次のとおりです。

P_i \eta_i^{n+1} = P_i \eta_i^n - \theta \Delta t \sum_{l=1}^{S_i} s_{i,l}\lambda_{j(i,l)} [\Delta\mathbf{Z}^\top \mathbf{U}]_{j(i,l)}^{n+1} - (1-\theta)\Delta t \sum_{l=1}^{S_i} s_{i,l}\lambda_{j(i,l)} [\Delta\mathbf{Z}^\top \mathbf{U}]_{j(i,l)}^{n}

$P_i$ はセル面積、 $\lambda_j$ は辺の長さ、 $s_{i,l}$ は外向きの符号です。

波動方程式への帰着と SPD の恩恵#

二式をそのまま置くと未知量は $N_z N_s + N_p$ 個です。運動量から $\mathbf{U}_j^{n+1}$ を $\eta^{n+1}$ の関数として解き、自由表面へ代入します。 $\mathbf{A}$ が正定値なら逆行列も正定値です。

代入の結果が離散波動方程式です。

P_i \eta_i^{n+1} - g \Delta t^2 \theta^2 \sum_{l=1}^{S_i} \frac{s_{i,l}\lambda_{j(i,l)}}{\delta_{j(i,l)}}\!\left[(\Delta\mathbf{Z})^\top \mathbf{A}^{-1}\Delta\mathbf{Z}\right]_{j(i,l)}^{n}\!\!\left(\eta^{n+1}_{i[j(i,l),2]} - \eta^{n+1}_{i[j(i,l),1]}\right) = R_i

未知量は $\eta^{n+1}$ だけです。大きさはセル数 $N_p$ です。行列は対称で、対角が強く、正定値です。PCG が合います。論文は、大きな格子でも反復が数十回で止まると書いています。

鉛直構造のない1Dでは $\mathbf{A}$ が1なので、コードの $k = g H \Delta t^2 \theta^2 / \Delta x^2$ が隣接係数になります。 $\Delta t$ が陽的限界の4倍で $c^2 = gH$ なら $k = 16\theta^2$ です。θ=0.6 なら $k=5.76$ です。内部の対角 $1+2k$ は12.52、非対角二つの和は11.52です。対角が勝つ幅はいつも1で、その1は貯蔵項 $P_i \eta$ が残したものです。刻みを大きくすると $k$ は対角と非対角の両方に付き、余剰は増えません。先のスライダーで CFL=6、θ=0.6 なら $k=12.96$ 、対角26.92、非対角の和25.92です。幅はやはり1です。より陰にしても行列はそれ以上安全になりません。死ぬのは振幅です。

$\eta$ が解ければ $\mathbf{U}_j^{n+1}$ は辺ごとの小さな三重対角に落ちます。鉛直速度は非圧縮の式から遡ります。

Python で 1D 半陰解法#

同じ分離で、 $N=200$ 、 $L=1$ 、 $H=1$ 、 $g=9.81$ 、θ=0.6 です。

import numpy as np
from scipy.linalg import solve_banded
 
class CasulliShallowWater1D:
    """1D 線形浅水方程式 — 半陰解法 θ 法"""
 
    def __init__(self, N=200, L=1.0, H=1.0, g=9.81, theta=0.6):
        self.N, self.L, self.H, self.g, self.theta = N, L, H, g, theta
        self.dx = L / N
        self.eta = np.zeros(N)            # セル中心の自由表面
        self.u = np.zeros(N + 1)          # 辺(face)の法線速度、両端 = 0
 
    def initial_bump(self, x0=0.3, sigma=0.05, amp=0.05):
        x = (np.arange(self.N) + 0.5) * self.dx
        self.eta = amp * np.exp(-((x - x0) ** 2) / (2 * sigma ** 2))
        self.u[:] = 0.0
 
    def step(self, dt):
        N, dx, g, H, th = self.N, self.dx, self.g, self.H, self.theta
        k = g * H * dt * dt * th * th / (dx * dx)   # 波動項係数
        alpha = H * dt / dx
 
        # G_j: 陽解部分を辺で集めた値(両端 = 0)
        G = np.zeros(N + 1)
        G[1:N] = self.u[1:N] - g * (dt / dx) * (1 - th) * (self.eta[1:] - self.eta[:-1])
 
        # 三重対角: -k η_{i-1} + (1+2k) η_i - k η_{i+1} = R_i
        ab = np.zeros((3, N))
        ab[0, 1:] = -k                         # 上対角
        ab[2, :-1] = -k                        # 下対角
        ab[1, :] = 1.0 + 2 * k
        ab[1, 0] = 1.0 + k                     # 左壁
        ab[1, -1] = 1.0 + k                    # 右壁
 
        rhs = self.eta.copy()
        rhs -= alpha * th * (G[1:] - G[:-1])
        rhs -= alpha * (1 - th) * (self.u[1:] - self.u[:-1])
 
        new_eta = solve_banded((1, 1), ab, rhs)
 
        new_u = self.u.copy()
        new_u[1:N] = self.u[1:N] - g * (dt / dx) * (
            th * (new_eta[1:] - new_eta[:-1]) + (1 - th) * (self.eta[1:] - self.eta[:-1])
        )
        new_u[0] = new_u[-1] = 0.0
        self.eta, self.u = new_eta, new_u
        return self.eta
 
if __name__ == "__main__":
    sim = CasulliShallowWater1D(N=200, theta=0.6)
    sim.initial_bump()
    c = np.sqrt(sim.g * sim.H)
    dt_explicit_limit = sim.dx / c
    dt = 4.0 * dt_explicit_limit                # 陽解法限界の 4 倍
    for _ in range(int(0.6 / dt)):
        sim.step(dt)
    print(f"max|η| = {np.max(np.abs(sim.eta)):.4f}  (安定維持)")

陽解法は $\Delta t = 4 \Delta t_\text{ex}$ ですぐ発散します。Casulli は同じ間隔で長さ0.6の区間を終えます。位相は遅れ、振幅は θ が1に近いほど余計に死にます。θ=0.5 は重みが半々なので2次です。0.6を超えると1次です。コードの初期値0.6は、すでにその1次側に足を置いています。

時間刻みを大きくしてみよう

水路長さ1、セル200なら $\Delta x = 0.005$ です。 $c=\sqrt{9.81}\approx 3.13$ m/s なので陽的限界は約0.0016秒、その4倍は約0.0064秒です。0.6をこの間隔で割って小数を捨てると93ステップになります。出力は max|η| が有限だという一行だけで、初期の隆起0.05がどれだけ減ったかは書きません。

その0.05は水深 $H=1$ の5%です。線形化 $\partial_t \eta + H\partial_x u = 0$ は、 $\eta$ が $H$ より十分小さいときに閉じます。だから波速は $\sqrt{gH}$ の一つの値で、三重対角の $k$ もセルごとに同じです。隆起を水深まで大きくすると $c$ は場所で変わり、この一定の $k$ では波動方程式が閉じません。

θ=0.5 に同じ4倍の刻みを入れると $k=16\times 0.25=4$ です。対角優位はやはり1です。安定限界という言葉は行列からではなく、時間の重み $(θ-1/2)$ から来ます。その値が0だと数値粘性は切れ、格子モードが減衰なしに往復します。

批判的考察

万能ではありません。

時間次数から値が違います。θ=0.5 は2次ですが安定境界にあり、重力波の波面が簡単に曲がります。0.6を超えると1次です。検潮記録の振幅が成果物なら、刻みを縮めるか、より高次の IMEX へ移します。

移流は陽的な項に残ります。単純な風上は数値粘性が強いです。Lagrange 軌跡を逆にたどって補間する EL を足す理由はそこにあります。格子一辺を越える $\Delta t$ でも安定しますが、補間次数と出発点の誤差で結果が揺れます。

静水圧 $\partial p/\partial z = -\rho g$ も前提です。砕波や落差のように鉛直加速度が大きい場所には、非静水圧の補正が要ります。

鈍角が一つ混ざると5分刻みは使えません#

5分と θ=0.6 の取引は、ステップ数を1500分の1にし、波高の1次減衰を受け入れることです。検潮記録の振幅が成果物なら、θ を1の側へ押しません。対角が勝つ幅はすでに1なので、より陰にしても行列はそれ以上安全になりません。死ぬのは波です。

この差分が $\partial \eta / \partial n$ になるのは、中心を結ぶ線が辺に垂直なときだけです。鈍角が一つ外心をセルの外へ押し出すと、その等式は終わります。辺に垂直でない成分が残ります。その成分を拡散フラックスと同じ形で戻せるかどうかは、非直交な面の拡散補正が問うている点です。

[論文レビュー] 嵐が来る前にシミュレーションを終わらせる — Casulli–Walters(2000) 非構造直交格子の浅水方程式