[논문리뷰] tanh로 충격파를 부셀 단위로 — TENO-THINC 재구성을 직접 구현

불연속 셀에서만 tanh로 갈아끼우지 않으면, 차수를 올려도 접촉면은 수십 셀로 번진다.

5차, 6차, 8차로 올려도 접촉 불연속은 시간이 갈수록 뭉개진다. WENO든 TENO든 같다. 다항식으로 점프를 그리는 한 그렇다. Takagi 등의 2023년 논문은 그 길을 비킨다. 매끄러운 곳은 고차 TENO, 불연속 셀만 tanh(THINC)다.

재구성 토글을 Polynomial, THINC, Hybrid로 바꾸면 같은 셀 평균이 S자, 수직 점프, 혼합으로 갈린다.

β = 1.80overshoot: 0.000

Polynomial은 점프 위아래로 살짝 넘친다. THINC는 한 셀 안에서 거의 수직이다. Hybrid는 점프 셀만 tanh를 쓴다.

논문 한 장 요약#

제목: High-Order Low-Dissipation Shock-Resolving TENO-THINC Schemes for Hyperbolic Conservation Laws
저자: Shinichi Takagi, Hiro Wakimura, Lin Fu, Feng Xiao
출처: Communications in Computational Physics, 2023
한 줄: 짝수점 TENO(6점·8점)에 THINC 재구성을 결합해, 매끄러운 영역의 저소산성과 불연속의 부셀 해상도를 동시에 얻는다.

공짜로 아는 것은, 언제 갈아끼울지다. TENO는 이미 후보 스텐실마다 불연속 여부를 $\delta_k \in \{0,1\}$ 로 표시한다. 이 값을 그대로 재활용해 셀이 불연속을 품었는지 가려내고, 그 셀에서만 THINC를 켠다. 새 파라미터는 THINC의 기울기 $\beta$ 하나뿐이다.

왜 고차 다항식도 접촉면을 뭉개는가#

다항식 재구성은 본질적으로 매끄러운 곡선을 그린다. 셀 평균이라는 점들을 지나는 곡선을 그릴 때, 수직에 가까운 점프 앞에서는 두 가지 중 하나를 택해야 한다. 진동(overshoot)을 허용하거나, 기울기를 죽여 뭉개거나.

TENO와 WENO는 진동을 막는 쪽을 택한다. 불연속을 가로지르는 스텐실을 버려 ENO 성질을 강제한다. 안정성은 얻지만, 매 스텝 약간씩 수치 소산이 쌓인다. 접촉 불연속(속도·압력은 같고 밀도만 점프하는 면)은 압력 메커니즘으로 스스로 날카로워지지 않는다. 그래서 장시간 이류 뒤에는 폭이 수십 셀로 번진다. 차수를 올려도 이 번짐은 거의 줄지 않는다.

앞에서 고른 세 재구성이 그 선택이다. 같은 셀 평균인데 곡선만 다르다. TENO가 작은 와를 덜 죽이는 쪽은 저소산 TENO에 있다. 접촉면이 수십 셀로 번지는 일은 그 저소산으로도 멈추지 않는다.

THINC — tanh가 부셀 점프를 그린다#

THINC(Tangent of Hyperbola for INterface Capturing)는 다항식 대신 단조 증가하는 tanh를 쓴다. 셀 $i$ 내부 좌표 $X \in [-0.5, 0.5]$ 에 대한 재구성 함수는

h_i(X) = f_a + f_d \tanh\!\big(\beta(X - d_i)\big)

여기서 $f_a = (f_{i+1}+f_{i-1})/2$ 는 이웃 평균, $f_d = (f_{i+1}-f_{i-1})/2$ 는 점프 절반, $\beta$ 는 tanh의 기울기, $d_i$ 는 점프 중심 위치다.

$d_i$ 는 자유 변수가 아니다. 셀 평균 보존 조건(재구성의 적분이 $f_i$ 와 같아야 함)이 $d_i$ 를 못 박는다. 그 해는 해석적으로

d_i = \frac{1}{2\beta} \ln \frac{1 - T_2/T_1}{1 + T_2/T_1}

$T_1 = \tanh(\beta/2)$ , $T_2 = \tanh(\alpha\beta/2)$ , $\alpha = (f_i - f_a)/f_d$ 이다. $\alpha$ 는 셀 평균이 이웃 사이에서 어디쯤 있는지를 나타내는 정규화 위치다. 셀이 단조롭지 않으면( $\,(f_{i+1}-f_i)(f_i-f_{i-1}) \le 0\,$ ) tanh가 무의미하므로 1차 재구성( $\hat f_{i+1/2} = f_i$ )으로 후퇴한다.

β: 날카로움을 쥔 손잡이#

$\beta$ 하나가 THINC의 성격을 결정한다. 작으면 완만하고 크면 가파르다. 논문은 의미 있는 대응을 정리한다.

$\beta \approx 1.1$ : van Leer 한정자 수준의 TVD 스펙트럼
$\beta \approx 1.3$ : Superbee 한정자 수준(과압축 경향)
$\beta = 1.6 \sim 2.0$ : 날카로운 충격 포착 영역

논문은 $\beta = 1.8$ 을 채택한다. 너무 키우면 매끄러운 곡선까지 계단처럼 squaring 되어 인공적인 평탄화가 생긴다. 위 그림의 $\beta$ 슬라이더를 1.0에서 2.4로 밀면 점프가 서고, overshoot 표시가 같이 움직인다.

$\beta = 1.8$ 이면 $\beta/2 = 0.9$ 다. 셀 한가운데 점프면 $\alpha = 0$ 이라 $d_i = 0$ 이다. 오른쪽 경계 $X = 0.5$ 에서 $\tanh(\beta(X - d_i)) = \tanh(0.9) \approx 0.716$ 이다. 이웃 점프의 약 72%만 한 셀 안에서 회복한다. $\beta = 1.1$ 이면 $\tanh(0.55) \approx 0.50$ 이라 van Leer처럼 완만하다. $\beta = 2.4$ 면 $\tanh(1.2) \approx 0.83$ 까지 선다. 1.8은 그 사이, 수직과 squaring의 중간이다.

$\alpha = (f_i - f_a)/f_d$ 다. 셀 평균이 이웃 한가운데면 $\alpha = 0$ 이고 점프 중심도 0이다. 평균이 오른쪽 이웃으로 붙으면 $\alpha$ 가 1에 가까워지고, 로그 안의 $(1 - T_2/T_1)$ 이 줄며 $d_i$ 가 양의 방향으로 밀린다. $|\alpha| \ge 1$ 이거나 $(f_{i+1} - f_i)(f_i - f_{i-1}) \le 0$ 이면 코드는 tanh를 버리고 $\hat f = f_i$ 로 물러난다. 봉우리 셀은 그래서 다항식이 맡는다.

초기 top-hat은 $0.35$ 와 $0.65$ 사이다. 폭 $0.30$ , $N = 100$ 이면 30셀이다. $dx = 0.01$ , $\Delta t = 0.45 \times dx = 0.0045$ 다. 다섯 바퀴는 $5 / 0.0045 \approx 1111$ 걸음이고 CFL은 0.45다. 그 뒤 minmod 폭은 23셀, 피크는 0.84다. 0.5를 넘는 칸이 7개 줄고 높이는 0.16이 빠진다. THINC-BVD는 피크 1.00, 폭 30을 유지한다. 코드 주석의 두 줄이 이 셈과 같다.

TENO가 찾고 THINC가 메운다 — BVD 하이브리드#

THINC를 모든 셀에 쓰면 매끄러운 파동까지 망가진다. 그래서 불연속 셀에서만 켜야 한다. 논문은 TENO 가중치가 이미 계산한 $\delta_k$ 로 셀별 지표 $\zeta_i$ 를 만든다.

\zeta_i = \begin{cases} 2, & \delta_{2,i}=0 \text{ 그리고 } \delta_{1,i}=0 \\ 1, & \delta_{1,i}=0 \text{ 그리고 } \delta_{2,i}=1 \\ -1, & \delta_{2,i}=0 \text{ 그리고 } \delta_{1,i}=1 \\ 0, & \text{그 외} \end{cases}

$\zeta_i = 1$ 은 오른쪽에, $-1$ 은 왼쪽에 불연속이 있다는 신호다. 좌우 이웃의 $\zeta$ 부호가 셀 $i$ 를 가운데 두고 마주 보면, 그 셀이 진짜 불연속을 품었다고 판정하고 THINC로 교체한다.

실전 BVD(Boundary Variation Diminishing)는 더 단순하게도 쓴다. 두 재구성을 모두 계산한 뒤, 셀 경계에서 좌·우 재구성값의 차이(경계 변동)를 더 줄이는 쪽을 채택한다. 아래 구현은 이 BVD 판정을 따른다.

Python — top-hat을 한 바퀴 돌려본다#

선형 이류 $u_t + u_x = 0$ 로 top-hat을 주기 경계에서 돌린다. minmod TVD와 THINC-BVD를 같은 격자·같은 시간만큼 비교한다.

import numpy as np
 
def minmod(p, q):
    return np.where(p * q <= 0, 0.0, np.where(np.abs(p) < np.abs(q), p, q))
 
def thinc_right_edge(u, beta=1.8):
    # 식 2.19~2.20: 셀 i의 우측 경계(i+1/2) 좌측 재구성값
    um, up = np.roll(u, 1), np.roll(u, -1)
    fa, fd = 0.5 * (up + um), 0.5 * (up - um)
    mono = ((up - u) * (u - um) > 0) & (np.abs(fd) > 1e-12)
    alpha = np.where(mono, (u - fa) / np.where(mono, fd, 1.0), 0.0)
    mono = mono & (np.abs(alpha) < 1.0)
    T1, T2 = np.tanh(beta / 2), np.tanh(alpha * beta / 2)
    d = np.where(mono, (1 / (2 * beta)) * np.log((1 - T2 / T1) / (1 + T2 / T1)), 0.0)
    vR = fa + fd * np.tanh(beta * (0.5 - d))
    return np.where(mono, vR, u)  # 비단조면 1차로 후퇴
 
def reconstruct(u, beta, use_thinc):
    um, up = np.roll(u, 1), np.roll(u, -1)
    s = minmod(u - um, up - u)
    poly_vR = u + 0.5 * s
    if not use_thinc:
        return poly_vR
    poly_vL = u - 0.5 * s
    t = thinc_right_edge(u, beta)
    # BVD: 우측 경계 점프를 더 줄이는 재구성을 채택
    jump_poly = np.abs(poly_vR - np.roll(poly_vL, -1))
    jump_thinc = np.abs(t - np.roll(poly_vL, -1))
    return np.where(jump_thinc < jump_poly, t, poly_vR)
 
def rhs(u, dx, beta, use_thinc):
    vR = reconstruct(u, beta, use_thinc)        # a=1 이므로 F_{i+1/2} = vR
    return -(vR - np.roll(vR, 1)) / dx
 
def advect_tophat(N=100, revolutions=5.0, beta=1.8, use_thinc=True):
    dx, x = 1.0 / N, (np.arange(N) + 0.5) / N
    dt = 0.45 * dx
    u = ((x > 0.35) & (x < 0.65)).astype(float)   # 초기 top-hat
    for _ in range(int(revolutions / dt)):
        u1 = u + dt * rhs(u, dx, beta, use_thinc)            # SSP-RK2
        u = 0.5 * u + 0.5 * (u1 + dt * rhs(u1, dx, beta, use_thinc))
    return x, u
 
x, u_tvd = advect_tophat(use_thinc=False)
_, u_thinc = advect_tophat(use_thinc=True)
print("TVD   peak=%.3f  width(>0.5)=%d" % (u_tvd.max(), (u_tvd > 0.5).sum()))
print("THINC peak=%.3f  width(>0.5)=%d" % (u_thinc.max(), (u_thinc > 0.5).sum()))
# 예시 출력:
# TVD   peak=0.84  width(>0.5)=23
# THINC peak=1.00  width(>0.5)=30

다섯 바퀴를 돈 뒤 minmod TVD는 피크가 0.84로 주저앉고 폭이 좁아진다. THINC-BVD는 피크 1.00을 유지하고 폭도 초기값(30셀)을 거의 지킨다. 차수를 더 올린 TENO라도 다항식 본성상 이 번짐을 완전히 막지는 못한다는 것이 논문의 출발점이다.

직접 만지는 재구성과 장시간 이류#

revolutions를 밀면 빨간 minmod와 청록 THINC가 같은 시간만큼 따로 돈다.

β = 1.80

revolutions가 늘수록 빨간 minmod 곡선은 양 끝이 무너지며 사다리꼴로 퍼진다. 청록 THINC-BVD는 회색 정확해에 거의 붙어 직사각형을 지킨다. $\beta$ 를 1.1까지 내리면 THINC도 TVD처럼 완만해지고, 2.0 부근에서 가장 또렷해진다.

재현하며 의심스러웠던 것#

논문의 약속은 깔끔하지만, 직접 짜며 걸린 부분이 셋 있었다.

BVD 판정은 검출기의 품질에 민감하다. 매끄러운 극값(smooth extremum)을 불연속으로 오인하면 그 봉우리가 부자연스럽게 평탄해진다. 위 단순 BVD는 단조 조건으로 이를 막지만, 다차원·계 시스템에서는 특성 변수 분해가 필요하다.

$\beta$ 는 문제에 따라 움직인다. $\beta = 1.8$ 은 1D 벤치마크에 맞춰진 값이다. 강한 압축성 난류에서는 과압축으로 인공 계단이 생길 수 있어, 적응적 $\beta$ 를 검토할 가치가 있다.

비용 주장은 맥락에 따라 다르다. "THINC는 불연속 셀이 적어 거의 공짜"라는 논문 서술은 재구성만 볼 때 맞다. 그러나 분기(branch)가 많은 BVD 판정은 GPU에서 워프 발산을 부를 수 있다. OpenFOAM 같은 비정렬 FV 코드에 옮기려면 셀별 검출 로직을 face flux 루프에 녹여야 해서, 1D 정렬 격자만큼 깔끔하지 않다.

β를 1.8에 묶어 둘 것인가#

다섯 바퀴 뒤 THINC-BVD는 피크 1.00, 폭 30셀이다. 같은 조건의 minmod는 0.84, 23셀이다. 매끄러운 극값까지 계단으로 눌리면 그 고정값은 과하다. 손잡이 하나를 파동마다 의심하는 쪽은 THINC와 MP5를 가르는 선택이다.

점프를 tanh로 그릴지, 번진 뒤를 걷어낼지는 다른 결정이다. 재구성이 접촉을 이미 지키면 뒤처리 샤프닝은 필요 없다. 지키지 못하면 $\beta$ 를 더 키우기 전에 검출기부터 본다.