Godunov·충격파와 entropy 선택

총량을 보존한다고 올바른 파를 고르는 것은 아니다

지역 상태에 따라 속도가 달라지는 교통 흐름을 생각해보자. 혼잡한 구간의 앞쪽에는 날카로운 전선이 생길 수 있고, 반대로 퍼지는 구간에서는 부채꼴 모양의 파가 열릴 수 있다. 수치 모형이 도로 전체의 총량을 보존하더라도 전선 속도를 틀리게 예측하거나, 퍼져야 하는 점프를 그대로 유지할 수 있다.

이 단원은 무점성 Burgers 수송으로 그 차이를 구체적으로 확인한다. Godunov와 미보정 Roe라는 두 공유 플럭스 방법으로 실제 셀 평균을 갱신한다. 둘 다 보존형이지만 특성 속도 0을 가로지르는 희박파에서는 면 플럭스가 다르다. Roe는 잘못된 정지 팽창 점프를 그대로 유지할 수 있다.

연속체 충격파의 적합성, 수치 면에서의 entropy 결손, 유한 시간 단계 뒤의 entropy 예산을 따로 살펴보자. 어느 한 수준의 검사에 통과했다고 나머지 둘까지 해결되지는 않는다. 특히 전역 entropy 상쇄는 잘못된 면 선택을 감출 수 있다.

무차원 비선형 보존법칙을 정하기

이 예제의 상태 q, 좌표 x, 시간 t는 무차원이다. 뒤에서 사용하는 이차량을 열역학적 엔트로피나 물리적 에너지라고 부르지 않는다. Burgers는 비선형 수송을 연구하는 스칼라 모형이며, 완전한 압축성 유체 시뮬레이션이 아니다.

$$q_t+\partial_x f(q)=0,\qquad f(q)=\frac{q^2}{2},\qquad f'(q)=q.$$

매끄러운 해에서는 수송 속도가 상태 자체와 같다. 큰 양의 값은 더 빠르게 오른쪽으로 이동하고 더 작은 음의 값은 더 빠르게 왼쪽으로 이동한다. 특성선이 모이면 충격파가, 벌어지면 희박파가 생긴다. 불연속이 생긴 뒤에는 q와 공간 미분을 점별로 곱한 식이 적분 보존법칙을 대신할 수 없다.

초기 조건은 일정한 왼쪽 상태에서 일정한 오른쪽 상태로 바뀌는 점프 하나다. 점프 위치는 계산 구간의 1/4, 1/2, 3/4 중 하나인 셀 경계에 맞춘다. 어느 셀에 임의의 혼합값을 표본으로 넣을지 고민하지 않고, 초기 셀 평균을 분명하게 정할 수 있다.

연속체에서 적합한 파를 선택하기

압축 점프는 왼쪽 상태가 오른쪽 상태보다 크다. 움직이는 전선 양쪽의 보존 수지를 맞추면 Rankine–Hugoniot 속도를 얻는다. 볼록한 수학적 entropy eta와 그 물리적 플럭스 psi도 정의하자.

$$s=\frac{q_L+q_R}{2},\qquad \eta(q)=\frac{q^2}{2},\qquad\psi(q)=\frac{q^3}{3},\qquad [\psi]-s[\eta]=-\frac{(q_L-q_R)^3}{12}.$$

대괄호는 오른쪽 값에서 왼쪽 값을 뺀다는 뜻이다. 압축 충격파에서는 마지막 식이 비양수여서 이 entropy 부등식과 맞는다. 점프를 뒤집으면 결손이 양수가 된다. 이 팽창 충격파는 보존법칙을 만족할 수는 있어도 적합한 해가 아니다.

왼쪽 상태가 더 작으면 적합한 해는 연속적인 희박파다. 양의 시간에서 부채꼴 내부의 상태는 초기 점프 위치를 기준으로 한 유사 변수와 같다. 양 끝은 각각 왼쪽 상태와 오른쪽 상태의 속도로 이동한다. 양쪽이 같으면 상수 해이고, 서로 반대인 같은 크기의 압축 상태는 정지 충격파를 만든다.

참조해도 각 셀에서 적분하기

계산은 유한 구간에서 하지만 정확한 참조는 전체 실수선의 Riemann 해다. 희박파는 유사 변수 xi로 정의하고, 충격파는 움직이는 불연속으로 정의한다.

$$\xi=\frac{x-x_0}{t},\qquad q_{\mathrm{rare}}(x,t)= \begin{cases}q_L,&\xi<q_L,\\ \xi,&q_L\le\xi\le q_R,\\ q_R,&\xi>q_R,\end{cases} \qquad q_{\mathrm{shock}}(x,t)=\begin{cases}q_L,&x<x_0+st,\\q_R,&x>x_0+st.\end{cases}$$

불연속 전선의 한 점에 어떤 값을 지정하는지는 셀 적분을 바꾸지 않는다. 수치 셀과 비교할 대상은 그 셀 전체에서 평균한 값이다. 충격파가 셀을 가로지르면 양쪽 상태가 차지한 길이의 비율로 평균한다. 희박파에서는 셀 중 부채꼴 안에 들어간 구간의 선형 프로필을 적분한다.

$$\bar q_i(t)=\frac1{\Delta x}\int_{x_i}^{x_{i+1}}q(x,t)\,dx,\qquad \mathcal E_1=\sum_i\Delta x\,|q_i-\bar q_i(t)|.$$

연속체 entropy 밀도 자체도 별도로 적분한다. q의 제곱을 2로 나눈 값의 평균과 평균 q를 제곱해 2로 나눈 값은 다르다. 충격파가 일부를 채운 셀이나 희박파 셀에서는 일반적으로 양의 평균화 간극이 생긴다. 이를 구분하지 않으면 연속체 entropy 예산을 다른 이산 이차량과 비교하는 셈이 된다.

구간별 상수인 수치 재구성과 연속 정확 프로필 사이의 절대 오차도 별도로 적분한다. 거친 셀의 평균은 거의 맞아도 셀 안의 날카로운 충격파를 전혀 해상하지 못할 수 있다. 따라서 셀 평균 오차가 격자 세분화에 따라 반드시 단조 감소하는 것은 아니다. 같은 물리 시간의 재구성 오차와 실제 충격 위치를 함께 확인해야 한다.

Godunov 면 플럭스를 선택하기

Godunov는 움직이지 않는 면에서 정확한 국소 Riemann 해가 주는 값을 사용한다. 희박파 전체가 오른쪽으로 가면 왼쪽 상태, 왼쪽으로 가면 오른쪽 상태, 속도 0을 가로지르면 0이 면에서 보인다. 충격파에서는 충격 속도의 부호에 따라 추적값을 고른다.

국소 상태	면의 선택	수치 플럭스
증가하며 둘 다 비음수	왼쪽 상태	왼쪽 상태 제곱 / 2
증가하며 둘 다 비양수	오른쪽 상태	오른쪽 상태 제곱 / 2
증가하며 0을 가로지름	sonic 값 0	0
감소하며 충격 속도가 양수	왼쪽 추적값	왼쪽 상태 제곱 / 2
감소하며 충격 속도가 음수	오른쪽 추적값	오른쪽 상태 제곱 / 2

속도가 0인 압축 충격파에서는 두 물리 플럭스가 같다. 구현에서는 플럭스 선택의 추적값을 왼쪽으로 선언한다. 이것은 불연속 전선의 점값이 하나로 유일하게 정해진다는 뜻이 아니다.

이 볼록 스칼라 플럭스에서는 증가 구간의 물리 플럭스 최솟값, 감소 구간의 양 끝 물리 플럭스 최댓값을 사용하는 독립적인 구성도 가능하다. 추적값을 고르는 제품 구현을 다른 계산 방식으로 확인하는 데 유용하다.

미보정 Roe의 반례를 드러내기

Roe는 국소 비선형 문제를 두 상태의 산술평균 속도로 수송하는 문제로 바꾼다. 중심 플럭스에 점성을 더한 이 공식도 보존형이지만, 특성 속도 0을 가로지르는 희박파에는 추가 처리가 필요하다.

$$\widehat q=\frac{a+b}{2},\qquad F_{\mathrm{Roe}}(a,b)=\frac{f(a)+f(b)}2-\frac{|a+b|}{4}(b-a),\qquad F_{\mathrm{Roe}}(-1,1)=\frac12,\quad F_{\mathrm{Godunov}}(-1,1)=0.$$

반례의 평균 속도는 0이다. Roe는 가상의 정지 불연속을 유지하지만 정확한 entropy 해는 −t에서 +t까지 퍼지는 부채꼴이다. 두 플럭스 모두 이웃 셀이 완벽하게 공유할 수 있으므로 보존성만으로 어느 파를 선택했는지 판별할 수 없다.

여기서는 미보정 Roe를 의도적으로 비교한다. 모든 Roe 기반 해법이 같은 실패를 한다는 뜻은 아니다. 실제 해법은 entropy fix를 적용할 수 있다. 이 실험의 목적은 빠진 조건 하나를 분리하고 진단이 그 실패를 잡는지 확인하는 것이다.

실제 상태를 갱신하고 경계를 선언하기

각 이전 시간의 셀 상태로 모든 면 플럭스를 구한 뒤 셀을 동시에 갱신한다. 두 해법은 각자의 상태 궤적을 실행한다. Godunov라는 이름을 Roe 실행 결과의 참조 표시에만 붙이는 것이 아니다.

$$r_i^n=\frac{F_{i-1/2}^n-F_{i+1/2}^n}{\Delta x},\qquad q_i^{n+1}=q_i^n+\Delta t\,r_i^n.$$

양 끝의 ghost 셀에는 해당 이전 시간의 정확한 전체선 Riemann 셀 평균을 넣는다. 이것은 외부 자료를 어떻게 공급하는지에 관한 선언이다. 수치적으로 완벽한 무반사 경계를 구현했다는 주장이 아니다. 내부 절단 오차는 여전히 면 플럭스에 영향을 줄 수 있다.

파가 계산 구간 밖으로 나간 뒤에도 이 경계 정의를 유지한다. 전체선 정확해는 외부 셀 평균과 내부 참조를 계속 제공한다. 희박파나 충격파가 경계를 지났는데도 초기 원방 상태를 그대로 고정하지 않는다. 시간 간격도 입력 그대로 실행한다. 불안정한 실험을 몰래 더 작은 시간 간격으로 바꾸지 않는다.

수치 경계 장부와 연속체 경계 장부를 따로 닫기

실제 갱신에 셀 폭을 곱해 합하면 내부 면은 이웃 셀과 상쇄되고 양 끝의 수치 경계 플럭스만 남는다.

$$\mathcal M_h^n=\sum_i\Delta x\,q_i^n,\qquad \mathcal M_h^{n+1}-\mathcal M_h^n =\Delta t(F_{\mathrm{left}}^n-F_{\mathrm{right}}^n).$$

수치 장부는 실제 계산한 그 플럭스를 누적한다. 정확한 연속체 경계 공급량을 대신 넣어 일치를 강제하지 않는다. 계산한 셀 변화도 따로 합산해 상쇄나 누적 오차가 관측될 수 있게 한다.

연속체 참조에는 별도의 시간 적분 경계 플럭스가 있다. 파가 양 끝에 도달하기 전에는 경계 추적값이 원방 상태다. 통과한 뒤에는 경계가 희박파 내부나 충격파 뒤에 놓일 수 있다. 독립 참조는 이런 통과 시간을 기준으로 시간 적분을 나눈다. 초기 시간에만 맞는 수지를 쓰지 않고 전체 구간에서 실제 적분 질량과 경계 공급량을 대조한다.

entropy 플럭스와 잠재함수를 구분하기

이차 entropy의 entropy 변수는 q다. 잠재함수는 q에 보존 플럭스를 곱한 뒤 psi를 뺀 값이며, 여기서는 물리적 entropy 플럭스의 절반이다. 둘을 혼동하면 수치 예산이 바뀐다.

$$\Psi(q)=qf(q)-\psi(q)=\frac{q^3}{6},\qquad \widehat\psi(a,b;F)=\frac{a+b}{2}F-\frac{a^3+b^3}{12},\qquad D(a,b;F)=\frac{b^3-a^3}{6}-(b-a)F.$$

F는 선택한 해법의 보존 플럭스다. 수치 entropy 플럭스는 두 상태로 정의한 양이며, 선택한 추적값의 물리적 entropy 플럭스와 같을 필요가 없다. 1에서 −1로 가는 정지 충격파의 수치 entropy 플럭스는 0이지만 양쪽 물리 플럭스는 1/3과 −1/3이다.

이 entropy에서 Godunov의 D는 모든 면에서 비음수다. 미보정 Roe의 −1에서 1로 가는 면에서는 D가 −2/3이고, 반대 방향의 압축 면에서는 +2/3이다. 두 셀 주기 패턴에는 두 면이 모두 있어서 전역 기여가 상쇄된다. 심지어 인접 셀의 예산에서도 상쇄된다. 그러므로 D를 면별로 살펴봐야 한다.

유한 시간 단계의 entropy 보정을 남겨두기

셀의 양쪽 면 기여를 모으면 정확한 반이산 항등식이 나온다. 실제 Euler 갱신을 전개하면 시간 이산화에 따른 양의 추가 항이 생긴다.

$$\Delta x\,\Delta\eta_i+\Delta t(\widehat\psi_{i+1/2}-\widehat\psi_{i-1/2}) =-\frac{\Delta t}{2}(D_{i-1/2}+D_{i+1/2}) +\frac{\Delta t^2\Delta x}{2}r_i^2.$$

우변이 비양수가 되는 충분한 국소 조건은 다음과 같다. 이 조건은 선언한 수치 entropy 플럭스와 이차 entropy에 대한 것이며, 모든 entropy나 보편적인 최적 CFL 상수를 말하는 것은 아니다.

$$\Delta t\,\Delta x\,r_i^2\le D_{i-1/2}+D_{i+1/2}.$$

현재 셀과 ghost 상태의 최대 절댓값에 기초한 실제 특성 Courant 지표도 따로 제공한다. 1보다 크면 통상적인 한 셀 전파 조건을 넘는다. 1보다 작다고 Roe의 transonic 면 선택이 고쳐지는 것은 아니다. 경계값에 정확히 걸치면 부동소수점 반올림이 엄격한 부등식 판단에 영향을 줄 수 있으므로 Boolean 표시와 실제 수치를 함께 확인하자.

유한 구간의 이차 entropy 총량은 반드시 감소하지 않는다. 경계의 entropy 공급이 내부 소산을 이길 수 있기 때문이다. 주기 합은 경계 공급을 없애지만, 양·음의 면 결손이 서로 상쇄되는 문제까지 없애지는 않는다.

주기 셀 한 단계를 손으로 계산하기

별도의 주기 손계산으로 폭이 1인 셀 네 개에 −1, −1, 1, 1을 넣고 시간 간격을 0.1로 두자. 증가 면은 0을 가로지르고, 닫는 면은 정지 압축 충격파다. 이는 한 단계 주기 예제이며 유한 구간의 전체선 참조 실험과는 구분한다.

Godunov의 면 플럭스는 1/2, 0, 1/2, 1/2다. 변화율은 0, 1/2, −1/2, 0이며 새 상태는 −1, −0.95, 0.95, 1이다. 총질량은 0으로 유지된다. 이차 entropy는 2에서 1.9025로 줄어든다. 반이산 항이 −0.1이고 양의 Euler 보정이 0.0025여서 총변화는 −0.0975다.

Roe는 모든 면에 1/2을 주므로 어느 셀도 바뀌지 않는다. 전역 entropy 변화는 0이지만 증가 면은 적합하지 않다. 모형이 사용하는 실제 갱신 함수로 이 계산을 재현할 수 있다.

import { burgersPeriodicStep } from '@/lib/learning-continuum-godunov';

const q = [-1, -1, 1, 1];
const good = burgersPeriodicStep(q, 1, 0.1, 'godunov');
const bad = burgersPeriodicStep(q, 1, 0.1, 'roe');
console.log(good.values.map(v => v.value)); // [-1, -0.95, 0.95, 1]
console.log(bad.values.map(v => v.value));  // [-1, -1, 1, 1]

비보존 갱신과 비교하고 작은 오차를 보존하기

양의 주기 상태에서는 q와 q의 미분을 곱한 형태의 한쪽 차분이 유용한 음성 대조다. q 제곱을 2로 나눈 플럭스의 미분을 공유 면으로 이산화한 알고리즘과 같지 않다.

$$q_i^{\mathrm{adv}}=q_i-\frac{\Delta t}{\Delta x}q_i(q_i-q_{i-1}),\qquad \Delta\mathcal M_{\mathrm{adv}}=-\frac{\Delta t}{2}\sum_i(q_i-q_{i-1})^2.$$

따라서 이 advective 갱신은 결과가 매끄러워 보여도 일반적으로 질량을 잃는다. 상태가 1, 1/2, 0, 1/4이고 시간 간격이 0.1이면 질량 변화는 −0.05625다. 같은 주기 스텐실의 공유 Godunov 갱신은 경계 공급과 대수적 질량 변화가 모두 0이다.

작은 오차에는 별도 처리가 필요하다. 정확한 파가 바로 이웃한 셀들을 벗어나기 전에는 정확한 셀 평균 변화가 시간에 대해 정확히 선형이다. 수치 계산의 초기 선형 응답과 누적 비선형 나머지를 분리해 저장하면, 전체 셀 값이나 주된 변화보다 수백 자릿수 작은 실제 수치 오차도 남길 수 있다.

부호·로그 표현은 일반 수치 범위 밖의 비영 곱을 보존한다. 큰 경계 항과 소산 항이 상쇄된 뒤 작은 이차 변화가 남는 경우에도 공유 면 일의 결합식을 합산 전에 대수적으로 평가한다. 개별 성분을 반올림해 더한 값은 진단으로 따로 남긴다. 작은 값을 잘라 0으로 만들지 않는다. 다만 산술은 여전히 유한 정밀도이며 셀 직접 합산이나 경계값 부근의 부호 판단에는 나쁜 조건이 생길 수 있다.

AI 모형 개발에 반례를 활용하기

학습한 보존 플럭스 모형은 압축 충격파, transonic 희박파, 상수 상태, 음의 수송, 경계를 통과하는 파에서 평가해야 한다. 파의 위치, 셀 평균 오차, 수치 경계 수지와 면별 entropy 결손을 확인하자. 작은 전역 잔차만으로는 국소적으로 잘못 선택한 파를 놓칠 수 있다.

압축 충격파에서는 정확한 전선 위치와 수치 셀 평균의 문턱 교차 위치를 비교한다. 퍼진 전선에는 전이 폭이 있고 불안정한 프로필에는 여러 교차가 생길 수 있다. 모든 실행에 정확한 전선 하나가 있는 것처럼 표시하지 말고 그 모호함을 보고해야 한다. 전선이 구간을 벗어난 뒤에는 내부 위치 추정 자체가 제공되지 않을 수 있다.

매끄러운 자료만 학습하면 보존형과 advective 형태의 차이가 드러나지 않을 수 있다. 면을 공유하는 구조는 내부 상쇄에 도움을 주지만 올바른 entropy 해, 경계 조건과 안정한 시간 간격을 자동으로 고르지는 않는다. 추가적인 수학 구조와 실패 사례에 맞춘 검증이 필요하다.

풀이가 있는 연습과 일차 자료

1. 왼쪽 2, 오른쪽 0의 파를 구하라. 압축 충격파이며 속도는 1이다. 시간 0.2에는 길이 0.2만큼 이동한다. entropy 결손은 −8/12, 즉 −2/3이다. 속도는 플럭스 점프를 상태 점프로 나누어 얻는다.
2. −1에서 1은 왜 다른가? 0을 가로지르는 증가 쌍이다. 정확한 희박파의 정지 면 상태는 0이어서 Godunov 플럭스도 0이다. Roe는 평균 속도 0을 사용해 1/2을 반환하고, 적합하지 않은 팽창 점프를 유지한다.
3. 네 셀의 entropy 변화를 확인하라. 두 셀만 크기 1에서 0.95로 바뀐다. 합친 이차 변화는 0.95의 제곱에서 1을 뺀 −0.0975다. 질량 변화는 상쇄된다. Roe의 상태가 바뀌지 않는 결과는 면별 entropy 검사에 통과하지 못한다.
4. 주기 entropy 변화가 0이면 모든 면이 적합한가? 아니다. 두 셀 Roe 패턴의 증가 면은 −2/3, 감소 면은 +2/3의 결손을 갖는다. 총합과 인접 합이 모두 상쇄된다. 합이나 셀 표시만 보지 말고 개별 면 결손을 확인해야 한다.
5. advective 질량 결손을 계산하라. 상태 (1, 1/2, 0, 1/4)의 이웃 점프 제곱 합은 9/8이다. −0.1/2를 곱하면 −9/160, 즉 −0.05625다. 공유 보존 플럭스는 같은 주기 스텐실에서 이런 내부 손실을 만들 수 없다.
6. 파가 유한 구간을 벗어나면 무엇이 바뀌는가? 전체선 정확한 참조는 유효하지만 경계 추적값이 초기 원방 상태와 달라질 수 있다. 이전 시간마다 정확한 외부 셀 평균을 공급하고, 통과 시간을 나누어 연속체 경계 플럭스를 독립 적분한다. 수치 질량 장부는 실제 수치 경계 플럭스와 대조한다.

일차 자료인 Clawpack Burgers 장은 충격파와 희박파 해를 설명한다. Burgers 근사 해법 장은 Roe의 transonic 실패와 entropy fix의 필요성을 다룬다. entropy 잠재함수, 수치 entropy 플럭스와 반이산·완전 이산 안정성의 구분은 Tadmor, Acta Numerica 2003의 틀을 직접 대수로 특수화했다. 문헌은 이론의 근거이며 실제 구현은 별도의 독립 고정밀 계산과 경계 적분 예산으로 검증한다.