ρ와 u는 자릿수까지 같은데 변형률만 45% 부푼다 — LBM 격자 세밀화의 비평형 재스케일

레벨 경계를 그대로 건너가도 되는 것은 ρ와 u뿐이다. 비평형 부분은 τ_f/(m·τ_c) 배로 다시 써야 한다.

ρ도 u도 자릿수까지 같았는데 응력만 부풀었다#

격자 볼츠만(LBM) 코드에 적응 격자 세밀화(AMR)를 붙이면 레벨 경계라는 새 물건이 생긴다. 경계 한 층에서 굵은 격자와 세밀 격자가 겹친다. 이 층에서는 값을 옮겨 적어야 한다.

옮겨 적기를 검증할 때 보통 밀도와 속도를 본다. 그런데 그것만 보면 버그가 통과한다. 분포함수 $f_i$ 를 세밀 레벨에 그대로 복사해도 $\rho$ 와 $\mathbf{u}$ 는 소수점 끝자리까지 같다. 같지 않은 것은 응력이다. 이 글은 그 차이가 정확히 몇 배인지, 왜 하필 그 배수인지를 다룬다.

결론부터 적으면 배수는 $\tilde\tau_f/(m\tilde\tau_c)$ 의 역수다. 흔한 설정에서 1.45배, 고레이놀즈 설정에서는 2배에 가깝다.

레벨을 넘어도 같아야 하는 값과 반드시 달라져야 하는 값#

레벨 사이 변환의 원칙은 하나다. 보존량은 같아야 한다. $\rho$ , $\mathbf{u}$ , $p$ 는 물리량이므로 격자를 바꾼다고 달라질 이유가 없다.

세밀화 비율을 $m = \Delta x_c/\Delta x_f$ 로 두자. 아래 첨자 $c$ 는 굵은 격자, $f$ 는 세밀 격자다. 음속 스케일링(acoustic scaling)을 쓰면 $\Delta t$ 도 같은 비율로 줄어든다.

\Delta x_f = \frac{\Delta x_c}{m}, \qquad \Delta t_f = \frac{\Delta t_c}{m}

이러면 격자 단위 속도 $\hat{u} = u\,\Delta t/\Delta x$ 가 두 레벨에서 같다. $\rho$ 도 같으니 평형분포 $f_i^{\mathrm{eq}}(\rho, \hat{u})$ 는 두 레벨에서 완전히 같은 값이다. 바꿔야 하는 것은 평형에서 벗어난 부분, 즉 $f_i^{\mathrm{neq}} = f_i - f_i^{\mathrm{eq}}$ 뿐이다.

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

tau_c0.80m2

Both lattices start from the same sine and advance over the same physical time. With tau rescaled the orange curve sits on the blue one; press tau copied and the fine level suddenly relaxes a different fluid. At step 0 the amplitudes are 1.0000 and 1.0000; the measured flux ratio is 0.0000 against the predicted r = 0.6875. Drag tau_c towards 0.51 and r falls to 1/2 — that is where copying a population costs you a factor of two.

같은 물리 영역·같은 물리 시간을 두 개의 실제 격자가 나란히 푼다. tau_c 슬라이더를 움직이면서 오른쪽 사다리의 measured q_f/q_c가 rescale r와 붙어 있는지 보라. tau copied 버튼을 누르면 세밀 레벨이 다른 점성의 유체가 되어 주황 곡선이 파란 곡선에서 떨어진다.

점성을 고정하면 τ가 따라 움직인다#

물리 점성은 레벨이 바뀐다고 달라지면 안 된다. 격자 단위 완화시간 $\tilde\tau$ 로 쓴 점성은 이렇다.

\nu = c_s^2\left(\tilde\tau - \tfrac{1}{2}\right)\frac{\Delta x^2}{\Delta t}

$c_s^2$ 는 격자 음속의 제곱, $\Delta x^2/\Delta t$ 는 확산계수의 단위다. $-\tfrac{1}{2}$ 는 사다리꼴 적분이 남긴 자리로, LBM 이산화가 남긴 Δt/2에서 따로 다뤘다.

음속 스케일링에서 $\Delta x^2/\Delta t$ 는 $1/m$ 배가 된다. $\nu_f = \nu_c$ 를 걸면 $\tilde\tau$ 가 끌려 나온다.

\tilde\tau_f = \frac{1}{2} + m\left(\tilde\tau_c - \frac{1}{2}\right)

$m=2$ , $\tilde\tau_c = 0.8$ 이면 $\tilde\tau_f = 1.1$ 이다. 두 배가 아니라는 점이 중요하다. $\tilde\tau$ 를 그냥 $m$ 배 하면 $\nu$ 가 어긋난다. $\tfrac{1}{2}$ 가 스케일링 대상이 아니기 때문이다.

비평형 부분에는 τ 말고 Δt도 붙는다#

Chapman–Enskog 전개의 1차 항을 Grad 근사로 쓰면 비평형 부분은 변형률에 비례한다.

f_i^{\mathrm{neq}} = -\frac{w_i \rho \tilde\tau}{c_s^2}\left(\mathbf{c}_i\mathbf{c}_i - c_s^2\mathbf{I}\right) : \hat{\mathbf{S}}

$w_i$ 는 가중치, $\mathbf{c}_i$ 는 격자 속도, $\hat{\mathbf{S}}$ 는 격자 단위 변형률이다. 물리 변형률 $\mathbf{S}$ 와의 관계는 $\hat{\mathbf{S}} = \mathbf{S}\,\Delta t$ 다. 같은 지점의 물리 변형률은 레벨과 무관하므로, $f^{\mathrm{neq}}$ 는 $\tilde\tau\,\Delta t$ 에 비례한다.

r \equiv \frac{f_i^{\mathrm{neq},f}}{f_i^{\mathrm{neq},c}} = \frac{\tilde\tau_f\,\Delta t_f}{\tilde\tau_c\,\Delta t_c} = \frac{\tilde\tau_f}{m\,\tilde\tau_c}

$r$ 이 재스케일 계수다. Dupuis와 Chopard가 정리한 형태와 같다. $\tilde\tau_c = 0.8$ , $m=2$ 이면 $r = 1.1/1.6 = 0.6875$ 다. 세밀 레벨의 비평형은 굵은 레벨의 69%다.

세 가지 전달 방식을 같은 표에 올린다#

레벨 경계에서 값을 넘기는 방법은 실무에서 셋으로 갈린다.

방식	넘기는 것	필요한 정보	비용	실패 지점
전체 복사	$f_i$ 그대로	없음	최저	$\rho,\mathbf{u}$ 는 맞고 응력이 $1/r$ 배로 부푼다
거시변수 + Grad 재구성	$\rho, \mathbf{u}$ 보간 후 $f^{\mathrm{neq}}$ 재계산	속도 구배	중간	구배를 유한차분으로 다시 떠야 한다
비평형 보간 + 스케일	$f^{\mathrm{eq}}$ 는 재계산, $f^{\mathrm{neq}}$ 는 $r$ 배	$\tilde\tau_c, \tilde\tau_f, m$	최저 수준	$r$ 의 $m$ 을 빠뜨리기 쉽다

두 번째와 세 번째는 결과가 같아야 한다. 셋 중 세 번째가 실무 표준인 이유는 구배를 다시 뜰 필요가 없어서다. $\tilde\tau$ 두 개와 $m$ 만 있으면 계수가 나온다.

공간 보간 자체는 평범하다. 굵은→세밀은 3차원에서 삼선형(trilinear) 보간, 세밀→굵은은 $2^d$ 개 셀의 평균이다. $d$ 는 차원 수다. 어려운 부분은 보간이 아니라 보간한 뒤 무엇에 $r$ 을 곱하느냐다.

Python으로 되찾은 변형률을 재봤다#

D2Q9 격자에서 물리 변형률을 하나 정해두고, 두 레벨의 $f^{\mathrm{neq}}$ 를 각각 만든다. 그다음 세밀 레벨의 $\tilde\tau_f$ 로 변형률을 되찾는다. 재스케일한 값과 복사한 값을 나란히 넣어 본다.

CS2 = 1.0 / 3.0
EX = [0, 1, 0, -1, 0, 1, -1, -1, 1]
EY = [0, 0, 1, 0, -1, 1, 1, -1, -1]
W = [4/9, 1/9, 1/9, 1/9, 1/9, 1/36, 1/36, 1/36, 1/36]
 
 
def equilibrium(rho, ux, uy):
    """D2Q9 평형분포. 레벨이 달라도 rho, u가 같으면 같은 값이 나온다."""
    out = []
    u2 = ux*ux + uy*uy
    for i in range(9):
        cu = EX[i]*ux + EY[i]*uy
        out.append(rho*W[i]*(1 + cu/CS2 + cu*cu/(2*CS2*CS2) - u2/(2*CS2)))
    return out
 
 
def grad_neq(rho, tau, sxy):
    """Grad 근사로 만든 비평형 부분. sxy는 그 레벨의 격자단위 전단변형률."""
    out = []
    for i in range(9):
        q_xy = EX[i]*EY[i]            # Q_i 의 xy 성분 (대각은 sxy와 짝이 없다)
        out.append(-(W[i]*rho*tau/CS2) * 2.0 * q_xy * sxy)
    return out
 
 
def recover_strain(fneq, rho, tau, dt):
    """비평형 모멘트에서 물리 단위 전단변형률을 되찾는다."""
    pi_xy = sum(fneq[i]*EX[i]*EY[i] for i in range(9))
    return -pi_xy / (2.0*rho*CS2*tau*dt)
 
 
def tau_on_level(tau_c, m):
    """점성을 고정했을 때 세밀 레벨이 가져야 할 완화시간."""
    return 0.5 + m*(tau_c - 0.5)
 
 
def nu_physical(tau, dx, dt):
    return CS2*(tau - 0.5)*dx*dx/dt
 
 
rho, ux, uy = 1.0, 0.05, 0.0
s_phys = 0.004          # 물리 전단변형률 — 레벨과 무관한 값
m = 2                   # 세밀화 비율
dx_c, dt_c = 1.0, 1.0
dx_f, dt_f = dx_c/m, dt_c/m
 
feq_c = equilibrium(rho, ux, uy)
feq_f = equilibrium(rho, ux, uy)
print("max |feq_c - feq_f| = %.3e" % max(abs(a-b) for a, b in zip(feq_c, feq_f)))
print()
print("%5s  %6s  %8s  %8s  %7s  %11s  %14s  %6s" % (
    "tau_c", "tau_f", "nu_c", "nu_f", "r_meas", "tf/(m*tc)", "copied/true", "err%"))
for tau_c in [0.51, 0.55, 0.60, 0.80, 1.20, 2.00]:
    tau_f = tau_on_level(tau_c, m)
    fneq_c = grad_neq(rho, tau_c, s_phys*dt_c)
    fneq_f = grad_neq(rho, tau_f, s_phys*dt_f)
    r = fneq_f[5]/fneq_c[5]
    s_ok = recover_strain(fneq_f, rho, tau_f, dt_f)
    s_bad = recover_strain(fneq_c, rho, tau_f, dt_f)
    print("%5.2f  %6.3f  %8.5f  %8.5f  %7.4f  %11.4f  %14.4f  %6.1f" % (
        tau_c, tau_f, nu_physical(tau_c, dx_c, dt_c), nu_physical(tau_f, dx_f, dt_f),
        r, tau_f/(m*tau_c), s_bad/s_ok, (s_bad/s_ok - 1)*100))
 
print()
worst = grad_neq(rho, 0.51, s_phys)
print("mass moment of f^neq     = %.3e" % sum(worst))
print("momentum moments of f^neq = %.3e, %.3e" % (sum(worst[i]*EX[i] for i in range(9)),
                                            sum(worst[i]*EY[i] for i in range(9))))

max |feq_c - feq_f| = 0.000e+00
 
tau_c   tau_f      nu_c      nu_f   r_meas    tf/(m*tc)     copied/true    err%
 0.51   0.520   0.00333   0.00333   0.5098       0.5098          1.9615    96.2
 0.55   0.600   0.01667   0.01667   0.5455       0.5455          1.8333    83.3
 0.60   0.700   0.03333   0.03333   0.5833       0.5833          1.7143    71.4
 0.80   1.100   0.10000   0.10000   0.6875       0.6875          1.4545    45.5
 1.20   1.900   0.23333   0.23333   0.7917       0.7917          1.2632    26.3
 2.00   3.500   0.50000   0.50000   0.8750       0.8750          1.1429    14.3
 
mass moment of f^neq     = 0.000e+00
momentum moments of f^neq = 0.000e+00, 0.000e+00

세 가지가 한 번에 나온다. 평형분포는 두 레벨에서 완전히 같다. $\nu_c$ 와 $\nu_f$ 열이 자릿수까지 일치한다. 측정한 $r$ 은 닫힌 형태 $\tilde\tau_f/(m\tilde\tau_c)$ 와 같다.

마지막 두 줄이 이 글의 제목이 나온 자리다. $f^{\mathrm{neq}}$ 의 0차·1차 모멘트가 정확히 0이다. 복사를 하든 재스케일을 하든 $\rho$ 와 운동량은 흔들리지 않는다. 틀리는 곳은 2차 모멘트 하나뿐이다.

τ가 0.5에 붙을수록 복사의 대가는 두 배로 간다#

오차 열을 위에서 아래로 읽으면 방향이 보인다. $\tilde\tau_c = 2.0$ 에서는 14%다. $\tilde\tau_c = 0.51$ 에서는 96%다. 닫힌 형태로 쓰면 $r$ 의 범위가 나온다.

r = \frac{1/2 + m(\tilde\tau_c - 1/2)}{m\,\tilde\tau_c}

$\tilde\tau_c \to 1/2$ 이면 $r \to 1/2$ 이고, $\tilde\tau_c \to \infty$ 이면 $r \to 1$ 이다. 점성이 크면 복사해도 티가 안 난다. 점성이 작으면 응력이 두 배가 된다.

문제는 AMR을 붙이는 이유가 대개 고레이놀즈라는 점이다. $\tilde\tau$ 를 0.5 근처까지 내려서 쓴다. 복사 버그가 가장 크게 터지는 영역이 곧 AMR을 쓰고 싶은 영역이다. 증상도 헷갈린다. 질량과 운동량은 보존되니 발산하지 않고, 레벨 경계선을 따라 얇은 소용돌이 층만 남는다.

$m$ 을 키우면 더 나빠진다. $m=4$ , $\tilde\tau_c = 0.8$ 이면 $\tilde\tau_f = 1.7$ 이다. 따라서 $r = 1.7/3.2 = 0.531$ 이고 복사 오차는 88%로 뛴다.

레벨 하나를 얹는 값은 m^(d+1)로 청구된다#

재스케일은 두 번만 일어난다. 굵은 레벨에서 세밀 레벨로 내려가는 explosion, 세밀 레벨에서 굵은 레벨로 올라오는 coalescence다. 그 사이는 평범한 collide-and-stream이다.

아래 시계에서 한 주기를 한 단계씩 걸어보자.

m2dwell42

Use next phase to walk the cycle one gate at a time: explosion, 2 fine sub-steps, coalescence. The two dashed lines are the only moments populations cross levels — everything between them is ordinary collide-and-stream. Raise m or switch to d = 3 and the work factor climbs as m^(d+1) = 8; currently at cycle 0, phase explosion.

next phase를 눌러 explosion → $m$ 번의 서브스텝 → coalescence 순서를 확인하라. 초록·보라 점선이 분포함수가 레벨을 건너는 유일한 두 순간이다. m과 d를 올리면 오른쪽 장부의 work factor가 $m^{d+1}$ 로 오른다.

이 지수가 AMR 설계의 실제 제약이다. 3차원에서 $m=2$ 면 패치 하나가 16배 비싸진다. 그래서 세밀 레벨을 어디에 얼마나 둘지가 스킴 선택보다 성능을 크게 가른다.

그래서 오버랩 층은 왜 한 겹이어야 하는가#

오버랩 층은 두 레벨의 노드가 모두 존재하는 셀 한 층이다. 왜 한 층인가. 스트리밍이 한 스텝에 한 칸만 움직이기 때문이다. 세밀 레벨이 한 서브스텝을 돌 때 바깥에서 들어오는 분포함수는 딱 한 칸 밖에서 온다. 그 한 칸만 채워주면 된다.

이 관점은 경계 노드가 잃어버리는 분포함수와 같다. 벽이든 레벨 경계든, 먼저 하는 일은 "스트리밍 후 어느 $f_i$ 가 비는지" 세는 것이다. 벽에서는 형상이 그 답을 정하고, 레벨 경계에서는 세밀화 비율 $m$ 이 정한다.

셀 중심(cell-centered) 격자를 쓰면 오버랩 층이 노드 중심(node-centered)보다 다루기 쉽다. 겹치는 노드가 없어서 소유권이 명확하고, 병렬 분할에서 통신 대상이 셀 목록 하나로 떨어진다. 대신 굵은→세밀 보간에서 위치가 반 칸씩 어긋나므로 보간 스텐실을 그에 맞게 잡아야 한다.

AMR 인터페이스를 새로 짜든 남의 코드를 읽든, 세 줄이면 진단이 끝난다. $\tilde\tau_f$ 가 $\tfrac{1}{2} + m(\tilde\tau_c - \tfrac{1}{2})$ 인가. $f^{\mathrm{neq}}$ 에 $r$ 을 곱하는가. $r$ 의 분모에 $m$ 이 있는가. 세 번째가 가장 자주 빠진다.