비뚤어진 격자에서 기울기가 틀어질 때 — Green

선형장은 스텐실이 비뚤어도 여섯 자리까지 맞고, 곡면은 먼 이웃을 n=1로 눌러야 오차가 줄어든다.

잔차가 이상하게 눕더니 발산했다. 로그를 뒤지니 벽 근처 한 셀에서 온도 기울기가 실제와 거의 직각으로 틀어져 있었다. 격자를 보니 그 셀만 유독 납작했다. 값은 멀쩡했다. 기울기를 만드는 방법이 격자 모양에 걸려 넘어진 것이다.

유한체적법(FVM)에서 셀 기울기는 확산 플럭스, 2차 재구성, 리미터의 입력이다. 셀 중심값만으로는 부족하다. 이웃 정보를 모아 $\nabla\phi_C$ 를 추정해야 한다. 비교는 Green–Gauss와 최소제곱, 두 가지다. 왜곡 격자에서 한쪽이 무너지는 이유와, 거리 가중이 줄이는 오차를 출력 숫자로 따라간다.

anisotropy를 0.2 근처로 내리면 스텐실이 납작해진다. weight exponent를 0에서 1로 올리면 주황 화살표가 초록 참값 쪽으로 다시 돈다.

angle error 16.8° · |∇φ| 0.64 vs 1.08 true — the reconstructed gradient is skewed

Stencil anisotropy: 0.35Field curvature κ: 1.40Weight exponent n: 0.0Neighbor clustering: 0.60

curvature κ를 같이 올리면 두 화살 사이가 벌어진다. 납작한 템플릿에서 가중이 방향을 되돌리는지만 보면 된다.

Green–Gauss: 면으로 둘러싸 평균한다#

발산정리를 셀에 그대로 적용한다. 부피 적분을 면 적분으로 바꾼다.

\nabla \phi_C = \frac{1}{V_C} \sum_{f} \phi_f \, \mathbf{n}_f \, A_f

$V_C$ 는 셀 부피, $\mathbf{n}_f$ 는 면 바깥 법선, $A_f$ 는 면 넓이, $\phi_f$ 는 면 위 값이다.

면 값 $\phi_f$ 는 직접 알 수 없다. 셀 중심값에서 보간한다. 가장 흔한 선형 보간은 이렇다.

\phi_f = g_f\, \phi_C + (1 - g_f)\, \phi_N

$g_f$ 는 면에서 이웃 중심 $N$ 까지의 거리 비로 정한 가중치다. 계산이 싸고 직관적이다. 구조 격자에서는 2차 정확도가 나온다.

왜곡 격자가 Green–Gauss를 흔든다#

문제는 $g_f$ 에 숨어 있다. 선형 보간은 면 중심과 두 셀 중심이 한 직선 위에 있다고 가정한다. 비정렬 격자에서는 이 가정이 거의 항상 깨진다.

면 중심이 두 셀 중심을 잇는 선에서 벗어나 있으면( $\text{skewness}$ ), 보간된 $\phi_f$ 가 편향된다. 셀이 납작하거나 이웃 크기가 크게 다르면 편향은 더 커진다. 이 오차는 반복하며 쌓인다. 결국 기울기가 엉뚱한 방향을 가리킨다.

보정 항을 붙이면 완화된다. 하지만 보정에는 다시 기울기가 필요하다. 순환이 생긴다. 반복 보정은 비싸고, 왜곡이 심하면 수렴이 느려진다.

최소제곱: 이웃에 평면을 맞춘다#

발상을 바꾼다. 면을 거치지 않는다. 셀 $C$ 주변 값 분포에 곧바로 평면 하나를 맞춘다.

이웃 $k$ 가 오프셋 $\mathbf{d}_k = \mathbf{x}_k - \mathbf{x}_C$ 에 있다면, 선형 근사는 $\phi_k \approx \phi_C + \nabla\phi_C \cdot \mathbf{d}_k$ 다. 모든 이웃에 대해 이 오차의 제곱합을 최소화한다.

\min_{\nabla\phi_C} \sum_k w_k \left( \phi_k - \phi_C - \nabla\phi_C \cdot \mathbf{d}_k \right)^2

$w_k$ 는 이웃별 가중치다. 미분해서 0으로 놓으면 정규방정식이 나온다.

\left( \sum_k w_k\, \mathbf{d}_k \mathbf{d}_k^{\mathsf T} \right) \nabla\phi_C = \sum_k w_k \left( \phi_k - \phi_C \right) \mathbf{d}_k

2차원이면 $2\times2$ 시스템이다. 풀어 쓰면 이렇다.

\begin{pmatrix} \sum w\, d_x^2 & \sum w\, d_x d_y \\[2pt] \sum w\, d_x d_y & \sum w\, d_y^2 \end{pmatrix} \begin{pmatrix} \partial_x\phi \\[2pt] \partial_y\phi \end{pmatrix} = \begin{pmatrix} \sum w\, \Delta\phi\, d_x \\[2pt] \sum w\, \Delta\phi\, d_y \end{pmatrix}

$\Delta\phi = \phi_k - \phi_C$ 다. 왼쪽 행렬은 이웃 배치만으로 정해진다. 셀마다 한 번 계산해 저장할 수 있다. 장점이 크다. 면 보간 가정이 없다. 왜곡에 둔감하다. 이웃이 3개 이상이면 항상 풀린다.

면 보간이 무너지는 지점은 비뚤어진 셀의 Green–Gauss와 최소제곱과 같은 자리다. 이웃이 세 개만 있어도 정규방정식은 풀리지만, 그 세 점이 한 줄이면 가로 방향의 답은 비어 있다.

거리로 가중한다#

가중치 $w_k$ 가 정확도를 좌우한다. 흔한 선택은 거리 역수다.

w_k = \frac{1}{\lvert \mathbf{d}_k \rvert^{\,n}}

$n=0$ 이면 모든 이웃을 똑같이 본다. 멀리 있는 이웃은 값 차이가 크고, 그만큼 곡률 오차도 크다. 그 이웃이 평면 맞춤을 잡아당긴다. $n$ 을 1로 올리면 먼 이웃의 발언권이 줄어든다. 국소 선형 가정이 잘 맞는 가까운 이웃에 무게가 실린다. 비점성 문제에서는 보통 $n=1$ 을 쓴다. 먼 이웃의 표를 깎아 곡률 오차를 줄이는 선택이다.

Python — 선형장은 정확히, 곡면은 가중으로#

최소제곱의 힘을 두 번에 나눠 확인한다. 먼저 선형장 $\phi = 1.3x - 0.7y + 2$ 의 기울기를 아무렇게나 흩뿌린 스텐실에서 재구성한다. 이론상 오차가 0이어야 한다. 다음으로 곡면 $\phi = e^{0.7x}\cos(0.9y)$ 에서 가중 없는 쪽과 $n=1$ 가중을 겨룬다.

import numpy as np
 
def weighted_lsq_gradient(offsets, dphi, n):
    # offsets: (K,2) neighbor offsets d_k; dphi: (K,) diffs phi_k - phi_C
    w = 1.0 / np.maximum(np.linalg.norm(offsets, axis=1), 1e-12) ** n
    M = (offsets[:, :, None] * offsets[:, None, :] * w[:, None, None]).sum(axis=0)
    b = (offsets * (dphi * w)[:, None]).sum(axis=0)
    return np.linalg.solve(M, b)
 
xc = np.array([0.4, -0.2])
 
# (A) 선형장 — 어떤 스텐실이든 최소제곱은 정확하다
lin = lambda x, y: 1.30 * x - 0.70 * y + 2.0
scatter = np.array([[0.45, 0.30], [-0.40, 0.12], [0.05, -0.48], [-0.33, -0.29],
                    [0.90, 0.15], [-0.11, 0.44], [0.28, -0.09]])   # 불규칙 스텐실
dlin = np.array([lin(*(xc + d)) - lin(*xc) for d in scatter])
gA = weighted_lsq_gradient(scatter, dlin, 1.0)
print(f"linear field   LSQ = ({gA[0]:+.6f}, {gA[1]:+.6f})   exact = (+1.300000, -0.700000)")
 
# (B) 곡면 — 거리 가중이 재구성 오차를 깎는다
phi  = lambda x, y: np.exp(0.7*x) * np.cos(0.9*y)
grad = lambda x, y: np.array([0.7*np.exp(0.7*x)*np.cos(0.9*y),
                             -0.9*np.exp(0.7*x)*np.sin(0.9*y)])
near = np.array([[0.10, 0.05], [-0.08, 0.09], [0.06, -0.10], [-0.09, -0.06], [0.11, 0.10]])
offs = np.vstack([near, [0.90, 0.70]])          # 멀리 떨어진 이웃 하나
dphi = np.array([phi(*(xc + d)) - phi(*xc) for d in offs])
gt   = grad(*xc)
for n in (0.0, 1.0):
    g = weighted_lsq_gradient(offs, dphi, n)
    err = np.linalg.norm(g - gt)
    print(f"curved field   n={n:.0f}  LSQ = ({g[0]:+.3f}, {g[1]:+.3f})   |error| = {err:.4f}")
print(f"curved field   true = ({gt[0]:+.3f}, {gt[1]:+.3f})")

출력은 이렇게 나온다.

linear field   LSQ = (+1.300000, -0.700000)   exact = (+1.300000, -0.700000)
curved field   n=0  LSQ = (+0.885, +0.200)   |error| = 0.0293
curved field   n=1  LSQ = (+0.891, +0.205)   |error| = 0.0222
curved field   true = (+0.911, +0.213)

선형장에서는 스텐실이 아무리 뒤틀려도 소수점 여섯 자리까지 정확하다. 출력이 보여 주는 (+1.300000, -0.700000)이 그 자리다. 선형 함수는 이웃이 어디에 있든 차분과 기울기의 내적이 딱 맞는다. 불규칙 스텐실에 (0.90, 0.15)가 끼어 있어도 오차는 생기지 않는다. 왜곡은 이 장에서 기울기를 밀지 못한다.

곡면은 다르다. 참값은 (+0.911, +0.213)이다. n=0의 (+0.885, +0.200)은 x가 0.026, y가 0.013 모자란다. n=1의 (+0.891, +0.205)는 모자람이 0.020과 0.008이다. 노름은 0.0293에서 0.0222로 내려간다. 차이 0.0071은 약 24%라서, 4분의 1이라는 말과 같다. 먼 이웃 (0.90, 0.70)의 성분은 가까운 다섯 점의 0.1 안팎보다 훨씬 크다. n=0이면 그 점이 같은 표로 평면을 끌어당긴다. n=1이면 거리가 표를 깎는다.

이웃이 한 줄로 서면 역행렬을 믿지 마라#

최소제곱이 왜곡에 강하다고 만능은 아니다. 왼쪽 행렬 $M = \sum w_k \mathbf{d}_k \mathbf{d}_k^{\mathsf T}$ 가 병들 수 있다. 이웃이 한 직선 위에 늘어서면 $M$ 이 특이에 가까워진다. 스텐실 방향으로는 기울기가 잘 정해지지만, 가로 방향으로는 거의 자유롭다. 조건수가 폭발한다.

Collinearity를 0.9 근처로 올리면 이웃이 한 줄로 모이고, cond(M)이 수십에서 무한대로 치솟는다.

cond(M) = 1.4 (λ_max 3.14 / λ_min 2.271) — well-conditioned

Collinearity: 0.20Neighbor count: 6

타원이 길고 얇게 늘어나는 축이 기울기를 못 정하는 축이다. 벽에 붙은 경계층 셀이 그 모양이다. 대각은 오프셋 제곱의 합이라, 이웃이 한 방향으로만 길면 다른 대각이 빈다. 경계 셀은 이웃이 유체 쪽으로 쏠린다. 벽 반대편 이웃이나 면 중심을 스텐실에 넣으면 빈 축이 조금 채워지고 조건수가 내려간다. 격자가 고정이면 $M$ 은 속도장이 바뀌어도 그대로다. 매 스텝 다시 세우지 말고 역행렬이나 QR을 저장한다. 거리 가중 $n$ 은 문제마다 다르다. 점성 항이 크면 $n=2$ 가 나을 때도 있다.

이동 최소제곱은 이 가중을 이웃마다 다시 짜는 쪽으로 한 걸음 더 간다. Green–Gauss는 싸고 스큐에 약하다. 가중 최소제곱은 스큐에 강하지만, 스텐실이 한 줄이면 행렬이 먼저 병든다. 방법을 고르기 전에 그 셀의 이웃이 평면에 퍼져 있는지를 본다.

비뚤어진 격자에서 기울기가 틀어질 때 — Green–Gauss와 가중 최소제곱 재구성