[論文レビュー] たった1ステージで2次精度 — Kapila多相流のためのGRPソルバ

界面で時間微分まで求めると、ステージを足さなくても時間2次が出ます。

時間2次が欲しければ Runge-Kutta を何段も回す、と覚えていました。この論文は1ステージでその次数を作ります。界面で値と時間微分を一緒に返す GRP（Generalized Riemann Problem）ソルバです。

論文: Tuowei Chen, Zhifang Du, Generalized Riemann Problem Method for the Kapila Model of Compressible Multiphase Flows, arXiv:2310.08241 (2023).

GRP のチェックを外すと勾配スライダーがロックされ、界面値は時間に凍ります。CFL を上げると、フラックスに入る半ステップの位置だけが前へ出ます。

advection speed a: 1.00left slope uₓ: 0.50CFL: 0.80GRP (off = Godunov)

勾配を 0 にするとシアンの点が止まり、Godunov と重なります。移流速度 $a$ の符号をひっくり返すと、値を持ってくるセルが変わります。橙の円は $\Delta t/2$ の時刻です。

1ステージで2次 — GRPが解く問題#

有限体積法の1ステップは界面フラックス $\mathbf{F}_{i+1/2}$ で決まります。Godunov 法は界面で Riemann 問題を解いて値を一つ得ます。その値は時間について定数です。時間2次が欲しければ、同じ空間演算を Runge-Kutta のように何度も繰り返すことになります。

GRP はフラックスを中間時刻 $t^{n+1/2}$ で評価します。

\mathbf{U}(\mathbf{x}_{i+1/2}, t^{n+1/2}) = \mathbf{U}^{n,*}_{i+1/2} + \frac{\Delta t}{2}\left(\frac{\partial \mathbf{U}}{\partial t}\right)^{n,*}_{i+1/2}

$\mathbf{U}^{n,*}$ は Riemann 解（界面値）、 $(\partial \mathbf{U}/\partial t)^{n,*}$ はその界面の瞬間の時間微分です。二つを Taylor 展開で足すと、界面値は $\Delta t/2$ だけ未来へ動きます。評価は1回です。

区分定数から区分線形へ — RiemannとGRPの違い#

Riemann 問題は界面の両側が区分定数です。GRP 問題は両側が区分線形です。初期不連続から出る波はもう直線ではなく、ソルバは Riemann 解と時間微分を一緒に返します。

スカラー移流 $u_t + a u_x = 0$ で、セル $i$ の勾配を $\sigma_i$ とします。 $a>0$ なら界面の Riemann 解は左セルからの外挿です。

U^{*}_{i+1/2} = u_i + \sigma_i \frac{\Delta x}{2}

$u_i$ はセル平均、 $\sigma_i$ はセル内勾配、 $\Delta x$ はセル幅です。先のシアンの点がこの式で、勾配を 0 にした還元もこの式から出ます。

Lax-Wendroff手続き — 時間微分を空間変化へ#

時間微分は支配方程式から取り出します。

\frac{\partial u}{\partial t} = -a \frac{\partial u}{\partial x} = -a\,\sigma_i

左の時間微分を右の空間微分に替えます。 $\sigma_i$ は再構成ですでに分かっています。半ステップの界面値は

U^{n+1/2}_{i+1/2} = U^{*}_{i+1/2} - a\,\sigma_i \frac{\Delta t}{2}

で、フラックスは $F = a\,U^{n+1/2}$ です。Euler や Kapila では $\partial u/\partial t = -A\,\partial u/\partial x$ に一般化され、 $A$ はヤコビアンです。論文はこの手続きで音波を線形化し、GRP を近似的に解きます。

Python — GRPスカラーソルバで収束次数を測る#

滑らかな初期値 $u_0=\sin(2\pi x)$ を1周期だけ移流させ、Godunov と GRP の $L_1$ 誤差の次数を比べます。領域は周期境界の線形移流です。

import numpy as np
 
# 1D 線形移流  u_t + a u_x = 0,  周期領域 [0,1]
# 1次 Godunov vs 単一ステージ2次 GRP(Lax-Wendroff)の比較
 
def cell_gradient(u, dx):
    # 滑らかな解の次数測定のため中心(無制限)勾配を使用
    return (np.roll(u, -1) - np.roll(u, 1)) / (2.0 * dx)
 
def godunov_step(u, a, dx, dt):
    # 区分定数: 風上界面値 (a > 0)
    flux = a * u                       # F_{i+1/2} = a u_i
    return u - dt / dx * (flux - np.roll(flux, 1))
 
def grp_step(u, a, dx, dt):
    g = cell_gradient(u, dx)           # 空間勾配 sigma_i
    u_star = u + g * dx / 2.0          # 界面Riemann解 U*
    u_half = u_star - a * g * dt / 2.0 # Lax-Wendroff半ステップ U^{n+1/2}
    flux = a * u_half                  # F_{i+1/2}^{n+1/2}
    return u - dt / dx * (flux - np.roll(flux, 1))
 
def l1_error(num, exact, dx):
    return np.sum(np.abs(num - exact)) * dx
 
def convergence_study(a=1.0, cfl=0.4, t_end=1.0):
    print(f"{'N':>5} {'Godunov':>16} {'GRP':>16}")
    prev = {}
    for N in [40, 80, 160, 320]:
        dx = 1.0 / N
        x = (np.arange(N) + 0.5) * dx
        u0 = np.sin(2 * np.pi * x)     # 滑らかな初期条件
        dt = cfl * dx / abs(a)
        nsteps = int(round(t_end / dt))
        dt = t_end / nsteps            # t_endに正確に到達
        ug = u0.copy(); ur = u0.copy()
        for _ in range(nsteps):
            ug = godunov_step(ug, a, dx, dt)
            ur = grp_step(ur, a, dx, dt)
        exact = np.sin(2 * np.pi * (x - a * t_end))
        eg = l1_error(ug, exact, dx)
        er = l1_error(ur, exact, dx)
        og = f"{np.log2(prev['g']/eg):.2f}" if 'g' in prev else "  - "
        orr = f"{np.log2(prev['r']/er):.2f}" if 'r' in prev else "  - "
        print(f"{N:>5} {eg:.2e}({og}) {er:.2e}({orr})")
        prev = {'g': eg, 'r': er}
 
convergence_study()

出力は次のとおりです。

    N          Godunov              GRP
   40 1.63e-01(  - ) 1.31e-03(  - )
   80 8.76e-02(0.90) 2.69e-04(2.28)
  160 4.54e-02(0.95) 6.32e-05(2.09)
  320 2.31e-02(0.97) 1.55e-05(2.03)

$\lambda\Delta t$ のスライダーを 1.6 より上げると、陽的 Euler が $[0,1]$ の帯を抜けるかどうかが変わります。初期の体積分率と平衡値を動かしても、Crank-Nicolson が帯の中に残るかを一緒に見ます。

stiffness λΔt: 1.60α₀ (initial): 0.85α_eq (target): 0.10

$\lambda\Delta t=3$ までシアンは平衡へ降り、赤は振動します。いまの表の次数 2.03 は滑らかな移流の数字で、このスライダーは、その半ステップの時間微分がソース積分に入るときの話です。

N=40 の誤差比は $1.63\times 10^{-1}/1.31\times 10^{-3}$ 、約 124 倍です。同じ格子で 100 倍小さい、という話はこの格子を指します。N=320 では $2.31\times 10^{-2}/1.55\times 10^{-5}$ 、約 1490 倍です。格子を 40 から 320 へ 8 倍細かくすると、Godunov の誤差は 7.1 倍しか減らず、GRP は 85 倍減ります。8 の二乗は 64 です。印刷された次数が 2.28 から 2.09、2.03 へ降りてくるあいだに、64 を少し超えています。Godunov 側の 0.90、0.95、0.97 も同じ向きに 1 へ付いていきます。

N=80 の比は $8.76\times 10^{-2}/2.69\times 10^{-4}$ 、約 326 倍です。N=160 は $4.54\times 10^{-2}/6.32\times 10^{-5}$ 、約 718 倍です。格子を2倍ずつ細かくするたびに誤差比は 124、326、718、1490 と増え、一段で2倍を少し超えます。

Kapilaモデルの硬い体積分率ソース#

Kapila の5方程式モデル（BN 7方程式の縮約）の体積分率方程式には、硬いソースが付きます。

\frac{\partial \alpha_1}{\partial t} + \mathbf{u}\cdot\nabla\alpha_1 = \alpha_1 K\,\nabla\cdot\mathbf{u}

$\alpha_1$ は相1の体積分率、 $K$ は Wood 音速から出る圧縮・膨張の係数です。この $K$ が体積分率に対してどう折れるかは、混合物の中で Wood 音速が落ちる計算と同じ場所です。ソースは $\alpha_1$ を $[0,1]$ の外へ押せます。陽的積分は剛性が大きくなると境界を越えて発散します。

GRP の時間微分は、界面値を新しい時刻 $t^{n+1}$ で知らせます。Gauss-Green で $\nabla\cdot\mathbf{u}$ のセル平均をその時刻に推定できます。ソースを Crank-Nicolson（半陰的）で積分し、体積分率の境界を守る時間刻みの制限まで立てます。先の $\lambda\Delta t$ 曲線がその制限の模型です。

再現して怪しかったこと

スカラー GRP は2次が出ましたが、Kapila 系はそれより粗いです。無制限の中心勾配は衝撃波で振動します。次数が出たのはテストが滑らかだったからです。不連続では minmod のような制限子が必要で、そのとき次数は局所的に1次へ落ちます。

Wood 音速は体積分率に対して非単調です。界面の数値拡散帯でマッハ数が揺れる原因です。GRP の時間微分は圧縮と膨張にはよく付いていきますが、その非単調そのものは消しません。論文もそこは正面から解くより、頑健さで避けています。

実務のコードはたいてい Strang 分割に Runge-Kutta を付けた構造なので、単一ステージの GRP をそのまま入れるのは難しいです。利点を残すにはフラックス関数を Lax-Wendroff 形に書き直します。移植の費用は残ります。線形化した音響近似の出所は、Ben-Artzi と Falcovitz の気体力学 GRP です。

ソースをフラックスの外に置いたまま次数を期待した失敗

表の 2.03 は、勾配が滑らかなときの数字です。衝撃波で制限子を入れると、その席は局所的に1次へ戻ります。次数を見てステージを消す前に、フラックス関数を書き直す予算があるかが先です。

ソースをフラックスに溶かさず外へ出すと、同じ剛性を別の時間積分が引き受けます。圧力を二つで解いてから一つに合わせる緩和は、その分離側で体積分率の境界を守るのか、それとも $\Delta t$ をまた細かく割るのか。この出力の 2.03 はその問いへの答えではありません。