壁では3個、角では5個 — LBM境界ノードが失う分布関数を数える

境界条件の実装はスキーム選びから始まりません。ノードごとに何個空くかを数えるところから始まります。

境界ノードは全体の1%、コードは半分#

格子ボルツマン法(LBM)のソルバを開くと、比率が妙です。衝突項は10行ほどです。ストリーミングは 5行です。境界条件は数百行あります。

計算量で見ると逆です。100×100の格子で境界ノードは400個ほど。全体の4%です。3次元に行くと1%を下回ります。演算の1%がコードの半分を占めています。

この不均衡には理由があります。境界条件そのものが難しいのではなく、ノードごとに解くべき問題の大きさが違うからです。その大きさを測る規則を先に決めれば、コードはまた短くなります。今回はその規則と、規則が決まったあとにデータ構造がどう決まるかを見ていきます。

空く場所を決めるのはスキームではなく形状です

ストリーミングは隣から値を引いてくる操作です。

f_k(\mathbf{x},\, t + \Delta t) = f_k^{\star}(\mathbf{x} - \mathbf{e}_k \Delta t,\; t)

ここで $f_k$ は方向 $k$ の分布関数、 $\mathbf{e}_k$ はその格子速度、星印は衝突直後の値です。値が来る場所は $\mathbf{x} - \mathbf{e}_k \Delta t$ 、すなわち上流の隣接ノードです。

その上流が固体なら、送るべき値がありません。そのリンクは空のまま到着します。ですからノード一つで空く分布関数の数は、とても素朴な量になります。そのノードの8近傍のうち固体であるセルの数と等しいのです。

bounce-back でも Zou–He でも、スキームはこの数を変えられません。数を決めるのは形状だけです。スキームが答えるのは次の問い、空いた場所を何で埋めるかだけです。

下の格子でノードを直接クリックしてみてください。

size3

Click nodes along the bottom wall: 3 red arrows every time. Click the inside corner where the step meets the floor — it jumps to 5, and the closure bar goes 3 short. The single node on top of the step corner is the opposite case: 1 unknown, and the 3 moment equations are one too many. Current pick: concave corner, 5 unknown.

床沿いを押していくと、赤い矢印はいつも3本です。段差が床と出会う内側の角では5本に跳ねます。段差の上の外側の角では1本に落ちます。変わったのは形状だけなのに、解くべき未知数の数が3倍以上開きます。

未知数の帳簿 — モーメント3本で覆えるもの#

空いた分布関数を埋めるには条件が要ります。使える条件は巨視量の定義だけです。

\rho = \sum_{k=0}^{8} f_k, \qquad \rho\, \mathbf{u} = \sum_{k=0}^{8} f_k\, \mathbf{e}_k

2次元ではこれが3本です。密度が1本、運動量が2本です。

未知数の側を数えましょう。空の分布関数が $|\mathcal{U}|$ 個あります。壁では普通、速度を与えて密度は知りません。ですから $\rho$ も未知数です。不足分はこう書けます。

d = |\mathcal{U}| + 1 - (D + 1), \qquad D = 2

$D$ は空間次元で、 $D+1$ が使えるモーメント式の本数です。平壁なら $|\mathcal{U}|=3$ なので $d=1$ 。式が1本足りません。Zou–He が非平衡 bounce-back を一つ足すのは、ちょうどここです。

f_k - f_k^{\mathrm{eq}} = f_{\bar{k}} - f_{\bar{k}}^{\mathrm{eq}}

$\bar{k}$ は $k$ の逆向きです。壁法線方向のリンク対にこの式を課すと帳簿が合います。詳しい導出は bounce-back と Zou–He を並べて比べた記事に書いてあります。

凹角では $|\mathcal{U}|=5$ です。 $d=3$ で、条件が3本足りません。平壁用に書いた閉包関係を一つそのまま持ってくると、二つが浮いたまま残ります。そこに残っている値は初期値か、前ステップの残りかすです。

これが「コードは回るのに角だけ値が変」の正体としてよくあるものです。発散はしません。静かに間違えます。

Python で格子を一枚なめてみました#

段差が一つある流路を作り、すべての流体ノードで空の方向を数えてみましょう。あとのデータ構造の話のために、キャッシュラインも一緒に数えます。

# D2Q9: 0 静止、1-4 軸方向、5-8 斜め方向
E = [(0, 0), (1, 0), (0, 1), (-1, 0), (0, -1), (1, 1), (-1, 1), (-1, -1), (1, -1)]
NX, NY = 24, 16
 
 
def solid_mask(nx, ny):
    """床に段差が一つある流路。"""
    m = [[False] * ny for _ in range(nx)]
    for i in range(nx):
        m[i][0] = True
        m[i][ny - 1] = True
    for i in range(8):
        for j in range(1, 5):
            m[i][j] = True
    return m
 
 
def unknown_dirs(m, i, j):
    """上流の隣接点 (i-ex, j-ey) が固体、または格子外である k。"""
    nx, ny = len(m), len(m[0])
    out = []
    for k in range(1, 9):
        si, sj = i - E[k][0], j - E[k][1]
        if not (0 <= si < nx and 0 <= sj < ny) or m[si][sj]:
            out.append(k)
    return out
 
 
def node_class(unk):
    axial = [k for k in unk if k <= 4]
    if len(axial) == 0:
        return "convex corner"
    if len(axial) == 1:
        return "flat wall"
    if len(axial) == 2:
        return "concave corner"
    return "slot / thin gap"
 
 
def scan_boundary(m):
    """分布関数を一つでも失う流体ノードすべてを、行優先の順で。"""
    ny = len(m[0])
    rows = []
    for i in range(len(m)):
        for j in range(ny):
            if m[i][j]:
                continue
            unk = unknown_dirs(m, i, j)
            if unk:
                rows.append((i * ny + j, node_class(unk), unk))
    return rows
 
 
def lines_touched(rows, n_nodes, layout):
    """空の分布関数を埋める際に読む64バイトライン(double 8個)の数。"""
    s = set()
    for lin, _, unk in rows:
        for k in unk:
            addr = k * n_nodes + lin if layout == "soa" else lin * 9 + k
            s.add(addr // 8)
    return len(s)
 
 
mask = solid_mask(NX, NY)
rows = scan_boundary(mask)
n_nodes = NX * NY
n_fluid = sum(1 for i in range(NX) for j in range(NY) if not mask[i][j])
 
print("lattice %dx%d   fluid %d   boundary %d (%.1f%% of fluid)"
      % (NX, NY, n_fluid, len(rows), 100.0 * len(rows) / n_fluid))
print()
print("%-16s %7s %6s %9s %9s" % ("class", "unk/node", "nodes", "unknowns", "closure"))
groups = {}
for lin, cls, unk in rows:
    groups.setdefault((cls, len(unk)), 0)
    groups[(cls, len(unk))] += 1
for (cls, n_unk) in sorted(groups, key=lambda g: (g[1], g[0])):
    n = groups[(cls, n_unk)]
    gap = n_unk + 1 - 3          # 空の分布関数 + rho 対 モーメント3本
    tag = "%+d" % gap if gap else "exact"
    print("%-16s %7d %6d %9d %9s" % (cls, n_unk, n, n_unk * n, tag))
print()
print("total unknown PDFs           %d" % sum(len(r[2]) for r in rows))
print("cache lines, SoA f[k][node]  %d" % lines_touched(rows, n_nodes, "soa"))
print("cache lines, AoS f[node][k]  %d" % lines_touched(rows, n_nodes, "aos"))

出力はこうなります。

lattice 24x16   fluid 304   boundary 72 (23.7% of fluid)
 
class            unk/node  nodes  unknowns   closure
convex corner          1      1         1        -1
flat wall              2      2         4     exact
flat wall              3     64       192        +1
concave corner         5      5        25        +3
 
total unknown PDFs           222
cache lines, SoA f[k][node]  153
cache lines, AoS f[node][k]  84

段差一つの形状でノードの種類が4つ出ました。平壁なのに空の方向が2つのノードも2個あります。段差の角のすぐ隣で、斜めリンクが一本生き残る場所です。長方形の箱だけを想定して書いたコードが実形状で崩れる地点は、こうして生まれます。

凸角では式が余ります

表で目に留まるのは最初の行です。凸角の不足分が $-1$ になっています。

空の分布関数は斜めリンク一本だけです。未知数はそれと $\rho$ 、合わせて2つ。モーメント式は3本です。式が1本余ります。

ここで3本すべてを課すと過剰決定になります。どの組み合わせを選んでも、残る1本は満たされません。無理に合わせると質量が漏れ始めます。

ですから凸角では普通、閉包関係を使いません。空いたリンク一本に bounce-back をかけて終わりです。方程式を解く代わりに、値を戻しておくわけです。

不足分の符号が処方を分けます。正なら条件を足す、0ならそのまま解く、負なら解くこと自体をやめる。一つのコードの中でこの3つが同時に出てきます。

方向性で分類してから配列が決まります

ここまで来ると、データ構造は自然に決まります。

ノードを二つの基準で分類します。一つ目は方向性です。どちら側の隣が欠けているか。2次元なら面が4つ、角が4つ、八つの区分です。二つ目は境界条件の種類です。壁か、速度入口か、圧力出口か。

この二軸の組み合わせごとに、空く方向の集合が固定されます。集合が固定されれば分岐が消えます。ループの中で if を使って方向を判定する代わりに、同じ処理を受けるノードを一かたまりに集め、そのかたまりを丸ごと回します。

そのためには、同じ区分のノードが配列上で連続に並ぶ必要があります。区分ごとのノード数を持つ配列を一つ、ノード番号を入れる配列(iNodeBC)を一つ用意します。前処理で一度だけ埋め、時間ループでは読むだけです。固定形状ならこのコストは全体で一度きりです。

第二段階は、各境界ノードに属する分布関数の番号をあらかじめ保存することです。このときノード一つの分布関数9個をメモリ上で隣接させておくと — 構造体配列(AoS, Array of Structure)配置 — 境界ループが引き込むメモリが減ります。

同じ未知数、違うメモリ

上のスクリプトが数えた二つの数字がその差です。SoA 配置で153ライン、AoS 配置で84ラインです。読む値の個数はどちらも222で同じです。違うのは配置だけです。

理由は、ストリーミングと境界ループのアクセスパターンが逆だからです。ストリーミングは方向 $k$ を一つ固定して格子全体をなめます。f[k][node] 配置が有利です。境界ループはノードを一つ固定して複数の方向をなめます。同じ配置では、そのノードの未知数は8つの方向ブロックに散らばっています。

下で配置を切り替えながら、同じスイープを回してみてください。

step h3speed9

Let one sweep finish on soa and read the line count, then hit aos — the same 0 nodes and the same unknowns, but the boxes light up in short runs instead of eight scattered bands. packed collapses the whole map into its top-left corner: the loop never reads a line it does not need. Now: 0 lines for 0 nodes.

soa で一周させてライン数を読み、aos を押して同じスイープをもう一度見てください。点灯するマスが、まばらな八つの帯から短いかたまりに変わります。packed は境界ノードを番号付けし直して並べた場合です。地図が左上の隅に折り畳まれます。

注意したいのは、これが全体の配置を変えろという話ではないことです。コード全体を AoS で回すとストリーミングと衝突が遅くなります。 MRT衝突をモーメント空間で扱った記事で見たとおり、衝突ループは方向ごとの連続アクセスを好みます。要点は、境界条件用の局所データ構造を別に持つことです。全ノードの数パーセントですから、コピーのコストもその程度です。

境界条件のバグがスキームのせいでないとき

新しい形状を載せたら壁の近くだけおかしい、というとき、手を付ける順番があります。

まずノードを数えます。上のスクリプトのように、種類ごとの個数と不足分を出力します。長方形の箱で flat wall しか出ていなかったのに、新しい形状で concave corner や slot / thin gap が現れたなら、その行がコードで処理されているかをまず確認します。

次は不足分の符号です。正の行に閉包関係が何本かかっているかを数えます。負の行でモーメントを強制していないかを見ます。

配置の話は最後です。値が合ってから見ることです。順番を変えると、速く間違った答えが出ます。