クアッドツリー AMR — 格子を必要な場所だけ育てる

τを0.05から0.4へ上げると裾は割れず、2:1を外すと面の流束が先に合わなくなる。

1984年、スタンフォードです。Marsha Berger と Joseph Oliger は衝撃波を一本追っていました。セルは100万です。意味のある仕事をしていたのは衝撃波付近の約1%だけで、残り99%は均一な自由流を同じ精度で計算し直していました。その年の Journal of Computational Physics が Berger–Oliger AMR です。同じ癖が、今も OpenFOAM の dynamicRefineFvMesh、AMReX、Chombo、p4est の中に残っています。

白い点をドラッグすると格子はその点の周りだけ割れ、threshold のスライダーが切る基準である傾き $\tau$ を変えます。

gradient threshold: 0.180max level: 5

drag the white dot to move the feature. higher threshold = coarser grid. cell count vs uniform refinement is the AMR savings.

$\tau$ を 0.05 まで下げると裾まで leaf が増え、0.4 を超えると粗い格子だけが残ります。add feature でもう一つの峰を置くと深いのは二つの核だけで、左上の数字は一様格子に対するセル数の比です。

均一格子はなぜ 90% も無駄になるのか#

3D の圧縮性 RANS を回したことがあるなら、この交換は身についています。一辺を2倍に粗くすると、1ステップのコストは8分の1になります。 $2^3=8$ だからです。その代わり衝撃波・せん断層・剥離の傾きまで一緒に潰れます。その細い帯に合わせて領域全体を細かくすると、自由流のセルまで小さくなります。メモリと時間が同時に膨らみます。

必要な帯だけ細かくすれば足ります。比は次の見積もりです。

\frac{N_\text{AMR}}{N_\text{uniform}} \approx \frac{V_\text{feature}}{V_\text{domain}} \cdot 4^{L} + (1 - \tfrac{V_\text{feature}}{V_\text{domain}})

$L$ は最大の細分化レベルです。 $V_\text{feature}$ は細かく解く体積、 $V_\text{domain}$ は領域全体です。 $4^{L}$ は、2次元で一段進むとセルが4倍になるという意味です。 $L=5$ なら $4^5=1024$ です。冒頭の1%を体積比として置くと、最も粗い格子に対するセル数は $0.01\times 1024+0.99=11.23$ 倍になります。同じレベルの一様格子は粗い格子の1024倍です。二つを割ると、AMR は一様な細格子の約1.1%です。衝撃波のように細かく解く体積が $h$ 程度まで薄いと、この割合はさらに下がります。一桁パーセントというのは、この割り算から出ます。3次元で一辺だけ半分にしてもセルは8倍です。その8倍を領域全体に掛けると費用は丸ごと増え、帯の中だけに掛けると8倍はそこで止まります。

h-, p-, r-: 適応化の三つの流派#

適応格子（adaptive mesh、解の分布に応じてセルサイズや次数が変わる格子）は、何を動かすかで分かれます。

h-adaptivity: セルサイズ $h$ を縮めます。いちばんよく使います。AMR と言えば普通はこちらです。
p-adaptivity: セルはそのままで、多項式次数 $p$ だけ上げます。DGM や SEM（spectral element method）でよく見ます。
r-adaptivity (relocation): セル数も次数も固定し、節点だけを勾配の大きい側へ動かします。メモリは一定です。格子の歪みは残ります。

ここから追うのは h-AMR、とくに 2D クアッドツリー です。3D では同じ考えがオクツリーになります。

クアッドツリー — セルが子を持つデータ構造

均一格子は配列一つで足ります。適応格子には木が必要です。セルを一つ refine すると子が四つ生まれ、子は親を指します。子がもう一度割れると孫です。

class QuadCell:
    __slots__ = ("cx", "cy", "size", "level", "children", "parent", "data")
 
    def __init__(self, cx, cy, size, level=0, parent=None):
        self.cx = cx           # セル中心 x
        self.cy = cy           # セル中心 y
        self.size = size       # 一辺の長さ
        self.level = level     # ルート=0、細分化のたびに +1
        self.children = None   # leaf なら None、細分化後は長さ 4 のリスト
        self.parent = parent
        self.data = 0.0        # 保存量(密度など)
 
    def is_leaf(self):
        return self.children is None
 
    def split(self):
        h = self.size * 0.5
        q = self.size * 0.25
        offsets = ((-q, -q), (+q, -q), (-q, +q), (+q, +q))
        self.children = [
            QuadCell(self.cx + dx, self.cy + dy, h, self.level + 1, self)
            for dx, dy in offsets
        ]
        for ch in self.children:
            ch.data = self.data    # 親値で初期化(保存性)

split で子の初期値を親と同じにするのは保存のためです。親の積分 $U \cdot V$ は、四つの子の和と等しくなければなりません。子一つの体積は親の4分の1です。値をそのまま写すと、四つの和は $4 \times (U \cdot V/4) = U \cdot V$ で親と一致します。そのままの複製は1次の補間です。実務では親の傾き（slope）を載せる MUSCL 型が、2次の prolongation としてよく付きます。傾きを載せても、四つの和は親の積分と一致したままです。そこがずれると、割ったセルは質量を作るか消します。

どこを切るか — 細分化基準

クアッドツリーは空の箱です。どのセルを切るかが細分化基準です。実務でよく出るのは五つです。

基準	形式	捕まえるもの
1 階微分	$h \lvert \nabla u \rvert$	衝撃波、せん断層
2 階微分	$h^2 \lvert \nabla^2 u \rvert$	曲率の大きな領域
Hessian Frobenius	$\sqrt{\sum H_{ij}^2}\, h^2$	異方性検出
Richardson 推定	$\lvert u_h - u_{2h} \rvert$	真の離散化誤差
物理しきい値	$\rho > \rho_, \lvert \omega \rvert > \omega_$	渦コア

OpenFOAM の dynamicRefineFvMesh は、ユーザが選んだ場にしきい値を掛けます。AMReX はセルごとの Boolean である tagcells() を空けておき、呼び出した側が埋めます。どこに印を付けるかが木の形を変えます。入れ物そのものより、こちらの方が結果を分けます。セル数を減らしながら同じ答えを残す印の付け方は、タグ付けと coarse-fine のリフラックスへ続きます。

ブラウザの格子は表の一行目だけを使います。セル中心の中心差分 $|\nabla u|$ に、セルサイズ $h$ を掛けます。積が $\tau$ を超えれば split、そうでなければ leaf です。

2:1 バランスと hanging node#

子と隣が会う面で問題が起きます。片側が大きなセル、もう片側が小さなセル二つだと、その面に hanging node（宙吊りの節点）が出ます。有限体積では、この節点がフラックスの収支を壊します。

通り道は二つです。

2:1 バランスを強制する: 隣り合う leaf のレベル差を1以下に抑えます。p4est、AMReX、Berger 系はこの方針です。
non-conforming 補間: hanging face に補間を別に置き、小さなセル二つの値から大きなセルの面値を作ります。DGM 系でよく使います。

強制するループは短いです。refine 候補を受け取り、隣にレベルが二つ以上低い leaf があれば、その隣もキューへ入れます。キューが空になるまで繰り返します。

def enforce_2to1(roots, refine_queue):
    """refine_queue 内の全セルが隣接と 2:1 を満たすように補強"""
    pending = list(refine_queue)
    queued = set(id(c) for c in pending)
    while pending:
        cell = pending.pop()
        for nb in face_neighbors(cell, roots):
            if nb.is_leaf() and cell.level - nb.level >= 1:
                if id(nb) not in queued:
                    pending.append(nb)
                    queued.add(id(nb))
    return [c for c in (refine_queue + pending) if c.is_leaf()]

face_neighbors は同じ面を共有する leaf を返します。木の探索なので $O(\log N)$ です。

Python で見る 1D 適応格子#

2D クアッドツリーを全部書くと長くなります。同じ考えを 1D 二分木 で一段だけ動かします。ガウスパルスを載せ、しきい値を超えた区間だけ二等分します。

import numpy as np
 
def gaussian_pulse(x, x0=0.6, sigma=0.04):
    return np.exp(-((x - x0) ** 2) / (2 * sigma ** 2))
 
def refine_indicator(cell_lo, cell_hi, fn, h_scale=1.0):
    """スケール済み 1 階微分 = セル両端値の差の絶対値"""
    u_lo = fn(cell_lo)
    u_hi = fn(cell_hi)
    return abs(u_hi - u_lo) * h_scale
 
def adaptive_bisect_1d(fn, tau=0.05, max_level=8):
    """しきい値基準でルート [0, 1] を再帰的に二等分"""
    leaves = []
    def recurse(lo, hi, level):
        ind = refine_indicator(lo, hi, fn)
        if level >= max_level or ind < tau:
            leaves.append((lo, hi, level))
            return
        mid = 0.5 * (lo + hi)
        recurse(lo, mid, level + 1)
        recurse(mid, hi, level + 1)
    recurse(0.0, 1.0, 0)
    return leaves
 
cells = adaptive_bisect_1d(gaussian_pulse, tau=0.05, max_level=7)
print(f"leaf cells: {len(cells)}")
print(f"均一レベル 7 比: {len(cells) / 2**7 * 100:.1f}%")
# 出力例: leaf cells: 41, 均一レベル 7 比: 32.0%

ガウス一つは領域の細い帯しか占めません。leaf は41枚、一様レベル7に対して32.0%です。パルスが二つなら深いのは二つの核だけで、あいだは粗いままです。

41枚と128セルのあいだ#

出力の41を $2^7=128$ で割ると 0.3203 です。コメントの32.0%は、この値をパーセントで小数一位まで切ったものです。128から41を引くと87個です。レベル7まで一様に敷いたときにできるセルの約68%を、その大きさでは作っていないことになります。コードの $\sigma=0.04$ は区間 $[0,1]$ の長さの4%です。傾きが $\tau=0.05$ を超える場所はその近傍に限られ、平らな自由流は両端の値の差がほぼ0なので、二等分は早く止まります。1D は一段で子が二つです。2D クアッドツリーは四つなので、同じ細い帯でもセル数の比は先の $4^{L}$ の側へより急になります。 $\tau$ を 0.4 まで上げると、0.05 と 0.4 のあいだにあった裾は、それ以上割れません。

並列の落とし穴 — load balancing と SFC#

分散メモリでは別の壁に当たります。あるランクが深い木ばかり持つと、そのランクが全体の wall clock を決めます。セル数のバランス（load balancing）が崩れた状態です。

木を1次元に並べて切れば足ります。空間充填曲線（space-filling curve, SFC） — Morton（Z-order）か Hilbert — でインデックスを作ると、空間で隣り合うセルが1次元でも近くに残ります。p4est の Morton インデックスは、座標のビットを交互に挟みます。

def morton_2d(ix, iy, level):
    """整数座標 (ix, iy) を Morton(Z-order)インデックスに"""
    code = 0
    for b in range(level):
        code |= ((ix >> b) & 1) << (2 * b)
        code |= ((iy >> b) & 1) << (2 * b + 1)
    return code
 
# 4x4 格子の Z-order
for iy in range(4):
    for ix in range(4):
        print(morton_2d(ix, iy, 2), end=" ")
    print()
# 0  1  4  5
# 2  3  6  7
# 8  9 12 13
# 10 11 14 15

Morton インデックスで並べてからランク数で切ると、隣のセルは同じランクに集まり、通信量は減ります。正方領域の角ではジャンプが出ます。Hilbert 曲線はそのジャンプを避けますが、インデックス計算はより高くつきます。実務では Morton が事実上の標準です。

印の場を外すと32%でも高くつく#

木の形はほぼ決まっています。子は四つで、Morton でランクを切れば隣は同じランクに集まります。4×4 で順番が 0, 1, 4, 5 から始まるのは、座標ビットを一段ずつ挟んだ結果です。費用が漏れるのはそのあとです。印を密度に掛けておき、誤差がせん断層にあると、レベル7まで降りても41枚の leaf は見当違いの帯に集まります。32%は $x_0=0.6$ 、 $\sigma=0.04$ に基準が合っていたときの割合です。基準の場が変わると式の $V_\text{feature}$ が変わり、粗い格子に対する11.23倍も一緒に動きます。

2:1 を外す選択もあります。隣のレベル差を1で止めると、hanging node は一面につき小さなセル二つまでです。その制限を外すと、面ごとに non-conforming 補間を別に書くことになります。一面の流束が親の積分 $U\cdot V$ からずれやすくなります。空間を割り終えても、refine を時間前進のどの行に置くかは残ります。その一行が誤差をどこまで増やすかは、LBM AMR の再帰時間前進に書いてあります。