셀의 18%로 같은 답을 낸다 — AMR 태깅 기준과 coarse-fine 리플럭싱

Löhner 센서와 버퍼 폭, 그리고 flux register가 막는 계면 질량 누출

5천만 셀짜리 해석에서, 답을 실제로 결정하는 셀은 몇 개인가? 충격파 한 장, 전단층 한 겹, 화염면 한 장. 세어 보면 대개 전체의 몇 퍼센트다. 나머지는 매끄러운 영역을 매끄럽게 계산하느라 도는 셀이다. 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다. 즉 격자를 두 배 세밀하게 만들 때 균일 격자는 셀이 8배 늘지만 AMR은 4배만 는다. 해상도를 올릴수록 이득이 커진다는 뜻이고, 이것이 AMR을 쓰는 유일한 이유다.

어디를 쪼갤지는 구배가 정하지 않는다#

가장 먼저 손이 가는 태깅 기준은 $|\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셀이 필요하다. 이 관계는 레벨과 무관하다 — subcycling을 쓰면 $\Delta t_\ell$ 과 $h_\ell$ 이 같은 비율로 줄기 때문이다.

아래 시뮬레이션에서 직접 조작해보자.

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를 0으로 내리고 regrid every를 30까지 밀면, 계면을 따라 빨간 셀이 자란다. 센서는 "여기 쪼개야 한다"고 말하는데 격자는 아직 안 따라온 셀들이다. n_buf를 2로 올리면 빨강이 사라지고 대신 leaf 셀 수가 올라간다. 버퍼를 넓히는 비용이 정확히 이 숫자다.

블록·셀·패치 — 자르는 단위가 정하는 것#

같은 태그 결과라도 어떤 단위로 쪼개느냐에 따라 셀 수와 코드 복잡도가 갈린다.

방식	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이 커터(hexahedron 하나를 8개로), refinementHistory가 되돌리기 위한 이력이다. patch 방식은 둘의 중간이고, 직사각형 안에서는 구조격자 루프를 그대로 돌릴 수 있어 벡터화에 유리하다.

계면 한 면에 답이 두 개 있다#

여기서부터가 진짜 문제다. 레벨 $\ell$ 과 $\ell+1$ 이 만나는 면을 보자. 굵은 셀 쪽에서는 그 면이 하나고, 가는 셀 쪽에서는 $r^{D-1}$ 개다( $r$ 은 refinement ratio, 보통 2). 게다가 subcycling을 쓰면 시간 스텝도 다르다.

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

굵은 격자가 한 스텝 갈 동안 가는 격자는 $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 안쪽 셀은 손대지 않는다 — 이미 가는 격자 값을 평균해 덮어썼기 때문이다. 이 한 줄이 전체 도메인의 질량을 기계 정밀도까지 되돌린다.

아래에서 리플럭싱을 켠 쪽과 끈 쪽을 나란히 돌려 보자.

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에 붙어 있다. pulse sigma를 줄여 펄스를 날카롭게 하면 빨간 이탈폭이 커지는데, 장부에 적히는 $\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다.

이제 리플럭싱. 2D 스칼라 이류에 두 레벨을 얹고, subcycling까지 넣은 뒤 전체 질량을 잰다.

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):                                   # subcycling: 굵은 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

60스텝에 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가 이긴다.

AMR을 켜기 전에 정할 세 가지#

하나. 태깅 기준은 무차원으로 만든다. 정규화된 2차 차분이면 임계값 하나가 모든 필드·모든 레벨에 통한다. 원시 구배를 쓰면 레벨마다 튜닝하게 된다.

둘. 버퍼 폭은 $\nu N_{\rm regrid}$ 로 계산해서 넣는다. 감으로 1셀 넣고 regrid를 20스텝마다 하면, 계산의 절반은 특징이 patch 밖으로 새어 나간 상태로 돈다.

셋. 리플럭싱은 옵션이 아니다. 보존형 스킴을 쓰고 있어도 coarse-fine 계면은 보존을 깬다. 연소나 다상유동처럼 질량이 곧 답인 문제에서는 flux register 없이 나온 수렴 그래프를 믿을 수 없다.