방향을 쪼개지 않고 이류를 푼다 — Bell-Colella-Glaz 비분할 업윈드

한 바퀴 뒤 질량 오차는 둘 다 반올림인데 봉우리는 0.349와 0.806으로 갈린다. 보존 합격이 원판 생존은 아니다.

1980년대 후반, 로런스 버클리 연구소의 John Bell과 Phillip Colella는 폭발 시뮬레이션에서 이상한 자국을 발견했다. 격자 축에 비스듬히 놓인 화염면이 시간이 갈수록 사선 방향으로 번지고 비틀렸다. 범인은 물리가 아니라 이산화였다. 한 방향씩 번갈아 푸는 방향분할(dimensional splitting)이 대각선 이류를 일그러뜨린 것이다. 해법은 Bell-Colella-Glaz(BCG) 비분할 2차 업윈드 스킴이다. 방향을 쪼개면 회전 원판이 마름모가 되고, 횡방향 미분 한 항이 그 구멍을 메운다.

CFL을 0에서 0.9로 밀면 주황 점이 셀 평균 쪽으로 붙고, limiter를 none으로 두면 점프 옆 직선이 이웃 값을 넘어선다.

CFL = 0.50

Cyan = limited PLM slope inside each cell · Amber dot = predicted left edge state at n+1/2 (advecting velocity u>0). Pick none near the jump to see the slope overshoot the neighbor values.

0.9에서는 기울기 가중치가 거의 사라진다. minmod로 되돌리면 그 넘침이 잘린다.

다루는 방정식은 가장 단순한 보존형 이류식이다. 발산이 0인(divergence-free) 속도장 $\mathbf{u}$ 가 스칼라 $s$ 를 실어 나른다.

\frac{\partial s}{\partial t} + \nabla\cdot(\mathbf{u}\,s) = 0

여기서 $s$ 는 농도·온도 같은 수동 스칼라, $\mathbf{u}$ 는 미리 주어진 속도장이다. 유한체적(FV)으로 셀 평균을 갱신하면 다음 형태가 된다.

s_{ij}^{n+1} = s_{ij}^{n} - \frac{\Delta t}{\Delta x}\left(F^x_{i+\frac12,j} - F^x_{i-\frac12,j}\right) - \frac{\Delta t}{\Delta y}\left(F^y_{i,j+\frac12} - F^y_{i,j-\frac12}\right)

$F^x = u\,s^{n+1/2}$ 는 face에서의 시간중심(time-centered) 플럭스다. 모든 정확도는 결국 face 위 값 $s^{n+1/2}$ 를 얼마나 잘 추정하느냐에 달려 있다.

한 방향씩 쪼개면 회전이 비틀린다#

방향분할은 2D 문제를 x 방향 1D 업데이트 다음 y 방향 1D 업데이트로 쪼갠다. 구현은 쉽다. 1D 솔버 하나면 끝이다. 문제는 대각선 흐름이다.

스칼라가 $45^\circ$ 로 흐르면, x를 먼저 풀고 y를 푸는 순서와 그 반대 순서가 다른 답을 낸다. 두 번의 1D 스텝 사이에 끼는 횡방향 결합(cross term)을 빠뜨리기 때문이다. 회전 유동에서는 이 오차가 매 스텝 누적되어 원판이 마름모로 일그러진다. Strang 분할로 순서를 번갈아도 2차 오차의 방향 비대칭만 줄어들 뿐, 근본 원인은 남는다.

BCG의 처방은 간단하다. x도 y도 동시에 푼다. 한 스텝 안에서 모든 face 값을 같은 시간 $t^{n+1/2}$ 로 예측하고, 그 예측에 횡방향 변화를 미리 반영한다.

시간을 공간으로 바꾸는 예측자#

핵심 아이디어는 face 값을 Taylor 전개로 시간·공간 양쪽으로 반 스텝 밀어내는 것이다. 셀 $i$ 의 오른쪽 face에서 좌측 상태를 보면,

s^{n+1/2}_{i+\frac12,\,L} = s^n_{ij} + \frac{\Delta x}{2}\,\partial_x s + \frac{\Delta t}{2}\,\partial_t s

여기서 시간 미분을 그대로 두면 쓸모가 없다. 이류 방정식 $\partial_t s = -u\,\partial_x s$ 를 대입해 시간 미분을 공간 미분으로 바꾼다(Cauchy-Kowalewski 트릭).

s^{n+1/2}_{i+\frac12,\,L} = s^n_{ij} + \frac{1}{2}\left(1 - \frac{\Delta t}{\Delta x}\,u_{i+\frac12,j}\right)\Delta_x s_{ij}

$\Delta_x s_{ij}$ 는 셀 $i$ 의 제한된 기울기(limited slope), $u_{i+1/2,j}$ 는 face의 이류 속도다. 괄호 안 $(1 - \mathrm{CFL})$ 항이 시간 전진을 담는다. CFL이 1에 가까우면 좌측 상태의 기여가 줄고, 0이면 단순한 공간 외삽이 된다.

오른쪽 셀에서 본 우측 상태도 대칭으로 만든다.

s^{n+1/2}_{i+\frac12,\,R} = s^n_{i+1,j} - \frac{1}{2}\left(1 + \frac{\Delta t}{\Delta x}\,u_{i+\frac12,j}\right)\Delta_x s_{i+1,j}

두 후보 중 하나를 고르는 것은 Godunov의 업윈드 규칙이다. face 속도의 부호가 바람을 가른다.

s^{n+1/2}_{i+\frac12} = \begin{cases} s_{i+\frac12,\,L}, & u_{i+\frac12} > 0 \\[2pt] s_{i+\frac12,\,R}, & u_{i+\frac12} < 0 \end{cases}

앞 그림에서 CFL을 0.9로 올리면 주황 점이 평균 쪽으로 붙고, limiter를 none으로 두면 직선이 이웃 값을 넘는다. 그 차이가 위 식의 시간 전진 인수다.

횡방향 미분이 모서리를 잇는다#

방향분할이 빠뜨린 결합을 BCG는 예측자 안에 직접 집어넣는다. x-face의 좌측 상태에 횡방향 플럭스 발산을 반 스텝만큼 빼준다.

s^{n+1/2}_{i+\frac12,\,L} \;\mathrel{-}=\; \frac{\Delta t}{2}\,\frac{\partial (v\,s)}{\partial y}\bigg|_{ij}

$v$ 는 y 방향 속도, $\partial_y(v s)$ 는 셀 $ij$ 에서의 횡방향 이류 변화다. 이 항을 추정하려면 먼저 y-face 위의 1D 예측 상태 $\hat s^y$ 를 구해 두어야 한다. 즉 BCG는 두 단계다. 먼저 각 방향의 1D 업윈드 상태를 예측하고, 그 예측을 횡방향 보정에 재사용한 뒤, 같은 업윈드 규칙으로 최종 face 값을 고른다. 이 "모서리 결합(corner coupling)"이 대각선 이류를 비틀림 없이 옮긴다.

기울기에 재갈을 물리다 — minmod와 superbee#

예측자의 정확도는 기울기 $\Delta_x s$ 에서 나온다. 중심차분 기울기를 그냥 쓰면 2차 정확도지만, 불연속 근처에서 진동(overshoot)을 만든다. 그래서 경사 제한자(slope limiter)로 재갈을 물린다. 가장 보수적인 minmod는 두 한쪽 차분 중 작은 쪽을, 부호가 다르면 0을 택한다.

\Delta_x s_{ij} = \operatorname{minmod}\!\left(s_{ij}-s_{i-1,j},\; s_{i+1,j}-s_{ij}\right)

superbee는 같은 단조성을 지키면서 기울기를 최대한 살려 접촉면을 더 날카롭게 끈다. 위 1D 데모에서 minmod와 superbee를 번갈아 보면, superbee 쪽 PLM 직선이 더 가파른 것을 확인할 수 있다.

같은 재갈을 플럭스에 물리는 1차원 계산은 Lax–Wendroff에 TVD 리미터를 얹는 방법이다. BCG는 리미터를 플럭스가 아니라 셀 기울기에 건다.

numpy로 돌리는 회전 원판#

이제 2D 비분할 BCG를 돌려본다. 토이 문제는 강체 회전(solid-body rotation)이다. 중심을 도는 발산-0 속도장에 스칼라 원판을 띄우고, 1차 업윈드와 BCG를 같은 격자에서 비교한다.

import numpy as np
 
N = 64
dx = 1.0 / N
omega = 2 * np.pi                      # 단위 시간당 1회전
 
def face_velocities():
    j = np.arange(N)
    i = np.arange(N)
    u = -omega * ((j + 0.5) * dx - 0.5)        # x-face 속도 (y에만 의존)
    u = np.broadcast_to(u[None, :], (N, N)).copy()
    v = omega * ((i + 0.5) * dx - 0.5)         # y-face 속도 (x에만 의존)
    v = np.broadcast_to(v[:, None], (N, N)).copy()
    return u, v
 
def minmod_pair(a, b):                          # 단조성 보존 기울기
    return np.where(a * b > 0, np.where(np.abs(a) < np.abs(b), a, b), 0.0)
 
def limited_slopes(s):
    slx = minmod_pair(s - np.roll(s, 1, 0), np.roll(s, -1, 0) - s)
    sly = minmod_pair(s - np.roll(s, 1, 1), np.roll(s, -1, 1) - s)
    return slx, sly
 
def donor_fluxes(s, u, v):                       # 1차 업윈드 (도너 셀)
    Fx = np.where(u > 0, u * np.roll(s, 1, 0), u * s)
    Fy = np.where(v > 0, v * np.roll(s, 1, 1), v * s)
    return Fx, Fy
 
def bcg_fluxes(s, u, v, dt):                     # 비분할 2차 예측자
    slx, sly = limited_slopes(s)
    # 1) 법선 방향 1D 예측 상태
    sLx = np.roll(s, 1, 0) + 0.5 * (1 - dt/dx * u) * np.roll(slx, 1, 0)
    sRx = s - 0.5 * (1 + dt/dx * u) * slx
    hatX = np.where(u > 0, sLx, np.where(u < 0, sRx, 0.5 * (sLx + sRx)))
    sBy = np.roll(s, 1, 1) + 0.5 * (1 - dt/dx * v) * np.roll(sly, 1, 1)
    sTy = s - 0.5 * (1 + dt/dx * v) * sly
    hatY = np.where(v > 0, sBy, np.where(v < 0, sTy, 0.5 * (sBy + sTy)))
    # 2) 횡방향(모서리) 보정 — 예측 상태를 재사용
    divVy = (np.roll(v * hatY, -1, 1) - v * hatY) / dx
    divUx = (np.roll(u * hatX, -1, 0) - u * hatX) / dx
    sLx -= 0.5 * dt * np.roll(divVy, 1, 0)
    sRx -= 0.5 * dt * divVy
    sBy -= 0.5 * dt * np.roll(divUx, 1, 1)
    sTy -= 0.5 * dt * divUx
    # 3) 같은 업윈드 규칙으로 최종 face 값 선택
    sX = np.where(u > 0, sLx, np.where(u < 0, sRx, 0.5 * (sLx + sRx)))
    sY = np.where(v > 0, sBy, np.where(v < 0, sTy, 0.5 * (sBy + sTy)))
    return u * sX, v * sY
 
def advance(s, u, v, dt, scheme):
    Fx, Fy = bcg_fluxes(s, u, v, dt) if scheme == 'bcg' else donor_fluxes(s, u, v)
    return s - dt * ((np.roll(Fx, -1, 0) - Fx) / dx + (np.roll(Fy, -1, 1) - Fy) / dx)
 
def make_disk():
    X, Y = np.meshgrid((np.arange(N)+0.5)*dx, (np.arange(N)+0.5)*dx, indexing='ij')
    return ((X - 0.5)**2 + (Y - 0.78)**2 < 0.13**2).astype(float)
 
u, v = face_velocities()
Umax = omega * 0.5 * np.sqrt(2)
dt = 0.6 * dx / Umax                    # CFL = 0.6
steps = int(round(1.0 / dt))            # 정확히 한 바퀴
 
for scheme in ('upwind1', 'bcg'):
    s = make_disk()
    mass0 = s.sum()
    for _ in range(steps):
        s = advance(s, u, v, dt, scheme)
    print(f"{scheme:8s}  peak={s.max():.3f}  "
          f"mass_err={abs(s.sum()-mass0)/mass0:.2e}")

출력은 다음과 비슷하다.

upwind1   peak=0.349  mass_err=3.1e-16
bcg       peak=0.806  mass_err=2.8e-16

스킴 토글을 1st-order upwind에서 BCG unsplit으로 넘기면, 봉우리가 0.349에서 0.806으로 가는 과정이 칸 위에서 다시 돈다.

CFL = 0.60

A scalar disk rotates rigidly about the center. Switch to 1st-order upwind and watch the disk melt into a ring within one revolution; BCG unsplit keeps its edge far longer.

1st-order upwind로 두면 원판이 한 바퀴도 못 돌아 고리처럼 퍼지고 가장자리가 흐려진다. BCG unsplit으로 바꾸면 경계가 훨씬 오래 살아남는다. CFL을 0.9로 올리면 BCG도 약간 거칠어지지만, 1차 업윈드와의 격차는 그대로다.

두 스킴 모두 질량은 기계 정밀도까지 보존한다. 보존형 플럭스 차분이라 당연하다. 갈리는 것은 봉우리(peak)다. 1차 업윈드는 한 바퀴 만에 원판을 절반 이하로 뭉개고, BCG는 원래 높이의 80%를 지킨다.

두 바퀴를 돈 뒤 남은 것#

코드에 적힌 수로 한 바퀴를 환산하면 이렇다. N=64라 간격은 1/64다. 모서리 속력은 각속도 2π에 0.5와 루트 2를 곱한 약 4.44다. CFL 0.6이면 시간 간격은 약 0.00211이고, 시간 1을 채우는 스텝은 474다. 원판 중심은 (0.5, 0.78)이라 회전 중심 (0.5, 0.5)에서 0.28 떨어져 있다. 그 자리의 속력은 2π 곱하기 0.28, 약 1.76이다. 시간 간격을 잡는 것은 이 속력이 아니라 모서리의 4.44다.

반지름 0.13의 지름 0.26은 격자 16.6칸이다. 면적을 칸 면적으로 나누면 약 217칸이 스칼라를 들고 있다. 기울기 가중치만 보면, CFL 0.6에서 왼쪽은 0.2, 오른쪽은 0.8이다. 절반 곱하기 (1-0.6)과 (1+0.6)이다. CFL 상한 0.8에서는 왼쪽 가중치가 0.1까지 준다. 출력의 봉우리 0.349와 0.806은 이 가중치와 횡방향 항을 474번 쌓은 높이다.

질량 오차 3.1e-16과 2.8e-16은 둘 다 반올림이다. 갈린 것은 보존이 아니라 잃어버린 높이 0.651과 0.194다. 1차가 잃은 높이는 BCG의 약 3.4배다. 횡방향 보정은 한 스텝에 시간 간격의 절반, 약 0.00106을 횡방향 발산에 곱해 뺀다. 그 항은 한 바퀴에 474번 들어간다. 수치 확산은 칸 수 64가 아니라, 한 칸에서 기울기를 얼마나 남기느냐에서 갈린다.

옮길 때 새는 곳#

CFL은 분할 스킴보다 빡빡하다. 비분할 2D는 대각선 정보 전파 때문에 안정 한계가 1차원보다 작다. 실무에서는 $\mathrm{CFL}\lesssim 0.8$ 로 잡고 시작한다.
속도장은 반드시 face에 둔다. 셀 중심 속도를 face로 보간하면 발산-0이 깨지고, 보존된 줄 알았던 질량이 샌다. MAC 배치(staggered)가 정답이다.
횡방향 항을 빼먹지 말 것. 법선 예측만 구현하면 그냥 차원분할 MUSCL이다. 코드가 돌긴 하지만 회전 비틀림은 그대로 남는다.
limiter는 극값에서 0이 된다. 봉우리 꼭대기에서 기울기가 0으로 깎이면서 국소 1차로 떨어진다. 이게 봉우리가 조금씩 낮아지는 이유다. 더 살리려면 극값 보존 limiter를 검토한다.

질량이 맞는데 원판만 죽는 선택#

봉우리 0.349는 한 바퀴 뒤 1차 업윈드가 남긴 높이다. BCG는 0.806이다. 질량 오차는 3.1e-16과 2.8e-16이라, 보존 로그만 보면 둘 다 통과다. 원판을 보면 한쪽만 남는다.

대각으로 도는 스칼라에서 1차 업윈드를 쓸지는 이 차이로 정하면 된다. 접촉면 높이를 버려도 되면 도너 셀로 충분하다. 높이를 남기려면 횡방향 항을 같은 스텝에 넣어야 한다. 격자가 비정렬로 바뀌면 minmod를 그대로 가져가는 순간 극값 처리가 달라진다. Barth–Jespersen 리미터가 그 가짜 극값을 이웃 최솟값과 최댓값으로 어떻게 자르는지가 다음 결정이다.