Motion in an Externally Generated Velocity Field
Osher & Fedkiw, Ch. 3
Reference: Stanley Osher and Ronald Fedkiw, Level Set Methods and Dynamic Implicit Surfaces, Springer, 2003, Ch. 3
앞의 두 챕터는 각각 “인터페이스를 어떻게 표현할 것인가”(implicit function, signed distance function)와 “ODE를 어떻게 시간 전진시킬 것인가”(Euler, Runge-Kutta)를 다뤘습니다. 이 챕터는 이 둘을 실제로 결합해서 인터페이스를 움직이는 첫 번째 구체적인 방법을 제시합니다. 레벨셋 방법의 이름 자체가 여기서 등장하는 방정식 하나(“level set equation”)에서 나왔을 만큼, 이 챕터는 레벨셋 방법 전체의 출발점입니다. 동시에 이 챕터는 왜 PDE를 단순한 중심차분으로 풀면 안 되는지, 왜 upwind 방식이 필요한지, 그리고 그 upwind 방식을 어떻게 5차 정확도까지 끌어올릴 수 있는지(HJ WENO)를 보여주는, 수치 스킴 설계의 축소판이기도 합니다.
Introduction: Lagrangian과 Eulerian, 두 가지 인터페이스 이동 방식
인터페이스 위의 각 점 $\vec{x}$에서 속도 $\vec{V}(\vec{x})$가 주어졌다고 합시다. 이 속도로 인터페이스를 움직이는 가장 직접적인 방법은, 인터페이스 위의 모든 점에 대해 상미분방정식
\[\frac{d\vec{x}}{dt} = \vec{V}(\vec{x})\]를 푸는 것입니다. 이것이 Lagrangian formulation입니다. 문제는 인터페이스 위에 점이 무한히 많다는 것이고, 그래서 실제로는 인터페이스를 선분(2차원)이나 삼각형(3차원)으로 이산화해서 그 꼭짓점들을 이동시키는 방식(front tracking)을 씁니다. 이 접근은 그럴듯해 보이지만 실무적으로 매우 까다롭습니다. 아주 단순한 속도장조차 경계 요소를 심하게 왜곡시킬 수 있고, 그러면 주기적으로 메쉬를 다듬어줘야 합니다. 더 심각한 것은 토폴로지가 바뀌는 경우(인터페이스가 합쳐지거나 갈라지는 경우)로, 이때는 요소들을 분리하고 다시 연결하는 “수술”이 필요한데, 이 수술 자체가 상당히 복잡한 별도의 알고리즘 문제가 됩니다.
레벨셋 방법은 이 문제를 완전히 다른 방식으로 우회합니다. 인터페이스를 점들의 집합으로 직접 추적하는 대신, 1장에서 다룬 implicit function $\phi$ 자체를 시간에 따라 진화시켜서 그 $\phi=0$ 등고선이 원하는 대로 움직이도록 만드는 것입니다. 이를 위한 방정식이 바로 convection(advection) 방정식입니다.
\[\phi_t + \vec{V}\cdot\nabla\phi = 0\]여기서 $\phi_t$는 시간에 대한 편미분이고, $\vec{V}\cdot\nabla\phi = u\phi_x+v\phi_y+w\phi_z$입니다. 이 식은 Eulerian formulation입니다. 점을 직접 추적(track)하는 대신, 고정된 격자 위에서 $\phi$ 값 자체가 인터페이스를 포착(capture)하도록 만드는 방식이기 때문입니다. 이 방정식은 Osher와 Sethian이 처음 도입했다고 해서 흔히 level set equation이라고 불리며, 연소(combustion) 분야에서는 화염 전면을 나타내는 $G(\vec{x})=0$ 등고선에 대해 똑같은 형태로 쓰이는 G-equation이라는 이름으로도 알려져 있습니다. 서로 다른 분야에서 독립적으로 같은 방정식에 도달했다는 사실은, 이 방정식이 “얇은 경계면의 움직임”이라는 문제에 대한 상당히 자연스러운 정식화라는 근거가 됩니다.
속도장을 어떻게 격자 전체에 정의할 것인가
여기서 실용적인 문제가 하나 생깁니다. $\vec{V}$가 인터페이스 위에서만 정의되어 있다면, Cartesian grid 위에서 이 PDE를 풀기가 애매해집니다. 그래서 보통 $\vec{V}$를 인터페이스 근방의 band, 혹은 도메인 전체로 확장해서 정의합니다. 그런데 이 확장 방식이 아무렇게나 되어서는 안 됩니다. 책이 드는 예시가 인상적인데, 만약 $\vec{V}$가 인터페이스 위에서만 $(1,0,0)$이고 나머지 영역에서는 정확히 0이라면, 참값은 “인터페이스가 오른쪽으로 속도 1로 이동”하는 것이지만, 격자점 대부분은 인터페이스 위에 있지 않으므로 $\vec{V}\cdot\nabla\phi$가 대부분의 격자점에서 그냥 0이 되어버립니다. 결국 수치적으로는 인터페이스가 거의 움직이지 않는 것으로 계산됩니다. $\vec{V}$를 연속함수로 제한해도 문제는 근본적으로 해결되지 않습니다. 인터페이스 주변 두께 $\epsilon$짜리 band 안에서만 속도가 유의미하게 변한다면, $\Delta x$가 $\epsilon$보다 충분히 작지 않은 한 여전히 대부분의 격자점이 이 band 바깥에 있어서 같은 문제가 재현됩니다.
이 문제를 해결하는 책의 결론은 흥미롭습니다. 인터페이스의 이동을 실제로 결정하는 것은 인터페이스에 접하는(tangential) 방향의 속도 변화뿐이고, 인터페이스에 수직인(normal) 방향으로 속도가 어떻게 변하든 그것은 인터페이스 이동과 무관합니다. 따라서 속도장의 변화를 최소화하면서도 정확한 인터페이스 운동을 보존하는 가장 자연스러운 선택은, 임의의 점 $\vec{x}$에서의 속도를 그 점에서 가장 가까운 인터페이스 점 $\vec x_C$에서의 속도 $\vec{V}(\vec x_C)$로 정의하는 것입니다. 이렇게 하면 인터페이스 위에서의 속도값은 그대로 보존하면서도, 인터페이스 법선 방향으로는 속도가 국소적으로 거의 일정해집니다. 책은 이 선택이 단순히 편의를 위한 트릭이 아니라는 점도 짚습니다. Zhao 등의 연구에 따르면, 이 “가장 가까운 인터페이스 점의 속도”를 사용해서 이류(advect)시키면 signed distance function이 signed distance function으로 잘 유지되는 경향이 있다는 것입니다. 이 속도 확장을 실제로 어떻게 계산하는지는 8장에서 다룬다고 예고합니다.
공간 이산화: Upwind Differencing
$\phi$와 $\vec{V}$가 격자 위에 정의되었다면, 이제 (3.2)를 실제로 풀어야 합니다. 가장 단순한 시간 이산화는 forward Euler 방법입니다.
\[\frac{\phi^{n+1}-\phi^n}{\Delta t} + \vec{V}^n\cdot\nabla\phi^n = 0\]문제는 공간 미분 $\nabla\phi^n$을 어떻게 근사하느냐입니다. 1장에서 소개한 forward/backward/central difference 중 아무거나 순진하게 쓰면 이 방법은 실패합니다. 왜 그런지 1차원 버전으로 살펴보면 이해가 쉽습니다.
\[\frac{\phi_i^{n+1}-\phi_i^n}{\Delta t} + u_i^n\phi_{x,i}^n = 0\]여기서 특성곡선법(method of characteristics)의 직관이 중요합니다. $u_i>0$이면 $\phi$ 값들이 왼쪽에서 오른쪽으로 흘러가고 있다는 뜻이므로, $x_i$에서 다음 시간 스텝의 값을 결정하는 정보는 $x_i$의 왼쪽에서 와야 합니다. 즉 $\phi_x$는 backward difference $D^-\phi$로 근사해야 합니다. 반대로 $u_i<0$이면 정보가 오른쪽에서 오므로 forward difference $D^+\phi$를 써야 합니다. 만약 이를 거꾸로, 즉 $u_i>0$일 때 $D^+\phi$를 쓰면 정작 $x_i$의 새로운 값을 결정해야 할 왼쪽의 정보를 아예 사용하지 않는 셈이 되어 근사가 원리적으로 틀리게 됩니다. 이렇게 속도의 부호에 따라 미분 근사 방향을 바꾸는 것을 upwind differencing(풍상차분)이라고 부릅니다. $u_i=0$이면 해당 항 자체가 사라지므로 근사할 필요가 없습니다. $D^-\phi$와 $D^+\phi$ 모두 1차 정확도이므로, 이 upwind 스킴 전체는 공간에 대해 1차 정확도를 갖습니다.
CFL 조건
forward Euler(시간)와 upwind(공간)를 결합한 이 방법은 $\Delta t\to 0$, $\Delta x\to 0$일 때 오차가 0으로 수렴한다는 의미에서 consistent합니다. 그런데 consistency만으로는 충분하지 않습니다. Lax-Richtmyer 동등성 정리(equivalence theorem)에 따르면, 선형 PDE에 대한 유한차분 근사가 실제로 참값에 수렴하려면 consistent할 뿐 아니라 stable해야 합니다. 안정성은 근사 과정에서 생기는 작은 오차가 시간이 지나도 증폭되지 않는다는 성질입니다.
이 안정성은 Courant-Friedreichs-Lewy 조건(CFL 조건)으로 강제됩니다. 직관적으로 말하면, 수치적으로 계산된 파동(정보)의 전파 속도가 실제 물리적 파동의 전파 속도보다 느려서는 안 된다는 조건입니다. 격자 위에서 한 스텝에 정보가 전달되는 속도는 $\Delta x/\Delta t$이고, 실제 물리적 속도는 $\lvert u\rvert$이므로, $\Delta x/\Delta t > \lvert u\rvert$가 필요합니다. 이를 스텝 크기에 대한 제약으로 다시 쓰면
\[\Delta t < \frac{\Delta x}{\max\{\lvert u\rvert\}}\]이고, 실무에서는 이를 CFL number $\alpha$로 다음과 같이 강제합니다.
\[\Delta t\left(\frac{\max\{\lvert u\rvert\}}{\Delta x}\right) = \alpha, \qquad 0<\alpha<1\]$\alpha=0.9$ 정도가 근사적으로 최적에 가까운 선택이고, $\alpha=0.5$는 조금 더 보수적인 선택으로 자주 쓰입니다. 다차원에서는
\[\Delta t\,\max\left\{\frac{\lvert u\rvert}{\Delta x}+\frac{\lvert v\rvert}{\Delta y}+\frac{\lvert w\rvert}{\Delta z}\right\} = \alpha\]또는
\[\Delta t\left(\frac{\max\{\lvert \vec{V}\rvert\}}{\min\{\Delta x,\Delta y,\Delta z\}}\right) = \alpha\]형태의 조건이 쓰입니다. 참고로 book은 central differencing을 그대로 쓰고 싶다면, forward Euler와 통상적인 CFL 조건($\Delta t\sim\Delta x$) 조합으로는 불안정하고, $\Delta t\sim(\Delta x)^2$ 수준의 훨씬 엄격한(=계산 비용이 큰) 조건이 필요하다고 지적합니다. 3차 정확도 TVD Runge-Kutta로 시간 이산화를 바꾸거나, 우변에 인공 점성(artificial viscosity) 항 $\mu\Delta\phi$($\mu\sim\Delta x$)를 추가하는 방법도 central differencing을 안정화시킬 수 있지만, 책은 이 세 대안보다 conservation law 수치해법에서 이미 검증된 upwind 계열 방법을 선호한다고 명시합니다.
더 높은 차수의 공간 정확도: Hamilton-Jacobi ENO
앞선 upwind 스킴은 $\phi_x^-$(backward, $u>0$일 때 사용)와 $\phi_x^+$(forward, $u<0$일 때 사용) 각각을 1차 정확도로만 근사했습니다. HJ ENO(Hamilton-Jacobi Essentially Non-Oscillatory) 방법은 어느 쪽을 쓸지 결정하는 upwind 원리는 그대로 유지하면서, $\phi_x^-$와 $\phi_x^+$ 자체의 근사 정확도를 끌어올립니다. 핵심 아이디어는 Harten 등이 conservation law를 위해 제안한 ENO 다항식 보간을 가져와서, “가장 매끄러운(smoothest)” 다항식 보간을 골라 미분하는 것입니다. Osher와 Sethian은 1차원 Hamilton-Jacobi 방정식이 conservation law의 적분이라는 사실을 이용해 이 ENO 아이디어를 Hamilton-Jacobi 방정식에 확장했습니다.
divided difference와 Newton 다항식 보간
이 방법은 Newton 다항식 보간의 표준 도구인 divided difference를 사용합니다. 0차 divided difference는 격자점 값 자체이고,
\[D_i^0\phi = \phi_i\]1차 divided difference는 격자점 사이 중점에서 정의됩니다.
\[D_{i+1/2}^1\phi = \frac{D_{i+1}^0\phi - D_i^0\phi}{\Delta x}\]이는 사실 익숙한 backward/forward difference와 정확히 같은 것입니다($D_{i-1/2}^1\phi=D_i^-\phi$, $D_{i+1/2}^1\phi=D_i^+\phi$). 2차 divided difference는 다시 격자점에서, 3차는 다시 중점에서 정의되며, 같은 패턴으로 반복됩니다. 이 divided difference들을 이용해서 $\phi(x)=Q_0(x)+Q_1(x)+Q_2(x)+Q_3(x)$ 형태의 다항식을 재구성하고, 이를 미분해서 $\phi_x(x_i)=Q_1’(x_i)+Q_2’(x_i)+Q_3’(x_i)$로 $\phi_x^-$와 $\phi_x^+$를 근사합니다($Q_0$는 상수항이라 미분하면 사라집니다). $Q_1’$ 항만 남기면($Q_2,Q_3$ 생략) 정확히 1차 upwind와 같아지므로, “1차 정확도 다항식 보간은 곧 1차 upwind differencing과 동일하다”는 것을 알 수 있습니다. 2차·3차 보정항($Q_2’, Q_3’$)을 추가하면 각각 2차, 3차 정확도가 됩니다.
어느 stencil을 선택할 것인가
여기서 흥미로운 설계 결정이 등장합니다. 2차 보정을 만들 때, 왼쪽 점을 하나 더 포함시킬지($D_k^2\phi$) 오른쪽 점을 하나 더 포함시킬지($D_{k+1}^2\phi$) 두 가지 선택지가 있습니다. ENO의 원칙은, 매끄럽고 완만하게 변하는 데이터는 divided difference 표에서 작은 값을 만들고, 불연속이거나 급격히 변하는 데이터는 큰 값을 만든다는 사실을 이용하는 것입니다. 즉 $\lvert D_k^2\phi\rvert$와 $\lvert D_{k+1}^2\phi\rvert$ 중 더 작은 쪽을 선택해서, 변화가 큰(불연속에 가까운) 영역을 되도록 보간에 포함시키지 않도록 합니다. 이렇게 고른 계수 $c$로
\[Q_2(x) = c(x-x_k)(x-x_{k+1})\]를 정의하고, 이를 미분해서 2차 보정항 $Q_2’(x_i)=c\,(2(i-k)-1)\Delta x$를 얻습니다. 여기서 정해진 $k^*$(더 작은 divided difference를 준 쪽의 인덱스)를 그대로 이어받아, 3차 보정 역시 $\lvert D_{k^*+1/2}^3\phi\rvert$와 $\lvert D_{k^*+3/2}^3\phi\rvert$를 비교해서 더 작은(=더 매끄러운) 쪽을 선택하는 방식으로 계속됩니다.
\[Q_3(x) = c^\*(x-x_{k^\*})(x-x_{k^\*+1})(x-x_{k^\*+2})\]이 선택 과정 전체를 요약하면, 매 단계마다 “덜 진동하는(가장 매끄러운) 다항식 보간을 선택한다”는 원칙을 반복 적용해서 stencil을 확장해 나가는 것이 HJ ENO의 본질입니다.
Hamilton-Jacobi WENO
HJ ENO는 3차 정확도를 위해 정확히 하나의 stencil(왼쪽/가운데/오른쪽 조합 중 하나)만을 선택합니다. Liu 등은 이 “셋 중 하나만 고르기”가 데이터가 매끄러운 영역에서는 지나치게 보수적(overkill)이라고 지적했습니다. 매끄러운 영역이라면 세 후보 stencil 모두 유용한 정보를 담고 있으니, 하나만 버리지 말고 가중 평균해서 쓰자는 것이 WENO(Weighted ENO)의 아이디어입니다. 이후 Jiang과 Shu가 이 가중치를 최적화해서 매끄러운 영역에서 5차 정확도를 달성하는 방법을 제안했고, Jiang과 Peng이 이를 Hamilton-Jacobi 방정식에 확장한 것이 HJ WENO입니다. 책은 이 방법이 3차 정확도 HJ ENO보다 오차를 한 자릿수(10배) 이상 줄여준다고 평가합니다.
$\phi_{x,i}^-$를 예로 들면, $v_1,\ldots,v_5$를 인접한 다섯 개의 backward difference 값이라 할 때 세 후보 근사는
\[\phi_x^1 = \frac{v_1}{3}-\frac{7v_2}{6}+\frac{11v_3}{6}, \qquad \phi_x^2 = -\frac{v_2}{6}+\frac{5v_3}{6}+\frac{v_4}{3}, \qquad \phi_x^3 = \frac{v_3}{3}+\frac{5v_4}{6}-\frac{v_5}{6}\]이고, HJ WENO 근사는 이 셋의 convex combination입니다.
\[\phi_x = \omega_1\phi_x^1 + \omega_2\phi_x^2 + \omega_3\phi_x^3, \qquad \omega_1+\omega_2+\omega_3=1,\ \ 0\le\omega_k\le1\]핵심은 매끄러운 영역에서 $\omega_1=0.1$, $\omega_2=0.6$, $\omega_3=0.3$을 쓰면 5차 정확도를 얻는다는 사실입니다. 다만 이는 데이터가 매끄러울 때만 유효한 최적값이고, 불연속 근처에서 이 고정 가중치를 그대로 쓰면 오히려 심각한 오차를 만듭니다. 그래서 데이터가 매끄럽지 않은 영역에서는 특정 stencil의 가중치를 0에 가깝게 낮춰서, 사실상 HJ ENO처럼 “덜 진동하는 stencil 하나만 선택”하는 것과 비슷하게 동작하도록 만들어야 합니다.
smoothness 지표와 가중치 계산
이 적응형 가중치는 각 stencil의 매끄러움을 측정하는 지표 $S_k$로부터 계산됩니다.
\[S_1 = \frac{13}{12}(v_1-2v_2+v_3)^2 + \frac{1}{4}(v_1-4v_2+3v_3)^2\]$S_2, S_3$도 인접한 $v$ 값들의 조합으로 비슷하게 정의됩니다(2차·1차 차분의 제곱합 형태). $S_k$가 클수록 그 stencil이 큰 변화(불연속에 가까운 데이터)를 포함하고 있다는 뜻입니다. 이 $S_k$로부터
\[\alpha_1 = \frac{0.1}{(S_1+\epsilon)^2}, \qquad \alpha_2 = \frac{0.6}{(S_2+\epsilon)^2}, \qquad \alpha_3 = \frac{0.3}{(S_3+\epsilon)^2}\]을 계산하고, 이를 정규화해서 $\omega_k=\alpha_k/(\alpha_1+\alpha_2+\alpha_3)$로 최종 가중치를 얻습니다. 여기서 $\epsilon$은 $S_k\approx 0$(완전히 매끄러운 경우)일 때 0으로 나누는 것을 막기 위한 작은 상수로, $\epsilon = 10^{-6}\max{v_1^2,\ldots,v_5^2}+10^{-99}$로 정의됩니다.
이 구조가 왜 잘 작동하는지 단계별로 보면 직관이 뚜렷합니다. 데이터가 충분히 매끄러워서 모든 $S_k$가 $\epsilon$보다 훨씬 작으면, $\alpha_k$는 거의 $0.1/\epsilon^2$, $0.6/\epsilon^2$, $0.3/\epsilon^2$의 비율을 그대로 유지하므로 정규화하면 정확히 $\omega_1\approx0.1$, $\omega_2\approx0.6$, $\omega_3\approx0.3$이 되어 최적의 5차 정확도가 나옵니다. 반대로 특정 stencil의 $S_k$가 다른 것들보다 훨씬 크면(그 stencil이 불연속을 가로지르고 있다는 신호), 그 $\alpha_k$는 상대적으로 작아지고 가중치 $\omega_k$도 작아져서, 그 stencil의 영향력이 자연스럽게 줄어듭니다. 세 $S_k$ 모두 크다면(데이터가 국소적으로 매우 험한 경우) 어떤 stencil도 특별히 신뢰할 수 없다는 뜻이고, 이는 HJ ENO도 마찬가지로 겪는 어려움이지만, 책은 이런 상황이 보통 공간·시간적으로 국소적이어서 상황이 지나가면 방법이 스스로 회복된다고 설명합니다.
한 가지 이론적으로 우아한 지점은, 가중치를 정확히 $0.1, 0.6, 0.3$이 아니라 $0.1+O((\Delta x)^2)$ 수준으로 살짝 흔들어도 5차 정확도가 그대로 유지된다는 사실입니다. 이는 세 가중치의 편차 계수 $C_1,C_2,C_3$가 $C_1+C_2+C_3=0$(가중치 합이 항상 1이어야 하므로)이라는 제약 때문에 편차로 인한 오차 항이 정확히 상쇄되기 때문입니다. 이 성질 덕분에 위에서 정의한 $S_k$ 기반의(다소 지저분한 형태의) 근사 가중치를 써도 매끄러운 영역에서의 정확도 차수가 손상되지 않는다는 것이 보장됩니다.
시간 이산화: TVD Runge-Kutta
HJ ENO/WENO는 공간 방향으로 최대 5차 정확도를 제공하지만, 지금까지 사용한 forward Euler 시간 이산화는 여전히 1차 정확도입니다. 책은 실무 경험상 레벨셋 방법이 공간 정확도에는 민감하지만 시간 절단 오차에는 상대적으로 덜 민감해서, 낮은 차수의 forward Euler로도 시간 방향은 충분히 버틸 수 있는 경우가 많다고 말합니다. 다만 더 높은 시간 정확도가 필요한 경우를 위해, Shu와 Osher는 TVD(total variation diminishing) Runge-Kutta 방법을 제안했습니다. 이는 method of lines 접근, 즉 공간 이산화로 PDE를 우선 거대한 ODE 시스템(2장에서 다룬 것과 같은 형태)으로 바꾼 뒤, 그 ODE를 시간 방향으로 얼마나 정확하게 전진시키느냐의 문제로 분리하는 방식에 기반합니다.
TVD RK의 설계 원리는 단순합니다. forward Euler 스텝 자체가 TVD(허위 진동을 만들지 않음)라고 가정하면, 이 Euler 스텝들을 연속으로 적용한 뒤 그 결과들을 계수가 모두 양수인 convex combination으로 평균 내는 것만으로도 전체 방법이 TVD로 유지됩니다. (다만 책은 HJ ENO/WENO를 upwind와 결합한 실제 스킴은 엄밀히는 TVD가 아니라 TVB(total variation bounded)에 가깝다고 솔직하게 짚습니다.)
2차 정확도 TVD RK(Heun의 predictor-corrector 방법, 수정 Euler 방법과 동일)는 2장에서 다룬 2차 Runge-Kutta 방법을 PDE의 공간 연산자에 그대로 적용한 것입니다. 먼저 Euler 스텝으로 $t^n+\Delta t$까지 전진하고,
\[\phi^{n+1} = \phi^n - \Delta t\,\vec{V}^n\cdot\nabla\phi^n\]이 결과에서 다시 한 번 Euler 스텝으로 $t^n+2\Delta t$까지 전진한 뒤,
\[\phi^{n+2} = \phi^{n+1} - \Delta t\,\vec{V}^{n+1}\cdot\nabla\phi^{n+1}\]마지막으로 초기값과 이 결과를 평균 냅니다.
\[\phi^{n+1} = \frac{1}{2}\phi^n + \frac{1}{2}\phi^{n+2}\]3차 정확도 TVD RK는 이 절차를 한 단계 더 반복합니다. 위와 같이 두 번의 Euler 스텝을 거친 뒤, $\frac34\phi^n+\frac14\phi^{n+2}$로 평균 내서 $t^n+\frac12\Delta t$ 시점의 근사값 $\phi^{n+1/2}$를 얻고, 여기서 다시 한 번 Euler 스텝을 밟은 뒤 $\frac13\phi^n+\frac23\phi^{n+3/2}$로 최종 평균을 내서 $t^n+\Delta t$ 시점의 3차 정확도 근사값을 얻습니다. 즉 3차 방법은 “Euler 스텝 세 번 + convex combination 두 번”이라는 조합으로, 순수하게 Euler 스텝이라는 하나의 building block을 반복 재사용해서 시간 정확도를 끌어올리는 구조입니다.
책은 4차 이상의 TVD RK 스킴도 존재하지만, 실무에서는 이 이상의 시간 정확도가 큰 이득을 주지 못한다고 지적합니다. HJ WENO 자체가 흥미로운 유동 영역에서는 이미 3차 정확도 HJ ENO와 비슷한 수준으로 정확도를 잃는 경우가 많기 때문에, 시간 정확도만 4차 이상으로 올려봐야 병목이 다른 곳(공간 근사)에 있다는 것입니다. 게다가 4차 이상의 TVD RK는 upwind와 downwind 차분을 모두 요구해서 공간 연산자 평가 비용이 두 배가 됩니다. 그래서 실무적으로는 3차 정확도 TVD RK가 시간 이산화의 실질적인 상한선처럼 쓰인다고 이해하면 됩니다.
Comment
이 챕터를 관통하는 하나의 원칙은 “정보가 오는 방향을 존중하라”는 것입니다. upwind differencing부터 HJ ENO의 stencil 선택, HJ WENO의 적응형 가중치까지, 모든 설계가 결국 “어느 방향의 데이터가 더 신뢰할 만한 정보를 담고 있는가”를 판단하고 그쪽에 더 무게를 싣는 방식으로 이루어져 있습니다. 이는 우연이 아니라, 이 챕터가 다루는 레벨셋 방정식이 본질적으로 쌍곡형(hyperbolic) 방정식이고, 쌍곡형 방정식은 정보가 특성곡선을 따라 유한한 속도로 전파된다는 성질을 갖기 때문입니다. 이 성질을 무시하고 중심차분처럼 양방향 정보를 동등하게 섞으면, 물리적으로 오지 않은 방향의 정보까지 억지로 섞어 넣는 셈이 되어 불안정해지는 것이 자연스러운 결과입니다.
또 하나 짚어둘 만한 점은, 이 챕터가 앞선 두 챕터의 도구들을 거의 그대로 재사용한다는 것입니다. HJ ENO는 1장의 divided difference·다항식 보간을, TVD RK는 2장의 Euler·Runge-Kutta 스텝을 building block으로 그대로 가져다 씁니다. 즉 이 챕터의 새로움은 완전히 새로운 계산 도구를 발명하는 데 있다기보다, 기존 도구들을 “쌍곡형 PDE에 안전하게 적용하려면 어떤 순서와 조합으로 써야 하는가”를 설계하는 데 있습니다. 이런 조합적 사고방식은 이후 챕터에서 다룰 mean curvature motion(포물형 방정식)이나 reinitialization 같은 다른 레벨셋 연산에서도 반복해서 등장할 것으로 보입니다.
다만 이 챕터가 다루지 않는 부분도 명확히 해둘 필요가 있습니다. 속도장을 인터페이스에서 도메인 전체로 확장하는 구체적인 알고리즘은 8장으로, $\phi$를 signed distance function으로 되돌리는 재초기화(reinitialization)는 이후 다른 챕터로 미뤄져 있습니다. 즉 이 챕터는 “속도장이 이미 잘 정의되어 있고 $\phi$가 signed distance function에 가깝다”는 전제 위에서 순수하게 PDE의 시공간 이산화만 다룬 것이고, 실제 레벨셋 시뮬레이션을 완성하려면 그 전제를 채워주는 나머지 챕터들이 필요합니다. 지금 단계에서 기억해야 할 핵심은 다섯 가지입니다: (1) 인터페이스는 점을 추적하지 않고 $\phi_t+\vec{V}\cdot\nabla\phi=0$을 풀어서 이류시킨다, (2) 공간 미분은 반드시 속도의 부호에 따라 upwind 방향으로 근사해야 한다, (3) 안정성을 위해 $\Delta t$는 CFL 조건 $\Delta t\lesssim\Delta x/\max\lvert\vec{V}\rvert$을 만족해야 한다, (4) HJ WENO는 세 ENO 후보의 매끄러움 기반 가중 평균으로 매끄러운 영역에서 5차 정확도를 낸다, (5) TVD RK는 Euler 스텝과 convex combination만으로 시간 정확도를 안전하게 끌어올린다.