[論文レビュー] 同じ点の色関数が経路によって0.544ずれた

カール拘束条件は初期条件ではなく毎ステップ守るべき恒等式であり、離散gradientとcurlが同じ4点を共有すれば丸め誤差の水準で自動的に保たれます。

同じ点で色関数が二つの値に読めた

界面場 $\mathbf{b}$ を400ステップ動かした後、格子中央のセルで色関数を復元してみました。 $\mathbf{b}$ を原点から積分すれば色関数が戻るはずです。定義が $\mathbf{b} = \nabla c$ だからです。 $x$ に先に進んで $y$ に曲がった値は一つの数でした。 $y$ に先に進んで $x$ に曲がった値は別の数でした。差は0.544です。色関数は0と1の間にある量です。

せん断の強さと Rusanov 散逸のスライダーを一緒に上げると、丸い界面は片方の格子だけで色を持ち始めます。

shear amplitude 1.00Rusanov dissipation 1.00

t = 0.000   max|curl b|: 0.00e+0 vs 0.00e+0   max|b|: 0.00 vs 0.00

Both panels start from the same interface field b = grad(phi) of one blob, and both are pushed by the same tanh shear layer. The left panel stays black at round-off level however hard the shear is driven; the right panel lights up within a fraction of a turnover. Raising the Rusanov dissipation slider makes the right panel worse, not better — it is the dissipation term itself that is not the gradient of any scalar.

左がスタッガード配置、右がコロケート配置です。左が黒いままのあいだ、右がいつ明るくなるかを見ると、経路差 0.544 がどちらに積もるかが分かります。

この記事は、粘性と表面張力を含む1階双曲型二相流モデルに対する Chiocchetti と Dumbser (2023) のスタッガード半陰的スキームを読み、そのうちカール拘束条件(curl involution)の部分だけを取り出して再現した記録です。なぜこの値がゼロでなければならないのか、なぜ通常の配置ではゼロにならないのか、そして何を変えれば丸め誤差まで落ちるのかを数値で確かめます。

界面場を独立変数に昇格させる理由

拡散界面(diffuse interface)モデルは、二相の境界を色関数 $c$ のなだらかな遷移で表します。表面張力の計算には $\nabla c$ とその発散が要ります。問題は微分の回数です。曲率は $c$ の2階微分なので、界面が数セル幅しかなければノイズがそのまま増幅されます。

1階双曲型(first order hyperbolic)への再定式化は、ここで微分を一段下げます。 $\mathbf{b} = \nabla c$ を独立変数に昇格させ、 $\mathbf{b}$ 自身の輸送方程式を解きます。 $c$ の輸送方程式 $\partial_t c + \mathbf{v}\cdot\nabla c = 0$ の両辺に $\nabla$ をかけるとこうなります。

\partial_t b_k + \partial_k (v_j b_j) = 0

$b_k$ は界面場の $k$ 成分、 $v_j$ は速度、 $v_j b_j = \mathbf{v}\cdot\mathbf{b}$ はスカラー一つです。これでシステム全体が1階微分だけを含みます。固有値が定義され、双曲型ソルバーがそのまま使えます。拘束が破れると、何も起きてはいけない圧力まで一緒に跳ねます。非保存項ひとつで圧力が21%上がった場所が同じ種類の帳簿です。

代償が一つあります。 $\mathbf{b}$ は $c$ とは独立に保存されます。二つはずれ得ます。

拘束条件は初期条件ではなく方程式の一部である

$\mathbf{b}$ が何らかのスカラーのgradientなら、回転は常にゼロです。

\nabla \times \mathbf{b} = \nabla \times \nabla c = \mathbf{0}

この関係を拘束条件(involution)と呼びます。輸送方程式自身がこの性質を保存します。 $\partial_t (\nabla \times \mathbf{b}) = -\nabla \times \nabla(v_j b_j) = \mathbf{0}$ だからです。連続方程式の水準では、 $t=0$ でゼロなら永遠にゼロです。

そのため初期条件の問題として扱われがちです。これがまさに落とし穴です。離散化された方程式では $\nabla_h \times \nabla_h$ がゼロである保証がありません。保証がなければ毎ステップ少しずつ漏れ、400ステップ後には冒頭の0.544になります。

カールがゼロでないと何が壊れるのか。Stokesの定理が答えます。閉じたループに沿った積分は、そのループが囲むカールの面積分に等しい。カールがあれば $\mathbf{b}$ はどの色関数のgradientでもありません。つまり保存された $\mathbf{b}$ に対応する界面がこの世に存在しません。

離散作用素がSchwarz恒等式を共有する場所#

論文の処方は単純です。離散gradientと離散curlが同じ4点を使うように配置することです。

色関数をセル中心に置き、 $b_x$ は $x$ 面 $(i+\frac12, j)$ に、 $b_y$ は $y$ 面 $(i, j+\frac12)$ に置きます。カールはノード $(i+\frac12, j+\frac12)$ で計算します。するとスカラー $\psi$ についてこう展開されます。

(\nabla_h \times \nabla_h \psi)_{i+\frac12,\, j+\frac12} = \frac{\psi_{i+1,j+1} - \psi_{i,j+1} - \psi_{i+1,j} + \psi_{i,j}}{h^2} - \frac{\psi_{i+1,j+1} - \psi_{i+1,j} - \psi_{i,j+1} + \psi_{i,j}}{h^2}

二つの分子は同じ4項を順番だけ変えて書いたものです。差は恒等的にゼロです。これが交差微分のSchwarz恒等式(mixed partials commute)を離散レベルにそのまま移した形です。近似ではありません。浮動小数点でも同じ4つの数を足し引きするので、丸め誤差までしか残りません。

この性質が時間積分に移る論理は一行です。毎ステップの更新量が $-\Delta t\, \nabla_h \psi$ の形なら、カールをかけた瞬間に消えます。つまり $\nabla_h \times \mathbf{b}$ はステップ数と無関係に初期値を保ちます。初期値がゼロなら最後までゼロです。

Pythonで同じせん断層に二つの離散化をかけた#

$64^2$ の周期格子に丸い界面を一つ置き、tanhのせん断層で押しました。片方は上のスタッガード作用素の組、もう片方は $\mathbf{b}$ をすべてセル中心に置いて成分ごとの Rusanovフラックスを使う配置です。どちらも同じ $\Delta t$ 、同じ速度場です。

import numpy as np
 
N, L = 64, 1.0
h = L / N
xc = (np.arange(N) + 0.5) * h
Xc, Yc = np.meshgrid(xc, xc, indexing="ij")
 
def shear_velocity(amp=1.0, delta=0.12):
    """Kelvin-Helmholtz型: tanhせん断層 + 正弦摂動"""
    return amp * np.tanh((Yc - 0.5) / delta), 0.35 * amp * np.sin(2 * np.pi * Xc)
 
def interface_potential(r0=0.22, w=0.05):
    r = np.sqrt((Xc - 0.5) ** 2 + (Yc - 0.5) ** 2)
    return 0.5 * (1.0 + np.tanh((r0 - r) / w))
 
def grad_stagger(phi):
    """b_x は x面 (i+1/2, j)、b_y は y面 (i, j+1/2)"""
    return (np.roll(phi, -1, 0) - phi) / h, (np.roll(phi, -1, 1) - phi) / h
 
def curl_node(bx, by):
    """ノード (i+1/2, j+1/2) のカール — grad_stagger と対になる作用素"""
    return (np.roll(by, -1, 0) - by) / h - (np.roll(bx, -1, 1) - bx) / h
 
def curl_centered(bx, by):
    """コロケート配置のカール: セル中心成分の中心差分"""
    return ((np.roll(by, -1, 0) - np.roll(by, 1, 0))
            - (np.roll(bx, -1, 1) - np.roll(bx, 1, 1))) / (2 * h)
 
def psi_cell(bx, by, vx, vy, upwind):
    """スカラー一つ psi = v . b をセル中心に集める"""
    if upwind:
        bxc = np.where(vx > 0, np.roll(bx, 1, 0), bx)
        byc = np.where(vy > 0, np.roll(by, 1, 1), by)
    else:
        bxc = 0.5 * (bx + np.roll(bx, 1, 0))
        byc = 0.5 * (by + np.roll(by, 1, 1))
    return vx * bxc + vy * byc
 
def update_compatible(bx, by, vx, vy, dt, upwind=True):
    """b <- b - dt * grad_h(psi): 加わる項がすべて離散gradient"""
    psi = psi_cell(bx, by, vx, vy, upwind)
    return (bx - dt * (np.roll(psi, -1, 0) - psi) / h,
            by - dt * (np.roll(psi, -1, 1) - psi) / h)
 
def update_collocated(bx, by, vx, vy, dt, lam):
    """成分ごとのRusanov: 散逸項はどのスカラーのgradientでもない"""
    psi = vx * bx + vy * by
    fx = 0.5 * (psi + np.roll(psi, -1, 0)) - 0.5 * lam * (np.roll(bx, -1, 0) - bx)
    fy = 0.5 * (psi + np.roll(psi, -1, 1)) - 0.5 * lam * (np.roll(by, -1, 1) - by)
    return (bx - dt * (fx - np.roll(fx, 1, 0)) / h,
            by - dt * (fy - np.roll(fy, 1, 1)) / h)
 
def path_gap(bx, by):
    """同じセルまで二つの経路で b を積分する。差が囲んだカール"""
    i1 = j1 = N // 2
    pa = bx[:i1, 0].sum() * h + by[i1, :j1].sum() * h   # x が先、次に y
    pb = by[0, :j1].sum() * h + bx[:i1, j1].sum() * h   # y が先、次に x
    return pa - pb
 
vx, vy = shear_velocity()
phi0 = interface_potential()
lam = float(np.max(np.abs(vx) + np.abs(vy)))
dt = 0.4 * h / lam
bs = list(grad_stagger(phi0))
bc = [np.gradient(phi0, h, axis=0, edge_order=2),
      np.gradient(phi0, h, axis=1, edge_order=2)]
 
print(f"N={N}  h={h:.4e}  dt={dt:.4e}  lambda={lam:.3f}")
print(f"{'t':>6}{'|b|max stag':>13}{'|curl|max stag':>16}"
      f"{'|b|max coll':>13}{'|curl|max coll':>16}")
for k in range(5):
    if k:
        for _ in range(100):
            bs = update_compatible(*bs, vx, vy, dt)
            bc = update_collocated(*bc, vx, vy, dt, lam)
    print(f"{k*100*dt:>6.3f}{np.max(np.hypot(*bs)):>13.3f}"
          f"{np.max(np.abs(curl_node(*bs))):>16.3e}"
          f"{np.max(np.hypot(*bc)):>13.3f}"
          f"{np.max(np.abs(curl_centered(*bc))):>16.3e}")
 
print(f"\npath-integral gap of phi at (0.5, 0.5):"
      f"  staggered {path_gap(*bs):+.3e}   collocated {path_gap(*bc):+.3e}")
 
bn = list(grad_stagger(phi0))
for _ in range(400):
    bn = update_compatible(*bn, vx, vy, dt, upwind=False)
print(f"staggered with NO upwinding: |b|max={np.max(np.hypot(*bn)):.3e} (unstable)"
      f"  |curl|max={np.max(np.abs(curl_node(*bn))):.3e}")

N=64  h=1.5625e-02  dt=4.6330e-03  lambda=1.349
     t  |b|max stag  |curl|max stag  |b|max coll  |curl|max coll
 0.000       10.167       0.000e+00        9.818       3.286e-04
 0.463       12.318       1.421e-12        6.538       5.545e+01
 0.927       18.063       1.748e-12        4.703       1.687e+02
 1.390        9.854       1.705e-12        4.124       2.269e+02
 1.853        5.746       2.132e-12        3.342       1.998e+02
 
path-integral gap of phi at (0.5, 0.5):  staggered -3.331e-16   collocated -5.442e-01
staggered with NO upwinding: |b|max=3.622e+03 (unstable)  |curl|max=1.091e-10

$\kappa$ スライダーを上げると二つの経路のバーが分かれ、色のついた面のカール積分がその差といっしょに増えます。

curl content κ 0.45loop size 0.52

A = 0.0000   B = 0.0000   gap = 0.000e+0

Two walkers start at the same corner and end at the same corner, one going x first and the other y first, each summing b along the way. Drag the curl content to zero and the two bars land on the identical number — the field has a potential, so the endpoint value is the potential and nothing else. Any nonzero curl splits the bars, and the gap always equals the shaded curl integral enclosed by the loop.

ゼロまで下げると二本のバーは同じ数で止まります。少し上げるだけで経路差がその面積分と揃います。

スタッガード側のカールは $10^{-12}$ に留まります。 $1/h^2 \approx 4096$ が掛かった倍精度の丸め誤差がその大きさです。コロケート側は時刻0.46で55、1.39で227まで上がります。 $|\mathbf{b}|$ の最大値が4程度なので、カールは場そのものより50倍以上大きい。残っているのは界面ではなく回転ノイズです。

$t=1.390$ のコロケート側でカールの最大は 226.9 です。同じ行の $|\mathbf{b}|$ 最大は 4.124 です。割ると 50 を超えます。場より回転の方が大きい。このカールを $(0.5, 0.5)$ まで二経路で積分した値が、経路差 $-5.442\times 10^{-1}$ です。色関数は 0 と 1 の間なのに、経路差はその幅の半分を超えます。

時刻の列も表にそのまま出ています。 $dt=4.6330\times 10^{-3}$ です。100 ステップが 0.463、400 ステップが 1.853 です。出力の最後の行と一致します。その間、コロケートのカールは $3.286\times 10^{-4}$ から 55.45、168.7、226.9 へ上がり、1.853 で 199.8 になります。 $|\mathbf{b}|$ が 3.342 まで削れたので、カールの絶対値も少し下がりました。消えたのではありません。スタッガード側は同じ 400 ステップで $1.421\times 10^{-12}$ から $2.132\times 10^{-12}$ です。 $1/h^2 \approx 4096$ を倍精度の丸めに掛けた大きさと揃います。

風上を切った行は別の失敗です。 $|\mathbf{b}|$ は $3.622\times 10^{3}$ まで育ちます。初期値 10.167 の約 360 倍です。カールは $1.091\times 10^{-10}$ に残ります。表の $\lambda=1.349$ と $h=1.5625\times 10^{-2}$ の積が、更新がスカラー一つの勾配でないときの一歩あたりの漏れの大きさです。400 ステップ分が積もると経路差 0.544 になります。恒等式は振幅を縛りません。散逸は別に入れます。

散逸を増やしてもカールが増えない条件

二本目の update_compatible は風上差分(upwind)を使います。散逸が入るということです。それでもカールは $10^{-12}$ のままです。ここがこのスキームの実質的な利点です。

理由は風上差分が入る場所にあります。風上判定はスカラー $\psi = \mathbf{v}\cdot\mathbf{b}$ をセル中心に集めるときにしか使われません。その後の更新は依然として $-\Delta t\, \nabla_h \psi$ の一塊です。 $\psi$ をどう作ろうと、gradientのカールはゼロです。

コロケート配置は逆です。Rusanov散逸項 $-\frac{\lambda}{2}(b_{x,i+1} - b_{x,i})$ は $b_x$ には $x$ 方向に、 $b_y$ には $y$ 方向にだけ付きます。二つの成分が異なる作用素を受けます。この組はどのスカラーのgradientでもないので、カールをかけると $\frac{\lambda h}{2}(\partial_x \partial_y^2 b_y - \partial_y \partial_x^2 b_x)$ が残ります。

制約は精度ではなく形です。更新量全体がスカラー一つの離散gradientでありさえすれば、リミッターでも高次再構成でも載せられます。格子が表せる項を先に数えて、スキームはその後に選びます。

カールがゼロでも計算は発散する

最後の出力行がその話をします。風上差分を切って中心平均だけにすると $|\mathbf{b}|$ が $3.6\times 10^3$ まで育ちます。初期値が10だったので360倍です。計算は発散しました。ところがカールは $1.09\times 10^{-10}$ です。依然として丸め誤差の水準です。

拘束条件の保存と安定性は別の性質です。スタッガード配置は「この場は依然として何らかの界面のgradientだ」を保証するだけで、その界面が有限の大きさに留まるとは約束しません。散逸は別に入れる必要があります。二つを混同すると、発散した計算を見て格子配置を疑うという誤った診断に至ります。

表面張力でこの性質が数値を変える場所

論文が扱うシステムは $\mathbf{b}$ から表面張力応力を直接作ります。応力テンソルはおよそ $\sigma(|\mathbf{b}|\,\mathbf{I} - \mathbf{b}\otimes\mathbf{b}/|\mathbf{b}|)$ の形です。 $\mathbf{b}$ に回転成分が混じると、この応力が存在しない界面の張力を作り出します。

症状は静止気泡で最初に見えます。理論上は速度がゼロであるべき問題で、寄生流(spurious current)が育ちます。カール成分は曲率を通じて応力に入り、応力は速度を作り、速度はさらに $\mathbf{b}$ を回転させます。フィードバックループが閉じると振幅は指数的に大きくなります。

論文は同じ格子の上で圧力を半陰的に解きます。音波だけを陰的に回し、残りは陽的に置きます。スタッガード配置はここで二度働きます。カール恒等式を与え、圧力方程式を対称正定値にして共役勾配法(CG)で解けるようにします。配置一つが二つの性質を同時に買ってくれます。

カールがゼロのまま発散したら格子を疑わない

$b_x$ と $b_y$ は別の面に住むので、 $\mathbf{v}\cdot\mathbf{b}$ を作るたびに平均が入ります。コードの psi_cell がその平均です。平均は高周波を削ります。界面が3〜4セルなら、この減衰を無視して配置を責めると界面幅が先に死にます。面とノードの境界ステンシルも別です。非構造格子へ移すと両立する作用素の組を作り直す必要があり、論文が直交格子に留まる理由はそこです。

判断は二つに分かれます。界面場や変形場を独立変数として持つなら、カール恒等式は選択肢ではありません。毎ステップ漏れ、漏れる速さは $\lambda h$ に比例します。 $\mathbf{b}$ を毎回 $c$ から計算し直すモデルなら、この漏れは起きません。どちらでも path_gap の十行で、保存したベクトル場に対応するスカラーがあるかを一度は数えます。その格子に圧力を陰的に載せる費用は、時間ステップを1001倍にした半陰的圧力と同じ帳簿です。

[論文レビュー] 同じ点の色関数が経路によって0.544ずれた — カール拘束条件とスタッガード格子