비균일 체적·공유 플럭스와 보존

크기가 다른 물통에서 시작하기

서로 연결된 물통 세 개를 생각해보자. 가운데 물통의 부피는 양옆 물통의 두 배다. 첫 물통에서 가운데 물통으로 물질 1그램이 이동하면, 내보낸 물통과 받은 물통에는 같은 질량이 반대 부호로 기록되어야 한다. 하지만 농도 변화는 같지 않다. 큰 물통에서는 같은 질량이 더 많은 부피에 희석되기 때문이다.

보존형 이산화의 핵심이 여기에 있다. 연결관의 전달량을 공유하는 것만으로는 부족하다. 질량을 농도 변화로 바꿀 때 사용하는 저장 체적도 맞아야 한다. 그래프 신경망에서 연결선의 메시지를 공유하더라도 노드별 체적을 무시하고 농도를 갱신하면 실제 질량을 만들거나 없앨 수 있다.

이 단원에서는 크기가 다른 셀이 원형으로 연결된 관을 사용한다. 셀 평균 농도를 실제로 시간 갱신하고, 의도적으로 잘못된 저장 규칙과 비교한다. 질량과 체적 가중 이차 노름을 따로 측정한다. 두 진단은 서로 다른 질문에 답한다. 질량을 보존하는 계산도 음의 농도를 만들거나 불안정해질 수 있다.

물리 모형과 단위를 먼저 정하기

관의 단면적 A는 일정하고 양수이며, 수송 속도 a는 양수나 음수인 상수다. 확산계수 D는 비음수다. 농도는 수동적이어서 그 값이 주어진 유속을 바꾸지 않는다. 오른쪽 끝은 첫 셀로 이어지는 주기 경계이며, 바깥 경계에서 들어오는 별도 공급원은 없다.

$$\partial_t c+a\partial_x c=D\partial_{xx}c,\qquad c(x+L,t)=c(x,t),\qquad A>0,\quad D\ge0.$$

길이와 위치의 단위는 m, 단면적은 m², 유속은 m/s, 확산계수는 m²/s, 농도는 kg/m³, 시간은 s다. 면 전체를 통과하는 적분 플럭스의 단위는 kg/s다. 농도에 셀 체적을 곱하면 kg 단위의 질량이 된다. 계산값을 비교하기 전에 이 단위 관계부터 확인하자.

이 실험에는 움직이는 격자, 반응, 가변 단면적, 벽 경계층, 운동량 방정식과의 결합이 없다. 이런 요소를 추가하면 그에 맞는 수지식을 다시 세워야 한다. 특히 셀 체적이 시간에 따라 바뀌면 저장량의 시간 미분에도 체적 변화가 들어간다. 현재의 고정 체적 공식을 그대로 사용할 수는 없다.

값을 넣기 전에 셀을 만들기

매개변수 공간을 N등분한 뒤 부드러운 주기 함수로 실제 경계 위치를 휜다. 변형 크기의 절댓값은 0.75 이하라서 사상이 증가한다. 양 끝은 정확히 0과 L로 지정한다. 셀 중심은 양 경계의 중점이며, 마지막 셀에서 첫 셀로 가는 중심 간 거리는 주기적으로 잇는다.

$$x_i=L\left[\frac{i}{N}+\frac{g}{2\pi}\sin\left(\frac{2\pi i}{N}\right)\right],\quad h_i=x_{i+1}-x_i,\quad V_i=Ah_i,\quad \gamma_i=\frac{DA}{d_i}.$$

여기서 d-i는 이웃한 두 중심 사이의 양의 거리다. 마지막 면에도 주기적으로 감싼 거리를 사용한다. 면 전도도 gamma의 단위는 m³/s다. 이 값에 양옆 셀의 농도 차이를 곱하면 확산으로 이동하는 질량 유량이 나온다.

구현에서는 일반 부동소수점 좌표로 저장한 격자를 계산 대상의 기하로 선언한다. 좌표 해상도보다 작은 변형은 저장된 격자를 바꾸지 못할 수 있다. 뒤에서 사용하는 부호·로그 표현은 작은 물리량의 곱을 보존하기 위한 것이며, 이미 반올림된 셀 좌표를 기호식이나 정확한 실수로 바꾸지는 않는다.

셀 평균은 한 점의 값이 아니라 적분이다

초기 농도는 평균 c-zero, 상대 대비 epsilon, 정수 모드 m, 주기 단위 위상 phi를 갖는 코사인 함수다. 대비는 −1에서 1 사이로 둔다. 셀에 넣는 값은 중심점에서 표본 하나를 뽑은 값이 아니라 셀 전체에 걸친 농도의 평균이다.

$$k=\frac{2\pi m}{L},\qquad c(x,0)=c_0[1+\epsilon\cos(kx+2\pi\phi)],\qquad q_i^0=c_0+c_0\epsilon\operatorname{sinc}(kh_i/2)\cos(kx_{c,i}+2\pi\phi), \quad\operatorname{sinc}(z)=\frac{\sin z}{z}.$$

넓은 셀에서는 진동의 더 넓은 구간을 평균하므로 진폭이 더 줄어든다. sinc 인자가 그 평균 효과를 나타낸다. 중심점 표본으로 대체하면 처음부터 다른 수치 문제를 풀게 된다. 높은 모드나 거친 격자에서 특히 차이가 크다. 격자를 바꿔 비교할 때도 각 격자의 실제 셀 적분으로 초기화해야 한다.

예를 들어 k와 셀 폭의 곱이 pi라면, 평균의 진동 성분에는 약 0.637인 2/pi가 곱해진다. 중심점 표본에는 이 감쇠가 없다. 이 차이는 첫 시간 단계 이전부터 존재한다. 이를 모두 시간 갱신 알고리즘의 오차로 해석하면 초기 조건 오차와 진화 오차를 섞게 된다.

각 면의 플럭스는 한 번 계산한다

면 i의 방향을 셀 i에서 다음 셀 j로 정한다. 상류 이류는 유속이 비음수이면 왼쪽 값을, 음수이면 오른쪽 값을 사용한다. 확산은 이산 농도 기울기의 반대 방향으로 물질을 옮긴다. 한 면에서 얻은 동일한 플럭스가 한쪽 셀에서는 빠지고 다른 쪽 셀에서는 더해진다.

$$F_i=aAq_{\mathrm{up},i}+\gamma_i(q_i-q_j),\qquad V_i r_i=F_{i-1}-F_i,\qquad q_i^{n+1}=q_i^n+\Delta t\,r_i^n.$$

왼쪽 면의 인덱스도 주기적으로 감싼다. 끝부분에 짝이 없는 별도 면이 생기지 않는다. 양의 플럭스는 정한 방향의 이동이고 음의 플럭스는 반대 방향의 이동이다. 유속이 음수가 되면 이류의 공급 셀이 바뀌지만 보존이라는 원칙은 바뀌지 않는다.

계산 순서는 현재 상태에서 모든 공유 면을 구하고, 모든 셀의 변화율을 조립한 다음, 모든 셀을 함께 갱신하는 것이다. 첫 셀을 먼저 바꾸고 그 새 값으로 다음 면을 계산하면 다른 알고리즘이 된다. 이 실험은 사용자가 정한 시간 간격과 정수 단계 수로 전진 Euler 갱신을 실제 수행한다.

질량 행렬로 보존성을 유도하기

행렬 B의 각 열은 한 면을 나타낸다. 출발 셀에는 +1, 도착 셀에는 −1을 넣는다. M은 실제 셀 체적을 대각 성분으로 갖는 질량 행렬이다. 닫는 면까지 포함해 B의 모든 열 합은 0이다.

$$M\dot q=-BF,\qquad \boldsymbol1^TB=0,\qquad \mathcal M=\boldsymbol1^TMq,\qquad \dot{\mathcal M}=0.$$

Euler 갱신에도 M을 곱하고 셀 전체를 합하면 정확한 산술에서 질량 변화가 0이다. 이 상쇄는 초기 코사인 함수에만 성립하는 우연이 아니라 행렬 구조의 성질이다. 학습한 면 플럭스를 사용하더라도 같은 면 값을 일관되게 공유하면 내부 전달량은 상쇄된다.

비균일 격자에서 농도의 단순 합에는 이런 물리적 의미가 없다. 작은 셀 하나와 큰 셀 하나를 똑같은 한 항으로 세기 때문이다. 실제 총량에는 체적 가중치가 필요하다. 그림으로 생각하면 농도는 막대 높이고, 질량은 높이와 실제 셀 폭과 단면적을 곱한 값이다. 높이만 더하면 폭에 관한 정보를 버리게 된다.

잘못된 저장 규칙도 독립적으로 실행하기

음성 대조는 공유 플럭스 공식을 유지하면서 모든 셀의 순유입량을 평균 체적 하나로 나눈다. 다음 단계에서는 이 잘못된 규칙으로 얻은 자체 상태를 사용해 플럭스를 다시 계산한다. 올바른 궤적에 잘못된 이름만 붙인 비교가 아니다.

$$V_{\mathrm{ref}}=\frac{1}{N}\sum_iV_i,\qquad \dot q_i^{\mathrm{wrong}}=\frac{F_{i-1}-F_i}{V_{\mathrm{ref}}},\qquad \dot{\mathcal M}_{\mathrm{wrong}} =-\sum_{i\to j}\frac{V_i-V_j}{V_{\mathrm{ref}}}F_i.$$

이 규칙은 가중하지 않은 농도 합은 보존하지만 일반적으로 실제 질량은 바꾼다. 결과가 매끄럽고 연결선 메시지가 상쇄되는 것처럼 보여도 저장 규칙이 틀릴 수 있다는 점에서 유용한 반례다. 문제는 면 값을 공유했는지가 아니라 유입 질량을 어떤 체적으로 나누었는지에 있다.

균일 격자에서는 모든 체적이 평균 체적과 같으므로 두 알고리즘이 일치한다. 이 대조를 반드시 포함하자. 여기서도 다르다면 저장 체적 외에 추가 오류가 있다는 뜻이다. 반대로 균일 격자에서의 일치만으로 불규칙한 셀을 올바르게 처리한다고 결론 내릴 수는 없다.

세 셀을 손으로 계산하기

셀 체적은 1, 2, 1 m³, 모든 면 전도도는 1 m³/s, 이류 속도는 0, 초기 농도는 1, 0, 0 kg/m³로 둔다. 정한 방향의 면 플럭스는 1, 0, −1 kg/s다. 마지막 면의 음의 플럭스는 첫 셀의 물질이 주기 경계를 건너 마지막 셀로 이동한다는 뜻이다.

항목	실제 체적 저장	평균 체적 저장
농도 변화율	−2, 1/2, 1	−3/2, 3/4, 3/4
0.1초 후 상태	0.8, 0.05, 0.1	0.85, 0.075, 0.075
0.1초 후 실제 질량	1 kg	1.075 kg
가중하지 않은 농도 합	0.95	1

올바른 갱신에서는 농도 합이 감소해도 수학적 질량은 정확히 1 kg이다. 이 감소를 누출로 해석하면 안 된다. 큰 가운데 셀은 같은 유입 질량에 대해 더 작은 농도 변화를 보인다. 잘못된 규칙은 농도 합을 그대로 유지하면서 0.075 kg을 만들어낸다.

이차 노름은 처음에 0.5 kg²/m³이고 실제 체적 규칙의 한 단계 뒤에는 0.3275 kg²/m³다. 이 손계산은 주기 격자 생성기와 별개로 대조한다. 격자 초기화가 맞다는 이유로 잘못된 갱신 규칙까지 검증된 것으로 간주하지 않기 위해서다.

import numpy as np
V = np.array([1., 2., 1.])
q = np.array([1., 0., 0.])
gamma = np.ones(3)
dt = 0.1
F = gamma * (q - np.roll(q, -1))  # zero advection
net = np.roll(F, 1) - F
correct = q + dt * net / V
wrong = q + dt * net / V.mean()
print(V @ correct, V @ wrong)  # 1.0, 1.075
print(0.5 * V @ (correct**2 - q**2))  # -0.1725

이차 노름의 예산을 유도하기

체적 가중 이차 노름을 정의하자. 농도의 크기와 변동을 측정하지만 단위가 J가 아니므로 물리적 운동에너지라고 부르면 안 된다. 수동 농도 모형에서의 수학적 이차량이다.

$$E=\frac12\sum_iV_iq_i^2,\qquad \dot E=-\sum_{i\to j}\gamma_i(q_i-q_j)^2 -\frac{|a|A}{2}\sum_{i\to j}(q_i-q_j)^2\le0.$$

확산 항은 각 셀 수지에 농도를 곱한 뒤 같은 면의 두 기여를 모으면 나온다. 면 전도도에 농도 점프의 제곱을 곱한 음의 항이 남는다. 일정 속도의 주기 상류 이류에서는 이웃한 제곱의 차이가 망원합으로 상쇄되고, 위 식의 추가 손실이 남는다.

첫 손실은 지정한 물리 확산에서 나온다. 두 번째는 일차 상류법의 수치 소산이다. D가 0이면 매끄러운 연속체 프로필은 이차 노름을 잃지 않고 평행 이동한다. 그러나 유한한 상류 격자에서는 대개 진폭이 줄어든다. 수치적 평활화와 물리 확산은 같은 현상으로 해석하지 않아야 한다.

시간 이산화가 결론을 바꾸는 이유

시간 연속 노름 미분이 비양수라는 사실만으로 유한한 Euler 단계의 노름 감소가 보장되지는 않는다. 갱신 뒤의 제곱을 전개하면 양의 추가 항이 생긴다. 이 항을 버리면 실제 실행한 알고리즘의 예산과 달라진다.

$$E(q+\Delta t r)-E(q) =\Delta t\,q^TMr+\frac{\Delta t^2}{2}r^TMr.$$

손계산 예제에서는 첫 항이 −0.2이고 양의 보정이 0.0275여서 총변화가 −0.1725다. 시간 간격을 키우면 양의 항은 제곱으로 커진다. 따라서 공유 면의 질량 장부는 닫혀 있어도 노름은 증폭될 수 있다.

비음 갱신 계수를 얻는 충분 조건은 셀마다 다르다. 유입 이웃의 계수는 비음수이고, 중심 계수도 비음수로 남아야 한다.

$$C_i=\frac{\Delta t}{V_i}\bigl(|a|A+\gamma_{i-1}+\gamma_i\bigr)\le1 \quad\text{for every cell }i.$$

이 조건에서는 새 값이 이전 이웃 값들의 볼록 결합이다. 상수 보존과 체적 가중 불변성을 함께 사용하면 Jensen 부등식으로 이차 노름의 수축을 얻는다. 가장 작은 셀이 시간 간격 제한을 정할 수 있다. 이 검사에서도 모든 체적을 평균값 하나로 바꾸면 위험한 셀을 놓칠 수 있다.

정확한 연속체 해와 비교하기

선택한 초기 조건의 연속체 해는 이동하면서 확산으로 감쇠하는 Fourier 모드다. 정확한 셀 평균에도 초기화와 같은 셀 적분을 적용한다.

$$\bar c_i(t)=c_0+c_0\epsilon e^{-Dk^2t} \operatorname{sinc}(kh_i/2)\cos(k(x_{c,i}-at)+2\pi\phi),\qquad \mathcal E_1=\sum_iV_i|q_i-\bar c_i(t)|.$$

이 참조는 주기적 상수 계수 문제에서 모든 시간에 유효하다. 중심점의 정확해도 아니고, 같은 이산 알고리즘을 다시 실행한 결과도 아니다. 절대 농도 오차를 체적으로 적분한 오차의 단위는 kg다. 최대 셀 평균 오차의 단위는 농도와 같다.

격자를 세분할 때는 모든 국소 제한을 만족하도록 시간 간격도 줄이고, 같은 물리적 최종 시간에서 비교하자. 단계 수만 맞추면 서로 다른 시간을 비교할 수 있다. 시간 간격을 고정한 채 셀만 늘리면 안정했던 거친 격자가 불안정한 미세 격자로 바뀔 수 있다. 특히 명시적 확산 제한은 셀 폭의 제곱에 비례해 강해진다.

작은 변화를 보존하고 반올림 진단을 드러내기

계산은 평균에서의 초기 편차와 누적 갱신량을 분리해 저장한다. 거의 같은 큰 농도 두 개를 뺀 뒤 작은 변화를 복원할 필요가 없다. 선형 연산자 R에 대해 누적 변화 z는 다음의 동등한 점화식을 따른다.

$$z^0=0,\qquad z^{n+1}=z^n+\Delta t\,[R(q^0-c_0\boldsymbol1)+Rz^n],\qquad q^n=c_0\boldsymbol1+(q^0-c_0\boldsymbol1)+z^n.$$

부호와 로그 크기를 사용하면 정확한 0과 일반 부동소수점으로 표현하기에 너무 작거나 큰 비영 값을 구분할 수 있다. 일반 수치가 제공되지 않는다는 이유로 그 값을 몰래 0으로 바꾸지 않는다. 확산계수와 면적, 유속과 면적 같은 작은 곱도 나눗셈이나 합산 이전에 이 표현으로 만든다.

그렇다고 모든 계산이 무한 정밀해지는 것은 아니다. 셀 기하는 반올림되어 있고 로그도 유한 정밀도다. 심한 불안정에서는 거대한 양·음 질량이 상쇄되어 직접 합산이 나쁜 조건을 가질 수 있다. 공유 면 질량 장부는 대수적 보존 항등식을 나타내고, 셀을 조립한 질량과 반올림된 총량의 직접 차감은 수치 합산 오차를 드러낸다. 이때 오차는 상쇄되는 질량들의 절댓값 규모에 비추어 평가해야 한다. 거대한 항의 차이로 얻은 작은 합이 정확하게 계산되었다고 가정하면 안 된다.

AI 모형에 연결하기

학습 모형이 방향 있는 각 면의 적분 플럭스 하나를 예측하고 마이너스 B를 곱하도록 설계할 수 있다. 여기에 실제 질량 행렬을 사용하면 전체 내부 전달량이 상쇄된다. 노드마다 농도 증분을 독립적으로 예측하면 추가 제약 없이는 같은 보장이 없다.

체적 가중치는 학습 손실에도 영향을 준다. 모든 셀을 동일하게 세면 작은 셀이 많은 영역이 같은 물리적 크기의 큰 셀 영역보다 손실에 더 큰 영향을 줄 수 있다. 체적 가중 오차는 물리 적분을 근사한다. 어떤 손실이 적합한지는 과학적 목표에 따라 정하되 가중치를 명시해야 한다.

보존은 여러 조건 중 하나다. 학습 플럭스는 상수장에서 일관성을 잃거나, 양성을 위반하거나, 시간 적분을 불안정하게 만들거나, 새 격자에서 일반화에 실패할 수 있다. 플럭스 식별에도 영공간이 있다. 이 고리의 모든 면에 같은 순환 플럭스를 더하면 B와의 곱은 바뀌지 않는다. 셀의 시간 변화만으로는 추가 정보 없이 그 순환 성분을 식별할 수 없다.

풀이가 있는 연습과 참고 자료

1. 같은 이동 질량이 왜 다른 농도 변화를 만드는가? 유입 질량을 받는 셀의 체적으로 나눈다. 0.1 kg은 1 m³ 셀에서 0.1 kg/m³, 2 m³ 셀에서 0.05 kg/m³의 변화를 만든다. 보존 대상은 질량이며 농도 증분의 동일성이 아니다.
2. 세 셀의 변화율을 다시 구하라. 면 플럭스 (1, 0, −1)의 순유입량은 (−2, 1, 1) kg/s다. 체적 (1, 2, 1)로 나누면 (−2, 1/2, 1)이다. 체적 가중 합은 −2 + 1 + 1 = 0이지만 단순 합은 −1/2다.
3. 잘못된 규칙의 질량 생성률을 구하라. 평균 체적은 4/3이다. 같은 순유입량을 이 값으로 나누면 (−3/2, 3/4, 3/4)다. 실제 체적 가중 합은 3/4 kg/s이고 0.1초에는 3/40 kg이 생성된다. 다음 단계의 플럭스는 잘못된 규칙의 자체 상태에서 다시 구해야 한다.
4. 노름 변화의 두 항을 분리하라. 같은 예제의 농도 점프 제곱 합은 2이므로 시간 연속 노름 변화율은 −2다. 체적 가중 변화율 제곱 합은 11/2다. 시간 간격 0.1에서 양의 보정은 0.01 곱하기 11/4, 즉 0.0275다. −0.2에 더하면 −0.1725다.
5. 질량 장부가 닫히면 계산이 안정한가? 아니다. 휜 격자에서 확산계수를 키우고 지나치게 큰 시간 간격을 그대로 실행하자. 공유 면은 여전히 상쇄되지만 중심 계수는 음수가 되고 농도가 음수가 되거나 양의 Euler 보정이 우세해질 수 있다. 실제 국소 계수와 결과 상태를 함께 확인해야 한다.
6. 셀 변화만으로 모든 면 플럭스를 식별할 수 있는가? 아니다. 고리의 모든 면에 같은 상수 플럭스를 더해도 어느 셀의 순유입량도 바뀌지 않는다. 자료는 플럭스 차이를 정하지만 순환 성분은 관측하지 못한다. 경계 측정, 구성 법칙 또는 추가 규약이 빠진 정보를 보완할 수 있다.

유한체적 갱신과 망원합 보존 논증의 기본 틀은 일차 자료인 Riemann Problems and Jupyter Solutions의 Approximate Solvers에서 확인할 수 있다. 유입·유출의 물리적 해석과 일정 속도의 이동 해는 같은 책의 Advection을 참고하자. 비균일 저장 반례, 이차 예산과 부호·로그 구현은 여기서 직접 특수화한 내용이며, 참고 문헌이 이 구현의 정확성을 인증한다는 뜻은 아니다.