벽거리를 O(N log N)에 푸는 법: Eikonal 방정식과 Fast Marching

벽거리는 전수 거리 쿼리가 아니라 도달 시간이다. 매 스텝 전 영역을 행진시킬지는 격자 변형이 정한다.

1996년, UC Berkeley의 James Sethian은 의료 영상에서 혈관 경계를 추적하는 알고리즘을 다듬고 있었다. 그 결과물인 Fast Marching Method는 지금 엉뚱한 곳에서 쓰인다. 항공기 날개 주위의 RANS 해석기가 매번 돌아갈 때, 바로 이 알고리즘이 각 격자점의 벽거리(wall distance)를 만들어낸다. Spalart-Allmaras 모델과 k-ω SST 모델은 격자마다 "가장 가까운 벽이 어디인가"라는 값을 입력으로 요구한다.

패턴은 단일 원통에서 시작한다. 속도 $F$ 를 올리면 같은 순간에 노란 띠가 더 멀리 나가고, 장애물을 슬릿으로 바꾸면 출발하는 벽의 모양이 바뀐다.

전파 속도 F =1.0프레임당 스텝40

진행 0%보라색 = 가까움, 노란색 = 멀리 · 벽에서부터 등거리 파면이 퍼진다

노란색 띠가 지금 업데이트되는 Considered 집합, 즉 "파면(front)" 자체다. 벽과 원통 경계에서 동시에 파면이 출발해, 서로 만나는 지점에서 거리가 확정된다. 두 원통 패턴에서는 두 파면의 충돌선이 바로 두 원통의 중간거리(medial axis) 근처에 형성된다. 슬릿 벽 패턴은 파면이 틈을 돌아 뒷면까지 도달하는 걸 보여준다. 이 "돌아 들어감"이 순진한 유클리드 거리와 갈라지는 지점이다. 실제 항공기·터빈 형상에서도 벽 뒤쪽의 난류 모델 입력은 바로 이 돌아 들어간 거리를 요구한다.

벽거리는 왜 난류 모델에 필요한가#

경계층 난류의 통계량은 벽에서의 거리 $y$ 에 강하게 의존한다. RANS(Reynolds-Averaged Navier-Stokes) 모델은 그 의존성을 대수식으로 끼워 넣는다. Spalart-Allmaras는 destruction 항에서 $d_w$ (벽거리)를 직접 사용한다. k-ω SST는 blending function $F_1$ , $F_2$ 안에 $d_w^2$ 가 들어간다. LES와의 혼합 모델(DES, DDES, IDDES)에서도 벽거리가 스위치 역할을 한다.

문제는 $d_w$ 가 메시의 기하량이라는 것이다. 메시가 움직이거나 변형되면 매 시간 스텝마다 다시 계산해야 한다. 회전 기계(컴프레서, 펌프)처럼 로터-스테이터 경계가 상대 운동을 하는 경우가 여기 해당한다.

순진한 방법의 함정#

바로 떠오르는 방법은 "모든 벽면과 모든 셀의 거리를 계산해 최솟값을 고른다"이다. 실제로 교과서 예제에서는 이걸 쓴다.

비용은 $O(N_{\text{cell}} \cdot N_{\text{wall}})$ 이다. 항공기 전체를 푼다고 해보자. 10M 셀, 0.5M 벽면이면 5조 번의 거리 계산이다. 이게 매 시간 스텝마다 돈다면 해석이 끝나지 않는다. KD-tree 같은 공간 자료구조로 $O(N_{\text{cell}} \log N_{\text{wall}})$ 까지 낮출 수 있지만, 복잡한 형상의 벽 근접 쿼리는 여전히 부담이다. 같은 전수 거리를 입자 접촉에서 자른 기록은 섹터 크기를 입자 하나에서 정하는 DEM에 있다. 벽거리와 대상은 다르고, 후보를 기하로 줄이는 이유는 같다.

Eikonal 방정식으로 생각을 바꾼다#

벽거리를 PDE의 해로 보자. "벽에서 0이고, 모든 방향으로 기울기 크기가 1"인 함수가 바로 거리장이다.

|\nabla d(\mathbf{x})| = \frac{1}{F(\mathbf{x})}, \qquad d|_{\Gamma_w} = 0

$d$ 는 벽까지의 "도달 시간", $F$ 는 전파 속도(기본 1), $\Gamma_w$ 는 벽 집합이다. $F=1$ 로 고정하면 $d$ 는 그대로 유클리드 거리다.

이 방정식은 쌍곡성(hyperbolic)이다. 정보가 특성선(벽에서 바깥으로 뻗는 직선)을 따라 한 방향으로만 흐른다. 이 단방향성이 효율적 알고리즘의 열쇠다. "이미 값이 정해진 점에서만 값이 들어온다"는 성질을 이용하면, 모든 점을 한 번씩만 방문해도 충분하다.

Fast Marching: Dijkstra의 연속판#

Sethian의 아이디어는 두 줄로 요약된다.

모든 격자점을 Accepted / Considered / Far 세 집합으로 나눈다.
Considered 중 도달 시간이 가장 작은 점을 Accepted로 확정하고, 그 이웃을 업데이트한다.

가장 작은 값을 빠르게 꺼내려면 min-heap이 자연스럽다. 그래프 최단 거리 알고리즘인 Dijkstra와 구조가 같다. 차이는 이웃 업데이트 방식이다. 그래프에서는 단순 덧셈이지만, Eikonal에서는 upwind 2차 방정식을 풀어 값을 얻는다.

2D 정렬 격자에서 업데이트 식은 다음과 같다.

(T - T_x)^2 + (T - T_y)^2 = \left(\frac{h}{F}\right)^2

$T_x$ , $T_y$ 는 각각 수평·수직 방향의 업윈드 이웃값(네 이웃 중 작은 쪽), $h$ 는 격자 간격이다. 두 근 중 $T \ge \max(T_x, T_y)$ 인 쪽만 채택한다. 이 조건을 인과성(causality)이라고 부른다. 값은 오직 "이미 정해진, 더 작은" 이웃에서만 들어와야 한다는 Eikonal 단방향성을 강제한다.

한 방향에서만 정보가 들어오는 경우(다른 이웃이 Far라면) 식은 단순해진다.

T = T_x + \frac{h}{F}

min-heap 삽입·삭제가 각각 $O(\log N)$ 이고, 각 점은 평균 상수 번 처리되므로 총 비용은 $O(N \log N)$ 이다.

Python 구현#

import heapq
import numpy as np
 
INF = np.inf
FAR, CONSIDERED, ACCEPTED = 0, 1, 2
 
 
def solve_eikonal_neighbor(T, j, i, h, F):
    """한 점에서 2차 upwind 해를 얻는다."""
    ny, nx = T.shape
    Tx = min(
        T[j, i - 1] if i > 0 else INF,
        T[j, i + 1] if i < nx - 1 else INF,
    )
    Ty = min(
        T[j - 1, i] if j > 0 else INF,
        T[j + 1, i] if j < ny - 1 else INF,
    )
    a, b = sorted([Tx, Ty])
    inv_f = h / F
    if not np.isfinite(b) or b >= a + inv_f:
        return a + inv_f
    diff = b - a
    disc = 2 * inv_f * inv_f - diff * diff
    if disc < 0:
        return a + inv_f
    return 0.5 * (a + b + np.sqrt(disc))
 
 
def march_wall_distance(walls, h=1.0, F=1.0):
    """2D 정렬 격자에서 Fast Marching으로 벽거리를 계산한다."""
    ny, nx = walls.shape
    T = np.full_like(walls, INF, dtype=np.float64)
    status = np.zeros_like(walls, dtype=np.uint8)
    heap = []
 
    T[walls] = 0.0
    status[walls] = ACCEPTED
 
    # 벽의 네 방향 이웃을 Considered로 초기화
    for j in range(ny):
        for i in range(nx):
            if not walls[j, i]:
                continue
            for dj, di in [(-1, 0), (1, 0), (0, -1), (0, 1)]:
                nj, ni = j + dj, i + di
                if 0 <= nj < ny and 0 <= ni < nx and status[nj, ni] != ACCEPTED:
                    new_t = solve_eikonal_neighbor(T, nj, ni, h, F)
                    if new_t < T[nj, ni]:
                        T[nj, ni] = new_t
                        status[nj, ni] = CONSIDERED
                        heapq.heappush(heap, (new_t, nj, ni))
 
    while heap:
        t, j, i = heapq.heappop(heap)
        if status[j, i] == ACCEPTED:
            continue  # 오래된 엔트리
        status[j, i] = ACCEPTED
        for dj, di in [(-1, 0), (1, 0), (0, -1), (0, 1)]:
            nj, ni = j + dj, i + di
            if not (0 <= nj < ny and 0 <= ni < nx):
                continue
            if status[nj, ni] == ACCEPTED:
                continue
            new_t = solve_eikonal_neighbor(T, nj, ni, h, F)
            if new_t < T[nj, ni]:
                T[nj, ni] = new_t
                status[nj, ni] = CONSIDERED
                heapq.heappush(heap, (new_t, nj, ni))
    return T
 
 
# 검증: 평면 채널
walls = np.zeros((61, 121), dtype=bool)
walls[0, :] = True
walls[-1, :] = True
 
d = march_wall_distance(walls, h=1.0)
 
# 중앙선(y=30)에서 벽거리는 30이어야 한다
print(f"중앙 벽거리 = {d[30, 60]:.4f} (정답 30.0)")
print(f"최대 상대 오차 = {np.max(np.abs(d[30, :] - 30)):.4f}")

평면 채널에서는 해석해 $d=\min(y, H-y)$ 가 깔끔하게 나온다. 위 코드의 중앙선 출력은 정확히 30.0에 거의 일치한다. 2D 첫 근사는 Upwind 1차 도함수라 코너 근처에서 약간의 오차가 있지만 실무에는 충분하다.

한쪽만 확정되면 칸마다 1이 더해진다#

채널 검증의 격자는 61×121, 간격은 $h=1$ , 속도는 $F=1$ 이다. 아래 벽은 인덱스 0, 위 벽은 60, 중앙선은 30이다. 한쪽 이웃만 Accepted이면 업데이트는 $T = T_x + h/F$ 로 줄어든다. $h/F=1$ 이라 벽에서 한 칸 떨어질 때마다 도달 시간이 1 늘어난다. 서른 칸이면 30이다. 위 벽에서 내려와도 30이다. 두 파면의 최솟값이 중앙의 30.0이다. 코드 주석이 정답으로 적은 숫자가 이 덧셈에서 나온다. 코너의 잔차는 판별식이 음수가 되어 1차 가지로 떨어지는 자리에서 생긴다.

10M 셀과 0.5M 벽면의 곱은 5조다. Fast Marching은 점마다 heap을 $O(\log N)$ 에 건드리고, 점은 Accepted로 한 번만 올린다. 셀이 $10^{7}$ 이면 $\log N$ 은 밑 10에서 7, 밑 2에서도 24 근처다. 연산은 수억을 크게 넘기기 어렵다. 5조와는 자릿수가 네 단계 이상 벌어진다. KD-tree도 $N_{\text{wall}}$ 의 로그를 셀 수에 다시 곱한다. 벽 쿼리가 행진보다 많이 남는다.

두 이웃의 차 $b-a$ 가 커서 $2(h/F)^2-(b-a)^2$ 가 음수면, 해는 작은 이웃에 $h/F$ 를 더한 1차 식으로 돌아간다. 채널처럼 한쪽이 아직 Far이면 이 가지가 처음부터 선택된다. 평면 벽의 거리장은 그래서 칸 수 그 자체다.

비정렬 격자로의 확장#

정렬 격자 위의 공식은 아름답지만, CFD의 절반은 비정렬 격자 위에 있다. Sethian & Vladimirsky (2000)는 이 문제를 풀었다. 한 점 $\mathbf{x}_0$ 의 값을 그 점을 포함하는 simplex(삼각형, 사면체)의 다른 꼭짓점 값으로부터 얻는다.

방향 벡터 $\mathbf{p}_r = \mathbf{x}_r - \mathbf{x}_0$ 를 행으로 모은 행렬 $P$ 와 방향도함수 $v_r = (T_0 - T_r)/|\mathbf{p}_r|$ 에서 정규방정식을 풀면, Eikonal $|\nabla T|^2 = 1/F^2$ 는 $T_0$ 에 대한 2차 방정식이 된다. 근 두 개 중 근사 그래디언트가 simplex 내부를 향하는 것만 채택한다. 이 시물렉스 인과성(simplex causality)이 정렬 격자의 단순 upwind 조건을 일반화한 것이다.

실무 체크리스트#

움직이는 메쉬: 전체 재계산 대신 변형이 큰 영역만 갱신하는 레이지 업데이트가 쓸모 있다. Fast Marching은 순차적이라 병렬화가 까다롭지만, narrow band 변형체는 구현하기 좋다.
Periodic BC: periodic 측에서 벽으로 값이 전파되도록 그래프를 이어 붙여야 한다. 안 그러면 주기 경계 근처 벽거리가 과대평가된다.
다중 재질: 로터-스테이터 해석에서는 상대 운동하는 두 벽 집합 모두 $\Gamma_w$ 에 포함시킨다. 터빈 팁 갭처럼 거리가 급변하는 곳은 격자가 충분히 촘촘해야 한다.
P-Poisson 대안: $|\nabla d|^p = 1$ ( $p=4 \sim 10$ )을 implicit하게 풀면 Fast Marching보다 느리지만 병렬 구현이 쉽다. 대형 병렬 솔버에서는 이쪽이 선호되기도 한다.
국소 음수 오차: 수치 discretization에 따라 매우 얇은 경계층에서 $d_w$ 가 음수로 나올 수 있다. 항상 $\max(d_w, 0)$ 로 클램프.

매 스텝 전 영역을 다시 행진시킬지, P-Poisson으로 넘길지#

움직이는 로터 격자에서 Fast Marching을 매 스텝 전 셀에 돌리면, 10M 셀의 heap 연산이 시간 전진마다 반복된다. 변형이 큰 띠만 갱신하지 않으면 $O(N \log N)$ 의 이득이 스텝 수만큼 다시 곱해진다.

병렬이 먼저면 $|\nabla d|^p = 1$ 쪽을 고른다. $p$ 는 4에서 10 사이다. 반복의 속도는 Jacobi에서 BiCGSTAB로 넘어가는 선형 솔버가 정한다. 행진의 인과성을 버릴 수 있을 때만 이 교체를 한다.

참고#

J. A. Sethian, A. Vladimirsky, "Fast methods for the Eikonal and related Hamilton–Jacobi equations on unstructured meshes", PNAS 97(11), 5699–5703, 2000.
P. Tucker, "Assessment of geometric multilevel convergence robustness and a wall distance method for flows with multiple internal boundaries", Applied Mathematical Modelling 22, 1998.