[논문 리뷰] WENO가 작은 와류를 뭉갤 때 — Fu(2019) TENO 저소산 유한체적 기법

같은 10%에서 χ는 1/4 근처에 남고, 다섯 점 스텐실은 꺼지지 않는다.

Stanford에서 Lin Fu가 난류 직접수치모사를 돌리다가 이상한 걸 봤다. 5차 WENO-JS는 충격을 잘 잡았다. 충격에서 먼 작은 와류까지 같이 흐려졌다. 차수는 5라고 적혀 있었다. 손에 남는 해상도는 3차 중심차분 쪽이었다.

$C_T$ 슬라이더를 $10^{-5}$ 에서 $10^{-3}$ 이나 $10^{-7}$ 로 옮기면 컷오프가 까다로워지거나 풀리고, 같은 범프가 주황과 청록으로 달리 깎인다.

N cells: 120

CFL: 0.40

TENO cutoff C_T: 1e-5

WENO-JSTENO-5L2 err — WENO: 0.0000 | TENO: 0.0000

한 바퀴 뒤에 사각 펄스 모서리는 둘 다 한두 칸인지 본다. $10^{-7}$ 에서 청록 펄스 옆에 잔진동이 생기는지도 본다.

논문 정보#

저자: Lin Fu (Stanford Turbulence Research Center)
학술지: Computer Physics Communications 235 (2019) 25–39
DOI: 10.1016/j.cpc.2018.10.009
키워드: TENO, High-order schemes, Shock-capturing, Low-dissipation, Finite-volume method

WENO-JS의 숨은 손실#

고전 WENO-JS(Jiang–Shu, 1996)는 작은 스텐실 $S_0, S_1, S_2$ 마다 평활도 $\beta_k$ 를 재고, 비선형 가중치로 섞는다.

w_k = \frac{\alpha_k}{\sum_j \alpha_j}, \qquad \alpha_k = \frac{d_k}{(\beta_k + \varepsilon)^2}

$d_k$ 는 5차를 되돌리는 최적 선형 가중치다. $\varepsilon$ 은 $10^{-6}$ 근처의, 0 나눗셈을 막는 상수다. 이 가중치는 연속이다. 스텐실이 아주 조금 흔들려도, $w_k=d_k$ 로 있어야 할 매끈한 구간에서 $w_k$ 가 움직이고, 더 낮은 차수 쪽에 무게가 실린다.

매끄러운 파동이 수십 스텝 가면 진폭이 깎이고 위상 오차가 쌓인다. 충격과 난류 와류가 같이 있는 LES에서는 이 소산이 치명적이다.

TENO의 두 아이디어: 강한 스케일 분리와 샤프 컷오프#

Fu는 두 군데를 바꿨다.

스텐실 구성을 바꿨다. 세 개의 3점 스텐실 $S_0, S_1, S_2$ 에 큰 5점 스텐실 $S_3$ 를 따로 둔다. $S_3$ 가 깨끗하면 그걸로 끝이다. 나머지는 쓰지 않는다.

가중치를 이진화했다. 매끈한지 판단은 컷오프 $C_T$ 가 한다.

\gamma_k = \left(1 + \frac{\tau_K}{\beta_k + \varepsilon}\right)^q, \qquad \chi_k = \frac{\gamma_k}{\sum_j \gamma_j}

\delta_k = \begin{cases} 0, & \chi_k < C_T \\ 1, & \text{otherwise} \end{cases}

$\tau_K=|\beta_0-\beta_2|$ 는 전체 기준이다. $q=6$ 으로 비율을 키운다. $\delta_k$ 는 0 아니면 1이다. 스텐실은 완전히 버리거나, 최적 가중치로 살아난다. 연속적인 전이가 없다. 매끈한 구간에서는 $S_3$ 의 5차 선형을 그대로 쓰고, 불연속 근처에서만 ENO처럼 스텐실을 고른다. 한 컷오프로 모든 파동을 맡길 수 있는지는 별개의 질문이다. 파동마다 재구성을 갈라야 한다는 논지는 그 시도를 의심한다.

수식: 부드러운 가중치에서 0/1 판정으로#

$S_3$ 가 매끈하면, $\delta_3=1$ 이면, 재구성은 5차 선형이다.

u_{i+1/2}^L = \frac{1}{60}\bigl(2\bar u_{i-2} - 13\bar u_{i-1} + 47\bar u_i + 27\bar u_{i+1} - 3\bar u_{i+2}\bigr)

$\bar u_j$ 는 셀 평균이다. 이 식은 완전히 선형이다. $C_T$ 가 0에 가까우면 중심차분처럼 소산이 작다. $\delta_3=0$ 이면 살아남은 작은 스텐실만 $d_k$ 로 섞는다.

u_{i+1/2}^L = \sum_{k=0}^{2} w_k\, u_{k,i+1/2}^L, \qquad w_k = \frac{d_k\, \delta_k}{\sum_j d_j\, \delta_j}

이게 ENO 성질이다. 불연속을 가로지르는 스텐실은 $\delta_k$ 가 0으로 꺼지고, 다른 스텐실로 넘어간다.

NumPy로 5차 스텐실 재구성#

실제 코드는 짧다. 다섯 점에서 평활도 셋과 후보 다항식 셋을 만든 뒤, TENO는 $C_T$ 로 끊는다. WENO-JS는 역제곱으로 섞는 쪽이 비교 대상이다. 이 함수 자체는 TENO 가지다.

import numpy as np
 
def teno_beta_indicators(um2, um1, u0, up1, up2):
    # 3점 스텐실 3개의 Jiang-Shu 평활도 지표
    b0 = (13/12)*(um2 - 2*um1 + u0)**2 + 0.25*(um2 - 4*um1 + 3*u0)**2
    b1 = (13/12)*(um1 - 2*u0 + up1)**2 + 0.25*(um1 - up1)**2
    b2 = (13/12)*(u0 - 2*up1 + up2)**2 + 0.25*(3*u0 - 4*up1 + up2)**2
    return b0, b1, b2
 
def teno_weighted_reconstruction(u, i, CT=1e-5, q=6):
    N = len(u)
    um2, um1, u0 = u[(i-2) % N], u[(i-1) % N], u[i]
    up1, up2 = u[(i+1) % N], u[(i+2) % N]
    b0, b1, b2 = teno_beta_indicators(um2, um1, u0, up1, up2)
    b3 = max(b0, b1, b2)
    tau = abs(b0 - b2)
    eps = 1e-40
    g = np.array([(1 + tau/(b + eps))**q for b in (b0, b1, b2, b3)])
    chi = g / g.sum()
    delta = (chi >= CT).astype(float)
    if delta[3] == 1:  # S3가 스무드 → 5차 선형 그대로
        return (2*um2 - 13*um1 + 47*u0 + 27*up1 - 3*up2) / 60
    d = np.array([1/10, 6/10, 3/10])
    alpha = d * delta[:3]
    if alpha.sum() < 1e-14:
        return u0  # 최후의 도너셀 폴백
    w = alpha / alpha.sum()
    p0 = (2*um2 - 7*um1 + 11*u0) / 6
    p1 = (-um1 + 5*u0 + 2*up1) / 6
    p2 = (2*u0 + 5*up1 - up2) / 6
    return w[0]*p0 + w[1]*p1 + w[2]*p2

일러 스텝을 붙이고 주기 경계로 돌리면, 부드러운 가우스 벌지는 TENO에서 거의 원형으로 돌아오고 WENO-JS에서는 30%쯤 낮아진다. 뼈대는 같고, 가중치 계산만 다르다.

가중치 식에 작은 흔들림을 넣는다. 이상 선형 가중치는 10분의 1, 10분의 6, 10분의 3이다. β가 셋 다 1이면 WENO의 w는 그 값으로 돌아온다. 세 번째만 10퍼센트 올려 β를 1, 1, 1.1로 둔다. 세 번째 알파는 0.3을 1.21로 나눈 약 0.248이다. 합으로 나누면 w는 약 0.105, 0.633, 0.262다. 이상 가중치와의 L1 차이는 약 0.077이다. 불연속이 아니다. 평활도 지표가 10퍼센트 움직인 것뿐이다.

같은 숫자를 TENO 쪽에 넣는다. 타우는 0.1이고 q는 6이다. 코드는 큰 스텐실 지표를 세 베타의 최댓값으로 두고 네 개를 함께 나눈다. 그러면 카이는 대략 0.256, 0.256, 0.244, 0.244다. 컷오프가 $10^{-5}$ 이면 네 자릿수 위다. 판정은 전부 1로 남고, 큰 스텐실이 살아 5차 선형식을 그대로 쓴다. 같은 10퍼센트에 WENO는 가중치가 흔들리고 TENO는 한 비트도 안 뒤집힌다.

본문 식의 엡실론은 $10^{-6}$ 근처다. 함수 안에서는 $10^{-40}$ 이다. 베타가 1 근처면 둘 다 결과에 안 나온다. 베타가 0에 붙는 스텐실에서만 분모가 갈린다. 컷오프와 별개로, 엡실론이 적힌 자리마다 다르다. 5차 식의 계수 2, −13, 47, 27, −3을 더하면 60이다. 60으로 나누면 상수장은 그대로다. 큰 스텐실이 살아 있을 때의 바닥이다.

C_T를 옮기면 남는 진폭#

앞의 이류는 같은 초기조건, 부드러운 사인 제곱 범프와 사각 펄스를 선형으로 민다. 주황이 WENO-JS, 청록이 TENO-5다. $C_T=10^{-5}$ 에서 한 바퀴 돌면 사각 펄스 모서리는 둘 다 한두 칸만 번진다. 왼쪽 사인 범프는 WENO-JS만 눈에 띄게 낮아진다. $C_T$ 를 $10^{-3}$ 으로 올리면 컷이 까다로워져 사각 펄스 근처에서 작은 스텐실로 자주 떨어지고, TENO도 소산이 는다. $10^{-7}$ 로 내리면 TENO는 대부분 5차 선형을 고집하고, 사각 펄스 옆에 잔진동이 보인다. $C_T$ 가 안정과 저소산의 저울이다.

비판적 고찰 — 튜닝의 함정#

논문이 자신 있게 말하지 않는 자리를 적어 둔다.

$C_T$ 는 문제에 붙는다. Fu는 $10^{-5}$ 를 권한다. 강한 충격과 아주 약한 음파가 같이 있으면 이 값에서 음파가 살짝 뜨고, $10^{-6}$ 으로 내리면 충격 근처에서 진동이 돈다. OpenFOAM이나 SU2에 옮긴 그룹들은 케이스마다 이 값을 다시 맞춘다고 적었다.

유한체적에서는 $S_3$ 의 평활도 $\beta_3$ 가 비싸다. 5점 다항식의 미분 적분이 여섯 쪽에 달한다. 유한차분으로 쓰면 훨씬 짧다. 이 논문의 저소산 Riemann 플럭스 전환까지 같이 쓰려면 유한체적 구조가 자연스럽다.

시간 적분과의 궁합은 비어 있다. SSP-RK3에서는 공간 정확도에 눌려 괜찮다. SSP-RK2나 음해법과 섞으면 $\delta_k$ 의 불연속이 Newton 수렴을 건드린다는 사후 보고가 있다.

실무에서 쓴다면 격자 수렴을 같은 $C_T$ 로 돌린다. 공간 차수가 실제로 5로 돌아오는지 본다. 돌아오지 않으면 $C_T$ 가 너무 커서 TENO가 상시 비선형 모드인 가능성이 높다.

C_T를 한 값에 고정하면 깨지는 검증#

$C_T=10^{-5}$ 를 고정한 채 격자를 줄이면, 매끈한 셀의 χ가 컷 아래로 떨어지는 순간 5차가 사라진다. 위에서 넣은 10퍼센트 흔들림은 컷을 넘지 못했다. 충격과 아주 약한 파가 같이 있으면 이야기가 달라진다.

그때 차수 검증이 5가 아니라 1 아래로 나오면, 먹은 것은 공간 연산자가 아니라 컷오프다. 같은 코드가 0.79차와 2.00차를 동시에 낸 검증은 케이스가 무엇을 재는지 묻는다. $C_T$ 를 고정한 TENO는 그 둘 중 어느 숫자를 내놓게 되는가.