セルの18%で同じ答えを出す — AMR タグ付け基準と coarse-fine リフラックス

バッファをCFLとregrid回数の積より狭くすると界面がpatchの外へ出ます。refluxを外すとdriftの符号が変わります。

5千万セルの解析で、答えを実際に決めているセルは何個でしょうか。衝撃波が一枚、せん断層が一層、火炎面が一枚。たいてい全体の数パーセントです。残りは、滑らかな領域を滑らかに計算するために回っているセルです。

n_buf を 0 に下げ、regrid every を 30 まで押し上げると、界面に沿って赤いセルが増えます。センサーだけが先に出て、格子は追いつきません。

threshold0.35n_buf2max level3regrid every14

leaf 0 / 1 · under-resolved 0

set n_buf to 0 and drag “regrid every” up to 30: red cells appear along the layer every cycle and clear the instant the mesh is rebuilt — that is the feature running out of its own patch. Push n_buf back to 2 and the red stops, but watch the leaf-cell count climb. Raising the threshold thins the patch the other way.

n_buf を 2 に上げると赤が消え、leaf セル数は上がります。バッファを広げる代金がこの数字です。

Adaptive mesh refinement(AMR・解が要求する場所にだけ格子を敷く手法)はこの比率を直接いじります。どこを割るかを決めるセンサー、そして割った後に界面から漏れる質量。二つの急所を数字で残します。

6万5千セルのうち、1万1千セル#

まず得られるものから数えます。Kelvin–Helmholtz せん断層(速度差のある二層の間で巻き上がる界面)を一枚だけ含む $256^2$ のドメインを取ります。均一格子なら 65,536 セルです。同じ界面を3レベルの AMR で包むと、leaf セルは 11,776 個。18% です。

この比率は偶然ではなく、次元が決めています。 $D$ 次元のドメインで界面は $D-1$ 次元です。均一格子のセル数が $N_{\rm uni} = (L/h)^D$ のとき、界面だけを覆うセル数は

N_{\rm AMR} \sim \left(\frac{L}{h}\right)^{D-1} = N_{\rm uni}^{(D-1)/D}

ここで $L$ はドメイン寸法、 $h$ は最細格子間隔です。2D で指数は 1/2、3D では 2/3。つまり格子を2倍細かくすると、均一格子はセルが8倍になるのに対し AMR は4倍で済みます。解像度を上げるほど得が大きくなります。AMR を使う理由はそこです。

今回の格子に式を当てると感覚がつかめます。均一の $256^2$ は 65,536 です。後で数える leaf の合計は 11,776 で、 $11776/65536=0.180$ です。2D でまっすぐな界面だけの見積もり $(L/h)^{1}=256$ は、均一格子の 0.39% です。巻き上がったせん断層にバッファと粗い leaf が乗って 18% になります。3D で間隔を半分にすると均一は $2^3=8$ 倍、界面側の AMR は $2^{2}=4$ 倍です。指数が $(D-1)/D$ なので 3D は 2/3 で、格子を細かくするほど均一比はさらに下がります。

どこを割るかは、勾配では決まらない

まず手が伸びるタグ付け基準は $|\nabla \phi| > \epsilon$ です。これが失敗する理由は単位にあります。圧力勾配は Pa/m、密度勾配は kg/m⁴。フィールドごと、問題ごと、さらにはレベルごとに $\epsilon$ を取り直すはめになります。

Löhner が1987年に出した基準は、この問題を正規化で消してしまいます。2階差分を1階差分の絶対値の和で割るのです。

E_i = \sqrt{ \sum_{d} \frac{\left(\phi_{i+e_d} - 2\phi_i + \phi_{i-e_d}\right)^2}{\Big(|\phi_{i+e_d}-\phi_i| + |\phi_i-\phi_{i-e_d}| + \varepsilon\big(|\phi_{i+e_d}| + 2|\phi_i| + |\phi_{i-e_d}|\big)\Big)^2} }

$e_d$ は $d$ 方向の隣接、 $\varepsilon$ はノイズフィルタ(通常 0.01〜0.05)です。分子と分母が同じ次元なので $E_i$ は無次元で、値はおおむね $[0,1]$ に収まります。閾値 0.3〜0.4 の一つが、圧力にも密度にも、レベル0にもレベル3にもそのまま通用します。

一点だけ注意があります。このセンサーは feature 検出器ではなく解像度検出器 です。界面がそのレベルですでに4〜5セルに解けていれば、2階差分が小さくなってセンサーは切れます。良い性質です — 必要な分だけ割って、勝手に止まってくれます。ただし「衝撃波は問答無用で最大レベルまで」を望むなら、このセンサーだけでは足りません。その場合は $\varepsilon$ のノイズフィルタ項を大きくするのではなく、別の物理基準(例: $\Delta p / p > 0.1$ )を OR で足します。

バッファは、次の regrid までかける保険#

センサーが立ったセルだけを割ると、次のステップで即座に破綻します。格子は毎ステップ作り直すわけではないからです。regrid は普通4〜20ステップに一度回り、その間も界面は移動し続けます。タグ付きセルを $n_{\rm buf}$ セル分だけ膨らませておく理由がこれです。

必要なバッファ幅は素直に計算できます。レベル $\ell$ で $N_{\rm regrid}$ ステップの間に特徴が移動する距離をセル数に換算すると

n_{\rm buf} \;\ge\; \frac{|u|_{\max}\,\Delta t_\ell\,N_{\rm regrid}}{h_\ell} \;=\; \nu\,N_{\rm regrid}

$\nu$ はそのレベルの CFL 数です。CFL 0.4 で10ステップごとに regrid するなら最低4セル要ります。この関係はレベルに依りません — サブサイクリングを使えば $\Delta t_\ell$ と $h_\ell$ が同じ比率で縮むからです。

下限を式にすると $n_{\rm buf}\ge 0.4\times 10=4$ です。コードのバッファ2は、この下限より狭いです。それより狭いと、次の regrid の前に界面がタグの外へ出ます。その赤いセルは、冒頭のダイヤルで n_buf を 0 にするとすぐ見えます。

ブロック・セル・パッチ — 切る単位が決めるもの

同じタグ結果でも、どの単位で割るかによってセル数とコードの複雑さが分かれます。

方式	refine 単位	データ構造	過剰 refine	代表実装
Block-based	固定サイズブロック( $8^3$ など)	ブロック octree	大きい	PARAMESH, FLASH
Cell-based	セル1個	セル単位 tree	なし	OpenFOAM `hexRef8`, RAGE
Patch-based	任意サイズの矩形 patch	box リスト	小さい	Chombo, BoxLib/AMReX

ブロック方式はデータ構造が最も単純で、キャッシュ局所性も良好です。その代わり、ブロック内にタグ付きセルが一つあるだけでブロック全体が割れます。セル方式は過剰がまったくない反面、隣接探索が毎回 tree 走査になります。OpenFOAM はここに属します — dynamicRefineFvMesh がエンジン、hexRef8 がカッター(六面体1個を8個に)、refinementHistory が巻き戻しのための履歴です。必要なマスだけを伸ばす四分木も単位はセルです。タグ一つでブロック全体は割りません。patch 方式は両者の中間で、矩形の内側では構造格子ループをそのまま回せるためベクトル化に有利です。

界面の一つの面に、答えが二つある

ここからが本題です。レベル $\ell$ と $\ell+1$ が接する面を見ます。粗いセル側ではその面は一つ、細かいセル側では $r^{D-1}$ 個です( $r$ は refinement ratio、通常2)。おまけにサブサイクリングを使うと時間ステップも違います。

\Delta t_\ell = \frac{\Delta t_0}{r^{\ell}}

粗い格子が1ステップ進む間に、細かい格子は $r$ 回進みます。ですからその一つの面を通して粗い格子が計算した flux $F^c$ と、細かい格子が $r^{D-1}$ 個の面で $r$ 回かけて実際に押し出した flux は、別の数字になります。二つの差が

\delta F_f = \frac{1}{r^{D}}\sum_{s=1}^{r}\sum_{f'\subset f} F^{s}_{f'} \;-\; F^{c}_{f}

であり、これが界面の作り出した質量です。無かった質量が生まれ、有った質量が消えます。保存形スキームを使っていても同じです — 保存性は一つの格子の内側でしか成り立たず、二つの格子が接する場所では保証されません。

帳簿に書いておいて、あとで一括清算する

解決策は単純です。粗いステップを始めるときに $F^c_f$ を帳簿(flux register)に書いておきます。細かい格子を $r$ 回回しながら、実際の flux を同じ帳簿に積み上げます。粗いステップが終わったら、差額 $\delta F_f$ を patch の外側の粗いセルに戻します。

\phi_c \;\leftarrow\; \phi_c \;\mp\; \frac{\Delta t_c}{h_c}\,\delta F_f

符号は、そのセルが面のどちら側にあるかで決まります。patch の内側のセルには手を付けません — すでに細かい格子の値を平均して上書きしているからです。この一行が、ドメイン全体の質量を機械精度まで戻します。

$r=2$ 、 $D=2$ なら $r^{D}=4$ です。細かい面 $r^{D-1}=2$ 個とサブサイクル $r=2$ 回をこの4で割り、粗い flux 一つと桁を合わせます。割り算が抜けると $\delta F$ の大きさがずれ、帳簿に書いて返しても質量は合いません。

コードで数えたセル数と質量ドリフト

まずタグ付けから。せん断層のスナップショット一枚に Löhner センサーをかけ、バッファを与えてからブロック単位でレベルを積み、leaf セルを数えます。

import numpy as np
 
N_EFF, BLOCK, MAX_LEVEL = 256, 4, 2   # 256^2 相当の解像度、ブロックあたり 4x4 セル、レベル 0〜2
EPS_L, THRESH, N_BUF = 0.02, 0.35, 2  # ノイズフィルタ / タグ閾値 / バッファ(セル単位)
 
 
def shear_layer(n, t=1.35):
    """Kelvin-Helmholtz 巻き上がりのスナップショット — ソルバなしでフィールドだけ作る。"""
    x = (np.arange(n) + 0.5) / n
    xx, yy = np.meshgrid(x, x, indexing='ij')
    warp = 0.06 * np.sin(2 * np.pi * xx + t) + 0.025 * np.sin(4 * np.pi * xx - 2 * t)
    return np.tanh((yy - 0.5 - warp) / 0.012)
 
 
def shift(a, d, ax):
    """隣接参照: x は周期、y は zero-gradient — 壁で偽のジャンプを作らない。"""
    return np.roll(a, d, ax) if ax == 0 else np.pad(a, 1, mode='edge')[1:-1, 1 + d:a.shape[1] + 1 + d]
 
 
def lohner_sensor(f):
    """正規化された2階差分。無次元なので閾値一つが全レベルに通用する。"""
    e2 = np.zeros_like(f)
    for ax in (0, 1):
        p, m = shift(f, -1, ax), shift(f, 1, ax)
        num = np.abs(p - 2.0 * f + m)
        den = np.abs(p - f) + np.abs(f - m) + EPS_L * (np.abs(p) + 2 * np.abs(f) + np.abs(m))
        e2 += (num / np.maximum(den, 1e-30)) ** 2
    return np.sqrt(e2)
 
 
def grow(mask, width):
    """バッファ: 次の regrid までに特徴が patch の外へ出ないようにする。"""
    for _ in range(width):
        out = mask.copy()
        for ax in (0, 1):
            out |= shift(mask, 1, ax) | shift(mask, -1, ax)
        mask = out
    return mask
 
 
def tag_blocks(f, thresh, n_buf):
    """セルセンサー -> セル単位バッファ -> タグ付きセルを一つでも含むブロックは refine。"""
    tagged = grow(lohner_sensor(f) > thresh, n_buf)
    nb = f.shape[0] // BLOCK
    return tagged.reshape(nb, BLOCK, nb, BLOCK).any(axis=(1, 3))
 
 
def leaf_cells(field):
    """粗いレベルから下へ辿り、実際に解かれるセルだけを数える。"""
    counts, live = [], None
    for lev in range(MAX_LEVEL + 1):
        n = N_EFF >> (MAX_LEVEL - lev)
        f = field.reshape(n, N_EFF // n, n, N_EFF // n).mean(axis=(1, 3))
        nb = n // BLOCK
        child = np.zeros((nb, nb), bool) if lev == MAX_LEVEL else tag_blocks(f, THRESH, N_BUF)
        live = np.ones((nb, nb), bool) if live is None else live
        counts.append(int((live & ~child).sum()) * BLOCK * BLOCK)
        live = np.kron(live & child, np.ones((2, 2), bool))
    return counts
 
 
counts = leaf_cells(shear_layer(N_EFF))
total, uniform = sum(counts), N_EFF * N_EFF
for lev, c in enumerate(counts):
    print(f'  level {lev}  h = 1/{N_EFF >> (MAX_LEVEL - lev):<3d}  leaf cells = {c:6d}')
print(f'  AMR total     = {total}')
print(f'  uniform 256^2 = {uniform}  ->  {100 * total / uniform:.1f} % of the cells')

出力はこうなります。

  level 0  h = 1/64   leaf cells =   3280
  level 1  h = 1/128  leaf cells =   1520
  level 2  h = 1/256  leaf cells =   6976
  AMR total     = 11776
  uniform 256^2 = 65536  ->  18.0 % of the cells

N_BUF を 0, 1, 2, 4 と変えると、合計は 9,616 → 10,864 → 11,776 → 14,224 と動きます。バッファ2セルの値段はセル2,160個、均一格子比で 3.3%p です。

0から2まで増えた 2,160 は 9,616 の約22%です。4セルまで行くと 14,224 で、バッファ0より 4,608 セル増えます。レベル別に見ると最細の 6,976 が全体 11,776 の59%です。レベル0の 3,280 は、滑らかな領域を粗い間隔で覆った分です。18% は界面だけの値ではなく、この三層の合計です。

次はリフラックスです。2D スカラー移流に二つのレベルを載せ、サブサイクリングまで入れてから全体質量を測ります。

import numpy as np
 
NC, R, CFL, NSTEP = 48, 2, 0.4, 60   # 粗い格子 / refinement ratio / CFL / 粗いステップ数
BOX = (12, 28, 16, 32)               # patch の角(粗い格子インデックス基準)
U, V = 1.0, 0.6
 
 
def upwind_faces(f, dx, dy):
    """周期ブロックの全 x・y 面での donor-cell flux。(Fx, Fy) を返す。"""
    fx = U * (f if U > 0 else np.roll(f, -1, 0))          # 面 i はセル i の左側
    fy = V * (f if V > 0 else np.roll(f, -1, 1))
    return np.roll(fx, 1, 0), np.roll(fy, 1, 1)
 
 
def march_block(f, fx, fy, dt, dx, dy):
    return f - dt / dx * (np.roll(fx, -1, 0) - fx) - dt / dy * (np.roll(fy, -1, 1) - fy)
 
 
def gaussian_patch(n, x0, y0, s):
    c = (np.arange(n) + 0.5) / n
    xx, yy = np.meshgrid(c, c, indexing='ij')
    return np.exp(-((xx - x0) ** 2 + (yy - y0) ** 2) / s ** 2)
 
 
def two_level_run(reflux):
    i0, i1, j0, j1 = BOX
    dx, dxf = 1.0 / NC, 1.0 / NC / R
    dt = CFL * dx / (abs(U) + abs(V))
    dtf = dt / R
 
    coarse = gaussian_patch(NC, 0.32, 0.42, 0.09)
    fine = np.kron(coarse[i0:i1, j0:j1], np.ones((R, R)))   # patch は粗い値と整合するように始める
 
    for _ in range(NSTEP):
        cfx, cfy = upwind_faces(coarse, dx, dx)
        # 粗い格子が patch 境界を通ると「思い込んだ」flux
        edge_c = {'lo_x': cfx[i0, j0:j1].copy(), 'hi_x': cfx[i1, j0:j1].copy(),
                  'lo_y': cfy[i0:i1, j0].copy(), 'hi_y': cfy[i0:i1, j1].copy()}
        coarse = march_block(coarse, cfx, cfy, dt, dx, dx)
 
        edge_f = {k: np.zeros_like(v) for k, v in edge_c.items()}
        for _ in range(R):                                   # サブサイクリング: 粗い1ステップにつき細かい R ステップ
            g = np.zeros((fine.shape[0] + 2, fine.shape[1] + 2))
            g[1:-1, 1:-1] = fine
            g[0, 1:-1] = np.repeat(coarse[i0 - 1, j0:j1], R)  # ghost: 粗い値を区分定数で注入
            g[-1, 1:-1] = np.repeat(coarse[i1, j0:j1], R)
            g[1:-1, 0] = np.repeat(coarse[i0:i1, j0 - 1], R)
            g[1:-1, -1] = np.repeat(coarse[i0:i1, j1], R)
            gfx, gfy = upwind_faces(g, dxf, dxf)
            fine = march_block(g, gfx, gfy, dtf, dxf, dxf)[1:-1, 1:-1]
            for k, s in (('lo_x', gfx[1, 1:-1]), ('hi_x', gfx[-1, 1:-1]),
                         ('lo_y', gfy[1:-1, 1]), ('hi_y', gfy[1:-1, -1])):
                edge_f[k] += s.reshape(-1, R).mean(axis=1) / R   # 面方向・時間方向の平均
 
        coarse[i0:i1, j0:j1] = fine.reshape(i1 - i0, R, j1 - j0, R).mean(axis=(1, 3))
 
        if reflux:                                            # 帳簿の清算
            coarse[i0 - 1, j0:j1] -= dt / dx * (edge_f['lo_x'] - edge_c['lo_x'])
            coarse[i1, j0:j1] += dt / dx * (edge_f['hi_x'] - edge_c['hi_x'])
            coarse[i0:i1, j0 - 1] -= dt / dx * (edge_f['lo_y'] - edge_c['lo_y'])
            coarse[i0:i1, j1] += dt / dx * (edge_f['hi_y'] - edge_c['hi_y'])
 
    mask = np.ones((NC, NC), bool)
    mask[i0:i1, j0:j1] = False
    return (coarse * mask).sum() * dx * dx + fine.sum() * dxf * dxf
 
 
m0 = gaussian_patch(NC, 0.32, 0.42, 0.09).sum() / NC ** 2
for tag, on in (('reflux off', False), ('reflux on ', True)):
    m = two_level_run(on)
    print(f'  {tag}:  mass = {m:.12f}   drift = {(m - m0) / m0:+.3e}')

  reflux off:  mass = 0.025702960114   drift = +1.006e-02
  reflux on :  mass = 0.025446894886   drift = +0.000e+00

pulse sigma を下げると赤い逸脱幅が大きくなり、帳簿の $\delta F$ も同じだけ大きくなります。

ratio r4CFL0.60pulse sigma0.050

step 0 · drift 0.00e+0 vs 0.00e+0

the pulse is harmless until it reaches the pink face. From that step on the red trace walks away from zero and never comes back, while the green one stays pinned at machine zero. Sharpen the pulse (sigma down) or raise r and the red excursion grows — the register grows with it, because it is exactly the same number with the opposite sign.

パルスがピンクの面に触れると、赤いトレースだけが0から離れ、緑は0に貼り付いたままです。リフラックスを切った側の drift +1.006e-02 がその赤い線です。

60ステップで1%です。出力の二つの質量を引くと約 $2.56\times 10^{-4}$ です。on 側の質量で割ると drift +1.006e-02 と同じ1%です。しかもこの値は単調に増えるわけでもありません。15ステップでは −1.3%、30ステップで −0.77%、60ステップで +1.0% でした。ブロブが patch を出入りするたびに符号が変わります。収束試験を回していてこんなグラフを見るとスキームを疑いたくなりますが、犯人はリフラックスを入れていない界面です。

ロードバランシング — 曲線一本でブロックを並ばせる

AMR を並列で回すと新しい問題が出ます。regrid のたびにブロック数がプロセッサごとに変わるのです。界面を抱えていたランクはセルが5倍になり、滑らかな領域しか持たないランクは変わりません。

標準的な処方は space-filling curve(SFC・多次元格子を一本線で辿る曲線)です。Morton(Z-order)や Hilbert 曲線で全 leaf ブロックに1次元インデックスを振り、その列をランク数だけ均等分割します。曲線の局所性のおかげで、隣り合うインデックスはおおむね空間的にも隣接するので、均等分割がそのまま通信量の小さい分割になります。p4est と AMReX がこの方式です。グラフ分割器(ParMETIS, Zoltan, Scotch)を使えば分割品質はさらに良くなりますが、regrid のたびに回し直す必要があるためコストは SFC よりずっと大きくなります。AMR は格子が頻繁に変わる側なので、たいてい SFC が勝ちます。

リフラックスを外すと、収束曲線がスキームを犯人にする

タグは無次元に置きます。正規化した2階差分なら、閾値一つがフィールドとレベルをまたぎます。生の勾配はレベルごとに取り直すことになります。バッファは $\nu N_{\rm regrid}$ で計算します。CFL 0.4 で10ステップなら4セルです。勘で1セルを入れ、regrid を20ステップごとにすると、ステップのかなりの部分は特徴が patch の外へ出たまま回ります。

リフラックスをオプションにすると、coarse-fine 界面が質量を作ります。60ステップの +1% は単調でもありません。符号の変わる drift を見て移流スキームを疑いやすいです。refine の順を一行ずらすと誤差が膨らんだ LBM の時間前進も同じ種類です。順序が違うと、セルを細かく割っても界面の答えは違います。