方向を分割せずに移流を解く — Bell-Colella-Glaz の非分割風上スキーム

1周後の質量誤差はどちらも丸めなのに、ピークは0.349と0.806に分かれます。保存の合格は円板が残ったことではありません。

1980年代後半、ローレンス・バークレー研究所の John Bell と Phillip Colella は、爆発シミュレーションに奇妙な傷跡を見つけました。格子軸に対して斜めに置いた火炎面が、時間とともに対角方向へにじみ、ねじれていったのです。原因は物理ではなく離散化でした。1方向ずつ順に解く方向分割(dimensional splitting)が、対角移流をゆがめていたのです。解決策は Bell-Colella-Glaz(BCG)の非分割2次風上スキームです。方向を分割すると回転円板はひし形になり、横方向微分が1本、その穴を埋めます。

CFLを0から0.9へ上げるとオレンジ点がセル平均の側へ寄り、limiterを none にするとジャンプ横の直線が隣の値を越えます。

CFL = 0.50

Cyan = limited PLM slope inside each cell · Amber dot = predicted left edge state at n+1/2 (advecting velocity u>0). Pick none near the jump to see the slope overshoot the neighbor values.

0.9では勾配の重みがほぼ消えます。minmodに戻すと、その飛び出しは切れます。

扱う方程式は、もっとも単純な保存形の移流式です。発散ゼロ(divergence-free)の速度場 $\mathbf{u}$ がスカラー $s$ を運びます。

\frac{\partial s}{\partial t} + \nabla\cdot(\mathbf{u}\,s) = 0

ここで $s$ は濃度や温度のような受動スカラー、 $\mathbf{u}$ は既知の速度場です。有限体積法(FV)でセル平均を更新すると、次の形になります。

s_{ij}^{n+1} = s_{ij}^{n} - \frac{\Delta t}{\Delta x}\left(F^x_{i+\frac12,j} - F^x_{i-\frac12,j}\right) - \frac{\Delta t}{\Delta y}\left(F^y_{i,j+\frac12} - F^y_{i,j-\frac12}\right)

$F^x = u\,s^{n+1/2}$ は face での時間中心(time-centered)フラックスです。精度はすべて、face 上の値 $s^{n+1/2}$ をどれだけうまく見積もれるかにかかっています。

1方向ずつ分割すると回転がねじれる#

方向分割は、2D 問題を x 方向の 1D 更新、続いて y 方向の 1D 更新へと分けます。実装は簡単で、1D ソルバが1つあれば済みます。問題は対角方向の流れです。

スカラーが $45^\circ$ で流れると、x を先に解いてから y を解く順序と、その逆順とで答えが変わります。どちらも2つの 1D ステップの間に挟まる横方向結合(cross term)を取りこぼすからです。回転流ではこの誤差が毎ステップ積み重なり、円板がひし形へとゆがみます。Strang 分割で順序を入れ替えても、2次誤差の方向非対称が消えるだけで、根本原因は残ります。

BCG の処方はシンプルです。x も y も同時に解く。1ステップの中で、すべての face 値を同じ時刻 $t^{n+1/2}$ で予測し、その予測に横方向の変化をあらかじめ織り込みます。

時間を空間に置き換える予測子

中心となる考え方は、Taylor 展開で face 値を時間・空間の両方へ半ステップ押し出すことです。セル $i$ の右 face における左状態を見ます。

s^{n+1/2}_{i+\frac12,\,L} = s^n_{ij} + \frac{\Delta x}{2}\,\partial_x s + \frac{\Delta t}{2}\,\partial_t s

時間微分をそのまま残しても役に立ちません。移流方程式 $\partial_t s = -u\,\partial_x s$ を代入し、時間微分を空間微分へ置き換えます(Cauchy-Kowalewski のトリック)。

s^{n+1/2}_{i+\frac12,\,L} = s^n_{ij} + \frac{1}{2}\left(1 - \frac{\Delta t}{\Delta x}\,u_{i+\frac12,j}\right)\Delta_x s_{ij}

$\Delta_x s_{ij}$ はセル $i$ の制限付き勾配(limited slope)、 $u_{i+1/2,j}$ は face の移流速度です。括弧内の $(1 - \mathrm{CFL})$ が時間前進を担います。CFL が 1 に近いほど左状態の寄与は小さくなり、0 なら単なる空間外挿になります。

隣のセルから見た右状態も対称に作ります。

s^{n+1/2}_{i+\frac12,\,R} = s^n_{i+1,j} - \frac{1}{2}\left(1 + \frac{\Delta t}{\Delta x}\,u_{i+\frac12,j}\right)\Delta_x s_{i+1,j}

2つの候補から1つを選ぶのが Godunov の風上規則です。face 速度の符号が風を決めます。

s^{n+1/2}_{i+\frac12} = \begin{cases} s_{i+\frac12,\,L}, & u_{i+\frac12} > 0 \\[2pt] s_{i+\frac12,\,R}, & u_{i+\frac12} < 0 \end{cases}

上の図でCFLを0.9にするとオレンジ点が平均へ寄り、none では直線が隣の値を越えます。その差が、直前の予測子にある時間前進の係数です。

横方向微分が角をつなぐ

方向分割が取りこぼす結合を、BCG は予測子の中へ直接組み込みます。x-face の左状態から、横方向フラックス発散を半ステップぶん差し引きます。

s^{n+1/2}_{i+\frac12,\,L} \;\mathrel{-}=\; \frac{\Delta t}{2}\,\frac{\partial (v\,s)}{\partial y}\bigg|_{ij}

$v$ は y 方向速度、 $\partial_y(v s)$ はセル $ij$ における横方向移流の変化です。この項を見積もるには、まず y-face 上の 1D 予測状態 $\hat s^y$ を求めておく必要があります。つまり BCG は2段階です。まず各方向の 1D 風上状態を予測し、その予測を横方向補正に再利用してから、同じ風上規則で最終の face 値を選びます。この「角の結合(corner coupling)」が、対角移流をねじれなく運びます。

勾配に手綱をかける — minmod と superbee#

予測子の精度は勾配 $\Delta_x s$ から生まれます。中心差分の勾配をそのまま使うと2次精度ですが、不連続付近で振動(overshoot)を起こします。そこで勾配制限子(slope limiter)で手綱をかけます。もっとも保守的な minmod は、2つの片側差分のうち小さい方を、符号が異なれば 0 を選びます。

\Delta_x s_{ij} = \operatorname{minmod}\!\left(s_{ij}-s_{i-1,j},\; s_{i+1,j}-s_{ij}\right)

superbee は同じ単調性を保ちながら勾配を最大限に立て、接触面をより鋭く引きます。上の 1D デモで minmod と superbee を切り替えると、superbee 側の PLM 直線がより急なのが確認できます。

同じ手綱をフラックスにかける1次元の計算は、Lax–WendroffにTVDリミッタを載せる方法です。BCGはリミッタをフラックスではなくセル勾配に掛けます。

numpy で回す回転円板#

ここで 2D 非分割 BCG を回します。トイ問題は剛体回転(solid-body rotation)です。中心を回る発散ゼロの速度場にスカラー円板を浮かべ、1次風上と BCG を同じ格子で比較します。

import numpy as np
 
N = 64
dx = 1.0 / N
omega = 2 * np.pi                      # 単位時間あたり1回転
 
def face_velocities():
    j = np.arange(N)
    i = np.arange(N)
    u = -omega * ((j + 0.5) * dx - 0.5)        # x-face 速度 (y のみに依存)
    u = np.broadcast_to(u[None, :], (N, N)).copy()
    v = omega * ((i + 0.5) * dx - 0.5)         # y-face 速度 (x のみに依存)
    v = np.broadcast_to(v[:, None], (N, N)).copy()
    return u, v
 
def minmod_pair(a, b):                          # 単調性を保つ勾配
    return np.where(a * b > 0, np.where(np.abs(a) < np.abs(b), a, b), 0.0)
 
def limited_slopes(s):
    slx = minmod_pair(s - np.roll(s, 1, 0), np.roll(s, -1, 0) - s)
    sly = minmod_pair(s - np.roll(s, 1, 1), np.roll(s, -1, 1) - s)
    return slx, sly
 
def donor_fluxes(s, u, v):                       # 1次風上 (ドナーセル)
    Fx = np.where(u > 0, u * np.roll(s, 1, 0), u * s)
    Fy = np.where(v > 0, v * np.roll(s, 1, 1), v * s)
    return Fx, Fy
 
def bcg_fluxes(s, u, v, dt):                     # 非分割2次予測子
    slx, sly = limited_slopes(s)
    # 1) 法線方向の 1D 予測状態
    sLx = np.roll(s, 1, 0) + 0.5 * (1 - dt/dx * u) * np.roll(slx, 1, 0)
    sRx = s - 0.5 * (1 + dt/dx * u) * slx
    hatX = np.where(u > 0, sLx, np.where(u < 0, sRx, 0.5 * (sLx + sRx)))
    sBy = np.roll(s, 1, 1) + 0.5 * (1 - dt/dx * v) * np.roll(sly, 1, 1)
    sTy = s - 0.5 * (1 + dt/dx * v) * sly
    hatY = np.where(v > 0, sBy, np.where(v < 0, sTy, 0.5 * (sBy + sTy)))
    # 2) 横方向(角)補正 — 予測状態を再利用
    divVy = (np.roll(v * hatY, -1, 1) - v * hatY) / dx
    divUx = (np.roll(u * hatX, -1, 0) - u * hatX) / dx
    sLx -= 0.5 * dt * np.roll(divVy, 1, 0)
    sRx -= 0.5 * dt * divVy
    sBy -= 0.5 * dt * np.roll(divUx, 1, 1)
    sTy -= 0.5 * dt * divUx
    # 3) 同じ風上規則で最終 face 値を選択
    sX = np.where(u > 0, sLx, np.where(u < 0, sRx, 0.5 * (sLx + sRx)))
    sY = np.where(v > 0, sBy, np.where(v < 0, sTy, 0.5 * (sBy + sTy)))
    return u * sX, v * sY
 
def advance(s, u, v, dt, scheme):
    Fx, Fy = bcg_fluxes(s, u, v, dt) if scheme == 'bcg' else donor_fluxes(s, u, v)
    return s - dt * ((np.roll(Fx, -1, 0) - Fx) / dx + (np.roll(Fy, -1, 1) - Fy) / dx)
 
def make_disk():
    X, Y = np.meshgrid((np.arange(N)+0.5)*dx, (np.arange(N)+0.5)*dx, indexing='ij')
    return ((X - 0.5)**2 + (Y - 0.78)**2 < 0.13**2).astype(float)
 
u, v = face_velocities()
Umax = omega * 0.5 * np.sqrt(2)
dt = 0.6 * dx / Umax                    # CFL = 0.6
steps = int(round(1.0 / dt))            # ちょうど1周
 
for scheme in ('upwind1', 'bcg'):
    s = make_disk()
    mass0 = s.sum()
    for _ in range(steps):
        s = advance(s, u, v, dt, scheme)
    print(f"{scheme:8s}  peak={s.max():.3f}  "
          f"mass_err={abs(s.sum()-mass0)/mass0:.2e}")

出力は次のようになります。

upwind1   peak=0.349  mass_err=3.1e-16
bcg       peak=0.806  mass_err=2.8e-16

スキームを 1st-order upwind から BCG unsplit へ切り替えると、ピークが0.349から0.806へ向かう過程が格子の上でやり直されます。出力と同じ数字です。

CFL = 0.60

A scalar disk rotates rigidly about the center. Switch to 1st-order upwind and watch the disk melt into a ring within one revolution; BCG unsplit keeps its edge far longer.

1st-order upwind のままでは、円板は1周し終わる前にリング状へ広がり、縁がぼやけます。BCG unsplit では境界がずっと長く残ります。CFLを0.9に上げるとBCGも少し粗くなりますが、1次との差は残ります。

どちらのスキームも質量を機械精度まで保存します。保存形のフラックス差分なので当然です。分かれるのはピークです。1次風上は1周で円板を半分以下に潰し、BCGは元の高さの80%を保ちます。

2周回った後に残るもの#

コードに書いてある数で1周の長さを換算するとこうなります。N=64なので間隔は1/64です。隅の速さは角速度2πに0.5とルート2を掛けた約4.44です。CFL 0.6なら時間刻みは約0.00211で、時刻1を埋めるステップは474です。円板の中心は(0.5, 0.78)なので、回転中心(0.5, 0.5)から0.28離れています。その場所の速さは2π×0.28、約1.76です。時間刻みを決めるのはこの速さではなく、隅の4.44です。

半径0.13の直径0.26は格子16.6個分です。面積をセル面積で割ると、約217セルがスカラーを持っています。勾配の重みだけ見ると、CFL 0.6で左は0.2、右は0.8です。(1-0.6)と(1+0.6)の半分です。上限0.8までCFLを上げると左の重みは0.1まで落ちます。出力のピーク0.349と0.806は、この重みと横方向項を474回積んだ高さです。

質量誤差3.1e-16と2.8e-16はどちらも丸めです。分かれるのは保存ではなく、失った高さ0.651と0.194です。1次が失う高さはBCGの約3.4倍です。横方向補正は1ステップで時間刻みの半分、約0.00106を横方向の発散に掛けて引きます。その項は1周に474回入ります。数値拡散はセル数64ではなく、1セルで勾配をどれだけ残すかで分かれます。

移植するときに漏れる場所

CFL は分割スキームより厳しい。 非分割 2D は対角方向の情報伝播のため、安定限界が 1 次元より小さくなります。実務では $\mathrm{CFL}\lesssim 0.8$ から始めます。
速度場は必ず face に置く。 セル中心速度を face へ補間すると発散ゼロが崩れ、保存したはずの質量が漏れます。MAC(スタッガード)配置が正解です。
横方向項を忘れない。 法線予測だけを実装すると、ただの次元分割 MUSCL です。コードは動きますが、回転のねじれは残ります。
limiter は極値で 0 になる。 ピークの頂点で勾配が 0 に削られ、局所的に 1 次へ落ちます。これがピークが少しずつ下がる理由です。さらに保ちたいなら極値保存 limiter を検討します。

質量は合うのに円板だけが死ぬ選択

ピーク0.349は、1周後に1次風上が残した高さです。BCGは0.806です。質量誤差は3.1e-16と2.8e-16なので、保存のログだけ見ればどちらも通ります。円板を見ると、残っているのは片方です。

対角に回るスカラーで1次風上を使うかは、この差で決めます。接触面の高さを捨ててよいならドナーセルで足ります。高さを残すなら、横方向項を同じステップに入れます。格子が非構造へ変わると、minmodをそのまま持っていく瞬間に極値の扱いが変わります。Barth–Jespersenリミッタが、その偽の極値を隣の最小値と最大値でどう切るかが次の判断です。