쿼드트리 AMR — 격자가 스스로 자라게 두는 법

τ를 0.4로 올리면 꼬리 셀은 더 이상 갈라지지 않고, 2:1을 빼면 그 면의 플럭스가 먼저 어긋난다.

1984년, Stanford. Marsha Berger와 Joseph Oliger는 충격파 한 줄을 쫓고 있었다. 셀은 100만 개. 그중 충격파 근처 1%만 의미 있는 일을 했다. 나머지 99%는 균일한 자유류(free-stream)를 같은 정밀도로 다시 풀고 있었다. 그해 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를 돌려본 사람은 이 교환을 안다. 한 변을 두 배로 굵게 하면 스텝당 비용은 8분의 1이 된다. $2^3=8$ 이라서다. 대신 충격파, 전단층(shear layer), 박리(separation)의 기울기도 같이 뭉개진다. 그 띠에 맞춰 도메인 전체를 잘게 쪼개면 자유류 셀까지 작아진다. 메모리와 시간이 같이 터진다.

필요한 띠만 잘게 두면 된다. 비율은 이 식이다.

\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$ 은 최대 refinement 레벨이다. $V_\text{feature}$ 는 잘게 풀 부피, $V_\text{domain}$ 은 도메인 전체다. $4^{L}$ 은 2차원에서 한 단계가 셀을 네 배로 늘린다는 뜻이다. $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, refine 한 번에 +1
        self.children = None   # leaf면 None, refine되면 길이 4 리스트
        self.parent = parent
        self.data = 0.0        # 보존량(density 등)
 
    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으로 자주 붙는다. 기울기를 얹어도 네 자식의 합은 부모 적분과 같아야 한다. 그 합이 어긋나면 갈라진 셀은 질량을 만들거나 지운다.

어디를 자를지 결정하기 — refinement criterion#

쿼드트리는 빈 통이다. 어느 셀을 자를지가 refinement criterion이다. 실무에서 자주 보이는 것은 다섯 가지다.

기준	수식	잡는 것
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 balance를 강제: 인접한 두 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 curve — 으로 인덱스를 만들면, 공간에서 붙은 셀이 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 curve는 그 점프를 피하지만 인덱스 계산이 더 비싸다. 실무에서는 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의 재귀 시간전진에 적혀 있다.