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

Berger–Oliger 이래 40년의 적응 격자 자료구조와 refinement 기준 정리

1984년, Stanford의 Marsha Berger와 Joseph Oliger는 충격파를 추적하는 일에 매달려 있었다. 그들의 코드는 셀 100만 개를 깔아두고도 충격파 근처 1%의 셀만 의미 있는 일을 했다. 나머지 99%는 균일한 자유류(free-stream)를 정밀하게 재계산하고 있었다. 두 사람이 그해 Journal of Computational Physics에 낸 논문이 Berger–Oliger AMR의 시작이었고, 같은 아이디어가 오늘 OpenFOAM의 dynamicRefineFvMesh, AMReX, Chombo, p4est 안에 살아 있다. 이 포스트는 그 핵심인 쿼드트리 자료구조와 refinement 기준을 따라가고, 1D 코드를 직접 돌린 뒤 브라우저에서 격자가 갈라지는 모습을 본다.

균일 격자가 90% 손해를 보는 이유#

3D 압축성 RANS를 한 번이라도 돌려본 사람은 안다. 격자를 절반으로 줄이면 스텝당 비용은 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}$ 은 전체 도메인 부피다. 충격파처럼 1차원 곡면 특이점이라면 $V_\text{feature}/V_\text{domain} \sim h$ 정도로 작아지므로, 레벨 5 균일 대비 셀 수가 한 자릿수 퍼센트로 떨어진다.

h-, p-, r-: 적응의 세 갈래#

적응 격자(adaptive mesh, 해의 분포에 따라 자체적으로 변하는 격자)는 어디를 바꾸느냐에 따라 셋으로 나뉜다.

h-adaptivity: 셀 크기 $h$ 를 줄인다. 가장 직관적이고 가장 흔히 쓰인다. AMR이라고 하면 보통 이 의미다.
p-adaptivity: 셀 크기는 그대로 두고 다항식 차수 $p$ 를 올린다. DGM·SEM(spectral element method) 계열에서 많이 본다.
r-adaptivity (relocation): 셀 개수와 차수를 그대로 두고 절점을 경사 큰 쪽으로 옮긴다. 메모리 변동이 없어 매력적이지만 격자 비틀림 관리가 까다롭다.

이 글은 가장 많이 쓰이는 h-AMR, 그중에서도 2D 쿼드트리(3D면 옥트리)에 집중한다.

쿼드트리: 격자가 자기 자식을 갖는 자료구조#

균일 격자는 1차원 배열만 있으면 된다. 적응 격자는 트리가 필요하다. 한 셀이 refine되면 네 개의 자식 셀이 생기고, 각 자식은 자기 부모를 가리킨다. 자식이 다시 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$ 가 네 자식의 합과 같아야 한다. 단순 복제는 정확도 1차의 보간이며, 실무에서는 부모의 미분(slope)을 함께 사용하는 MUSCL 형태 보간이 자주 쓰인다.

어디를 자를지 결정하기 — 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()를 호출자가 채우게 한다. 무엇을 기준으로 두느냐가 AMR 결과를 결정한다 — 격자 자체보다 더 중요하다.

이 포스트의 인터랙티브 viz는 가장 단순한 1번(스케일된 1차 미분)을 쓴다. 셀 중심에서 중앙 차분으로 $|\nabla u|$ 를 구하고 셀 크기 $h$ 를 곱한다. 이 곱이 threshold $\tau$ 를 넘으면 split, 아니면 leaf.

2:1 균형과 hanging node#

자식 셀과 이웃이 만나는 면에서 문제가 생긴다. 한쪽은 큰 셀, 다른 쪽은 작은 셀 두 개로 나뉘어 있다면 그 면에는 hanging node(매달린 절점)가 생긴다. 유한체적에서 hanging node는 플럭스 평형을 깨뜨린다.

해결은 두 가지다.

2:1 balance를 강제: 인접한 두 leaf의 레벨 차이를 1 이하로 제한한다. p4est, AMReX, Berger 계열이 모두 이 정책을 쓴다.
non-conforming 보간으로 처리: hanging face에 별도의 보간식을 두고 작은 셀 둘의 값으로 큰 셀 면값을 만든다. DGM 계열에서 일반적이다.

2:1 balance를 강제하는 알고리즘은 짧다. refine 후보 리스트를 받아, 각 후보의 이웃 중 레벨이 두 단계 이상 작은 것이 있으면 그 이웃도 함께 refine 큐에 넣는다. 이를 큐가 빌 때까지 반복.

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%

가우시안 한 개는 도메인의 좁은 영역만 차지하므로 균일 격자 대비 30% 안팎으로 셀 수가 떨어진다. 펄스가 둘이면 두 영역만 깊게 들어가고 그 사이는 거칠게 남는다 — 바로 AMR이 약속한 그림이다.

직접 격자가 갈라지는 걸 보자#

아래 시뮬레이션에서 직접 조작해보자. 흰 점은 가우시안 특이점이고 — 클릭해서 끌면 격자 분할이 따라온다.

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.

threshold를 0.05까지 내리면 가우시안 꼬리까지 격자가 따라가고 leaf 수가 빠르게 늘어난다. 0.4 위로 올리면 거의 균일 거친 격자만 남는다. add feature로 둘째 특이점을 추가하면 두 코어만 깊게 풀고 그 사이는 거칠게 남는 것을 확인할 수 있다 — 동일한 균일 격자 대비 셀 수 비율이 viz 좌상단에 표시된다.

병렬화의 함정 — load balancing과 SFC#

분산 메모리에서 AMR은 또 다른 함정을 만난다. 한 랭크가 깊은 트리만 잡으면 그 랭크가 전체 시뮬레이션의 wall clock을 결정한다. 셀 수의 균형(load balancing)이 깨진 것이다.

해법은 트리 노드를 1차원으로 줄지어 자르는 것이다. 공간 채움 곡선(space-filling curve, SFC) — Morton(Z-order)이나 Hilbert curve — 으로 셀 인덱스를 만들면 트리의 공간적 인접성이 1D 인접성으로 보존된다. 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 인덱스로 셀을 정렬한 뒤 랭크 수만큼 잘라 분배하면, 인접한 셀은 같은 랭크에 모이고 통신량이 최소화된다. 단점은 정사각 도메인의 모서리에서 점프(jump)가 생긴다는 것 — Hilbert curve는 그 점프를 피하지만 인덱스 계산이 더 비싸다. 실무에서는 Morton이 사실상 표준이다.

다시 펴들 때 기억할 것#

AMR의 가치는 격자 데이터 구조가 아니라 refinement criterion에서 나온다. 데이터 구조는 트리·SFC가 거의 표준이지만 무엇을 잡을지는 문제마다 다르다.
2:1 balance를 강제하면 hanging node 처리가 단순해지고, 강제하지 않으면 non-conforming 보간을 모든 면에 짜야 한다.
병렬에서는 트리 자체보다 load balancing이 wall clock을 결정한다. Morton/Hilbert SFC + 동적 재분배가 사실상의 표준이다.